††thanks: These authors contribute equally.††thanks: These authors contribute
equally.

# Spin-charge separation and unconventional superconductivity in t-J model on
honeycomb lattice

Jian-Jian Miao  Department of Physics, The Chinese University of Hong Kong,
Sha Tin, New Territories, Hong Kong, China  Zheng-Yuan Yue  Department of
Physics, The Chinese University of Hong Kong, Sha Tin, New Territories, Hong
Kong, China  Hao Zhang  School of Physics and Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  Wei-Qiang Chen<EMAIL_ADDRESS>Department of Physics and Shenzhen Institute for Quantum Science and
Engineering, Southern University of Science and Technology, Shenzhen 518055,
China  Shenzhen Key Laboratory of Advanced Quantum Functional Materials and
Devices, Southern University of Science and Technology, Shenzhen 518055, China
Zheng-Cheng Gu<EMAIL_ADDRESS>Department of Physics, The Chinese
University of Hong Kong, Sha Tin, New Territories, Hong Kong, China

###### Abstract

The physical nature of doped Mott-insulator has been intensively studied for
more than three decades. It is well known that the single band Hubbard model
or $t$-$J$ model on the bipartite lattice is the simplest model to describe a
doped Mott insulator. Unfortunately, the key mechanism of superconductivity in
these toy models is still under debate. In this paper, we propose a new
mechanism for the $d+id$-wave superconductivity (SC) that occurs in the small-
doping region of the honeycomb lattice $t$-$J$ model based on the Grassmann
tensor product state numerical simulation and spin-charge separation
formulation. Moreover, in the presence of anti-ferromagnetic order, a
continuum effective field theory for holons is developed near half-filling. It
reveals the competition between repulsive and attractive holon interactions
induced by spinon fluctuations and gauge fluctuations, respectively. At a
large value of $t/J$, the repulsive interaction dominates, leading to the non-
Fermi liquid like behavior; while in a moderate range of $t/J$, the attractive
interaction dominates, leading to the SC order. Possible experimental
detection of spin-charge separation phenomena is also discussed.

##  I Introduction

Since the discovery of cuprates Bednorz and Müller (1986), the mechanism of
high-$T_{c}$ superconductivity (SC) remains one of the long-standing hardcore
problems in condensed matter physics. The $t$-$J$ model is a well-accepted
microscopic model which potentially captures the essential physics of
$\mathrm{CuO_{2}}$ layers in high-$T_{c}$ cuprates, and can be derived from
the large-$U$ limit of three-band Hubbard model for the $\mathrm{CuO_{2}}$
layers Zhang and Rice (1988). The $t$-$J$ model Hamiltonian reads:

$\displaystyle H=-t{\sum_{\langle
ij\rangle,\sigma}}(\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}+h.c.)+J{\sum_{\langle
ij\rangle}}\left(\mathbf{S}_{i}\cdot\mathbf{S}_{j}-\frac{1}{4}n_{i}n_{j}\right),$
(1)

where
$n_{i}=\sum_{\sigma}c_{i\sigma}^{\dagger}c_{i\sigma}\equiv\sum_{\sigma}n_{i\sigma}$
is the electron density operator,
$\hat{c}_{i\sigma}=(1-n_{i,-\sigma})c_{i\sigma}$ is the electron annihilation
operator defined in no-double-occupancy subspace, and
$\mathbf{S}_{i}=(1/2)\sum_{\alpha\beta}c^{\dagger}_{i\alpha}\boldsymbol{\sigma}_{\alpha\beta}c_{i\beta}$
is the spin-1/2 operator. The $t$-term permits the motion of holes while the
$J$-term is the superexchange interaction.

Even though having been studied for decades, the global phase diagram of the
$t$-$J$ model on a square lattice is still controversial. Recently, the study
of the $t$-$J$ model on a honeycomb lattice attracts a lot of interest since
these two types of lattices share similar features – both of them are
bipartite lattices that stabilize the anti-ferromagnetic (AFM) order at half-
filling. Remarkably, Grassmann tensor product state methods discovered the
emergence of uniform $d+id$-wave SC order at very small doping, while the
stripe and the $d$-wave SC orders coexist at relatively large doping in the
honeycomb lattice $t$-$J$ model Gu  _et al._ (2013); Xu  _et al._ (2022). In
the weak coupling limit, the uniform $d+id$ superconductivity may also occur
in the doped graphene systems revealed by various numerical methods, such as
renormalized mean-field theory Black-Schaffer  _et al._ (2014), quantum Monte
Carlo Pathak  _et al._ (2010); Ma  _et al._ (2011); Ying and Wessel (2018);
Jiang  _et al._ (2014); Wang  _et al._ (2012), renormalization group Wang  _et
al._ (2012); Kiesel  _et al._ (2012); Nandkishore  _et al._ (2012), and
dynamical cluster approximation Xu  _et al._ (2016). Other pairing symmetries,
such as $s$ and $p+ip$ waves, are also discovered Faye  _et al._ (2015); Uchoa
and Castro Neto (2007); Gu  _et al._ (2020). Apparently, the phase diagram of
the $t$-$J$ model on a honeycomb lattice can provide deep insights into the
key mechanism of high-$T_{c}$ superconductivity on a square lattice, which is
closely related to realistic experimental materials.

Three decades ago, P. W. Anderson proposed the resonating valence bond (RVB)
picture to account for the high-$T_{c}$ superconductivity Anderson (1987), and
the corresponding spin-charge separation scenario is systematically studied in
terms of slave-boson formalism for the $t$-$J$ model Lee  _et al._ (2006).
Considering the AFM order as a characteristic feature in the small doping
region, the slave-fermion formalism is a more appealing approach for small
doping; the condensation of Schwinger bosons leads to a long-range magnetic
order Auerbach (1998). Nevertheless, previous slave-fermion studies of the
square lattice $t$-$J$ model suggested the spiral order at finite doping
instead of the AFM order Jayaprakash  _et al._ (1989), which contradicts
numerical and experimental results. However, as there are two sites per unit
cell on honeycomb lattice, mean-field solutions with AFM order is still
possible at finite doping. We believe that such kind of translation invariant
mean-field theory might be crucial for the emergence of uniform $d+id$-wave SC
order.

In this paper, we first use the state-of-the-art Grassmann tensor product
state numerical method to obtain the global phase diagram of $t$-$J$ model on
honeycomb lattice at finite doping. We adopt the slave-fermion mean-field
approach and it gives rise to an AFM order on a honeycomb lattice that is
consistent with our numerical results. At small doping, we find coexisting
short-range AFM and FM correlations, and holon mobility. The mutual dependence
of ferromagnetic (FM) correlation and holon mobility is also reproduced, which
is consistent with the physics of the Nagaoka state in the infinite-$U$
Hubbard model near half-filling Nagaoka (1966). We further employ the
functional field integral formalism to derive the effective interacting holon
theory. The effective holon Hamiltonian contains repulsive interaction
generated by exchanging the spinons and attractive interaction via gauge
fluctuations. For a moderate $t/J$, the attractive interaction overcomes the
repulsive one and leads to the SC orders. We stress that the attractive
interaction is only possible in the systems with spin-charge separation. We
hope that such a mechanism can shed light on the emergence of SC orders in the
$t$-$J$ model on the square lattice as well and help us to understand the
origin of high temperature superconductivity in cuprates.

The rest of the paper is organized as follows: Sec. II presents the Grassmann
tensor product state numerical results, including the phase diagram at finite
doping of the $t$-$J$ model on a honeycomb lattice. The slave-fermion mean-
field theory of the $t$-$J$ model on honeycomb lattice is presented in Sec.
III. In Sec. IV, we derive the low-energy effective interacting holon theory
by including the spinon-holon coupling and gauge fluctuations. Finally, we
summarize the results and discuss the potential experimental detection of
spin-charge separation phenomena in Sec. V.

Figure 1: The honeycomb lattice. The two triangular sub-lattices (with basis
vectors $\\{a_{n}\\}_{n=1,2}$) are named $A,B$ and colored in green, red
respectively. A unit cell is shown by the shaded area.
$\\{\delta_{n}\\}_{n=1,2,3}$ are nearest neighbor vectors (with
$\delta_{1}=a\hat{x}$, $a$ being the nearest neighbor distance).
$\\{\eta_{n}\\}_{n=1,2,3}$ are vectors to the unit cells containing the
nearest neighbors, which are related to $\delta_{n}$ by
$\eta_{n}=\delta_{n}-a\hat{x}$ (in particular $\eta_{1}=0$). Figure 2: Ground
state energy versus doping at $t/J=20$.

Figure 3: Staggered magnetization, amplitudes of singlet and triplet SC order
parameters versus doping at $t/J=20$.

##  II Grassmann Tensor Product State Numerical results

We first use the state-of-the-art Grassmann tensor product state method to
investigate the ground state phase diagram of the $t$-$J$ model on a honeycomb
lattice. The ground state wave function is obtained via the so-called simple
update (SU) scheme, and we also choose different bond dimensions $D=10,12,14$
to examine the $D$ dependence of the ground state energy and physical
measurement such as the staggered magnetization $m=\sqrt{\langle
S^{x}_{i}\rangle^{2}+\langle S^{y}_{i}\rangle^{2}+\langle
S^{z}_{i}\rangle^{2}}$ (with $\langle{\mathbf{S}_{i\in
A}}\rangle=-\langle{\mathbf{S}_{i\in B}}\rangle$; $A,B$ are the two sub-
lattices of the honeycomb lattice), singlet superconductivity (SC) order
parameters $\Delta^{n}_{s}=\langle
c_{i\uparrow}c_{i+\delta_{n}\downarrow}-c_{i\downarrow}c_{i+\delta_{n}\uparrow}\rangle/\sqrt{2}$
and triplet SC order parameters
$\boldsymbol{\Delta}^{n}_{t}=\langle\sum_{\alpha,\beta}c_{i\alpha}(i\sigma^{y}\boldsymbol{\sigma})_{\alpha\beta}c_{i+\delta_{n},\beta}\rangle/\sqrt{2}$
($\\{\delta_{n}\\}_{n=1,2,3}$ are the 3 nearest neighbor vectors of a site
$i\in A$, as shown in Fig. 1).

The calculated singlet SC order parameters show $d+id$ pairing, characterized
by:
$\Delta^{2}_{s}/\Delta^{1}_{s}\simeq\Delta^{3}_{s}/\Delta^{2}_{s}\simeq\Delta^{1}_{s}/\Delta^{3}_{s}\simeq
e^{2\pi i/3}$. Similar to a previous study of $t/J=3$ caseXu  _et al._ (2022),
we find the ground state energy shows very good convergence for different
choices of $D$, e.g., for the case with $t/J=20$ in Fig. 2 . Nevertheless,
very different from the $t/J=3$ case, we find that both magnetization and SC
order parameters have a rather strong $D$ dependence, especially for larger
doping. The calculated magnetization, magnitudes of singlet and triplet SC
order parameters (averaged over the 3 nearest neighbors) at $t/J=20$ are shown
in Fig. 3. After proper extrapolation to the infinite $D$ limit (see more
details Appendix A), we estimate that the stagger magnetization will vanish
for $\delta>0.03$ and the singlet SC order parameter will vanish for
$\delta>0.05$. The triplet SC order parameter has a rather strong $D$
dependence even at very small doping, and our simple extrapolation suggests it
should be very close to zero at any doping.

By carefully analyzing the data at other $t/J$ ratios, we obtain the
approximate phase diagram at finite doping as shown in Fig. 4. Exactly at
half-filling, the Hamiltonian of $t$-$J$ model reduces to the AFM Heisenberg
model. Upon doping, the AFM order is suppressed and finally vanishes, while
the $d+id$ SC orders emerge. The singlet SC order dominates and exists
persistently to relatively large doping (for each fixed $t/J$). Interestingly,
for large enough $t/J$ at each fixed doping, both AFM and SC orders eventually
vanish. We conjecture that the system should enter a potentially non-Fermi
liquid phase, and we shall explore the physical nature of such an exotic phase
in the next two sections. The raw data determining the phase diagram are given
in Appendix A.

Figure 4: The phase diagram of $t$-$J$ model on honeycomb lattice at finite
doping. AFM, SC, and NFL denote antiferromagnetic order, superconducting order
and non-Fermi liquid phase respectively.

##  III Slave-fermion approach

###  III.1 Mean-field theory

Figure 5:  Solution to the mean-field equations for doping $0\leq\delta\leq
0.2$ on a system of $64\times 64$ unit cells at low temperature $\beta=800$.
(a) The values of $A,B$ and $-D$ at fixed $t/J=10$. A second-order phase
transition (indicated by a vertical dotted line) happens at $\delta=0.09$. (b)
The values of $A,B$ and $-D$ at fixed doping $\delta=0.05$. Phase transitions
(indicated by horizontal dotted lines) happens at $t/J\approx 0.3$ (first-
order) and $18.7$ (second-order). (c) The mean-field phase diagram.

To gain more physical understanding for the result from numerical simulations,
we first adopt the slave-fermion mean-field theory to analytically study the
$t$-$J$ model on the honeycomb lattice. In the slave-fermion representation,
the electron operator can be decomposed into fermionic holons and bosonic
spinons in the no-double-occupancy subspace as:

$c^{\dagger}_{i\sigma}=b^{\dagger}_{i\sigma}h_{i},$ (2)

where $b_{i\sigma}$ are Schwinger-boson (spinon) operators and $h_{i}$ are
slave-fermion (holon) operators, with no-double-occupancy constraint at each
site $i$:

${\textstyle\sum_{\sigma}}b_{i\sigma}^{\dagger}b_{i\sigma}+h^{\dagger}_{i}h_{i}=1.$
(3)

In the no-double-occupancy subspace, the spin and electron density operators
on each site can be expressed in spinons only:

$\mathbf{S}_{i}=\frac{1}{2}\sum_{\alpha,\beta}b^{\dagger}_{i\alpha}\boldsymbol{\sigma}_{\alpha\beta}b_{i\beta},\quad
n_{i}=\hat{n}^{b}_{i}\equiv\sum_{\sigma}b^{\dagger}_{i\sigma}b_{i\sigma}.$ (4)

Then we can re-express the $t$-$J$ model Eq. (1) as:

$\displaystyle
H=t\sum_{\langle{ij}\rangle,\sigma}(h^{\dagger}_{j}h_{i}b^{\dagger}_{j\sigma}b_{i\sigma}+h.c.)-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}b^{\dagger}_{j,-\sigma}b^{\dagger}_{i\sigma}b_{i\sigma^{\prime}}b_{j,-\sigma^{\prime}}.$
(5)

Through out the whole paper, we label spin-up as $+1$, and spin-down as $-1$.
In the summation over bonds, we choose $i$ to be on sub-lattice $A$. To obtain
the mean-field theory, the no-double-occupancy constraint is only imposed _on
average_ :

$\langle{h^{\dagger}_{i}h_{i}}\rangle=\delta,\quad{\textstyle\sum_{\sigma}}\langle{b^{\dagger}_{i\sigma}b_{i\sigma}}\rangle=1-\delta\quad(0\leq\delta\leq
1),$ (6)

which is done by introducing chemical potentials $\lambda_{i}$ for spinons and
$\mu_{i}$ for holons at each site. We further define the spinon pairing (RVB)
operators and spinon hopping operators on each bond $\langle{ij}\rangle$ as:

$\hat{A}_{ij}=\frac{1}{2}\sum_{\sigma}\sigma
b_{i\sigma}b_{j,-\sigma},\quad\hat{B}_{ij}=\frac{1}{2}\sum_{\sigma}b_{i\sigma}^{\dagger}b_{j\sigma}.$
(7)

Both are invariant under global $SU(2)$ spin rotations. The mean-field ansatz
reads:

$A_{ij}=\langle{\hat{A}_{ij}}\rangle,\quad
B_{ij}=\langle{\hat{B}_{ij}}\rangle,\quad
D_{ij}=\langle{h^{\dagger}_{i}h_{j}}\rangle.$ (8)

Here $\langle{\,\cdot\,}\rangle$ refers to the mean-field average, and $j$ is
the nearest neighbor of $i$. All these parameters are in general complex
numbers. We choose $i\in A$ to avoid redundancy, since $A_{ji}=-A_{ij}$,
$B_{ji}=B^{*}_{ij}$, and $D_{ji}=D^{*}_{ij}$. Physically, $A$ and $B$
represent the short-range AFM and FM correlations respectively, and $D$
represents the mobility of holons Jayaprakash  _et al._ (1989). In addition,
we further require that:

$\displaystyle\langle{b_{i\sigma}b_{j\sigma^{\prime}}}\rangle$
$\displaystyle=\sigma A_{ij}\delta_{\sigma,-\sigma^{\prime}},$ (9)
$\displaystyle\langle{b^{\dagger}_{i\sigma}b_{j\sigma^{\prime}}}\rangle$
$\displaystyle=B_{ij}\delta_{\sigma\sigma^{\prime}},$
$\displaystyle\langle{b^{\dagger}_{i\sigma}b_{i\sigma^{\prime}}}\rangle$
$\displaystyle=(1/2)(1-\delta)\delta_{\sigma\sigma^{\prime}}.$

The mean-field Hamiltonian becomes (see details in Appendix B.1):

$\displaystyle H_{\text{MF}}$ $\displaystyle=H_{h}+H_{b}+H_{0},$ (10)
$\displaystyle H_{h}$
$\displaystyle=2t\sum_{\langle{ij}\rangle}(h^{\dagger}_{i}h_{j}B^{*}_{ij}+h.c.)-\sum_{i}\mu_{i}h^{\dagger}_{i}h_{i},$
(11) $\displaystyle H_{b}$
$\displaystyle=2t\sum_{\langle{ij}\rangle}(D_{ij}\hat{B}^{\dagger}_{ij}+h.c.)+\sum_{i,\sigma}\lambda_{i}b^{\dagger}_{i\sigma}b_{i\sigma}$
$\displaystyle\quad+J\sum_{\langle{ij}\rangle}(B^{*}_{ij}\hat{B}_{ij}-2A^{*}_{ij}\hat{A}_{ij}+h.c.),$
(12) $\displaystyle H_{0}$
$\displaystyle=-2t\sum_{\langle{ij}\rangle}(D_{ij}B^{*}_{ij}+h.c.)-J\sum_{\langle{ij}\rangle}(|B_{ij}|^{2}-2|A_{ij}|^{2})$
$\displaystyle\quad+\sum_{i}[\mu_{i}\delta-\lambda_{i}(1-\delta)]-\frac{3}{4}NJ(1-\delta)^{2},$
(13)

where $N$ is the number of _unit cells_. To keep the translation symmetry and
$C_{3}$ rotation symmetry on the honeycomb lattice, we choose the _real and
uniform_ ansatz:

$\displaystyle A_{ij}=A,\ \ B_{ij}=B,\ \ D_{ij}=D,\ \ \mu_{i}=\mu,\ \
\lambda_{i}=\lambda.$ (14)

Then we perform Fourier transform

$\displaystyle h_{i\in s}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{k}h^{s}_{k}e^{ik\cdot R_{i}},$
$\displaystyle b_{i\in s,\sigma}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{k}b^{s}_{k\sigma}e^{ik\cdot R_{i}}.$

where $s=A,B$ labels the sub-lattice, and $R_{i}$ is the position of the unit
cell containing site $i$. Defining

$\displaystyle h_{k}=\begin{bmatrix}h^{A}_{k}\\\\[3.00003pt]
h^{B}_{k}\end{bmatrix},\quad
b_{k\sigma}=\begin{bmatrix}b^{A}_{k\sigma}\\\\[3.00003pt]
b^{B}_{k\sigma}\end{bmatrix},$ (15) $\displaystyle r=tB,\quad p=tD+JB/2,\quad
q=JA,$ (16)

the mean-field Hamiltonian in momentum space is:

$\displaystyle H_{h}$ $\displaystyle=\sum_{k}h^{\dagger}_{k}H^{h}_{k}h_{k},$
(17) $\displaystyle H_{b}$
$\displaystyle=\sum_{k}\begin{bmatrix}b^{\dagger}_{k\uparrow}&b_{-k\downarrow}\end{bmatrix}H^{b}_{k}\begin{bmatrix}b_{k\uparrow}\\\
b^{\dagger}_{-k\downarrow}\end{bmatrix}-2N\lambda,$ (18)

where the Hermitian matrices $H^{h}_{k},H^{b}_{k}$ are defined as:

$\displaystyle H^{h}_{k}$ $\displaystyle=\begin{bmatrix}-\mu&\xi^{f}_{k}\\\
\xi^{f*}_{k}&-\mu\end{bmatrix},$ (19) $\displaystyle H^{b}_{k}$
$\displaystyle=\begin{bmatrix}\lambda&\xi^{b}_{k}&&-\Delta_{k}\\\
\xi^{b*}_{k}&\lambda&\Delta^{*}_{k}&\\\ &\Delta_{k}&\lambda&\xi^{b}_{k}\\\
-\Delta^{*}_{k}&&\xi^{b*}_{k}&\lambda\end{bmatrix},$ (20)
$\displaystyle\xi^{f}_{k}$
$\displaystyle=2r\gamma_{k},\quad\xi^{b}_{k}=p\gamma_{k},\quad\Delta_{k}=q\gamma_{k}.$
(21)

Here $\gamma_{k}=\sum_{\eta}e^{ik\cdot\eta}$ and $\eta$ sums over vectors to
the 3 unit cells containing nearest neighbor sites (see Fig. 1). To
diagonalize Eq. (17), we define quasi-holons $f^{s}_{k}$ (fermions) by a
unitary transformation:

$\begin{bmatrix}h^{A}_{k}\\\\[1.00006pt]
h^{B}_{k}\end{bmatrix}=U_{k}\begin{bmatrix}f^{A}_{k}\\\\[1.00006pt]
f^{B}_{k}\end{bmatrix},\quad
U_{k}\equiv\frac{1}{\sqrt{2}}\begin{bmatrix}s_{r}&e^{i\theta_{k}}\\\
e^{-i\theta_{k}}&-s_{r}\end{bmatrix},$ (22)

where $\theta_{k}$ is the phase of $\gamma_{k}$ (i.e.
$\gamma_{k}=|\gamma_{k}|e^{i\theta_{k}}$), and
$s_{r}\equiv\operatorname{sgn}(r)$, with the definition
$\operatorname{sgn}(0)=1$. To diagonalize Eq. (18), we define quasi-spinons
$\beta^{s}_{k\sigma}$ (bosons) by a Bogoliubov transformation

$\begin{gathered}\begin{bmatrix}b_{k\uparrow}\\\
b^{\dagger}_{-k\downarrow}\end{bmatrix}=W_{k}\begin{bmatrix}\beta_{k\uparrow}\\\
\beta^{\dagger}_{-k\downarrow}\end{bmatrix},\\\
W_{k}=\frac{1}{\sqrt{2}}{\small\begin{bmatrix}u_{k}s_{q}e^{i\theta_{k}}&u_{k}s_{q}e^{i\theta_{k}}&v_{k}s_{q}e^{i\theta_{k}}&v_{k}s_{q}e^{i\theta_{k}}\\\
u_{k}s_{p}s_{q}&-u_{k}s_{p}s_{q}&-v_{k}s_{p}s_{q}&v_{k}s_{p}s_{q}\\\
-v_{k}s_{p}e^{i\theta_{k}}&v_{k}s_{p}e^{i\theta_{k}}&u_{k}s_{p}e^{i\theta_{k}}&-u_{k}s_{p}e^{i\theta_{k}}\\\
v_{k}&v_{k}&u_{k}&u_{k}\end{bmatrix}},\end{gathered}$ (23)

where $s_{p}=\operatorname{sgn}(p)$, $s_{q}=\operatorname{sgn}(q)$ and

$\displaystyle
u_{k}=\bigg{(}\frac{\lambda+\chi_{k}}{2\chi_{k}}\bigg{)}^{1/2},\quad
v_{k}=\bigg{(}\frac{\lambda-\chi_{k}}{2\chi_{k}}\bigg{)}^{1/2},$ (24)
$\displaystyle\chi_{k}=\sqrt{\lambda^{2}-|\Delta_{k}|^{2}}.$ (25)

The diagonalized mean-field Hamiltonian reads:

$\displaystyle H_{\text{MF}}$
$\displaystyle=\sum_{k,s}E^{f}_{ks}f^{s\dagger}_{k}f^{s}_{k}+\sum_{k,s,\sigma}E^{b}_{ks}\beta^{s\dagger}_{k\sigma}\beta^{s}_{k\sigma}$
$\displaystyle\quad+H_{0}-2N\lambda+\sum_{k,s}E^{b}_{ks},$ (26)

with quasi-holon $f^{s}_{k}$ and quasi-spinon $\beta^{s}_{k\sigma}$
dispersion:

$\displaystyle E_{k\pm}^{f}$ $\displaystyle=\pm|\xi^{f}_{k}|-\mu,$ (27)
$\displaystyle E_{k\pm}^{b}$ $\displaystyle=\pm|\xi^{b}_{k}|+\chi_{k}.$ (28)

The self-consistency equations are derived by minimizing the mean-field free
energy with respect to $A$, $B$, $D$, $\mu$ and $\lambda$:

$\displaystyle\delta$
$\displaystyle=\frac{1}{2N}\sum_{k}[n_{f}(E^{f}_{k+})+n_{f}(E^{f}_{k-})],$
(29) $\displaystyle 1-\delta$
$\displaystyle=\frac{1}{N}\sum_{k}\Big{\\{}\frac{\lambda}{\chi_{k}}[1+n_{b}(E^{b}_{k+})+n_{b}(E^{b}_{k-})]-1\Big{\\}},$
(30) $\displaystyle A$ $\displaystyle=\frac{q}{2\alpha
N}\sum_{k}\frac{|\gamma_{k}|^{2}}{\chi_{k}}[1+n_{b}(E^{b}_{k+})+n_{b}(E^{b}_{k-})],$
(31) $\displaystyle B$ $\displaystyle=\frac{\operatorname{sgn}(p)}{2\alpha
N}\sum_{k}|\gamma_{k}|[n_{b}(E^{b}_{k+})-n_{b}(E^{b}_{k-})],$ (32)
$\displaystyle D$ $\displaystyle=\frac{\operatorname{sgn}(r)}{2\alpha
N}\sum_{k}|\gamma_{k}|[n_{f}(E^{f}_{k+})-n_{f}(E^{f}_{k-})].$ (33)

Here $\alpha=3$ is the coordination number; $n_{b}$ and $n_{f}$ are usual
boson and fermion distributions defined as:

$n_{b}(x)=\frac{1}{e^{\beta x}-1},\quad n_{f}(x)=\frac{1}{e^{\beta x}+1}.$
(34)

Fig. 5 shows the solution of the self-consistency equations on a honeycomb
lattice of $64\times 64$ unit cells with periodic boundary condition,
calculated at $\beta=800$, which is large enough to approximate the zero
temperature behavior. Similar to previous slave-fermion mean-field theory on
the square lattice Jayaprakash  _et al._ (1989); Yoshioka (1989), we obtain 3
phases in the small-doping region. When $t/J$ is small or at zero doping,
there is an AFM phase in which $A>0$ and $B,D=0$. As $t/J$ increases above a
certain value at finite doping, the system enters a phase with $A,B,D\neq 0$
via a first-order transition, where short-range FM and AFM correlations
coexist. Note that the minimum of the spinon dispersion $E^{b}_{k-}$ is always
at $k=0$; thus at zero temperature (when spinon condensation occurs), the
spiral magnetic order is avoided, in contrast to the mean-field theory on the
square lattice Jayaprakash  _et al._ (1989). The behavior of $A,B,D$ in this
phase is shown in Fig. 5 (a) and (b). In the upper-right part of the phase
diagram there is an FM phase with $A=0$ and $B,D\neq 0$, which is separated
from the $A,B,D\neq 0$ phase by a second-order transition. This phase appears
because of a large holon mobility, similar to the physics of Nagaoka states in
the Hubbard model with infinite $U$ and near half-filling. However, the FM
phase is not observed from numerical simulations at large but finite $t/J$. We
believe that the long-range FM order in this phase is an artifact of the mean-
field theory, and the FM order is expected to be eliminated by gauge
fluctuations induced by no-double-occupancy constraints. Below we shall
discuss more details for the gauge structure for mean-field Hamiltonian.

###  III.2 Gauge structure of the mean-field Hamiltonian

Since the mean field theory does not impose the no-double-occupancy constraint
on each site exactly, we must consider the gauge fluctuations to restore such
a nontrivial constraint for the low energy subspace. In the slave-fermion
representation Eq. (2), the physical electron operator $c_{i\sigma}$ is
invariant under the local $U(1)$ gauge transformation, i.e., for any site
$i\in A,B$:

$h_{i}\to h_{i}e^{-i\theta_{i}},\quad b_{i\sigma}\to
b_{i\sigma}e^{-i\theta_{i}}.$ (35)

After projection into the no-double-occupancy physical space, each of the
three phases from the mean field theory is characterized by its invariant
gauge group (IGG) Wen (2002); Wang (2010), which consists of local gauge
transformations of the form Eq. (35) that leave the ansatz unchanged. The FM
phase ($A=0$ and $B,D\neq 0$) has the usual $U(1)$ gauge structure Eq. (35),
with IGG $=U(1)$. However, the AFM phase ($A\neq 0$ and $B,D=0$) has a
_staggered_ $U(1)$ gauge structure

$h_{i\in s}\rightarrow h_{i}e^{-is\theta_{i}},\quad b_{i\in
s,\sigma}\rightarrow b_{i\sigma}e^{-is\theta_{i}}.$ (36)

where we interpret $s=+1$ on sub-lattice $A$, and $s=-1$ on $B$. This opposite
sign on two sub-lattices indicate opposite gauge charges on $A$ and $B$ sub-
lattices. In the next section, we will show that such a staggered $U(1)$ gauge
structure eventually leads to an effective attraction between holons on
different sub-lattices, and is crucial for the emergence of superconductivity.

Finally, the $A,B,D\neq 0$ phase has a $\mathbb{Z}_{2}$ gauge structure

$h_{i}\to s_{i}h_{i},\quad b_{i\sigma}\to s_{i}b_{i\sigma}.$ (37)

where $s_{i}=\pm 1$ is a sign factor depending on the site $i$, and the
corresponding IGG is $\mathbb{Z}_{2}$. In this phase, besides the holon and
spinon quasi-particle excitations, there is another excitation corresponding
to the $\mathbb{Z}_{2}$ flux of $\mathbb{Z}_{2}$ gauge fields, dubbed as vison
Senthil and Fisher (2000). The $\mathbb{Z}_{2}$ gauge fluctuation corresponds
to phase fluctuations of order parameters. It is coupled to the holon and
spinon in a $\mathbb{Z}_{2}$ gauge invariant way and has a finite gap(if we
assume the spinon does not condense), which does not confine the holon and
spinon; hence it does not change the qualitative results obtained in this
section. The existence of the vison excitation is an indication of
fractionalization in topological order Senthil and Fisher (2001), and the
corresponding topological ground state degeneracy can be detected by DMRG
calculations on cylinder or torus.

##  IV Mechanism of superconductivity

To understand the key mechanism of SC order in the $t$-$J$ model, we need to
investigate the effective interactions among holons. We begin with the AFM
phase with $A\neq 0$ and $B,D=0$, and restrict the discussion to the small
doping region ($\delta\ll 1$). After we establish a low energy effective field
theory to understand the interactions among holons in this phase, we shall try
to understand the whole phase diagram later.

###  IV.1 Repulsive holon interaction at large t / J

The first step beyond the mean-field theory is to add the spinon-holon
interaction from Eq. (5), and the total Hamiltonian reads:

$\displaystyle H$ $\displaystyle=H_{h}+H_{b}+H_{\text{int}},$ (38)
$\displaystyle H_{\text{int}}$
$\displaystyle=t\sum_{\langle{ij}\rangle,\sigma}(h^{\dagger}_{j}h_{i}b_{i\sigma}^{\dagger}b_{j\sigma}+h.c.).$
(39)

To study the effective interaction between holons induced by the spinons, we
shall use the coherent state path integral formalism. The spinon (or holon)
operators are replaced by the corresponding complex (or Grassmann) numbers,
with the following notations:

$(b_{i\sigma},b_{i\sigma}^{\dagger})\to(b_{i\sigma},\bar{b}_{i\sigma}),\quad(h_{i},h^{\dagger}_{i})\to(h_{i},\bar{h}_{i}).$

The partition function of fermion-boson coupled theory in the functional field
integral representation is (number terms are dropped)

$\displaystyle Z$ $\displaystyle=\int
D\bar{h}\,Dh\,D\bar{b}\,Db\,e^{-S[\bar{h},h;\bar{b},b]},$ (40) $\displaystyle
S$ $\displaystyle=S_{h}+S_{b}+S_{\text{int}},$ (41) $\displaystyle S_{h}$
$\displaystyle=\int_{0}^{\beta}d\tau\Big{[}\sum_{i}\bar{h}_{i}\partial_{\tau}h_{i}+H_{h}(\tau)\Big{]},$
(42) $\displaystyle S_{b}$
$\displaystyle=\int_{0}^{\beta}d\tau\Big{[}\sum_{i,\sigma}\bar{b}_{i\sigma}\partial_{\tau}b_{i\sigma}+H_{b}(\tau)\Big{]},$
(43) $\displaystyle S_{\text{int}}$
$\displaystyle=\int_{0}^{\beta}d\tau\,H_{\text{int}}(\tau).$ (44)

Now $h$ and $b$ are Grassmann and complex variables, denoting the holon and
spinon degrees of freedom in the coherent state representation. We then
integrate out spinons to derive the effective holon Hamiltonian:

$\displaystyle H^{\text{eff}}_{h}$
$\displaystyle=H_{h}+H^{\text{eff}}_{\text{int}},$ (45) $\displaystyle
H^{\text{eff}}_{\text{int}}$
$\displaystyle=\frac{1}{N}{\sum_{\\{k\\}}}^{\prime}V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}h^{B\dagger}_{k_{1}}h^{A\dagger}_{k_{2}}h^{A}_{k^{\prime}_{1}}h^{B}_{k^{\prime}_{2}},$
(46)

where ${\sum_{\\{k\\}}}^{\prime}$ is summation over all momentum with
conservation ($k_{1}+k_{2}=k^{\prime}_{1}+k^{\prime}_{2}$). In the small
doping limit, $B,D\approx 0$, and the interaction strength is (see Appendix C)

$\displaystyle
V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}=-\frac{t^{2}}{2N}\sum_{p}\bigg{\\{}\frac{1}{\chi_{p}\chi_{p^{\prime}}(\chi_{p}+\chi_{p^{\prime}})}$
$\displaystyle\times\Big{[}J_{b}^{2}(\gamma^{*}_{p}\gamma^{*}_{p^{\prime}}\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}+p^{\prime}}+\gamma_{p}\gamma_{p^{\prime}}\gamma^{*}_{k^{\prime}_{1}+p}\gamma_{k_{2}-p})$
$\displaystyle\quad+(\chi_{p}\chi_{p^{\prime}}-\lambda^{2})(\gamma^{*}_{k^{\prime}_{1}+p}\gamma_{k_{2}+p^{\prime}}+\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}-p})\Big{]}\bigg{\\}},$
(47)

where $p^{\prime}=p+k_{1}-k^{\prime}_{1}$, and $J_{b}=JA$. We have already
take the static limit approximation to extract the frequency-independent
effective interaction, which is mediated by exchanging two spinons with
momentum $p$ and $p^{\prime}$. As the minimum of the spinon dispersion is at
$k=0$, the dominant contributions to the interaction coefficient in the low
energy come from $p\sim p^{\prime}\sim 0$, and we replace
$\gamma_{k^{\prime}_{1}-p}\to\gamma_{k^{\prime}_{1}}$ and
$\gamma_{k_{2}+p^{\prime}}\to\gamma_{k_{2}}$ in Eq. (47), leading to

$\displaystyle
V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}\approx-t^{2}\gamma_{-k^{\prime}_{1}}\gamma_{k_{2}}f(k_{1}-k^{\prime}_{1}),$
(48) $\displaystyle
f(k)\equiv\frac{1}{N}\sum_{p}\left.\frac{J_{b}^{2}\operatorname{Re}(\gamma_{p}\gamma_{p^{\prime}})+(\chi_{p}\chi_{p^{\prime}}-\lambda^{2})}{\chi_{p}\chi_{p^{\prime}}(\chi_{p}+\chi_{p^{\prime}})}\right|_{p^{\prime}=p+k}.$
(49)

In the small doping limit, the holon momenta $k_{1},k^{\prime}_{1}$ are both
expected to be small. Let us then calculate:

$\displaystyle f(0)$
$\displaystyle=\frac{1}{N}\sum_{p}\frac{1}{2\chi_{p}^{3}}[J_{b}^{2}\operatorname{Re}(\gamma_{p}^{2})+(\chi_{p}^{2}-\lambda^{2})]$
$\displaystyle=-\frac{1}{N}\sum_{p}\frac{J_{b}^{2}}{\chi_{p}^{3}}(\operatorname{Im}\gamma_{p})^{2}<0,$
(50)

which is _negative_. Assuming a weak momentum dependence of $f(k)$, we can
approximate $f(k)$ by a _negative_ constant $C/J_{b}$ (with $C<0$), an inverse
Fourier transform of Eq. (46) leads to a simple interaction of holons in real
space:

$\displaystyle H^{\text{eff}}_{\text{int}}$
$\displaystyle=-\frac{t^{2}C}{J_{b}N}{\sum_{\\{k\\}}}^{\prime}\gamma_{-k^{\prime}_{1}}\gamma_{k_{2}}h^{B\dagger}_{k_{1}}h^{A\dagger}_{k_{2}}h^{A}_{k^{\prime}_{1}}h^{B}_{k^{\prime}_{2}}$
$\displaystyle=-\frac{t^{2}C}{J_{b}}\sum_{j\in
B}\sum_{\delta,\delta^{\prime}}h^{\dagger}_{j}h^{\dagger}_{j-\delta^{\prime}}h_{j-\delta}h_{j},$
(51)

where $\delta,\delta^{\prime}$ sum over the nearest neighbor vectors
$\\{\delta_{i}\\}_{i=1}^{3}$ (see Fig. 1). In particular,
$H^{\text{eff}}_{\text{int}}$ contains the following nearest-neighbor
interaction terms (when $\delta=\delta^{\prime}$)

$-\frac{t^{2}C}{J_{b}}\sum_{j\in
B}\sum_{\delta}h^{\dagger}_{j}h^{\dagger}_{j-\delta}h_{j-\delta}h_{j}=-\frac{t^{2}C}{J_{b}}\sum_{\langle{ij}\rangle}h^{\dagger}_{i}h_{i}h^{\dagger}_{j}h_{j}.$
(52)

$C<0$ implies that the nearest neighbor interaction between holons induced by
exchanging spinons is _repulsive_.

###  IV.2 U(1) gauge fluctuations in the t-J Model

Next we impose the no-double-occupancy constraint on each site exactly. This
is done by a delta-function in the partition function:

$\displaystyle\prod_{j}\delta\Big{(}\bar{h}_{j}h_{j}+\sum_{\sigma}\bar{b}_{j\sigma}b_{j\sigma}-1\Big{)}$
$\displaystyle=\int
D\lambda\,\exp\bigg{[}{-i\sum_{j}\lambda_{j}\Big{(}\bar{h}_{j}h_{j}+\sum_{\sigma}\bar{b}_{j\sigma}b_{j\sigma}-1\Big{)}}\bigg{]},$
(53)

where $\\{\lambda_{j}\\}$ are real Lagrange multipliers, and
$D\lambda\equiv\prod_{j}d\lambda_{j}$. With the Hamiltonian Eq. (5), the
coherent state path integral representation of the partition function is:

$\displaystyle Z=\int
D\bar{b}\,Db\,D\bar{h}\,Dh\,D\lambda\,e^{-\int_{0}^{\beta}d\tau\,L},$
$\displaystyle
L=\sum_{j}\Big{[}\bar{h}_{j}(\partial_{\tau}+i\lambda_{j})h_{j}+\sum_{\sigma}\bar{b}_{j\sigma}(\partial_{\tau}+i\lambda_{j})b_{j\sigma}-i\lambda_{j}\Big{]}$
$\displaystyle+t\sum_{\langle{ij}\rangle,\sigma}(\bar{h}_{j}h_{i}\bar{b}_{i\sigma}b_{j\sigma}+h.c.)-\frac{J}{2}\sum_{\begin{subarray}{c}\langle{ij}\rangle\\\
\sigma,\sigma^{\prime}\end{subarray}}\sigma\sigma^{\prime}\bar{b}_{j,-\sigma}\bar{b}_{i\sigma}b_{i\sigma^{\prime}}b_{j,-\sigma^{\prime}}.$

The $J$-term is then decoupled by a Hubbard-Stratonovich transformation in the
AFM channel (the FM channel is unimportant in the small-doping region).
Introducing auxiliary complex fields $\Delta_{ij}$ ($i\in A$, $j\in B$) and
their conjugate $\bar{\Delta}_{ij}$ on each bond, we arrive at

$\displaystyle Z$ $\displaystyle=\int
D\bar{b}\,Db\,D\bar{h}\,Dh\,D\bar{\Delta}\,D\Delta\,D\lambda\,e^{-\int_{0}^{\beta}d\tau\,L},$
(54) $\displaystyle L$
$\displaystyle=\sum_{j}\Big{[}\bar{h}_{j}(\partial_{\tau}+i\lambda_{j})h_{j}+\sum_{\sigma}\bar{b}_{j\sigma}(\partial_{\tau}+i\lambda_{j})b_{j\sigma}\Big{]}$
$\displaystyle\quad-i\sum_{j}\lambda_{j}+t\sum_{\langle{ij}\rangle,\sigma}(\bar{h}_{j}h_{i}\bar{b}_{i\sigma}b_{j\sigma}+h.c.)$
$\displaystyle\quad+\frac{J}{2}\sum_{\langle{ij}\rangle}\Big{[}|\Delta_{ij}|^{2}-\sum_{\sigma}\sigma(\bar{\Delta}_{ij}b_{i\sigma}b_{j,-\sigma}+h.c.)\Big{]},$
(55)

where $D\Delta\equiv\prod_{\langle{ij}\rangle}d\Delta_{ij}$ (with $i\in A$),
and $|\Delta_{ij}|^{2}\equiv\bar{\Delta}_{ij}\Delta_{ij}$. The Lagrangian Eq.
(55) is invariant (up to a total $\tau$-derivative term) under the staggered
$U(1)$ gauge transformation:

$\begin{aligned} b_{i\in s,\sigma}&\to b_{i\sigma}e^{-is\theta_{i}},\\\
h_{i\in s}&\to h_{i}e^{-is\theta_{i}},\end{aligned}\quad\begin{aligned}
\Delta_{ij}&\to\Delta_{ij}e^{-i(\theta_{i}-\theta_{j})},\\\ \lambda_{i\in
s}&\to\lambda_{i}+s\,\partial_{\tau}\theta_{i}.\end{aligned}$ (56)

The transformation can be reformulated using $U(1)$ gauge fields: define the
vector potential $\mathbf{A}$ (on bonds) and the scalar potential $A_{\tau}$
(on sites)

$\theta_{i}-\theta_{j}=x_{ij}\cdot\mathbf{A}(i,j),\quad
A_{\tau}(i)=\partial_{\tau}\theta_{i}.$ (57)

Here $x_{ij}$ is the vector from site $j$ to site $i$. The gauge
transformation of $\Delta,\lambda$ are then

$\displaystyle\Delta_{ij}$
$\displaystyle\to\Delta_{ij}e^{-ix_{ij}\cdot\mathbf{A}(i,j)},$ (58)
$\displaystyle\lambda_{i\in s}$ $\displaystyle\to\lambda_{i}+sA_{\tau}(i).$

In the partition function, we keep configurations of $\Delta_{ij}$,
$\bar{\Delta}_{ij}$, $\lambda_{i}$ at the saddle point with fluctuations of
the gauge field:

$\displaystyle\Delta_{i,i+\delta}$ $\displaystyle\to\Delta
e^{i\delta\cdot\mathbf{A}(i,i+\delta)}$ $\displaystyle(i\in A),$ (59)
$\displaystyle\bar{\Delta}_{i,i+\delta}$ $\displaystyle\to\Delta
e^{-i\delta\cdot\mathbf{A}(i,i+\delta)}$ $\displaystyle(i\in A),$
$\displaystyle\lambda_{i\in s}$ $\displaystyle\to-i\lambda+sA_{\tau}(i),$

where $\Delta=2A$ and $\lambda$ are real solutions from the mean-field theory.
The fluctuation of the norm of $\Delta_{ij}$ is ignored. Then

$\displaystyle Z$ $\displaystyle=\int
D\bar{b}\,Db\,D\bar{h}\,Dh\,DA\,e^{-\int_{0}^{\beta}d\tau\,L},$ (60)
$\displaystyle L$ $\displaystyle=L_{h}+L_{\text{heis}}+L_{t}+L_{0},$ (61)

with each part in the Lagrangian given by

$\displaystyle\ \ \ L_{h}=\sum_{s}\sum_{i\in
s}\Big{[}\bar{h}_{i}(\partial_{\tau}+isA_{\tau}(i))h_{i}+\lambda\bar{h}_{i}h_{i}\Big{]},$
(62) $\displaystyle L_{\text{heis}}=\sum_{s}\sum_{i\in
s}\Big{[}\bar{b}_{i\sigma}(\partial_{\tau}+isA_{\tau}(i))b_{i\sigma}+\lambda\bar{b}_{i\sigma}b_{i\sigma}\Big{]}$
$\displaystyle-\frac{J\Delta}{2}\sum_{i\in
A}\sum_{\delta,\sigma}\sigma\Big{[}e^{i\delta\cdot\mathbf{A}(i,i+\delta)}\bar{b}_{i\sigma}\bar{b}_{i+\delta,-\sigma}+h.c.\Big{]},$
(63) $\displaystyle\ \ \ \ L_{t}=t\sum_{i\in
A}\sum_{\delta,\sigma}\bar{h}_{i+\delta}h_{i}\bar{b}_{i\sigma}b_{i+\delta,\sigma}+h.c.,$
(64) $\displaystyle\ \ \ \ L_{0}=-\sum_{s}\sum_{i\in
s}(\lambda+isA_{\tau}(i))+\frac{\alpha JN}{2}\Delta^{2}.$ (65)

Here $s$ sums over the two sub-lattices $A,B$. We also used the fact that the
nearest neighbors of a site $i\in A$ are at positions
$\\{i+\delta_{l}\\}_{l=1}^{3}$ on the honeycomb lattice.

We then take the continuum limit by setting the nearest neighbor distance
$a\to 0$ and system size $\to\infty$. A summation over sites in one sub-
lattice is replaced by an integral over space:

$\sum_{i\in A}\text{ or }\sum_{j\in
B}\to\int\frac{d^{2}x}{2\Omega},\quad\Omega=\frac{3\sqrt{3}}{4}a^{2}.$ (66)

Here $a$ is the nearest neighbor distance, and $2\Omega$ is the area of a unit
cell. The continuum fields are

$\displaystyle h_{i\in s}$ $\displaystyle\to h^{s}(x),$ $\displaystyle b_{i\in
s,\sigma}$ $\displaystyle\to b^{s}_{\sigma}(x),$ (67) $\displaystyle
A_{\tau}(i)$ $\displaystyle\to A_{\tau}(x),$
$\displaystyle\mathbf{A}(i,i+\delta)$
$\displaystyle\to\mathbf{A}(x+\tfrac{\delta}{2}).$

where $x$ is the position of the site $i$. For convenience, we redefine
spinon-holon coupling $\alpha t\to t$, and introduce $Q\equiv\alpha
J\Delta/8$. The continuum Lagrangian is (see details in Appendix D)

$\displaystyle
L_{h}=\int\frac{d^{2}x}{2\Omega}\sum_{s}\Big{[}\bar{h}^{s}(\partial_{\tau}+isA_{\tau})h^{s}+\lambda\bar{h}^{s}h^{s}\Big{]},$
(68) $\displaystyle
L_{\text{heis}}=\int\frac{d^{2}x}{2\Omega}\sum_{s,\sigma}\Big{[}\bar{b}^{s}_{\sigma}(\partial_{\tau}+isA_{\tau})b^{s}_{\sigma}+\lambda\bar{b}^{s}_{\sigma}b^{s}_{\sigma}\Big{]}$
$\displaystyle\
+Qa^{2}\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\sigma\Big{[}(\nabla-i\mathbf{A})\bar{b}^{A}_{\sigma}\cdot(\nabla+i\mathbf{A})\bar{b}^{B}_{-\sigma}+h.c.\Big{]}$
$\displaystyle\
-4Q\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\sigma\Big{[}\bar{b}^{A}_{\sigma}\bar{b}^{B}_{-\sigma}+h.c.\Big{]},$
(69) $\displaystyle
L_{t}=t\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\Big{[}\bar{h}^{B}h^{A}\bar{b}^{A}_{\sigma}b^{B}_{\sigma}+h.c.\Big{]},$
(70) $\displaystyle
L_{0}=\int\frac{d^{2}x}{2\Omega}\Big{[}-2\lambda+\frac{\alpha
J}{2}\Delta^{2}\Big{]}.$ (71)

Define $A_{\mu}=(A_{\tau},\mathbf{A})$ and
$\partial_{\mu}=(\partial_{\tau},\nabla)$; the gauge transformation in
continuum limit is

$\displaystyle b^{s}_{\sigma}(x)$ $\displaystyle\to
b^{s}_{\sigma}(x)e^{-is\theta(x)},$ (72) $\displaystyle h^{s}(x)$
$\displaystyle\to h^{s}(x)e^{-is\theta(x)},$ $\displaystyle A_{\mu}(x)$
$\displaystyle\to A_{\mu}(x)+\partial_{\mu}\theta(x).$

###  IV.3 Effective theory of holons

Following Ref. 27, we combine the Schwinger bosons into a low energy field $z$
and a high energy field $\pi$:

$z_{\sigma}=(b^{A}_{\sigma}+\sigma\bar{b}^{B}_{-\sigma})/2,\quad\pi_{\sigma}=(b^{A}_{\sigma}-\sigma\bar{b}^{B}_{-\sigma})/2.$
(73)

From Eq. (72), the gauge transformation of $z,\pi$ fields are

$z_{\sigma}\to
z_{\sigma}e^{-i\theta},\quad\pi_{\sigma}\to\pi_{\sigma}e^{-i\theta}.$ (74)

We then express the Lagrangian $L$ in terms of $z,\pi$ fields. The transformed
$L_{\text{heis}}$ (Eq. (69)) is

$\displaystyle
L_{\text{heis}}=\int\frac{d^{2}x}{\Omega}\sum_{\sigma}\bar{z}_{\sigma}[(\lambda-4Q)-Qa^{2}(\nabla+i\mathbf{A})^{2}]z_{\sigma}$
$\displaystyle+\int\frac{d^{2}x}{\Omega}\sum_{\sigma}[\bar{\pi}_{\sigma}(\partial_{\tau}+iA_{\tau})z_{\sigma}-\pi_{\sigma}(\partial_{\tau}-iA_{\tau})\bar{z}_{\sigma}]$
$\displaystyle+\int\frac{d^{2}x}{\Omega}\sum_{\sigma}\bar{\pi}_{\sigma}[(\lambda+4Q)+Qa^{2}(\nabla+i\mathbf{A})^{2}]\pi_{\sigma},$
(75)

and the spinon-holon interaction $L_{t}$ (Eq. (70)) becomes

$L_{t}=(-t)\int\frac{d^{2}x}{\Omega}\sum_{\sigma}\sigma(\bar{\pi}_{\sigma}\bar{h}^{B}h^{A}\bar{z}_{-\sigma}+h.c.).$
(76)

To derive the effective theory of holons, we first integrate over the high-
energy field $\pi$. The low-energy effective Lagrangian of $h,z,A$ fields is

$\displaystyle L_{hzA}=L_{0}+L_{h}$
$\displaystyle+\int\frac{d^{2}x}{\Omega}\sum_{\sigma}\Big{\\{}(\lambda-4Q)|z_{\sigma}|^{2}+Qa^{2}|(\nabla+i\mathbf{A})z_{\sigma}|^{2}$
$\displaystyle\quad+\frac{1}{\lambda+4Q}\big{[}|(\partial_{\tau}+iA_{\tau})z_{\sigma}|^{2}-t^{2}\bar{h}^{A}h^{B}\bar{h}^{B}h^{A}z_{\sigma}\bar{z}_{\sigma}$
$\displaystyle\qquad+t\sigma(\bar{h}^{A}h^{B}z_{-\sigma}\partial_{\tau}z_{\sigma}-\bar{h}^{B}h^{A}\bar{z}_{-\sigma}\partial_{\tau}\bar{z}_{\sigma})\big{]}\Big{\\}}.$
(77)

After scaling $\tau\to\tau/c$ (therefore $A_{\tau}\to cA_{\tau}$) and defining

$\displaystyle\Gamma$ $\displaystyle=\sqrt{\lambda^{2}-16Q^{2}},$
$\displaystyle m$ $\displaystyle=\Gamma/c,$ (78) $\displaystyle c$
$\displaystyle=\sqrt{Q(\lambda+4Q)}a,$ $\displaystyle g^{-1}$
$\displaystyle=Qa^{2}/(c\Omega),$

we obtain the effective action

$\displaystyle S_{hzA}$ $\displaystyle=S_{h}+S_{zA}+S_{t},\quad
S_{t}=S_{t}^{1}+S_{t}^{2},$ (79) $\displaystyle S_{h}$
$\displaystyle=\frac{1}{\Omega}\int
d^{3}x\sum_{s}\Big{[}\bar{h}^{s}(\partial_{\tau}+isA_{\tau})h^{s}+\frac{\lambda}{c}\bar{h}^{s}h^{s}\Big{]},$
(80) $\displaystyle S_{zA}$ $\displaystyle=\frac{1}{g}\int
d^{3}x\sum_{\sigma}\Big{[}m^{2}|z_{\sigma}|^{2}+|(\partial_{\mu}+iA_{\mu})z_{\sigma}|^{2}\Big{]},$
(81) $\displaystyle S_{t}$ $\displaystyle=\frac{t}{cg}\int
d^{3}x\sum_{\sigma}\sigma(\bar{h}^{A}h^{B}z_{-\sigma}\partial_{\tau}z_{\sigma}-\bar{h}^{B}h^{A}\bar{z}_{-\sigma}\partial_{\tau}\bar{z}_{\sigma})$
$\displaystyle\quad+\frac{t^{2}}{c^{2}g}\int
d^{3}x\,\bar{h}^{A}h^{A}\bar{h}^{B}h^{B}\sum_{\sigma}\bar{z}_{\sigma}z_{\sigma}.$
(82)

Here $\int d^{3}x\equiv\int_{0}^{c\beta}d\tau\int d^{2}x$, and constant terms
from $L_{0}$ are dropped. Note that at half-filling, the holons are
eliminated, and the effective action reduces to Eq. (81) only, which is the
$\mathrm{\mathbb{C}P^{1}}$ model Read and Sachdev (1989). Let us re-express
$z$ fields as

$\displaystyle z_{\sigma}(\tau,x)$
$\displaystyle=\sqrt{\rho_{\sigma}(\tau,x)}\exp[i\phi_{\sigma}(\tau,x)],$ (83)
$\displaystyle\bar{z}_{\sigma}(\tau,x)$
$\displaystyle=\sqrt{\rho_{\sigma}(\tau,x)}\exp[-i\phi_{\sigma}(\tau,x)].$

At zero temperature, the spinon (and thus the $z$-field) condenses. At
sufficiently low temperature, we can replace both $\rho_{\uparrow}$ and
$\rho_{\downarrow}$ by a number $\rho_{0}/2$, where $\rho_{0}$ is of the same
order as the density $n_{0}$ of condensed spinons; the fluctuations of the
$\rho_{\sigma}$ fields will be ignored. Then

$\displaystyle\sum_{\sigma}|(\partial_{\mu}+iA_{\mu})z_{\sigma}|^{2}$
$\displaystyle=\sum_{\sigma}\frac{1}{4\rho_{\sigma}}\Big{[}(\partial_{\mu}\rho_{\sigma})^{2}+4\rho_{\sigma}^{2}(A_{\mu}+\partial_{\mu}\phi_{\sigma})^{2}\Big{]}$
$\displaystyle\approx\sum_{\sigma}\rho_{\sigma}A_{\mu}^{2}\approx\rho_{0}A_{\mu}^{2}.$
(84)

In the last line, we dropped the derivative $\partial_{\mu}\rho_{\sigma}$. By
a gauge transformation, the derivative of $\phi$ is absorbed into $A_{\mu}$,
and $A_{\mu}$ acquires a mass $\sqrt{\rho_{0}}$ via the Anderson-Higgs
mechanism. For analysis near the mean-field saddle point, we can assume a weak
$\tau$-dependence of $z$ and drop the $\tau$-derivatives of $z$ in Eq. (82).
Then

$\displaystyle S_{zA}$ $\displaystyle\approx\frac{\rho_{0}}{g}\int
d^{3}x\,(m^{2}+A_{\mu}^{2}),$ (85) $\displaystyle S_{t}$
$\displaystyle\approx\frac{t^{2}\rho_{0}}{c^{2}g}\int
d^{3}x\,\bar{h}^{A}h^{A}\bar{h}^{B}h^{B},$ (86)

yielding an effecive action independent of $\phi$. The integration over $\phi$
simply gives an (infinite) constant factor, and is dropped thereafter (this
removes the gauge redundancy). We reorganize terms in $S_{hzA}$ according to
the order of $A$:

$\displaystyle S_{hzA}$ $\displaystyle=S_{0}+S_{1}+S_{2},$ (87) $\displaystyle
S_{0}$ $\displaystyle=S_{t}+\int d^{3}x\,m^{2}g^{-1}\rho_{0}$
$\displaystyle\quad+\Omega^{-1}\int
d^{3}x\sum_{s}\Big{(}\bar{h}^{s}\partial_{\tau}h^{s}+\frac{\lambda}{c}n^{h}_{s}\Big{)},$
(88) $\displaystyle S_{1}$ $\displaystyle=i\Omega^{-1}\int
d^{3}x\,A_{\tau}(n^{h}_{A}-n^{h}_{B}),$ (89) $\displaystyle S_{2}$
$\displaystyle=g^{-1}\int d^{3}x\,\rho_{0}A_{\mu}^{2}.$ (90)

Here $n^{h}_{s}\equiv\bar{h}^{s}h^{s}$ ($s=A,B$). Finally we integrate over
$A$:

$\displaystyle Z_{hzA}\cong\int D\bar{h}\,Dh\,DA\,e^{-S_{hzA}}$
$\displaystyle=\int D\bar{h}\,Dh\,e^{-S_{0}}\underbrace{\int
DA\,e^{-(S_{1}+S_{2})}}_{\equiv Z^{\prime}}$ $\displaystyle
Z^{\prime}\equiv\underbrace{\int D\mathbf{A}\exp\bigg{\\{}\int
d^{3}x\,(-g^{-1}\rho_{0}\mathbf{A}^{2})\bigg{\\}}}_{\text{a constant factor}}$
$\displaystyle\times\int DA_{\tau}\exp\bigg{\\{}\int
d^{3}x\,\Big{[}{-\frac{1}{g}\rho_{0}A_{\tau}^{2}-\frac{i}{\Omega}(n^{h}_{A}-n^{h}_{B})A_{\tau}}\Big{]}\bigg{\\}}$
$\displaystyle\propto\exp\bigg{\\{}\int
d^{3}x\,\frac{g}{4\rho_{0}}\Big{[}\frac{i}{\Omega}(n^{h}_{A}-n^{h}_{B})\Big{]}^{2}\bigg{\\}}$
$\displaystyle=\exp\bigg{\\{}-\int
d^{3}x\,\frac{g}{4\rho_{0}\Omega^{2}}(n^{h}_{A}-n^{h}_{B})^{2}\bigg{\\}},$

Then the effective action of the holons is:

$S_{\text{eff}}=S_{0}+\int
d^{3}x\,\frac{g}{4\rho_{0}\Omega^{2}}(n^{h}_{A}-n^{h}_{B})^{2}.$ (91)

We identify terms representing the interaction between $A$-holons and
$B$-holons (recall we redefined $\alpha t\to t$ earlier, now we restore the
original $t$):

$H_{\text{int}}^{AB}=\bigg{[}\frac{\alpha^{2}t^{2}\rho_{0}}{c^{2}g}-\frac{g}{2\rho_{0}\Omega^{2}}\bigg{]}\int
d^{2}x\,n^{h}_{A}n^{h}_{B}.$ (92)

The first term (inherited from Eq. (86)) is repulsion, produced by the
exchange of two spinons between holons, in agreement with the previous
calculation (in Sec. IV.1) on a discrete lattice; the second term is
attraction, produced by the gauge fluctuation.

The most important implication of Eq. (92) is that there exists a critical
$t/J$ below which the attraction due to gauge fluctuation should overcome the
repulsion due to spinon exchange, leading to a superconducting order. For a
large value of $t/J$, the repulsion dominates, hence the absence of SC order.
This critical value of $t/J$ is estimated from

$\frac{\alpha^{2}t^{2}\rho_{0}}{c^{2}g}\lesssim\frac{g}{2\rho_{0}\Omega^{2}}.$
(93)

Note that $\lambda\simeq 4Q$ at low temperature and small doping. Thus

$\frac{t}{J}\lesssim\frac{\Delta}{\sqrt{2}\rho_{0}}.$ (94)

From the mean-field calculation at half-filling and zero temperature (see
Appendix B.2), we obtain $\Delta\approx 1.2$, and $\rho_{0}\sim n_{0}\approx
0.5$, leading to $t/J\lesssim 1.7$. However, from the renormalization group
perspective (by integrating over high-energy holons) Nagaosa (1999), the
effective attractive between holons will renormalize to a bigger value at low
energy. Thus the actual range of $t/J$ that allows superconductivity is much
larger than our naive estimation. The net attractive interaction leads to
($p+ip$)-wave pairing of the spinless holons Cheng  _et al._ (2010). Together
with the singlet-pairing of the spinons, we obtain the ($d+id$)-wave pairing
of the electrons.

For the $A,B,D\neq 0$ region, the order parameter $B$ can be regarded as an
additional Higgs field which further breaks down the staggered $U(1)$ gauge
field into a $\mathbb{Z}_{2}$ gauge field, and we believe that such an
additional gauge field will not affect the attractive interaction in the
presence of AFM long-range order. In fact, a nonzero B is also crucial for the
holon mobility, e.g., for the introduction of holon kinetic term, and
superconductivity is only possible in this phase. Nevertheless, at very small
doping, our estimation for the effective holon interactions in previous
sections is still valid since the correction induced by B is always
proportional to holon concentration $\delta$. Physically, as long as we start
from the short-range RVB phase $A\neq 0$, the staggered $U(1)$ gauge
fluctuations will always induce attractive interactions between holons
belonging to different sub-lattices, and this naturally explains why $d+id$ SC
order still survive even in the absence of AFM long-range order in this phase.
We stress that the spin-charge separation scenario plays an essential role for
the emergence of SC order here, especially for the region without AFM order,
and there is no other conventional picture can explain such a novel $d+id$ Sc
oder emerged in a doped-Mott insulator!

Finally, for the region $A=0$ and $B,D\neq 0$, the usual uniform $U(1)$ gauge
fluctuations will kill the FM long-range order at finite $t/J$ and induce
repulsive interactions among holons. Such a deconfined phase of $U(1)$ gauge
field does not support superconductivity and would be a natural candidate for
the non-Fermi liquid phase emerged in the large $t/J$ limit.

##  V Conclusion and discussion

In conclusion, we investigate the global phase diagram of the $t$-$J$ model on
a honeycomb lattice at small doping region using the Grassmann tensor network
numerical method. Furthermore, the slave-fermion mean-field theory is employed
to account for the global phase diagram. The effective interacting holon
theory is rigorously derived in the presence of the AFM background. The novel
phenomena of spin-charge separation and attractive interaction among holons
naturally explain the emergence of $d+id$ SC orders. Remarkably, the predicted
phase boundary for $d+id$ SC orders via self-consistent mean-field theory
approach and effective field theory analysis is intrinsically close to the
results from Grassmann tensor product numerical simulations. We stress that
the spin-charge separation mechanism occurs only in strongly correlated
systems and is essentially different from the weak coupling mechanism in BCS
theory.

We also would like to point out that the attractive interaction between holons
belonging to different sub-lattices is induced by gauge fluctuations, and the
_staggered_ $U(1)$ gauge structure is crucial to the attractive interaction,
since in this case the two sub-lattices $A$ and $B$ carry opposite gauge
charges. Although our derivation is in the region for $A\neq 0$ and $B,D=0$,
we believe that this mechanism is still valid in the region $A,B,D\neq 0$, and
it explains the emergence of SC orders even in the absence of long-range AFM
order, as long as the system is still dominated by short-range AFM spin
fluctuations.

It is well known that in the conventional BCS theory, the existence of zero-
energy bound fermion state in the vortex of superconductor depends on the
nontrivial topological class of the bulk. According to the topological
classification, the $d+id$ superconductivity discovered on the honeycomb
lattice breaks the time reversal symmetry, but has the $SU(2)$ spin-rotation
symmetries. Thus it belongs to the class $C$ of ten Altland-Zirnbauer classes
Chiu  _et al._ (2016); Sato and Ando (2017). Althugh the class $C$ has no
zero-energy states in the vortex from the perspective of noninteracting
topological classification, we argue that the unconventional superconductivity
in the $t$-$J$ model on a honeycomb lattice emerges from the spin-charge
separation is a different case. The bosonic spinons carrying the spin degrees
of freedom condense in the AFM region. The $p+ip$ superconductivity of the
spinless holons belongs to the class $D$ and there will be a zero-energy bound
state in the vortex, similar to the Majorana zero mode in conventional $p+ip$
superconductor Kopnin and Salomaa (1991); Volovik (1999); Read and Green
(2000), which serves as a solid evidence for the spin-charge separation
scenario. Experimentally, the doped spin $1/2$ honeycomb lattice
antiferromagnet $\mathrm{InV_{1/3}Cu_{2/3}O_{3}}$ Möller  _et al._ (2008)
would be an appealing candidate to examine the emergence of topological
Majorana zero mode in its vortex core.

## Acknowledgment

This work is supported by General Research Fund Grant No. 14302021 and
NSFC/RGC Joint Research Scheme No. N-CUHK427/18 from Research Grants Council,
and Direct Grant No. 4053416 from the Chinese University of Hong Kong. WQC is
supported by the National Key R&D Program of China (Grants No.
2022YFA1403700), NSFC (Grants No. 11861161001), the Science, Technology and
Innovation Commission of Shenzhen Municipality (No. ZDSYS20190902092905285),
Guangdong Basic and Applied Basic Research Foundation under Grant No.
2020B1515120100, and Center for Computational Science and Engineering at
Southern University of Science and Technology.

##  Appendix A Details of the Grassmann tensor numerical calculation

Figure 6: Staggered magnetization, amplitudes of singlet and triplet SC order
parameters versus doping at $t/J=3,5,10$.

Figure 7: Staggered magnetization, amplitudes of singlet and triplet SC order
parameters versus doping at $t/J=15,25,30$.

We use state-of-the-art Grassmann tensor product numerical method to obtain
data in Sec. II. The standard form of GTPS is used as our variational wave
function. Translation-invariant ansatz is assumed and is specified by just two
different Grassmann tensors $\mathbf{T}_{A}$ and $\mathbf{T}_{B}$ on sub-
lattices A and B of each unit cell:

$\displaystyle\Psi({m_{i}},{m_{j}})=\text{tTr}\int\Pi_{\langle
ij\rangle}\mathbf{g}_{aa^{\prime}}\Pi_{i\in
A}\mathbf{T}^{m_{i}}_{A;abc}\Pi_{j\in
B}T^{m_{j}}_{B;a^{\prime}b^{\prime}c^{\prime}}.$ (95)

with

$\displaystyle{\textbf{T}}^{m_{i}}_{A;abc}$ $\displaystyle=$
$\displaystyle{T}^{m_{i}}_{A;abc}\theta_{\alpha}^{P^{f}(a)}\theta_{\beta}^{P^{f}(b)}\theta_{\gamma}^{P^{f}(c)},$
$\displaystyle{\textbf{T}}^{m_{j}}_{B,a^{\prime}b^{\prime}c^{\prime}}$
$\displaystyle=$
$\displaystyle{T}^{m_{j}}_{B;a^{\prime}b^{\prime}c^{\prime}}\theta_{\alpha^{\prime}}^{P^{f}(a^{\prime})}\theta_{\beta^{\prime}}^{P^{f}(b^{\prime})}\theta_{\gamma^{\prime}}^{P^{f}(c^{\prime})},$
$\displaystyle\textbf{g}_{aa^{\prime}}$ $\displaystyle=$
$\displaystyle\delta_{aa^{\prime}}{\mathrm{d}\theta}_{\alpha}^{P^{f}(a)}{\mathrm{d}\theta}_{\alpha^{\prime}}^{P^{f}(a^{\prime})}.$
(96)

More details about this ansatz can be found in the appendix of Ref. 16. The
variational ground state is obtained from the simple update imaginary time
evolution algorithm Gu (2013). The desired doping is achieved by adjusting
chemical potential. Then we measure the physical quantities by the
renormalization group algorithms (GTERG/wGTERG) in Refs. 37; 36. The total
system size is up to $2\times 3^{6}$ sites and all calculations are performed
with periodic boundary conditions. We use three different virtual bond
dimensions $D=10,12,14$ of the GTPS. Different dopings are obtained by tuning
the chemical potential.

The data of $D=\infty$ is obtained by the following approximate extrapolation
method. For each fixed $D$, we get a series of data
$(\delta_{i}^{D},O_{i}^{D})$, where $\delta_{i}^{D}$ is the doping and
$O_{i}^{D}$ is the physical quantity that we are interested in. We choose some
equally spaced doping points $\tilde{\delta}_{i}^{D}$. For each doping point
$\tilde{\delta}_{i}$ and each $D$, we use the linear interpolation to get the
approximated $\tilde{O}_{i}^{D}$ corresponding to $\tilde{\delta}_{i}$ using
two nearby data points we have, ($\delta_{i,L}^{D}$,$O_{i,L}^{D}$) and
($\delta_{i,R}^{D}$,$O_{i,R}^{D}$),

$\frac{\tilde{O}_{i}^{D}-O_{i,L}^{D}}{\tilde{\delta}_{i}-\delta_{i,L}^{D}}=\frac{O_{i,R}^{D}-O_{i,L}^{D}}{\delta_{i,R}^{D}-\delta_{i,L}^{D}}.$
(97)

We then use least squares linear fit using data points
$(D^{-1},\tilde{O}_{i}^{D})$ and then find $\tilde{O}_{i}^{\infty}$ by setting
$D^{-1}=0$. The linear fit $O=k_{i}D^{-1}+b_{i}$ minimizes the squared error

$E_{i}=\sum_{D}|k_{i}D^{-1}+b_{i}-\tilde{O}_{i}^{D}|^{2}.$ (98)

To reduce fluctuation, we use only points with good fitting (R-squared value
greater than 0.7). We also set negative $O_{i}^{\infty}$ as 0 due to physical
consideration. $O_{i}^{\infty}$ at zero doping is set to be 0 for
superconducting order parameter and $O_{i}^{14}$ of smallest doping we have.
Here we provide figures of other values of $t/J$ (Figs. 6, 7) below in
addition to the $t/J=20$ results in the main text (Fig. 3).

##  Appendix B Details of the mean-field theory

###  B.1 Mean-field Hamiltonian of the t-J model

In this appendix, we derive the mean-field Hamiltonian Eqs. (10) to (13). The
mean-field decoupling is done by ignoring higher order fluctuations around the
mean-field average Bruus and Flensberg (2004): for any two operators $A$ and
$B$,

$AB\approx A\langle{B}\rangle+\langle{A}\rangle
B-\langle{A}\rangle\langle{B}\rangle.$ (99)

The $t$-terms in Eq. (5) are decoupled to

$\displaystyle H_{t}$ $\displaystyle\equiv
2t\sum_{\langle{ij}\rangle}(h^{\dagger}_{i}h_{j}\hat{B}^{\dagger}_{ij}+h.c.)$
$\displaystyle\approx
2t\sum_{\langle{ij}\rangle}(D_{ij}\hat{B}^{\dagger}_{ij}+h^{\dagger}_{i}h_{j}B^{*}_{ij}+h.c.)$
$\displaystyle\quad-2t\sum_{\langle{ij}\rangle}(D_{ij}B^{*}_{ij}+h.c.).$ (100)

By Wick’s theorem, the decoupling of $J$-terms in Eq. (5)

$H_{J}=-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}b^{\dagger}_{j,-\sigma}b^{\dagger}_{i\sigma}b_{i\sigma^{\prime}}b_{j,-\sigma^{\prime}}$

involves contribution from 3 channels:

$\displaystyle H_{J}$ $\displaystyle\approx
H_{J}^{(A)}+H_{J}^{(B)}+H_{J}^{(\delta)},$ $\displaystyle H_{J}^{(A)}$
$\displaystyle=-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}\Big{[}b^{\dagger}_{j,-\sigma}b^{\dagger}_{i\sigma}\langle{b_{i\sigma^{\prime}}b_{j,-\sigma^{\prime}}}\rangle+\langle{b^{\dagger}_{j,-\sigma}b^{\dagger}_{i\sigma}}\rangle
b_{i\sigma^{\prime}}b_{j,-\sigma^{\prime}}-\langle{b^{\dagger}_{j,-\sigma}b^{\dagger}_{i\sigma}}\rangle\langle{b_{i\sigma^{\prime}}b_{j,-\sigma^{\prime}}}\rangle\Big{]},$
$\displaystyle H_{J}^{(B)}$
$\displaystyle=-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}\Big{[}b^{\dagger}_{j,-\sigma}b_{i\sigma^{\prime}}\langle{b^{\dagger}_{i\sigma}b_{j,-\sigma^{\prime}}}\rangle+\langle{b^{\dagger}_{j,-\sigma}b_{i\sigma^{\prime}}}\rangle
b^{\dagger}_{i\sigma}b_{j,-\sigma^{\prime}}-\langle{b^{\dagger}_{j,-\sigma}b_{i\sigma^{\prime}}}\rangle\langle{b^{\dagger}_{i\sigma}b_{j,-\sigma^{\prime}}}\rangle\Big{]},$
$\displaystyle H_{J}^{(n)}$
$\displaystyle=-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}\Big{[}b^{\dagger}_{j,-\sigma}b_{j,-\sigma^{\prime}}\langle{b^{\dagger}_{i\sigma}b_{i\sigma^{\prime}}}\rangle+\langle{b^{\dagger}_{j,-\sigma}b_{j,-\sigma^{\prime}}}\rangle
b^{\dagger}_{i\sigma}b_{i\sigma^{\prime}}-\langle{b^{\dagger}_{j,-\sigma}b_{j,-\sigma^{\prime}}}\rangle\langle{b^{\dagger}_{i\sigma}b_{i\sigma^{\prime}}}\rangle\Big{]}.$

With the mean-field ansatz Eqs. (8) and (9), we calculate the $A$-channel
terms:

$\displaystyle H_{J}^{(A)}$
$\displaystyle=-2J\sum_{\langle{ij}\rangle}(A^{*}_{ij}\hat{A}_{ij}+h.c.-|A_{ij}|^{2}),$
(101)

the $B$-channel terms:

$\displaystyle H_{J}^{(B)}$
$\displaystyle=-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}\Big{[}B^{*}_{ij}b^{\dagger}_{i\sigma}b_{j,-\sigma^{\prime}}+h.c.-|B_{ij}|^{2}\Big{]}\delta_{\sigma,-\sigma^{\prime}}$
$\displaystyle=\frac{J}{2}\sum_{\langle{ij}\rangle}\Big{[}B^{*}_{ij}\sum_{\sigma}b^{\dagger}_{i\sigma}b_{j\sigma}+h.c.-2|B_{ij}|^{2}\Big{]}=J\sum_{\langle{ij}\rangle}(B^{*}_{ij}\hat{B}_{ij}+h.c.-|B_{ij}|^{2}),$
(102)

and the $n$-channel terms:

$\displaystyle H_{J}^{(n)}$
$\displaystyle=-\frac{J}{2}\sum_{\langle{ij}\rangle}\sum_{\sigma,\sigma^{\prime}}\sigma\sigma^{\prime}\Big{[}\frac{1-\delta}{2}(b^{\dagger}_{j,-\sigma}b_{j,-\sigma}+b^{\dagger}_{i\sigma}b_{i\sigma})-\frac{(1-\delta)^{2}}{4}\Big{]}\delta_{\sigma\sigma^{\prime}}$
$\displaystyle=-\frac{J}{2}\sum_{\langle{ij}\rangle}\Big{[}\frac{1-\delta}{2}(\hat{n}^{b}_{i}+\hat{n}^{b}_{j})-\frac{(1-\delta)^{2}}{2}\Big{]}=-\frac{1}{4}\alpha
J(1-\delta)\sum_{i}\hat{n}^{b}_{i}+\frac{1}{4}\alpha NJ(1-\delta)^{2}.$ (103)

Here $\alpha=3$ is the coordination number. Finally, the mean-field
Hamiltonian is

$\displaystyle H_{\text{MF}}$ $\displaystyle=H_{t}+H_{J}+H_{\delta},$ (104)
$\displaystyle H_{\delta}$
$\displaystyle={\textstyle\sum_{i}}\left[\lambda_{i}(\hat{n}^{b}_{i}-1+\delta)-\mu_{i}(\hat{n}^{h}_{i}-\delta)\right].$
(105)

Here $H_{\delta}$ is added to impose the no-double-occupancy constraint Eq.
(6), in which

$\hat{n}^{b}_{i}={\textstyle\sum_{\sigma}}b^{\dagger}_{i\sigma}b_{i\sigma},\quad\hat{n}^{h}_{i}=h^{\dagger}_{i}h_{i}.$

Finally, we redefine $\lambda_{i}$ by a constant shift

$\displaystyle\lambda_{i}-\frac{1}{4}\alpha J(1-\delta)\to\lambda_{i},$

and separate $H_{\text{MF}}$ to a holon part $H_{h}$, a spinon part $H_{b}$
and a number term $H_{0}$, leading to Eqs. (10) to (13) in the main text.

###  B.2 Spinon condensation at zero temperature

At zero temperature, Bose condensation of spinons will occur at $k=0$, as
$\lambda$ decreases to the critical value $\lambda_{c}$ determined from
$\min_{k}E^{b}_{k-}=E^{b}_{0-}=0$, yielding

$\lambda_{c}=\sqrt{p^{2}+q^{2}}\ \alpha.$ (106)

We need to modify the mean-field equations Eqs. (30) to (32) by including the
density of condensed spinons at $k=0$:

$n_{0}\equiv\frac{1}{N}\sum_{\sigma}\langle{b^{s\dagger}_{0\sigma}b^{s}_{0\sigma}}\rangle\quad(s=A,B),$
(107)

which is the same on the two sub-lattices. In the $B,D=0$ (AFM) phase,
$\lim_{\beta\to\infty}\chi_{0}\to 0$, and $n_{0}$ is equal to the limit

$n_{0}=\lim_{\beta\to\infty}\frac{1}{N}\frac{\lambda}{\chi_{0}}[1+2n_{b}(\chi_{0})].$
(108)

In phases with $B,D\neq 0$, $\chi_{0}\neq 0$. Then $n_{b}(E^{b}_{0-})$ becomes
a macroscopic number $N_{0}\gg 1$, which can be related to $n_{0}$ by

$n_{0}=\frac{1}{N}\frac{\lambda}{\chi_{0}}(1+N_{0})\approx\frac{\lambda}{\chi_{0}}\frac{N_{0}}{N}.$
(109)

Then the modified self-consistency equations at $T=0$ are

$\displaystyle 1-\delta$ $\displaystyle=n_{0}+\frac{1}{N}\sum_{k\neq
0}\Big{[}\frac{\lambda}{\chi_{k}}-1\Big{]},$ (110) $\displaystyle A$
$\displaystyle=\frac{q\alpha}{2\lambda}n_{0}+\frac{q}{2\alpha N}\sum_{k\neq
0}\frac{|\gamma_{k}|^{2}}{\chi_{k}},$ (111) $\displaystyle B$
$\displaystyle=-\frac{p\alpha}{2\lambda}n_{0}(1-\delta_{\chi_{0}}),$ (112)

where $\delta_{x}=1$ if $x=0$ and $0$ otherwise. In particular, at half-
filling ($\delta=0$), we get $A=0.605$ and $n_{0}=0.484$ as $N\to\infty$. In
the FM phase ($A=0$), we get maximum condensation $n_{0}=1-\delta$.

##  Appendix C Effective holon interaction due to spinons

This appendix derives Eq. (47), the effective holon interaction by exchanging
two spinons. At very small doping, we take the approximation $B,D\approx 0$.
Then the spinon part of the Hamiltonian Eq. (20) is simplified to

$H^{b}_{k}=\begin{bmatrix}\lambda&&&-\Delta_{k}\\\ &\lambda&\Delta^{*}_{k}&\\\
&\Delta_{k}&\lambda&\\\ -\Delta^{*}_{k}&&&\lambda\end{bmatrix},$ (113)

where $\Delta_{k}=J_{b}\gamma_{k}$ (with $J_{b}=JA$). We then Fourier
transform the holons and spinons to momentum-frequency representation:

$\displaystyle h^{s}_{k}(\tau)$
$\displaystyle=\frac{1}{\sqrt{\beta}}\sum_{k,\omega}e^{i\omega\tau}h^{s}_{k\omega},$
(114) $\displaystyle b^{s}_{k\sigma}(\tau)$
$\displaystyle=\frac{1}{\sqrt{\beta}}\sum_{k,\nu}e^{i\nu\tau}b^{s}_{k\nu\sigma},$
(115)

where $\omega,\nu$ sum over discrete fermion and boson Matsubara frequencies,
respectively. The spinon part of the action is

$S_{b}=\sum_{k,k^{\prime}}\sum_{\nu^{\prime},\nu^{\prime}}\left[\bar{b}_{k\nu\uparrow},b_{-k,-\nu\downarrow}\right](G_{0})^{-1}_{k\nu,k^{\prime}\nu^{\prime}}\begin{bmatrix}b_{k^{\prime}\nu^{\prime}\uparrow}\\\\[1.99997pt]
\bar{b}_{-k^{\prime},-\nu^{\prime}\downarrow}\end{bmatrix},$ (116)

where we introduced the shorthand notation

$b_{k\nu\sigma}\equiv\begin{bmatrix}b^{A}_{k\nu\sigma}\\\\[3.00003pt]
b^{B}_{k\nu\sigma}\end{bmatrix},$

and the matrix $G_{0}$ is (with
$\chi_{k}=\sqrt{\lambda^{2}-|\Delta_{k}|^{2}}$):

$\displaystyle(G_{0})_{k\nu,k^{\prime}\nu^{\prime}}=(\mathcal{G}_{0})_{k\nu}\delta_{k-k^{\prime}}\delta_{\nu-\nu^{\prime}},$
(117)
$\displaystyle(\mathcal{G}_{0})_{k\nu}\equiv\frac{1}{\nu^{2}+\chi_{k}^{2}}\begin{bmatrix}\lambda-i\nu&&&\Delta_{k}\\\
&\lambda-i\nu&-\Delta^{*}_{k}&\\\ &-\Delta_{k}&\lambda+i\nu&\\\
\Delta^{*}_{k}&&&\lambda+i\nu\\\ \end{bmatrix}.$ (118)

Action terms corresponding to spinon-holon interaction Eq. (39) is Fourier
transformed to (with $\gamma_{k}=\sum_{\eta}e^{ik\cdot\eta}$):

$\displaystyle
S_{\text{int}}=\sum_{k,k^{\prime}}\sum_{\nu^{\prime},\nu^{\prime}}\left[\bar{b}_{k\nu\uparrow},b_{-k,-\nu\downarrow}\right]T_{k\nu,k^{\prime}\nu^{\prime}}\begin{bmatrix}b_{k^{\prime}\nu^{\prime}\uparrow}\\\\[1.99997pt]
\bar{b}_{-k^{\prime},-\nu^{\prime}\downarrow}\end{bmatrix},$ (119)
$\displaystyle T_{k\nu,k^{\prime}\nu^{\prime}}=\begin{bmatrix}&t^{BA}_{1}&&\\\
t^{AB}_{1}&&&\\\ &&&t^{AB}_{2}\\\
&&t^{BA}_{2}&\end{bmatrix}_{k\nu,k^{\prime}\nu^{\prime}}.$ (120)

The matrix elements of $T$ are

$\displaystyle t^{AB}_{1,k\nu,k^{\prime}\nu^{\prime}}$
$\displaystyle=\frac{t}{\beta
N}\sum_{q}{\sum_{\omega,\omega^{\prime}}}^{\prime}\bar{h}^{A}_{k^{\prime}+q,\omega}h^{B}_{k+q,\omega^{\prime}}\gamma_{q},$
(121) $\displaystyle t^{AB}_{2,k\nu,k^{\prime}\nu^{\prime}}$
$\displaystyle=\frac{t}{\beta
N}\sum_{q}{\sum_{\omega,\omega^{\prime}}}^{\prime}\bar{h}^{A}_{k^{\prime}+q,\omega}h^{B}_{k+q,\omega^{\prime}}\gamma_{q+k+k^{\prime}},$
(122) $\displaystyle t^{BA}_{1,k\nu,k^{\prime}\nu^{\prime}}$
$\displaystyle=\frac{t}{\beta
N}\sum_{q}{\sum_{\omega,\omega^{\prime}}}^{\prime}\bar{h}^{B}_{k^{\prime}+q,\omega^{\prime}}h^{A}_{k+q,\omega}\gamma^{*}_{q},$
(123) $\displaystyle t^{BA}_{2,k\nu,k^{\prime}\nu^{\prime}}$
$\displaystyle=\frac{t}{\beta
N}\sum_{q}{\sum_{\omega,\omega^{\prime}}}^{\prime}\bar{h}^{B}_{k^{\prime}+q,\omega^{\prime}}h^{A}_{k+q,\omega}\gamma^{*}_{q+k+k^{\prime}},$
(124)

where

${\sum_{\omega,\omega^{\prime}}}^{\prime}\equiv\sum_{\omega,\omega^{\prime}}\delta_{(\nu+\omega)-(\nu^{\prime}+\omega^{\prime})}.$
(125)

To obtain an effective theory for the holons, we integrate out the spinons.
Collecting terms involving spinons in the action, we get

$S_{b}+S_{\text{int}}=\sum_{k,k^{\prime}}\sum_{\nu^{\prime},\nu^{\prime}}\left[\bar{b}_{k\nu\uparrow},b_{-k,-\nu\downarrow}\right]G^{-1}_{k\nu,k^{\prime}\nu^{\prime}}\begin{bmatrix}b_{k^{\prime}\nu^{\prime}\uparrow}\\\\[1.99997pt]
\bar{b}_{-k^{\prime},-\nu^{\prime}\downarrow}\end{bmatrix},$

where the matrix $G^{-1}$ is given by

$G^{-1}_{k\nu,k^{\prime}\nu^{\prime}}=(G_{0}^{-1}+T)_{k\nu,k^{\prime}\nu^{\prime}}.$
(126)

The effective action for holons $S^{\text{eff}}_{h}[\bar{h},h]$ is defined by

$\displaystyle Z$ $\displaystyle=\int
D\bar{h}\,Dh\,e^{-S^{\text{eff}}_{h}[\bar{h},h]},$ $\displaystyle
e^{-S^{\text{eff}}_{h}}$ $\displaystyle\equiv\int
D\bar{b}\,Db\,e^{-S}=e^{-S_{h}}\int
D\bar{b}\,Db\,e^{-(S_{b}+S_{\text{int}})}.$

Performing Gaussian integration over spinons, we get

$\displaystyle\int D\bar{b}\,Db\,e^{-(S_{b}+S_{\text{int}})}$
$\displaystyle\propto(\det G^{-1})^{-1}=\exp(-\ln\det G^{-1}).$ (127)

Thus the effective holon action is given by

$S^{\text{eff}}_{h}=S_{h}+\ln\det G^{-1}.$ (128)

The log-determinant can be expanded as a power series of $t$; using $\ln\det
A=\operatorname{tr}\ln A$ and $\ln(1+x)=\sum_{n=1}^{\infty}(-1)^{n+1}x^{n}/n$,
we obtain

$\displaystyle\ln\det G^{-1}$
$\displaystyle=\operatorname{tr}\ln(G_{0}^{-1}+T)$
$\displaystyle=\operatorname{tr}[\ln G_{0}^{-1}+\ln(1+G_{0}T)]$
$\displaystyle=\operatorname{tr}\ln
G_{0}^{-1}+\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\operatorname{tr}(G_{0}T)^{n}.$
(129)

The term $\operatorname{tr}\ln G_{0}^{-1}$ is a constant independent of $h$
and can be omitted. Then the effective holon action is

$\displaystyle S^{\text{eff}}_{h}[\bar{h},h]$
$\displaystyle=S_{h}^{0}+\sum_{n=1}^{\infty}S_{h}^{n},$ (130) $\displaystyle
S_{h}^{n}$
$\displaystyle\equiv\frac{(-1)^{n+1}}{n}\operatorname{tr}(G_{0}T)^{n}.$ (131)

We keep up to second order terms. The first order term
$\operatorname{tr}(G_{0}T)=0$. The second order term is

$\displaystyle S_{h}^{2}$
$\displaystyle\equiv-\frac{1}{2}\operatorname{tr}(G_{0}T)^{2}=-\frac{1}{2}\sum_{k,\nu}\sum_{k_{1},\nu_{1}}\sum_{k_{2},\nu_{2}}\sum_{k^{\prime},\nu^{\prime}}\operatorname{tr}[(G_{0})_{k\nu,k_{1}\nu_{1}}T_{k_{1}\nu_{1},k_{2}\nu_{2}}(G_{0})_{k_{2}\nu_{2},k^{\prime}\nu^{\prime}}T_{k^{\prime}\nu^{\prime},k\nu}]$
$\displaystyle=\frac{1}{2}\sum_{k,\nu}\sum_{k^{\prime},\nu^{\prime}}\frac{1}{(\nu^{2}+\chi_{k}^{2})(\nu^{\prime
2}+\chi_{k^{\prime}}^{2})}\Big{[}\Delta_{k}\Delta_{k^{\prime}}(t^{AB}_{1}t^{\prime
BA}_{2}+t^{BA}_{2}t^{\prime
AB}_{1})+\Delta^{*}_{k}\Delta^{*}_{k^{\prime}}(t^{BA}_{1}t^{\prime
AB}_{2}+t^{AB}_{2}t^{\prime BA}_{1})$
$\displaystyle\qquad-(i\nu-\lambda)(i\nu^{\prime}-\lambda)(t^{AB}_{1}t^{\prime
BA}_{1}+t^{BA}_{1}t^{\prime
AB}_{1})-(i\nu+\lambda)(i\nu+\lambda^{\prime})(t^{AB}_{2}t^{\prime
BA}_{2}+t^{BA}_{2}t^{\prime AB}_{2})\Big{]}.$

Here we simply write
$t^{ss^{\prime}}_{i}=t^{ss^{\prime}}_{i,k\nu,k^{\prime}\nu^{\prime}}$,
$t^{\prime ss^{\prime}}_{i}=t^{ss^{\prime}}_{i,k^{\prime}\nu^{\prime},k\nu}$
(where $s,s^{\prime}=A,B$ and $i=1,2$). Note that each term is invariant under
the exchange of $(k,\nu)\leftrightarrow(k^{\prime},\nu^{\prime})$. Then

$\displaystyle S_{h}^{2}$
$\displaystyle=\sum_{k,\nu}\sum_{k^{\prime},\nu^{\prime}}\frac{1}{(\nu^{2}+\chi_{k}^{2})(\nu^{\prime
2}+\chi_{k^{\prime}}^{2})}\Big{[}\Delta_{k}\Delta_{k^{\prime}}t^{BA}_{2}t^{\prime
AB}_{1}+\Delta^{*}_{k}\Delta^{*}_{k^{\prime}}t^{BA}_{1}t^{\prime AB}_{2}$
$\displaystyle\qquad-(i\nu-\lambda)(i\nu^{\prime}-\lambda)t^{BA}_{1}t^{\prime
AB}_{1}-(i\nu+\lambda)(i\nu^{\prime}+\lambda)t^{BA}_{2}t^{\prime
AB}_{2}\Big{]}.$

Evaluating summation over $\nu^{\prime}$ first after substituting in the
definition of $t_{i}^{ss^{\prime}}$ (Eqs. (121) to (124)), we obtain:

$\displaystyle S_{h}^{2}$ $\displaystyle=\bigg{(}\frac{t}{\beta
N}\bigg{)}^{2}\sum_{k,\nu}\sum_{k^{\prime},\nu^{\prime}}\sum_{q_{1},q_{2}}\sum_{\omega_{1},\omega^{\prime}_{1}}\sum_{\omega_{2},\omega^{\prime}_{2}}\frac{\delta_{(\nu+\omega_{1})-(\nu^{\prime}+\omega^{\prime}_{1})}\delta_{(\nu^{\prime}+\omega_{2})-(\nu+\omega^{\prime}_{2})}}{(\nu^{2}+\chi_{k}^{2})(\nu^{\prime
2}+\chi_{k^{\prime}}^{2})}\bar{h}^{B}_{k^{\prime}+q_{1},\omega^{\prime}_{1}}h^{A}_{k+q_{1},\omega_{1}}\bar{h}^{A}_{k+q_{2},\omega_{2}}h^{B}_{k^{\prime}+q_{2},\omega^{\prime}_{2}}$
$\displaystyle\quad\times\Big{[}\Delta_{k}\Delta_{k^{\prime}}\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}}+\Delta^{*}_{k}\Delta^{*}_{k^{\prime}}\gamma^{*}_{q_{1}}\gamma_{q_{2}+k+k^{\prime}}-(i\nu-\lambda)(i\nu^{\prime}-\lambda)\gamma^{*}_{q_{1}}\gamma_{q_{2}}-(i\nu+\lambda)(i\nu^{\prime}+\lambda)\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}+k+k^{\prime}}\Big{]}$
$\displaystyle=\bigg{(}\frac{t}{\beta
N}\bigg{)}^{2}\sum_{k,k^{\prime},\nu}\sum_{q_{1},q_{2}}{\sum_{\\{\omega\\}}}^{\prime}\frac{1}{(\nu^{2}+\chi_{k}^{2})[(\nu+\omega_{1}-\omega^{\prime}_{1})^{2}+\chi_{k^{\prime}}^{2}]}\bar{h}^{B}_{k^{\prime}+q_{1},\omega^{\prime}_{1}}h^{A}_{k+q_{1},\omega_{1}}\bar{h}^{A}_{k+q_{2},\omega_{2}}h^{B}_{k^{\prime}+q_{2},\omega^{\prime}_{2}}$
$\displaystyle\quad\times\Big{[}\Delta_{k}\Delta_{k^{\prime}}\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}}+\Delta^{*}_{k}\Delta^{*}_{k^{\prime}}\gamma^{*}_{q_{1}}\gamma_{q_{2}+k+k^{\prime}}$
$\displaystyle\qquad-(i\nu-\lambda)[i(\nu+\omega_{1}-\omega^{\prime}_{1})-\lambda]\gamma^{*}_{q_{1}}\gamma_{q_{2}}-(i\nu+\lambda)[i(\nu+\omega_{1}-\omega^{\prime}_{1})+\lambda]\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}+k+k^{\prime}}\Big{]},$
(132)

where

${\sum_{\\{\omega\\}}}^{\prime}\equiv\sum_{\omega_{1},\omega^{\prime}_{1}}\sum_{\omega_{2},\omega^{\prime}_{2}}\delta_{(\omega_{1}+\omega_{2})-(\omega^{\prime}_{1}+\omega^{\prime}_{2})}$

We keep only the _instantaneous_ holon interaction by taking the _static
limit_ , i.e. setting $\omega_{1}=\omega^{\prime}_{1}$ in the coefficient of
$\bar{h}^{B}h^{A}\bar{h}^{A}h^{B}$. Then Eq. 132 reduces to

$\displaystyle S_{h}^{2}$ $\displaystyle=\bigg{(}\frac{t}{\beta
N}\bigg{)}^{2}\sum_{k,k^{\prime},\nu}\sum_{q_{1},q_{2}}{\sum_{\\{\omega\\}}}^{\prime}\frac{\bar{h}^{B}_{k^{\prime}+q_{1},\omega^{\prime}_{1}}h^{A}_{k+q_{1},\omega_{1}}\bar{h}^{A}_{k+q_{2},\omega_{2}}h^{B}_{k^{\prime}+q_{2},\omega^{\prime}_{2}}}{(\nu^{2}+\chi_{k}^{2})(\nu^{2}+\chi_{k^{\prime}}^{2})}$
$\displaystyle\times\Big{[}\Delta_{k}\Delta_{k^{\prime}}\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}}+\Delta^{*}_{k}\Delta^{*}_{k^{\prime}}\gamma^{*}_{q_{1}}\gamma_{q_{2}+k+k^{\prime}}-(i\nu-\lambda)^{2}\gamma^{*}_{q_{1}}\gamma_{q_{2}}-(i\nu+\lambda)^{2}\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}+k+k^{\prime}}\Big{]}.$
(133)

In the zero-temperature limit, the summation over $\nu$ can be replaced by an
integral:

$\displaystyle\frac{1}{\beta}\sum_{\nu}$
$\displaystyle\to\int\frac{d\nu}{2\pi},$
$\displaystyle\int\frac{d\nu}{2\pi}\frac{1}{(\nu^{2}+\chi_{k}^{2})(\nu^{2}+\chi_{k^{\prime}}^{2})}$
$\displaystyle=\frac{1}{2\chi_{k}\chi_{k^{\prime}}(\chi_{k}+\chi_{k^{\prime}})},$
$\displaystyle\int\frac{d\nu}{2\pi}\frac{(i\nu\pm\lambda)^{2}}{(\nu^{2}+\chi_{k}^{2})(\nu^{2}+\chi_{k^{\prime}}^{2})}$
$\displaystyle=\frac{\lambda^{2}-\chi_{k}\chi_{k^{\prime}}}{2\chi_{k}\chi_{k^{\prime}}(\chi_{k}+\chi_{k^{\prime}})}.$

We can then inverse Fourier transform $\omega$ to $\tau$:

$\displaystyle\frac{1}{\beta}{\sum_{\\{\omega\\}}}^{\prime}\bar{h}^{B}_{k^{\prime}+q_{1},\omega^{\prime}_{1}}\bar{h}^{A}_{k+q_{2},\omega_{2}}h^{A}_{k+q_{1},\omega_{1}}h^{B}_{k^{\prime}+q_{2},\omega^{\prime}_{2}}$
$\displaystyle=\int
d\tau\,\bar{h}^{B}_{k^{\prime}+q_{1},\tau}\bar{h}^{A}_{k+q_{2},\tau}h^{A}_{k+q_{1},\tau}h^{B}_{k^{\prime}+q_{2},\tau}.$
(134)

Therefore (the order of $h$ is changed)

$\displaystyle
S_{h}^{2}=-\frac{t^{2}}{2N^{2}}\sum_{k,k^{\prime}}\sum_{q_{1},q_{2}}\frac{\bar{h}^{B}_{k^{\prime}+q_{1}}\bar{h}^{A}_{k+q_{2}}h^{A}_{k+q_{1}}h^{B}_{k^{\prime}+q_{2}}}{\chi_{k}\chi_{k^{\prime}}(\chi_{k}+\chi_{k^{\prime}})}$
$\displaystyle\ \
\times\Big{[}\Delta_{k}\Delta_{k^{\prime}}\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}}+\Delta^{*}_{k}\Delta^{*}_{k^{\prime}}\gamma^{*}_{q_{1}}\gamma_{q_{2}+k+k^{\prime}}$
$\displaystyle\quad+(\chi_{k}\chi_{k^{\prime}}-\lambda^{2})(\gamma^{*}_{q_{1}}\gamma_{q_{2}}+\gamma^{*}_{q_{1}+k+k^{\prime}}\gamma_{q_{2}+k+k^{\prime}})\Big{]}.$
(135)

Redefining summation variables:

$\left.\begin{aligned} k_{1}&=k+q_{1}\\\ k_{2}&=k^{\prime}+q_{2}\\\
q&=k^{\prime}-k\\\ p&=-k^{\prime}\end{aligned}\right\\}\ \Rightarrow\
\left\\{\begin{aligned} k&=-p-q\\\ k^{\prime}&=-p\\\ q_{1}&=k_{1}+p+q\\\
q_{2}&=k_{2}+p\end{aligned}\right.,$

we obtain

$\displaystyle
S_{h}^{2}=-\frac{t^{2}}{2N^{2}}\sum_{k_{1},k_{2}}\sum_{p,q}\frac{\bar{h}^{B}_{k_{1}+q}\bar{h}^{A}_{k_{2}-q}h^{A}_{k_{1}}h^{B}_{k_{2}}}{\chi_{p+q}\chi_{p}(\chi_{p+q}+\chi_{p})}$
$\displaystyle\
\times\Big{[}\Delta^{*}_{p+q}\Delta^{*}_{p}\gamma^{*}_{k_{1}-p}\gamma_{k_{2}+p}+\Delta_{p+q}\Delta_{p}\gamma^{*}_{k_{1}+p+q}\gamma_{k_{2}-p-q}$
$\displaystyle\quad+(\chi_{p+q}\chi_{p}-\lambda^{2})(\gamma^{*}_{k_{1}+p+q}\gamma_{k_{2}+p}+\gamma^{*}_{k_{1}-p}\gamma_{k_{2}-p-q})\Big{]}.$
(136)

Here we used

$\chi_{-k}=\chi_{k},\quad\gamma_{-k}=\gamma^{*}_{k},\quad\Delta_{-k}=\Delta^{*}_{k}.$

From this effective action, we read off terms in the Hamiltonian that
represent interaction between holons (rename $k_{1},k_{2}$ to $k,k^{\prime}$):

$\displaystyle
H^{\text{eff}}_{\text{int}}=\frac{1}{N}\sum_{k,k^{\prime},q}V_{kk^{\prime}q}h^{B\dagger}_{k+q}h^{A\dagger}_{k^{\prime}-q}h^{A}_{k}h^{B}_{k^{\prime}},$
(137) $\displaystyle
V_{kk^{\prime}q}\equiv-\frac{t^{2}}{2N}\sum_{p}\frac{1}{\chi_{p+q}\chi_{p}(\chi_{p+q}+\chi_{p})}$
$\displaystyle\ \
\times\Big{[}\Delta^{*}_{p+q}\Delta^{*}_{p}\gamma^{*}_{k-p}\gamma_{k^{\prime}+p}+\Delta_{p+q}\Delta_{p}\gamma^{*}_{k+p+q}\gamma_{k^{\prime}-p-q}$
$\displaystyle\quad+(\chi_{p+q}\chi_{p}-\lambda^{2})(\gamma^{*}_{k+p+q}\gamma_{k^{\prime}+p}+\gamma^{*}_{k-p}\gamma_{k^{\prime}-p-q})\Big{]}.$
(138)

This effective holon interaction is due to the exchange of two spinons of
momenta $p$ and $p+q$. Let us call the holon momenta as

$k+q=k_{1},\quad k=k^{\prime}_{1},\quad k^{\prime}=k^{\prime}_{2},\quad
k^{\prime}_{2}+k^{\prime}_{1}-k_{1}=k_{2}.$

Then we get an alternative expression of $H^{\text{eff}}_{\text{int}}$:

$\displaystyle
H^{\text{eff}}_{\text{int}}=\frac{1}{N}\sum_{\\{k,k^{\prime}\\}}V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}h^{B\dagger}_{k_{1}}h^{A\dagger}_{k_{2}}h^{A}_{k^{\prime}_{1}}h^{B}_{k^{\prime}_{2}},$
(139) $\displaystyle
V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}\equiv-\frac{t^{2}}{2N}\sum_{p}\frac{\delta_{(k_{1}+k_{2})-(k^{\prime}_{1}+k^{\prime}_{2})}}{\chi_{p}\chi_{p+k_{1}-k^{\prime}_{1}}(\chi_{p}+\chi_{p+k_{1}-k^{\prime}_{1}})}$
$\displaystyle\times\Big{[}\Delta^{*}_{p}\Delta^{*}_{p+k_{1}-k^{\prime}_{1}}\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k^{\prime}_{2}+p}+\Delta_{p}\Delta_{p+k_{1}-k^{\prime}_{1}}\gamma^{*}_{k_{1}+p}\gamma_{k_{2}-p}$
$\displaystyle\
+(\chi_{p}\chi_{p+k_{1}-k^{\prime}_{1}}-\lambda^{2})(\gamma^{*}_{k_{1}+p}\gamma_{k^{\prime}_{2}+p}+\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}-p})\Big{]}.$
(140)

Define $p^{\prime}\equiv p+k_{1}-k^{\prime}_{1}$. We then rewrite the
expression by keeping $k^{\prime}_{1},k_{2}$ (the momenta of $h^{A}$) only:

$\displaystyle
V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}\equiv-\frac{t^{2}}{2N}\sum_{p}\bigg{\\{}\frac{\delta_{(k_{1}+k_{2})-(k^{\prime}_{1}+k^{\prime}_{2})}}{\chi_{p}\chi_{p^{\prime}}(\chi_{p}+\chi_{p^{\prime}})}$
$\displaystyle\times\Big{[}\Delta^{*}_{p}\Delta^{*}_{p^{\prime}}\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}+p^{\prime}}+\Delta_{p}\Delta_{p^{\prime}}\gamma^{*}_{k^{\prime}_{1}+p^{\prime}}\gamma_{k_{2}-p}$
$\displaystyle\
+(\chi_{p}\chi_{p^{\prime}}-\lambda^{2})(\gamma^{*}_{k^{\prime}_{1}+p^{\prime}}\gamma_{k_{2}+p^{\prime}}+\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}-p})\Big{]}\bigg{\\}}.$

Finally, with

${\sum_{\\{k\\}}}^{\prime}\equiv\sum_{k_{1},k_{2}}\sum_{k^{\prime}_{1},k^{\prime}_{2}}\delta_{(k_{1}+k_{2})-(k^{\prime}_{1}+k^{\prime}_{2})}$

we get (substituting in $\Delta_{k}=J_{b}\gamma_{k}$)

$\displaystyle
H^{\text{eff}}_{\text{int}}=\frac{1}{N}{\sum_{\\{k\\}}}^{\prime}V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}h^{B\dagger}_{k_{1}}h^{A\dagger}_{k_{2}}h^{A}_{k^{\prime}_{1}}h^{B}_{k^{\prime}_{2}}.$
(141) $\displaystyle
V_{k_{1}k_{2}k^{\prime}_{1}k^{\prime}_{2}}=-\frac{t^{2}}{2N}\sum_{p}\bigg{\\{}\frac{1}{\chi_{p}\chi_{p^{\prime}}(\chi_{p}+\chi_{p^{\prime}})}$
$\displaystyle\ \
\times\Big{[}J_{b}^{2}(\gamma^{*}_{p}\gamma^{*}_{p^{\prime}}\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}+p^{\prime}}+\gamma_{p}\gamma_{p^{\prime}}\gamma^{*}_{k^{\prime}_{1}+p}\gamma_{k_{2}-p})$
$\displaystyle\quad+(\chi_{p}\chi_{p^{\prime}}-\lambda^{2})(\gamma^{*}_{k^{\prime}_{1}+p}\gamma_{k_{2}+p^{\prime}}+\gamma^{*}_{k^{\prime}_{1}-p}\gamma_{k_{2}-p})\Big{]}\bigg{\\}},$
(142)

which is Eq. (47).

##  Appendix D Continuum limit of the _t-J_ model Lagrangian

This section derives the Lagrangian Eqs. (68) to (71) in the continuum limit.
The following identities on the honeycomb lattice are useful:

$\displaystyle\textstyle\sum_{\delta}\delta\cdot\mathbf{A}$ $\displaystyle=0,$
$\displaystyle\textstyle\sum_{\delta}(\delta\cdot\nabla)^{2}f$
$\displaystyle=\textstyle(\alpha/2)a^{2}\nabla^{2}f,$
$\displaystyle\textstyle\sum_{\delta}(\delta\cdot\mathbf{A})(\delta\cdot\mathbf{B})$
$\displaystyle=\textstyle(\alpha/2)a^{2}\mathbf{A}\cdot\mathbf{B},$
$\displaystyle\textstyle\sum_{\delta}(\delta\cdot\mathbf{A})^{2}$
$\displaystyle=\textstyle(\alpha/2)a^{2}\mathbf{A}^{2}.$

Here $a$ is the nearest neighbor distance, $\delta$ sums over nearest
neighbors of a site on sub-lattice $A$ (see Fig. 1), $f,g$ are any scalar
fields and $\mathbf{A},\mathbf{B}$ are any vector fields. The continuum limit
of $L_{h}$ (Eq. (62)) and $L_{0}$ (Eq. (65)) are

$\displaystyle L_{h}$
$\displaystyle=\int\frac{d^{2}x}{2\Omega}\sum_{s}\Big{[}\bar{h}^{s}(\partial_{\tau}+isA_{\tau})h^{s}+\lambda\bar{h}^{s}h^{s}\Big{]},$
(143) $\displaystyle L_{0}$
$\displaystyle=\int\frac{d^{2}x}{2\Omega}\bigg{[}-\sum_{s}(\lambda+isA_{\tau})+\frac{\alpha
J}{2}\Delta^{2}\bigg{]}$
$\displaystyle=\int\frac{d^{2}x}{2\Omega}\bigg{[}-2\lambda+\frac{\alpha
J}{2}\Delta^{2}\bigg{]}.$ (144)

The continuum limit of the spinon-holon interaction $L_{t}$ (Eq. (64)) is

$\displaystyle L_{t}\to t\int\frac{d^{2}x}{2\Omega}\sum_{\delta,\sigma}$
$\displaystyle\quad\Big{[}\bar{h}^{B}(x+\delta)h^{A}(x)\bar{b}^{A}_{\sigma}(x)b^{B}_{\sigma}(x+\delta)+h.c.\Big{]}$
$\displaystyle\approx
t\int\frac{d^{2}x}{2\Omega}\sum_{\delta,\sigma}\Big{\\{}[(1+\delta\cdot\nabla+\tfrac{1}{2}(\delta\cdot\nabla)^{2})\bar{h}^{B}]h^{A}$
$\displaystyle\quad\times\bar{b}^{A}_{\sigma}[(1+\delta\cdot\nabla+\tfrac{1}{2}(\delta\cdot\nabla)^{2})b^{B}_{\sigma}]+h.c.\Big{\\}}.$

We expand this in powers of the nearest neighbor distance $a$. The first order
terms in $a$ turn out to vanish. The second order terms are not important in
describing the spinon-holon interaction. We then only keep the zeroth order
terms:

$\displaystyle L_{t}$ $\displaystyle\approx
t\int\frac{d^{2}x}{2\Omega}\sum_{\delta,\sigma}(\bar{h}^{B}h^{A}\bar{b}^{A}_{\sigma}b^{B}_{\sigma}+h.c.)$
$\displaystyle=\alpha
t\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}(\bar{h}^{B}h^{A}\bar{b}^{A}_{\sigma}b^{B}_{\sigma}+h.c.).$
(145)

Finally we calculate the continuum limit of $L_{\text{heis}}$ (Eq. (63)). Let
us separate it to two parts:

$\displaystyle L_{\text{heis}}$ $\displaystyle=L_{b}^{0}+L_{J},$ (146)
$\displaystyle L_{b}^{0}$ $\displaystyle=\sum_{s,\delta}\sum_{i\in
s}\Big{[}\bar{b}_{i\sigma}(\partial_{\tau}+isA_{\tau}(i))b_{i\sigma}+\lambda\bar{b}_{i\sigma}b_{i\sigma}\Big{]},$
(147) $\displaystyle L_{J}$
$\displaystyle=-\frac{J\Delta}{2}\sum_{\begin{subarray}{c}i\in A\\\
\delta,\sigma\end{subarray}}\sigma\Big{[}e^{i\delta\cdot\mathbf{A}(i,i+\delta)}\bar{b}_{i\sigma}\bar{b}_{i+\delta,-\sigma}+h.c.\Big{]}.$
(148)

The continuum limit of $L_{b}^{0}$ is

$L_{b}^{0}\to\int\frac{d^{2}x}{2\Omega}\sum_{s,\sigma}\Big{[}\bar{b}^{s}_{\sigma}(\partial_{\tau}+isA_{\tau})b^{s}_{\sigma}+\lambda\bar{b}^{s}_{\sigma}b^{s}_{\sigma}\Big{]}.$
(149)

Next,

$\displaystyle
L_{J}\to-\frac{J\Delta}{2}\int\frac{d^{2}x}{2\Omega}\sum_{\delta,\sigma}\sigma$
$\displaystyle\times\Big{[}e^{i\delta\cdot\mathbf{A}(x)}\bar{b}^{A}_{\sigma}(x-\tfrac{\delta}{2})\bar{b}^{B}_{-\sigma}(x+\tfrac{\delta}{2})+h.c.\Big{]}$
$\displaystyle\approx-\frac{J\Delta}{2}\int\frac{d^{2}x}{2\Omega}\sum_{\delta,\sigma}\sigma\bigg{\\{}[1+i\delta\cdot\mathbf{A}+\tfrac{1}{2}(i\delta\cdot\mathbf{A})^{2}]$
$\displaystyle\quad\times\Big{[}(1-\tfrac{\delta}{2}\cdot\nabla+\tfrac{1}{2}(\tfrac{\delta}{2}\cdot\nabla)^{2})\bar{b}^{A}_{\sigma}$
$\displaystyle\qquad\cdot(1+\tfrac{\delta}{2}\cdot\nabla+\tfrac{1}{2}(\tfrac{\delta}{2}\cdot\nabla)^{2})\bar{b}^{B}_{-\sigma}\Big{]}+h.c.\bigg{\\}}.$
(150)

We keep terms up to second order in $a$. The zeroth order terms in $a$ are

$\displaystyle L_{J}^{0}$
$\displaystyle=-4Q\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\sigma(\bar{b}^{A}_{\sigma}\bar{b}^{B}_{-\sigma}+h.c.).$
(151)

Here $Q\equiv\alpha J\Delta/8$. The first order terms $L_{J}^{1}=0$. Using the
identity

$\displaystyle\int
d^{2}x\,\Big{[}(\nabla^{2}\bar{b}^{A}_{\sigma})\bar{b}^{B}_{-\sigma}+\bar{b}^{A}_{\sigma}(\nabla^{2}\bar{b}^{B}_{-\sigma})\Big{]}$
$\displaystyle=\int
d^{2}x\,\Big{[}\nabla^{2}(\bar{b}^{A}_{\sigma}\bar{b}^{B}_{-\sigma})-2\nabla\bar{b}^{A}_{\sigma}\cdot\nabla\bar{b}^{B}_{-\sigma}\Big{]}$
$\displaystyle=-2\int
d^{2}x\,\nabla\bar{b}^{A}_{\sigma}\cdot\nabla\bar{b}^{B}_{-\sigma},$ (152)

we calculate the second order terms in $a$:

$\displaystyle
L_{J}^{2}=-\frac{J\Delta}{2}\int\frac{d^{2}x}{2\Omega}\sum_{\delta,\sigma}\sigma\bigg{\\{}\Big{[}\tfrac{1}{2}[(-\tfrac{\delta}{2}\cdot\nabla)^{2}\bar{b}^{A}_{\sigma}]\bar{b}^{B}_{-\sigma}+\bar{b}^{A}_{\sigma}\tfrac{1}{2}[(\tfrac{\delta}{2}\cdot\nabla)^{2}\bar{b}^{B}_{-\sigma}]+(-\tfrac{\delta}{2}\cdot\nabla)\bar{b}^{A}_{\sigma}(\tfrac{\delta}{2}\cdot\nabla)\bar{b}^{B}_{-\sigma}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad+(i\delta\cdot\mathbf{A})\Big{[}[(-\tfrac{\delta}{2}\cdot\nabla)\bar{b}^{A}_{\sigma}]\bar{b}^{B}_{-\sigma}+\bar{b}^{A}_{\sigma}[(\tfrac{\delta}{2}\cdot\nabla)\bar{b}^{B}_{-\sigma}]\Big{]}+\tfrac{1}{2}(i\delta\cdot\mathbf{A})^{2}\bar{b}^{A}_{\sigma}\bar{b}^{B}_{-\sigma}\bigg{\\}}+h.c.$
$\displaystyle=Qa^{2}\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\sigma\Big{\\{}\tfrac{-1}{4}\Big{[}(\nabla^{2}\bar{b}^{A}_{\sigma})\bar{b}^{B}_{-\sigma}+\bar{b}^{A}_{\sigma}(\nabla^{2}\bar{b}^{B}_{-\sigma})\Big{]}+\tfrac{1}{2}\nabla\bar{b}^{A}_{\sigma}\cdot\nabla\bar{b}^{B}_{-\sigma}+i\mathbf{A}\cdot\Big{[}(\nabla\bar{b}^{A}_{\sigma})\bar{b}^{B}_{-\sigma}-\bar{b}^{A}_{\sigma}(\nabla\bar{b}^{B}_{-\sigma})\Big{]}+\mathbf{A}^{2}\bar{b}^{A}_{\sigma}\bar{b}^{B}_{-\sigma}\Big{\\}}+h.c.$
$\displaystyle=Qa^{2}\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\sigma\Big{\\{}\nabla\bar{b}^{A}_{\sigma}\cdot\nabla\bar{b}^{B}_{-\sigma}+i\mathbf{A}\cdot\Big{[}(\nabla\bar{b}^{A}_{\sigma})\bar{b}^{B}_{-\sigma}-\bar{b}^{A}_{\sigma}(\nabla\bar{b}^{B}_{-\sigma})\Big{]}+\mathbf{A}^{2}\bar{b}^{A}_{\sigma}\bar{b}^{B}_{-\sigma}\Big{\\}}+h.c.$
$\displaystyle=Qa^{2}\int\frac{d^{2}x}{2\Omega}\sum_{\sigma}\sigma\Big{[}(\nabla-i\mathbf{A})\bar{b}^{A}_{\sigma}\cdot(\nabla+i\mathbf{A})\bar{b}^{B}_{-\sigma}+h.c.\Big{]}.$
(153)

## References

  * Bednorz and Müller (1986) J George Bednorz and K Alex Müller, “Possible high $T_{c}$ superconductivity in the Ba-La-Cu-O system,” Zeitschrift für Physik B Condensed Matter 64, 189–193 (1986).
  * Zhang and Rice (1988) F. C. Zhang and T. M. Rice, “Effective hamiltonian for the superconducting Cu oxides,” Phys. Rev. B 37, 3759–3761 (1988).
  * Gu  _et al._ (2013) Zheng-Cheng Gu, Hong-Chen Jiang, D. N. Sheng, Hong Yao, Leon Balents, and Xiao-Gang Wen, “Time-reversal symmetry breaking superconducting ground state in the doped mott insulator on the honeycomb lattice,” Phys. Rev. B 88, 155112 (2013).
  * Xu  _et al._ (2022) Zheng-Tao Xu, Zheng-Cheng Gu, and Shuo Yang, “Competing orders in the honeycomb lattice $t$-$J$ model,” arXiv preprint arXiv:2208.13681 (2022).
  * Black-Schaffer  _et al._ (2014) Annica M. Black-Schaffer, Wei Wu, and Karyn Le Hur, “Chiral $d$-wave superconductivity on the honeycomb lattice close to the mott state,” Phys. Rev. B 90, 054521 (2014).
  * Pathak  _et al._ (2010) Sandeep Pathak, Vijay B. Shenoy, and G. Baskaran, “Possible high-temperature superconducting state with a $d+id$ pairing symmetry in doped graphene,” Phys. Rev. B 81, 085431 (2010).
  * Ma  _et al._ (2011) Tianxing Ma, Zhongbing Huang, Feiming Hu, and Hai-Qing Lin, “Pairing in graphene: A quantum monte carlo study,” Phys. Rev. B 84, 121410 (2011).
  * Ying and Wessel (2018) Tao Ying and Stefan Wessel, “Pairing and chiral spin density wave instabilities on the honeycomb lattice: A comparative quantum monte carlo study,” Phys. Rev. B 97, 075127 (2018).
  * Jiang  _et al._ (2014) Shenghan Jiang, Andrej Mesaros, and Ying Ran, “Chiral spin-density wave, spin-charge-chern liquid, and $d+id$ superconductivity in $1/4$-doped correlated electronic systems on the honeycomb lattice,” Phys. Rev. X 4, 031040 (2014).
  * Wang  _et al._ (2012) Wan-Sheng Wang, Yuan-Yuan Xiang, Qiang-Hua Wang, Fa Wang, Fan Yang, and Dung-Hai Lee, “Functional renormalization group and variational monte carlo studies of the electronic instabilities in graphene near $\frac{1}{4}$ doping,” Phys. Rev. B 85, 035414 (2012).
  * Kiesel  _et al._ (2012) Maximilian L. Kiesel, Christian Platt, Werner Hanke, Dmitry A. Abanin, and Ronny Thomale, “Competing many-body instabilities and unconventional superconductivity in graphene,” Phys. Rev. B 86, 020507 (2012).
  * Nandkishore  _et al._ (2012) Rahul Nandkishore, Leonid S Levitov, and Andrey V Chubukov, “Chiral superconductivity from repulsive interactions in doped graphene,” Nature Physics 8, 158–163 (2012).
  * Xu  _et al._ (2016) Xiao Yan Xu, Stefan Wessel, and Zi Yang Meng, “Competing pairing channels in the doped honeycomb lattice hubbard model,” Phys. Rev. B 94, 115105 (2016).
  * Faye  _et al._ (2015) J. P. L. Faye, P. Sahebsara, and D. Sénéchal, “Chiral triplet superconductivity on the graphene lattice,” Phys. Rev. B 92, 085121 (2015).
  * Uchoa and Castro Neto (2007) Bruno Uchoa and A. H. Castro Neto, “Superconducting states of pure and doped graphene,” Phys. Rev. Lett. 98, 146801 (2007).
  * Gu  _et al._ (2020) Zheng-Cheng Gu, Hong-Chen Jiang, and G. Baskaran, “Emergence of $p+ip$ superconductivity in two-dimensional doped dirac systems,” Phys. Rev. B 101, 205147 (2020).
  * Anderson (1987) Philip W Anderson, “The resonating valence bond state in la2cuo4 and superconductivity,” science 235, 1196–1198 (1987).
  * Lee  _et al._ (2006) Patrick A. Lee, Naoto Nagaosa, and Xiao-Gang Wen, “Doping a mott insulator: Physics of high-temperature superconductivity,” Rev. Mod. Phys. 78, 17–85 (2006).
  * Auerbach (1998) Assa Auerbach,  _Interacting electrons and quantum magnetism_ (Springer Science & Business Media, 1998).
  * Jayaprakash  _et al._ (1989) C. Jayaprakash, H. R. Krishnamurthy, and Sanjoy Sarker, “Mean-field theory for the $t$-$J$ model,” Phys. Rev. B 40, 2610–2613 (1989).
  * Nagaoka (1966) Yosuke Nagaoka, “Ferromagnetism in a narrow, almost half-filled $s$ band,” Phys. Rev. 147, 392–405 (1966).
  * Yoshioka (1989) Daijiro Yoshioka, “Slave-fermion mean field theory of the hubbard model,” in  _Strong Correlation and Superconductivity_ (Springer, 1989) pp. 124–131.
  * Wen (2002) Xiao-Gang Wen, “Quantum orders and symmetric spin liquids,” Phys. Rev. B 65, 165113 (2002).
  * Wang (2010) Fa Wang, “Schwinger boson mean field theories of spin liquid states on a honeycomb lattice: Projective symmetry group analysis and critical field theory,” Phys. Rev. B 82, 024419 (2010).
  * Senthil and Fisher (2000) T. Senthil and Matthew P. A. Fisher, “${Z}_{2}$ gauge theory of electron fractionalization in strongly correlated systems,” Phys. Rev. B 62, 7850–7881 (2000).
  * Senthil and Fisher (2001) T. Senthil and Matthew P. A. Fisher, “Fractionalization in the cuprates: Detecting the topological order,” Phys. Rev. Lett. 86, 292–295 (2001).
  * Read and Sachdev (1989) N. Read and Subir Sachdev, “Valence-bond and spin-peierls ground states of low-dimensional quantum antiferromagnets,” Phys. Rev. Lett. 62, 1694–1697 (1989).
  * Nagaosa (1999) Naoto Nagaosa,  _Quantum field theory in condensed matter physics_ (Springer Science & Business Media, 1999).
  * Cheng  _et al._ (2010) Meng Cheng, Kai Sun, Victor Galitski, and S. Das Sarma, “Stable topological superconductivity in a family of two-dimensional fermion models,” Phys. Rev. B 81, 024504 (2010).
  * Chiu  _et al._ (2016) Ching-Kai Chiu, Jeffrey C. Y. Teo, Andreas P. Schnyder, and Shinsei Ryu, “Classification of topological quantum matter with symmetries,” Rev. Mod. Phys. 88, 035005 (2016).
  * Sato and Ando (2017) Masatoshi Sato and Yoichi Ando, “Topological superconductors: a review,” Reports on Progress in Physics 80, 076501 (2017).
  * Kopnin and Salomaa (1991) N. B. Kopnin and M. M. Salomaa, “Mutual friction in superfluid ${}^{3}\mathrm{He}$: Effects of bound states in the vortex core,” Phys. Rev. B 44, 9667–9677 (1991).
  * Volovik (1999) GE Volovik, “Fermion zero modes on vortices in chiral superconductors,” Journal of Experimental and Theoretical Physics Letters 70, 609–614 (1999).
  * Read and Green (2000) N. Read and Dmitry Green, “Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum hall effect,” Phys. Rev. B 61, 10267–10297 (2000).
  * Möller  _et al._ (2008) A. Möller, U. Löw, T. Taetz, M. Kriener, G. André, F. Damay, O. Heyer, M. Braden, and J. A. Mydosh, “Structural domain and finite-size effects of the antiferromagnetic ${S}=1/2$ honeycomb lattice in $\mathrm{InCu_{2/3}V_{1/3}O_{3}}$,” Phys. Rev. B 78, 024420 (2008).
  * Gu (2013) Zheng-Cheng Gu, “Efficient simulation of grassmann tensor product states,” Phys. Rev. B 88, 115139 (2013).
  * Gu  _et al._ (2010) Zheng-Cheng Gu, Frank Verstraete, and Xiao-Gang Wen, “Grassmann tensor network states and its renormalization for strongly correlated fermionic and bosonic states,” arXiv preprint arXiv:1004.2563 (2010).
  * Bruus and Flensberg (2004) Henrik Bruus and Karsten Flensberg,  _Many-body quantum theory in condensed matter physics: an introduction_ (OUP Oxford, 2004).