# Doping a Mott insulator with excitons in moiré bilayer: fractional
superfluid, neutral Fermi surface and Mott transition

Ya-Hui Zhang  Department of Physics and Astronomy, Johns Hopkins University,
Baltimore, Maryland 21218, USA

###### Abstract

In this paper we explore possible phases arising from doping neutral excitons
into a Mott insulator in the context of moiré bilayers. We consider two moiré
layers coupled together only through inter-layer repulsion and there is a
U(2)$\times$ U(2) symmetry. The densities of the two layers can be tuned to be
$n_{t}=x,n_{b}=1-x$ with $n_{t}+n_{b}=1$. $x=0$ limit is a layer polarized
Mott insulator and small $x$ regime can be reached by doping inter-layer
excitons at density $x$. Charge gap can remain finite at small $x$, as is
demonstrated experimentally in the WSe2-hBN-WSe2/WS2 system. To capture the
intertwinement of the spin and exciton degree of freedom, we propose a four-
flavor spin model. In addition to the obvious possibility of single exciton
condensation phase, we also identified more exotic physics with
fractionalization: (I) We define spin gap $\Delta_{t},\Delta_{b}$ for the two
layers respectively. As long as the spin gap at either layer is finite, single
exciton condensation is impossible and we can only have paired exciton
condensation. If both spin gaps are finite, it can be a fractional exciton
superfluid with paired exciton condensation coexisting with $Z_{2}$ spin
liquid. Numerical evidences for such a phase will be provided. (II) If the
layer polarized Mott insulator at $x=0$ is in a U(1) spin liquid with spinon
Fermi surface, the natural phase at $x>0$ hosts neutral Fermi surface formed
by fermionic excitons. There are metallic counter-flow transport and also
Friedel oscillations in layer polarization in this exotic phase. Our numerical
simulation in one dimension observes an analog of this neutral Fermi surface
phase. (III) There could be a metal-insulator transition (MIT) by tuning $x$
in the universality class of a bandwidth tuned MIT. We provide one theory for
such a continuous Mott transition and predicted a universal drag resistivity.

###### pacs:

Valid PACS appear here

##  I Introduction

Recently the idea of moiré bilayerZhang and Vishwanath (2020); Zhang  _et al._
(2021a) was introduced by us to simulate “spin” physics and realize quantum
spin liquid (QSL)Savary and Balents (2016); Zhou  _et al._ (2017); Knolle and
Moessner (2019) phases using the layer pseudospin. Similar to the quantum Hall
bilayerEisenstein (2014); Li  _et al._ (2017); Liu  _et al._ (2017), we
consider a bilayer system where two layers are coupled through inter-layer
repulsion without inter-layer tunneling (see Fig. 1). Thus the total charge of
each layer is separately conserved. If we label the charge of the two layers
as $Q_{t}$ and $Q_{b}$ respectively, we can define two conserved quantum
numbers: $Q=Q_{t}+Q_{b}$ and $P_{z}=\frac{1}{2}(Q_{t}-Q_{b})$. We can view $Q$
as the total charge, then $P_{z}$ is a pseudospin and its transport can be
conveniently probed by electrically counter-flow measurementZhang and
Vishwanath (2020). We also need the inter-layer distance $d$ to be much
smaller than the intra-layer lattice constant $a_{M}$ so that the inter-layer
repulsion is comparable to the intra-layer repulsion. To achieve that, it is
natural to build up a bilayer using moire superlattices based on grapheneCao
_et al._ (2018); Chen  _et al._ (2019a, b); Yankowitz  _et al._ (2019); Chen
_et al._ (2019a, 2020); Lu  _et al._ (2019); Cao  _et al._ (2020); Liu  _et
al._ (2020); Shen  _et al._ (2020); Polshyn  _et al._ (2020); Chen  _et al._
(2021); Sharpe  _et al._ (2019); Serlin  _et al._ (2020) or transition metal
dichalcogenides (TMD)Tang  _et al._ (2020); Regan  _et al._ (2020); Wang  _et
al._ (2020); Gu  _et al._ (2021); Zhang  _et al._ (2021b); Li  _et al._
(2021); Ghiotto  _et al._ (2021). Here we will focus on TMD moiré superlattice
in which Mott physics and anti-ferromagnetic spin coupling have been well
studied experimentally in single moiré layerTang  _et al._ (2020); Regan  _et
al._ (2020); Wang  _et al._ (2020); Li  _et al._ (2021); Ghiotto  _et al._
(2021). Previous theoretical studies have also shown that an extended spin 1/2
Hubbard model is a good description of the single moiré layerWu  _et al._
(2018, 2019); Pan  _et al._ (2020a, b); Pan and Sarma (2021, 2022); Zhang  _et
al._ (2020). Especially for TMD hetero-bilayer there is a good SU(2) spin
rotation symmetry because valley contrasting flux is negligibleWu  _et al._
(2018).

Unlike the quantum Hall bilayer, the spin within each moiré layer is still
active and we have both the layer pseudospin and the real spin. In our
previous workZhang and Vishwanath (2020); Zhang  _et al._ (2021a), we consider
the ideal limit where the real spin and the layer pseudospin together form a
SU(4) spin. Although various interesting phases including a chiral spin liquid
was predicted theoretically in the SU(4) symmetric point, we note that this
ideal limit requires the interlayer distance to be almost zero, which is in
tension with the requirement of zero inter-layer tunneling. Hence the
realization of the ideal SU(4) model in realistic systems is a challenge.
Twisted AB stacked TMD homobilayer has been experimentally realizedXu  _et
al._ (2022), where inter-layer tunneling is suppressed by spin conservation
instead of hBN barrier. However, it is not clear now whether the inter-layer
tunneling in that system is completely zero, which is required for counterflow
transport. Experimentally it is much easier to completely suppress the inter-
layer tunneling by putting a thick hexagon boron nitride (hBN) between the two
layersShimazaki  _et al._ (2020); Gu  _et al._ (2021); Zhang  _et al._
(2021b). In Ref. Gu  _et al._ , 2021; Zhang  _et al._ , 2021b, a TMD moiré
layer based on WSe2/WS2 on the bottom is separated from a moiré-less TMD
monolayer WSe2 on the top by a hBN barrier. Let us label the density of the
two layers as $n_{t},n_{b}$ per moiré unit cell. A correlated insulator with
charge gap is found for the filling $n_{t}=x,n_{b}=1-x$ in the range
$x\in(0,x_{c})$ with the critical exciton density $x_{c}$ as large as
$0.5$Zhang  _et al._ (2021b). The SU(4) symmetry is strongly broken in this
case and we have only U(2)$\times$ U(2) symmetry corresponding to the separate
conservation of charge and spin in the two layers. This moiré+monolayer setup
is easier to fabricate than the moiré+moiré system which we previously
proposedZhang and Vishwanath (2020); Zhang  _et al._ (2021a). However, it is
not clear whether interesting spin liquid phase can still be realized in this
anisotropic case with lower symmetry. The main purpose of this work is to
propose a theoretical framework to analyze the spin and exciton dynamics and
demonstrate the existence of exotic phases with fractionalization even in this
simpler system.

We first want to highlight the essential difference between the system
considered here and the traditional excitonic insulatorJérome  _et al._
(1967); Zittartz (1967); Halperin and Rice (1968); Comte and Nozieres (1982).
In the conventional case, a bosonic exciton operator is well defined on top of
a band insulator. In contrast, the vacuum of the inter-layer excitons in our
moiré bilayer is a layer polarized Mott insulator at the $n_{t}=0,n_{b}=1$
limit. Unlike a band insulator, the spin configuration of a Mott insulator is
not uniquely frozen and the spin remains as an active degree of freedom at low
energy. In the usual excitonic insulator, one can make a particle hole
transformation relative to the band insulator, so a hole operator is well
defined. Under such a particle hole transformation, the band insulator becomes
the vacuum of the exciton. Then the active degree of freedoms are just small
densities of electron and holes on top of the band insulator, which then form
excitons. In contrast, the spin configuration of a Mott insulator is not
frozen and Mott insulator can not be viewed as a vacuum. Excitons are not the
only active degree of freedom at finite $x$ and they need to interact with the
spin 1/2 moments. It is clear that we need a more sophisticated model to
capture the exciton-spin dynamics than a simple boson model usually assumed
for excitonic insulator. If one applies a strong magnetic field to polarize
the spin, then the Mott insulator becomes a band insulator and we are in the
conventional excitonic insulator. In this paper we will be mainly focused on
the zero magnetic field case where spin is not dead and more interesting
physics than single exciton condensation may arise.

We propose a four-flavor spin model to capture the low energy spin-exciton
dynamics. This four dimensional Hilbert space can be constructed by the tensor
product of two spin $1/2$ corresponding to $\vec{P}$ and $\vec{S}$. Here
$\vec{S}$ represents the real spin and $\vec{P}$ represents the layer
pseudospin. $P_{z}=\frac{1}{2}(n_{t}-n_{b})$ is the layer charge. The physical
operators can be constructed as $P_{\mu}S_{\nu}$ with $\mu,\nu=0,1,2,3$.
$P_{0}S_{0}$ is the identity operator and there are $15$ non-trivial operators
representing the neutral spin-exciton degree of freedom. They correspond to
the $15$ generators of SU(4) group. Actually at leading order the low energy
model is just a four-flavor Heisenberg spin model similar to the SU(4) spin
model studied in Ref. Zhang  _et al._ , 2021a. Now we need to include
anisotropy terms coming from finite inter-layer distance $d>0$ and imbalanced
filling $n_{t}<<n_{b}$, but the structure of the spin model remains similar to
the SU(4) model, though the symmetry is now reduced to $(SU(2)_{t}\times
SU(2)_{b}\times U(1)_{l})/Z_{4}$111$Z_{4}$ here corresponds to
$c_{i;a\sigma}\rightarrow e^{i\frac{2\pi n}{4}}c_{i;a\sigma},n=0,1,2,3$, which
is shared by the global charge transformation and the spin transformation such
as generated by $e^{i\frac{1}{2}P_{z}2\pi n}e^{iS_{t;z}2\pi n}$.. Here
$SU(2)_{t}$ is generated by $\vec{S}_{t}=(P_{z}+\frac{1}{2})\vec{S}$,
$SU(2)_{b}$ is generated by $\vec{S}_{b}=(\frac{1}{2}-P_{z})\vec{S}$ and
$U(1)_{l}$ is generated by $P_{z}$. Physically $\vec{S}_{t},\vec{S}_{b}$ are
the spin operators of the two layers respectively and $U(1)_{l}$ corresponds
to the layer polarization or a dipole charge.

This spin model can be naturally derived from a four-flavor Hubbard model
suitable for moiré+moiré case. In the strong interaction limit, the system is
in a Mott insulator at total filling $\nu_{T}=n_{t}+n_{b}=1$. But even for the
moiré+monolayer case, we argue that the same model works for the small x limit
because the top layer can inherit the same moiré superlattice from the bottom
layer through inter-layer repulsion. Actually at leading order, the effective
exciton-spin model below the charge gap is constrained to be in the form we
propose by the $U(2)\times U(2)$ symmetry of the system. We note a simple
Hubbard model is not sufficientPan and Sarma (2022) to capture the long range
Coulom interaction. However, at the filling $\nu_{T}=1$ below the charge gap,
the effect of longer range interaction is to renormalize the values of the
parameters in our four-flavor spin model. In this paper we will treat these
couplings as phenomenological parameters and explore the phase diagram.
Derivation of these parameters for a specific realistic material will be left
to future work. We also emphasize that our model only works at energy scale
well below the charge gap. When the charge gap is small, additional ring
exchange terms beyond bilinear terms are needed. When the charge gap
approaches zero, a bosonic model is not sufficient anymore and one needs a
fermionic model to also capture the possibility of electron-hole gas. We also
do not include disorder in our study, which may lead to Anderson localization
at small $x$Ahn and Sarma (2022). Future experiments are needed to address the
question whether the insulator seen in Ref. Gu  _et al._ , 2021; Zhang  _et
al._ , 2021b is actually an Anderson insulator of electrons or excitons.

We study the four-flavor spin model on triangular lattice using either
Schwinger boson mean field or Abrikosov fermion mean field theory. Deep inside
the Mott insulator, Schwinger boson theory is convenient to capture the
$120^{\circ}$ magnetic order of $\vec{S}_{b}$ in the layer polarized limit
with $x=0$. If the $120^{\circ}$ order at the bottom layer is strong, the
physics is reduced to a spin $1/2$ boson model and we obtain either an exciton
condensation phase with $\vec{S}_{t}$ ordered or a paired exciton condensation
phase with $\vec{S}_{t}$ gapped. In paired exciton condensation phase, only a
Cooper pair of excitons (bi-exciton) is condensed and the single exciton is
gapped. In a more interesting case, we find the $120^{\circ}$ order in the
bottom layer can also be depleted by reducing the density $n_{b}$, after which
we get a paired exciton condensation phase coexisting with $Z_{2}$ spin
liquid. We will provide numerical evidences for spin gaps through density
matrix renormalization group (DMRG) simulationWhite (1992, 1993); McCulloch
(2008). We actually find that there are two different $Z_{2}$ spin liquids
with $0$ or $\pi$ flux for the Schwinger bosons. Both $0$ and $\pi$ flux
$Z_{2}$ spin liquids have been proposed before for a single layer spin $1/2$
model on triangular latticeWang and Vishwanath (2006). But they are usually
unstable to Schwinger boson condensation because the boson density is too
large. Here because we can reduce the density of the boson in the bottom layer
to $n_{b}=1-x$, it is possible to stabilize the $Z_{2}$ topological order. The
$Z_{2}$ fractionalization coexists with the paired exciton condensation, and
thus the phases can be called as fractional superfluids. The possibility of
$Z_{2}$ spin liquids in this system opens the exciting future direction to
search for superconductivity upon doping the bilayer Mott insulator, following
Anderson’s resonating-valence-bond (RVB) theoryAnderson (1987).

For the weak Mott regime, magnetic order may not exist even in the layer
polarized limitLi  _et al._ (2021). In the single moiré layer based on
MoTe2/WSe2, the weak Mott insulator was shown to be more consistent with a
U(1) spin liquid phase with spinon Fermi surfaceLi  _et al._ (2021) given the
large spin susceptibility at lowest temperature. Smoking gun evidence for such
a neutral fermi surface has remained elusive because of the lack of the probe
of neutral excitations. We propose to add a monolayer TMD separated from
MoTe2/WSe2 by a hBN similar to the setup of Ref. Gu  _et al._ , 2021; Zhang
_et al._ , 2021b. If the layer polarized limit is indeed a spinon Fermi
surface state, then we should use the four-flavor Abrikosov fermion
$f_{i;a\sigma}$ to attack the four-flavor spin model. In the $x=0$ limit, the
band for $f_{i;t\sigma}$ is empty while $f_{i;b\sigma}$ forms the spinon Fermi
surface. Upon increasing $x$, there is a one component to two components
Lifshitz transition and we reach a U(1) spin liquid phase where both layers
host a spinon Fermi surface. Physically $f_{t;\sigma}$ can be interpreted as a
fermionic exciton formed by binding electron $c_{i;t\sigma}^{\dagger}$ with a
holon in the bottom layer. This phase will support metallic counter-flow
transport $\rho_{\text{counterflow}}(T)\sim T^{\alpha}$ with some non-Fermi-
Liquid (NFL) exponent $\alpha$ due to gauge fluctuation. The exact value of
$\alpha$ is still not well established theoreticallyLee and Nagaosa (1992);
Maslov  _et al._ (2011); Hartnoll  _et al._ (2014); Lee (2018, 2021) and we
hope the future experiment in moiré bilayer could provide experimental
constraint on the theory of spinon Fermi surface phase in 2+1d. The neutral
Fermi surface could also be revealed by Friedel oscillation in terms of the
layer polarization, which may be measured electrically.

In the experiments of Ref. Gu  _et al._ , 2021; Zhang  _et al._ , 2021b, the
charge gap closes beyond a critical doping $x_{c}$ for the exciton density(
see Fig. 1(b)). It is thus interesting to study the metal insulator transition
(MIT) tuned by $x$. We will argue that such a transition is in the class of
bandwidth tuned MIT. The easiest way to understand the transition is from
tuning the inter-layer distance $d$. When $d$ is large, the two layers
decouple and must be in metallic phases with Fermi surfaces. At fixed x, when
we decrease $d$ to be smaller than a critical distance $d_{c}$, the system
becomes insulating because of the inter-layer repulsion
$U^{\prime}n_{t}n_{b}$, which is analogous to the Hubbard U term. Therefore
the distance tuned MIT is similar to the MIT tuned by $U/t$ in the simple
spinful Hubbard model. It is easy to imagine that $d_{c}$ decreases with $x$
and then one can drives the same MIT by tuning $x$ at fixed $d$. We also
provide a theory of this MIT by generalizing the theory of Ref. Senthil, 2008.
We assume the Mott insulator side hosts spinon Fermi surfaces and the MIT is
driven by the condensation of the slave boson. The theory is essentially the
same as in the single layer case, but in the moiré bilayer case one can drive
currents in the two layers separately. Hence we can get a $2\times 2$ tensor
for the resistivity $\rho_{xx}$. In the slave boson theory of MIT, there is a
universal drag resistivity $\rho_{xx;12}$ at order $h/e^{2}$, which may be
easier to measure than the jump of the residual resistivity proposed in the
single layer caseSenthil (2008). Drag resistivity essentially measures the
pseudospin-charge separation and thus may also be very useful to characterize
the metallic state upon doping the bilayer Mott insulator by tuning
$n_{t}+n_{b}$ away from $1$.

The rest of the paper is structured in the following way. In Sec. II we
introduce the four-flavor spin model and its derivation. We clarify the
relation ship between exciton and layer pseudospin language in Sec. III. In
Sec. IV we introduce the framework to deal with the model using either
schwinger boson or Abrikosov fermion parton construction. In Sec. V we discuss
various possible exciton superfluid phases coexisting with magnetic order or
Z2 spin liquids. This includes a spin polarized BEC of excitons, a paired
superfluid of excitons with a spin gap in one of the two layers. Most
interestingly we show that it is possible to deplete the magnetic order of the
Mott insulator by doping excitons, leading to a Z2 spin liquid coexisting with
exciton superfluid. A phase diagram in $x-J^{\prime}_{b}$ space will be
provided, with $x$ as the exciton density and $J^{\prime}_{b}$ as the next
nearest neighbor coupling in the Mott layer. The phase diagram hosts two
different Z2 spin liquids corresponding to zero and $\pi$ flux ansatz of
Schwinger boson, which connect to the $120^{\circ}$ ordered phase and stripe
ordered antiferromagnetic phase respectively. In Sec.VI we discuss the
possibility of fermionic exciton and neutral Fermi surface and in Sec. VII we
show a critical theory of continuous mott insulator transition tuned by the
exciton density $x$.

Figure 1: Illustration of a moiré bilayer. Two layers are separated by an
insulating barrier and each layer has the same superlattice with lattice
constant $a_{M}$.

##  II Anisotropic SU(4) spin model

We consider a moiré bilayer with two moiré superlattice layers separated by an
insulating barrier, as shown in Fig. 1. The two moiré superlattices are
assumed to be aligned with each other and share the same triangular lattice.
The inter-layer tunneling is forbidden by insulating barrier. In the following
we are going to derive a four-flavor exciton-spin model for filling
$\nu_{T}=1$ at the energy scale well below the charge gap. We will justify
this model using a simple four-flavor Hubbard model. But we note at leading
order the form of this four-flavor spin model is constrained by the
U(2)$\times$ U(2) symmetry of the system and we expect the same model even for
more complicated models, for example with longer range Coulomb interaction.

We consider the following four-flavor Hubbard modelZhang  _et al._ (2021a)

$\displaystyle H$
$\displaystyle=-\sum_{a=t,b}t_{a}\sum_{ij}\sum_{\sigma=\uparrow,\downarrow}c^{\dagger}_{i;a\sigma}c_{j;a\sigma}-\frac{D}{2}\sum_{i}(n_{i;t}-n_{i;b})$
$\displaystyle+\frac{1}{2}\sum_{i}\sum_{a}U_{a}n_{i;a}(n_{i;a}-1)+U^{\prime}\sum_{i}n_{i;t}n_{i;b}$
$\displaystyle+\sum_{a}V_{a}\sum_{\langle
ij\rangle}n_{i;a}n_{j;a}+V^{\prime}\sum_{\langle
ij\rangle}(n_{i;t}n_{j;b}+n_{i;b}n_{j;t})$ (1)

where $a=t,b$ labels the layer index and $\sigma=\uparrow,\downarrow$ labels
the spin. In TMD, spin is locked to valley due to a strong ising spin-orbit
coupling. Here we focus on two narrow moiré bands from the two valleys on top
of the original valence band. Within this two flavor space, there is a good
SU(2) spin rotation symmetryWu  _et al._ (2018). .
$n_{i;a}=\sum_{\sigma}c^{\dagger}_{i;a\sigma}c_{i;a\sigma}$ is the density at
layer $a$ on site $i$. $U_{a}$ and $U^{\prime}$ are intra-layer and inter-
layer on-site Hubbard repulsion. $V_{a}$ and $V^{\prime}$ are intra-layer and
inter-layer nearest neighbor repulsion. $D$ is a displacement field. Note here
we assume that the two layers can have different hopping and interaction, so
one layer can be made less correlated, for example with $U_{b}>>U_{t}$. The
model has exact U(4) symmetry if $t_{t}=t_{b}$, $U_{t}=U_{b}=U^{\prime}$,
$V_{t}=V_{b}=V^{\prime}$ and $D=0$. We will work in the case with significant
anisotropic terms breaking the U(4) symmetry. Even in the anisotropic case, we
still have a large symmetry $U(2)\times U(2)$, where each U(2) corresponds to
the charge and spin conservation at one of the two layers. Note the spin
$\vec{S}_{i;a}=\frac{1}{2}c^{\dagger}_{i;a\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}c_{i;a\sigma^{\prime}}$
at each layer $a=t,b$ is separately conserved, the same as the charge
$n_{i;a}=\sum_{\sigma}c^{\dagger}_{i;a\sigma}c_{i;a\sigma}$. In the above we
only included on-site and nearest neighbor repulsion. For filling $\nu_{T}=1$,
long ranger Coulomb interaction should only renormalizes the values of the
parameters in the low energy spin model we are going to study.

If we assume $U_{t},U_{b},U^{\prime}$ are all much larger than the hopping
$t_{t},t_{b}$, we can have a Mott insulator at total filling $n_{t}+n_{b}=1$,
as illustrated in Fig. 2(a). This is true in WSe2-WS2-WSe2 system or twisted
AB stacked WSe2 bilayerZhang  _et al._ (2021a). However, we will also consider
the case that the top layer is less correlated and has a larger $t_{t}$. This
is true in the moiré+monolayer system such as in the recent experimentsGu  _et
al._ (2021); Zhang  _et al._ (2021b). In this case we expect that the Mott gap
survives only when the density of the top layer is smaller than a critical
value, as illustrated in Fig. 2. In this paper we mainly focus on the case
with $n_{t}=x,n_{b}=1-x$ with small $x$, so there is a charge gap and the
difference between the moiré+moiré and moiré+monolayer system does not matter
too much.

Figure 2: Illustration of the bilayer Mott insulator at total filling
$n_{t}+n_{b}=1$. The red color labels the region with a Mott charge gap. The
blue color labels the region without charge gap. (a) If
$U_{t},U_{b},U^{\prime}>>t_{t},t_{b}$, the whole line of $n_{t}+n_{b}=1$ is in
a Mott insulator. This applies to moiré bilayer formed in WSe2-WS2-WSe2 or
twisted AB stacked WSe2 bilayer where each layer feels a strong moiré
superlattice potential. (b) If the top layer is less correlated with a large
$t_{t}$, then the charge gap is finite only in the region with $n_{t}$ smaller
than a critical value $x_{c}$. This is true in moiré+monolayer setup such as
in the WSe2-hBN-WSe2/WS2 systemGu  _et al._ (2021); Zhang  _et al._ (2021b).

Within the region with a finite charge gap, the low energy physics is governed
by an anisotropic SU(4) model obtained from the standard $t/U$ expansion:

$\displaystyle H_{S}$ $\displaystyle=J_{t}\sum_{\langle
ij\rangle}\vec{S}_{t}(i)\cdot\vec{S}_{t}(j)+J_{b}\sum_{\langle
ij\rangle}\vec{S}_{b}(i)\cdot\vec{S}_{b}(j)+J^{\prime}_{b}\sum_{\langle\langle
ij\rangle\rangle}\vec{S}_{b}(i)\cdot\vec{S}_{b}(j)+\frac{1}{2}J_{pz}\sum_{\langle
ij\rangle}P_{z}(i)P_{z}(j)$ $\displaystyle+\frac{1}{2}J_{p}\sum_{\langle
ij\rangle}\big{(}P_{x}(i)P_{x}(j)+P_{y}(i)P_{y}(j)\big{)}\left(4\vec{S}(i)\cdot\vec{S}(j)+S_{0}(i)S_{0}(j)\right)$
(2)

where $J_{t}=\frac{4t_{t}^{2}}{U_{t}}$, $J_{b}=\frac{4t_{b}^{2}}{U_{b}}$,
$J_{p}=\frac{4t_{t}t_{b}}{U^{\prime}}$ and $J_{pz}=2\delta
V-\frac{1}{2}(J_{t}+J_{b})+(\frac{4t_{t}^{2}}{U^{\prime}}+\frac{4t_{b}^{2}}{U^{\prime}})$,
where $\delta V=V_{t}+V_{b}-2V^{\prime}$. In the above we ignored the constant
term and the term linear to $P_{z}(i)$, which is just a chemical potential. We
also include a term $J^{\prime}_{b}$ coming from the next nearest neighbor
hopping in the bottom layer. Such a term will be included in the discussion of
Sec. V.3, but can be ignored for now.

In the $SU(4)$ symmetric limit with $J_{t}=J_{b}=J_{p}=J_{pz}=J$, we recover
the SU(4) spin model

$\displaystyle
H_{S}=\frac{J}{8}\left(4\vec{P}(i)\cdot\vec{P}(j)+P_{0}(i)P_{0}(j)\right)\left(4\vec{S}(i)\cdot\vec{S}(j)+S_{0}(i)S_{0}(j)\right)$
(3)

In the above we keep up to the order of $t^{2}/U$, but higher order terms can
be generated in the standard way. In the case where the top layer is not
correlated, we have large value of $\frac{t_{t}}{U_{t}}$ and
$\frac{t_{t}}{U^{\prime}}$, so $J_{t}$ can be significantly larger than
$J_{b}$ and $J_{p}$. We may also worry that the $t/U$ expansion fails because
$\frac{t_{t}}{U_{t}}$ or $\frac{t_{t}}{U^{\prime}}$ are too large. However,
the term $\frac{t_{t}^{2}}{U_{t}}$ exists only when there are two nearby
particles in the top layer and is associated with a factor
$n_{t}(i)n_{t}(j)\sim x^{2}$. Similarly the terms
$\frac{t_{t}t_{b}}{U^{\prime}}$ and $\frac{t_{t}^{2}}{U^{\prime}}$ are
associated with a factor $n_{t}(i)n_{b}(j)\sim x$ in the $x<<1$ limit.
Therefore the terms including $t_{t}$ are suppressed at the dilute exciton
limit and the $t/U$ expansion is valid in the $x<<1$ limit. When we increase
the exciton density, the expansion will fail eventually, which indicates a
Mott transition with the charge gap closes. This Mott transition will be
discussed in Sec. VII.

###  II.1 Moiré+monolayer case

In the $x<<1$ limit, we can have a Mott insulator even if the top layer is not
strongly correlated. Actually the top layer even does not need to have a moiré
superlattice. The point is that the doped electrons in the top layer inherit
the moiré potential of the bottom layer and are trapped in the triangular
moiré superlattice sites due to the attraction from the hole in the bottom
layer. The only assumption we need to make is that the exciton is tightly
bound with the 2D position of the electron-hole pair vertically aligned. At
low energy the electron in the top layer can not hop away from the hole in the
bottom layer, which will break the exciton and costs energy. With this
assumption, we only need to keep bosonic degree of freedom in the low energy.
As we will argue below, the low energy physics is captured by a similar spin-
layer four-flavor model as in Eq. 2.

When the top layer is moiré-less, in principle we need to keep
$O\big{(}(\frac{a_{M}}{a})^{2}\big{)}$ states for the top layer within each
moiré unit cell of the bottom layer. But because of the inter-layer coulomb
interaction, it is reasonable to imagine that there are only finite number of
states needed to be kept. For example, if there is a hole at the moiré site
$i$ in the bottom layer, then the relevant state for the top layer is state
$c^{\dagger}_{i;t\sigma}\ket{0}=\int
d^{2}\mathbf{x}\Phi(\mathbf{x}-\mathbf{R_{i}})f^{\dagger}_{\sigma}(\mathbf{x})\ket{0}$.
Here $f^{\dagger}_{\sigma}(\mathbf{x})$ is the microscopic electron creation
operator on the top layer, which lives in the original lattice with small
lattice constant $a$. Note we still have two flavors
$\sigma=\uparrow,\downarrow$ coming from the spin-valley locking.
$\Phi(\mathbf{x}-\mathbf{R_{i}})$ is the wavefunction for the lowest exciton
state when the hole is localized at moiré site $i$. $c^{\dagger}_{i;t\sigma}$
can be viewed as a coarse graining operator for the top layer living on the
moiré site $i$. At low energy, it is reasonable to assume that we only need to
keep the exciton created by $c^{\dagger}_{i;t\sigma}c_{i;b\sigma^{\prime}}$.
Note here $c_{i;b\sigma^{\prime}}$ is the annihilation operator of electron at
the moiré site $i$ of the bottom layer. In low energy we can safely ignored
other states in the top layer and keep only four states at each site:
$\ket{1}_{i}=c^{\dagger}_{i;t\uparrow}\ket{0}$,
$\ket{2}=c^{\dagger}_{i;t\downarrow}\ket{0}$,
$\ket{3}_{i}=c^{\dagger}_{i;b\uparrow}\ket{0}$,
$\ket{4}=c^{\dagger}_{i;b\downarrow}\ket{0}$. This is exactly the same Hilbert
space defined in Eq. 2. The symmetry of this system is still the same as the
moiré bilayer case, with spin of the two layers separately conserved. The
pseudospin $P_{z}$ is still conserved. Within this restricted Hilbert space,
the low energy effective Hamiltonian is constrained by symmetry to be in the
same form as in Eq. 2. Thus we conclude that Eq. 2 is also the correct
effective model for the moiré+monolayer case.

The parameters can of course be different from the case with both layers to be
moiré layers. We will treat the couplings $J_{t},J_{b},J_{p},J_{pz}$ as
phenomenological parameters. $J_{b}$ is the spin-spin coupling of the bottom
layer and is fixed by the physics of the moiré layer itself and we know
$J_{b}>0$. $J_{t}$ is the spin-spin coupling of the to layer. $J_{p}$ is
basically the hopping of the exciton. $J_{pz}$ is the nearest neighbor
repulsion between two excitons.

Figure 3: Illustration of a toy model for the moiré-monolayer case. The bottom
layer has a moiré superlattice labeled by the solid line. The top layer is
just a monolayer TMD without moiré superlattice. However, it feels interaction
from the bottom moiré superlattice. We assume that only three orbitals in the
top layer need to be kept per moiré unit cell. They live on the A,B,C
sublattice, which together form a triangular lattice with lattice constant
$\frac{a_{M}}{\sqrt{3}}$ denoted by the dashed line. Only the A sublattice is
on top of the moiré site in the bottom layer. If we fix $n_{b}=1$ and increase
$n_{t}$, the doped electron in the top layer will be repelled to B,C
sublattice and form a honeycomb lattice. However, if we increase $n_{t}$ while
keeping $n_{b}+n_{t}=1$, the doped electron in the top layer will be bound to
the hole in the bottom layer at the moiré site $A$. In this case, states at
$B$ and $C$ sublattice are penalized, but they can assist super-exchange
through virtual hopping.

Next we offer a simple toy model to show that we still expect $J_{t}>0$ and
$J_{p}>0$. As argued before, we only need to keep one coarse graining orbital
on top of the moiré site for the top layer. This is the A sublattice in Fig.
3. But there are some other orbitals which may be important in the virtual
process to generate the $J_{p}$ and $J_{t}$ terms. For example, it is also
important to keep the B,C sublattice in Fig. 3. With these three orbitals, we
can write down an effective model for the moiré-monolayer system:

$\displaystyle H$
$\displaystyle=-t_{b}\sum_{ij}c^{\dagger}_{i;b\sigma}c_{j;b\sigma}+\frac{U_{b}}{2}\sum_{i}n_{i;b}^{2}$
$\displaystyle-
t_{t}\sum_{\tilde{i}\tilde{j}}c^{\dagger}_{\tilde{i};t\sigma}c_{\tilde{j};t\sigma}+\frac{U_{t}}{2}\sum_{\tilde{i}}n_{\tilde{i};t}^{2}+V_{t}\sum_{\langle\tilde{i}\tilde{j}\rangle}n_{\tilde{i};t}n_{\tilde{j};t}$
$\displaystyle+U^{\prime}\sum_{i\in A}n_{i;t}n_{i;b}+V^{\prime}\sum_{\langle
i\tilde{j}\rangle}n_{i;b}n_{\tilde{j};t}$ (4)

where $i,j\in A$ live on the bottom layer in the big lattice formed by A
sites. $\tilde{i},\tilde{j}$ live on the small lattice formed by A,B,C sites
together for the top layer. Note $V_{t}$ is between the nearest neighbor of
the smaller triangular lattice denoted by the dashed line. $V^{\prime}$ is the
inter-layer repulsion between the moiré site $i$ in the bottom layer with the
nearest neighbor B or C site in the top layer. $U^{\prime}$ is the on-site
inter-layer repulsion. Generically we expect $U^{\prime}>V^{\prime}$,
$V_{t}>V^{\prime}$. In the above we ignored the intra-layer and inter-layer
interaction between A sites.

Let us fix our density to be $n_{t}=x,n_{b}=1-x$. In the low energy, we only
need to keep the states with one electron in either layer at site $A$. We can
not occupy both layers of the same A site, which is forbidden by $U^{\prime}$.
The state at B and C is also penalized. For example, if we move one electron
in the top layer from site A to site B, there will be energy cost of
$2V^{\prime}$, assuming nearest neighbor A sites of B are occupied in the
bottom layer (otherwise the energy cost is $2V_{t}$ or $V_{t}+V^{\prime}$,
which is even larger). With this constraint, the low energy model must be a
spin-layer four-flavor model only on A site.

Next we try to derive the parameters in Eq. 2. First, the $J_{b}$ term is not
changed and should still be $J_{b}=\frac{4t_{b}^{2}}{U_{b}}$. For $J_{t}$, we
now need to consider a four-step process through the intermediate B or C
sites. This gives a term
$\frac{4t_{t}^{4}}{U_{t}(V_{t}+V^{\prime})^{2}}(4\vec{S}_{t}(i)\cdot\vec{S}_{t}(j)-n_{t}(i)n_{t}(j))+\frac{2t_{t}^{4}}{V^{\prime}(V_{t}+V^{\prime})^{2}}(4\vec{S}_{t}(i)\cdot\vec{S}_{t}(j)+n_{t}(i)n_{t}(j))$.
In total we get:
$J_{t}=\frac{4t_{t}^{4}}{U_{t}(V_{t}+V^{\prime})^{2}}+\frac{2t_{t}^{4}}{V^{\prime}(V_{t}+V^{\prime})^{2}}$.
The $n_{t}(i)n_{t}(j)$ term will modify the $J_{pz}$ term. The most important
term is the hopping term $J_{p}$. This is generated by a three-step process:
electron in the top layer hops from site i to site j through an intermediate B
or C site; the electron in the bottom layer hops from site j to site i. We
find that
$J_{p}=\frac{4t_{t}^{2}t_{b}}{U^{\prime}V^{\prime}}+\frac{t_{t}^{2}t_{b}}{V^{\prime
2}}$. As for the repulsion term $J_{pz}$, there will be contributions from the
super-exchange. But the dominant term is still from the anisotropy part of the
nearest neighbor repulsion $\delta V$, which is ignored in Eq. 4. Thus we can
just use $J_{pz}=2\delta V$. For the moiré-monolayer system, the asymmetry of
the two layers are larger. Hence we expect a larger $J_{pz}$ term compared to
the moiré bilayer case. However, when the exciton density is low, the
repulsion term should not be too important.

Based on the toy model, we can show that the low energy Hamiltonian is still
in the form of Eq. 2 with $J_{p}>0,J_{t}>0,J_{pz}>0$. In this simple model we
only keep three states for the top layer. Strictly speaking, there are
$O((\frac{a_{M}}{a})^{2})\sim O(1000)$ orbitals within each moiré unit cell.
Therefore we should keep more orbitals at the sublattice A,B,C, but they will
have higher on-site energy $\Delta$ compared to those we are keeping now. We
can also get super-exchange coupling by hopping to these high energy orbitals
in the virtual process, but this does not change the physics qualitatively.

One additional implication of the toy model in Eq. 4 is the following: if we
fix $n_{b}=1$ at each site and only dope the top layer, then the electron in
the top layer will hop in the honeycomb lattice formed by B and C sublattice,
while the state at site A is pushed to higher energy by the $U^{\prime}$ term.
This may be one interesting way to realize a honeycomb lattice Hubbard model,
which we leave to future analysis.

##  III Dictionary between exciton and layer-pseudospin language

Conventionally the physics in coulomb coupled double layer is thought to be
associated with the inter-layer excitons. One may wonder why we use the
notation “layer pseudo-spin” instead of exciton in this paper. In this section
we want to show that the layer pseudo-spin is mathematically equivalent to
inter-layer exciton in the spinless case and is more convenient in the spinful
case.

###  III.1 Spinless case

We consider the simple case with spinless electron first. With a strong Zeeman
field, the lattice Hubbard model in Eq. 1 reduces to

$\displaystyle H$
$\displaystyle=-\sum_{a=t,b}t_{a}\sum_{ij}c^{\dagger}_{i;a}c_{j;a}-D\sum_{i}P_{z}(i)+U^{\prime}\sum_{i}n_{i;t}n_{i;b}$
$\displaystyle+V\sum_{\langle ij\rangle}n_{i}n_{j}+\delta V\sum_{\langle
ij\rangle}P_{z}(i)P_{z}(j)+\frac{V_{t}-V_{b}}{2}\sum_{ij}n_{i}P_{z}(j)$ (5)

where $V=\frac{V_{t}+V_{b}}{2}+V^{\prime}$ and $\delta
V=V_{t}+V_{b}-2V^{\prime}$. $P_{z}(i)=\frac{1}{2}(n_{t}(i)-n_{b}(i))$. In the
above we absorbed term linear to $P_{z}$ into the D term.

If we set $t_{a}=t_{b}=t$, $V_{t}=V_{b}=V^{\prime}$, we can see that the above
model just reduces to the standard spin 1/2 Hubbard model if we interpret
$a=t,b$ as the spin index. The $D$ term is basically a Zeeman field in spin
language. If we fix our density to be $n_{t}=x,n_{b}=1-x$, then the system is
in an insulating state when $U^{\prime}/t>>1$. The insulator can be called
inter-layer excitonic insulator because electron and hole in the two layers
are bound. Alternatively in the Hubbard model language the insulator is
clearly a Mott insulator, with a possible partial spin polarization
$P_{x}=x-\frac{1}{2}$. If we decrease $U^{\prime}/t$, there can be a metal-
insulator transition. In exciton language, the metal side can be viewed as an
electron-hole liquid and the Mott transition is associated with the exciton
dissociation. The situation is a little different if $t_{t}>t_{b}$. In this
case the Mott transition may happen if increasing $x$. Nevertheless, the
insulating side is still captured by a spin $1/2$ model, though with
anisotropy term breaking the SU(2) spin rotation down to $U(1)$.

Deep inside the insulating phase, the spin model Eq. 2 reduces to

$H_{S}=J_{p}\sum_{\langle
ij\rangle}\big{(}P_{x}(i)P_{x}(j)+P_{y}(i)P_{y}(j)\big{)}+J_{z}\sum_{\langle
ij\rangle}P_{z}(i)P_{z}(j)$ (6)

where $J_{z}=\frac{1}{2}J_{pz}+\frac{1}{4}(J_{t}+J_{b})$. Note in the above we
rewrite every operator in the form $P_{\mu}S_{\nu}$ and then replace
$S_{z}=\frac{1}{2},S_{0}=1,S_{x,y}=0$. We have also ignored the constant term
and the term linear to $P_{z}$.

Eq. 6 is just a spin $1/2$ model with only XY symmetry if we view $\vec{P}$ as
a spin $1/2$ degree of freedom. It is well-known that such a model can be
mapped to a hard-core bosonic model:

$H_{b}=t_{b}\sum_{\langle ij\rangle}(b^{\dagger}_{i}b_{j}+h.c.)+V\sum_{\langle
ij\rangle}n^{b}_{i}n^{b}_{j}$ (7)

where $t_{b}=\frac{1}{2}J_{p}$ and $V=J_{z}$. Again we ignored the constant
term and the term linear to $n=P_{z}+\frac{1}{2}$. We have used the mapping
$P^{\dagger}_{i}\rightarrow b_{i}^{\dagger},P^{-}_{i}\rightarrow
b_{i},\frac{1}{2}+P_{z}\rightarrow n^{b}_{i}$. In the hard-core boson
language, we label $\ket{P_{z}=-\frac{1}{2}}$ as $\ket{0}$ with $n^{b}=0$. We
label $\ket{P_{z}=\frac{1}{2}}$ as $\ket{1}$ with $n^{b}=1$. We list a
dictionary between the layer pseudo-spin and the exciton language in Table. 1.
We also list several possible spin phases as examples. In the spinless case
one can use the exciton language and write down the low energy Hamiltonian as
a boson model, which is mathematically equivalent to certain spin $1/2$ model.
On triangular lattice, the ground state is known to be in a $120^{\circ}$
order, or an exciton supersolid in the exciton language.

| pseudo-spin language | exciton language
---|---|---
$c^{\dagger}_{t}\ket{0}$ | $\ket{P_{z}=\frac{1}{2}}$ | $b^{\dagger}\ket{0}$
$c^{\dagger}_{b}\ket{0}$ | $\ket{P_{z}=-\frac{1}{2}}$ | $\ket{0}$
$\frac{1}{2}(n_{t}-n_{b})$ | $P_{z}$ | $n_{b}-\frac{1}{2}$
$c^{\dagger}_{b}c_{t}$ | $P^{-}$ | $b$
$c^{\dagger}_{t}c_{b}$ | $P^{\dagger}$ | $b^{\dagger}$
Phase I | XY FM | exciton condensation
Phase II | 120∘ Neel order | exciton supersolid
Phase III | valence bond solid | crystallized insulator
Phase IV | chiral spin liquid | exciton FQHE
Phase V | spinon Fermi surface | exciton metal
phase VI | Fermi liquid | electron hole liquid
Counter-flow | spin conductivity | exciton conductivity
Table 1: A dictionary between the layer pseudo-spin language and the exciton
language in the spinless case. $b^{\dagger}$ creates exciton starting from the
bottom layer polarized limit. XY FM denotes the XY ferromagnetism with
$\langle S_{x}\rangle\neq 0$. FQHE denotes fractional quantum Hall state. The
120∘ order (supersolid phase) is known to be the ground state of Eq. 6 (Eq.
7). But we also list some other phases which may be stabilized by additional
ring exchange terms. Counter-flow measures the pseudo-spin transport or
equivalently the exciton transport, which can distinguish these different
phases unambiguously.

###  III.2 Spinful case

Next we discuss the case without any magnetic field and the real spin
$\vec{S}$ must be included in the low energy dynamics. In this case, one can
not capture the full dynamics with a simple boson Hubbard model. We are forced
to keep the spin operator $\vec{S}$ in the Hamiltonian. We can still use the
mapping $P^{\dagger}\rightarrow b^{\dagger}$ to write Eq. 2 as a boson-spin
model:

$\displaystyle H_{S}$ $\displaystyle=J_{t}\sum_{\langle
ij\rangle}\vec{S}_{t}(i)\cdot\vec{S}_{t}(j)+J_{b}\sum_{\langle
ij\rangle}\vec{S}_{b}(i)\cdot\vec{S}_{b}(j)$
$\displaystyle~{}~{}~{}+\frac{1}{2}J_{pz}\sum_{\langle
ij\rangle}n^{b}_{i}n^{b}_{j}$
$\displaystyle~{}~{}~{}+\frac{1}{4}J_{p}\sum_{\langle
ij\rangle}\big{(}b^{\dagger}_{i}b_{j}+b^{\dagger}_{j}b_{i}\big{)}\left(4\vec{S}(i)\cdot\vec{S}(j)+S_{0}(i)S_{0}(j)\right)$
(8)

The above Hamiltonian is equivalent to Eq. 2. However, the connection to the
SU(4) symmetric limit is not explicit. Also the model looks complicated and is
not a simple boson Hubbard model anymore. Besides, now the inter-layer exciton
creation operator $c^{\dagger}_{b\sigma}c_{t\sigma^{\prime}}$ must carry spin
index $\sigma,\sigma^{\prime}=\uparrow,\downarrow$. One can see that there are
several different flavors of exciton. Both $b^{\dagger}$ and
$b^{\dagger}\vec{S}$ are physical exciton creation operators. Therefore,
exciton is not a very convenient representation in the spinful case. In the
following we will use the $\vec{P}$ instead and take the anisotropic SU(4)
model in Eq. 2. The “exciton-spin” physics is fully captured in this model.

##  IV Exciton fractionalization: bosonic and fermionic parton theory

Figure 4: (a)A conventional exciton formed by a $S_{t}=\frac{1}{2}$ electron
in the top layer and a $S_{b}=\frac{1}{2}$ hole in the bottom layer. The
exciton carries the dipole charge $P_{z}=1$ and $S_{t}=S_{b}=\frac{1}{2}$. It
carries two spin indexes:
$\Phi^{\dagger}_{\sigma\sigma^{\prime}}=c^{\dagger}_{t\sigma}c_{b\sigma^{\prime}}$.
In the model of Eq. 2,the exciton creation operator
$\Phi^{\dagger}_{\sigma\sigma^{\prime}}$ can be identified as $P^{\dagger}$
and $P^{\dagger}\vec{S}$. (b) Fractionalization of exciton. The hole operator
in the bottom layer is divided into a holon and a spinon, which interact with
each other through a gauge field (blue wave line). Then an electron in the top
layer and a holon bind into a fractionalized exciton (surrounded by the pink
dashed square), which carries $P_{z}=\frac{1}{2}$, $S_{t}=\frac{1}{2}$ and
$S_{b}=0$. Meanwhile the spinon (surrounded by the dashed blue square) carries
$P_{z}=-\frac{1}{2}$, $S_{t}=0$ and $S_{b}=\frac{1}{2}$. The dashed pink
square ($f^{\dagger}_{t\sigma}$ or $b^{\dagger}_{t\sigma}$) and the dashed
blue square ($f_{b\sigma}$ or $b_{b\sigma}$) can be viewed as one half of the
conventional exciton and they can be either fermionic or bosonic. Both the
conventional exciton condensation phase and more exotic phase with neutral
Fermi surface can be conveniently described with these fractionalized degree
of freedom.

Usually we think the inter-layer exciton is bosonic because it is formed by
electron hole pair. In our case, the exciton creation operator
$\Phi^{\dagger}_{\sigma\sigma^{\prime}}=c^{\dagger}_{t\sigma}c_{b\sigma^{\prime}}$
carries two spin indexes and there are four of them forming a vector
representation of $SO(4)\cong(SU(2)\times SU(2))/Z_{2}$. They can be
identified as the $P^{\dagger}$ and $P^{\dagger}\vec{S}$ operators in Eq. 2. A
phase with single exciton condensation corresponds to
$\langle\Phi^{\dagger}_{\sigma\sigma^{\prime}}\rangle\neq 0$, which needs to
also break the spin rotation symmetry of both layers. This means that single
exciton condensation needs to be accompanied by magnetic orders in both
layers. Description of the various phases directly using these exciton
operators $\Phi^{\dagger}_{\sigma\sigma^{\prime}}$ is too complicated. It is
more convenient to use fractionalized operators which correspond to one half
of the conventional exciton operator and carry a spin $1/2$ in one of the two
layers.

The fractionalization of the exciton can be understood in the following way.
When we dope holes into the Mott insulator of the bottom layer, it is known
that a useful picture is through spin-charge separationLee  _et al._ (2006):
holon carries charge and neutral spinon carries a spin $1/2$. Then we can
imagine an “exciton” formed by the electron in the top layer with the holon in
the bottom layer. This is the object inside the pink square in Fig. 4(b). The
holon is a fractionalized degree of freedom, and its statistics can be either
bosonic or fermionic depending on the exact spin state of the Mott insulator.
As a result, this “exciton” formed by electron-holon pair can be either
bosonic or fermionic. As shown in Fig. 4(b), we can think that the
conventional exciton in Fig. 4(a) is divided into two parts: the pink square
represents the exciton-holon pair and the blue square represents the spinon in
the bottom layer. There are two possibilities: (I) If the spinon in the bottom
layer is fermionic and labeled as $f_{i;b\sigma}$, then the hole operator is
$c_{i;b\sigma}=\varphi_{i}f_{i;b\sigma}$ with $\varphi_{i}$ as a bosonic
holon. In this case, the pink square is labeled as
$f^{\dagger}_{i;t\sigma}=c^{\dagger}_{i;t\sigma}\varphi_{i}$. The conventional
exciton operator can be written as
$\Phi^{\dagger}_{i;\sigma\sigma^{\prime}}=f^{\dagger}_{i;t\sigma}f_{i;b\sigma^{\prime}}$.
$f^{\dagger}_{i;\sigma}=c^{\dagger}_{i;t\sigma}\varphi_{i}$ and
$f_{i;b\sigma^{\prime}}$ correspond to the pink and blue square in Fig. 4(b).
(II) The bottom layer hole operator is written as
$c_{i;b\sigma}=h_{i}b_{i;\sigma}$. $h_{i}$ is a fermionic holon and
$b_{i;\sigma}$ is a bosonic spinon. In this case we have
$\Phi^{\dagger}_{i;\sigma\sigma^{\prime}}=b^{\dagger}_{i;t\sigma}b_{i;b\sigma^{\prime}}$.
$b^{\dagger}_{i;\sigma}=c^{\dagger}_{i;t\sigma}h_{i}$ and
$b_{i;b\sigma^{\prime}}$ correspond to the pink and blue square in Fig. 4(b).

The parton operators $f_{i;a\sigma}$ are just a four-flavor generalization of
the familiar Abrikosov fermionAbrikosov (1965) parton theory of spin $1/2$
system. Similarly $b_{i;a\sigma}$ is the Schwinger boson partonArovas and
Auerbach (1988); Read and Sachdev (1991). We can just view the system in a
Mott insulator of a four-flavor Hubbard model ad defined in Eq. 1. Inside the
Mott insulator, we can just restrict ourselves to the four-flavor spin model
defined in Eq. 2. We can then simply use either a four-flavor Schwinger boson
or four-flavor Abrikosov fermion to represent the anisotropic SU(4) spin.

In the Schwinger boson theory, the neutral particle-hole operators are
represented as:

$\displaystyle\vec{S}_{i;a}$
$\displaystyle=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}}b^{\dagger}_{i;a\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}b_{i;a\sigma^{\prime}}$
$\displaystyle P^{\dagger}_{i}$
$\displaystyle=\sum_{\sigma}b^{\dagger}_{i;t\sigma}b_{i;b\sigma}$
$\displaystyle P^{-}_{i}$
$\displaystyle=\sum_{\sigma}b^{\dagger}_{i;b\sigma}b_{i;t\sigma}$
$\displaystyle P_{i;z}$
$\displaystyle=\frac{1}{2}\sum_{\sigma}(b^{\dagger}_{i;t\sigma}b_{i;t\sigma}-b^{\dagger}_{i;b\sigma}b_{i;b\sigma})$
$\displaystyle P^{\dagger}_{i}\vec{S}_{i}$
$\displaystyle=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}b^{\dagger}_{i;t\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}b_{i;b\sigma^{\prime}}$
$\displaystyle P^{-}_{i}\vec{S}_{i}$
$\displaystyle=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}b^{\dagger}_{i;b\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}b_{i;t\sigma^{\prime}}$
(9)

with the constraint
$\sum_{a=t,b}\sum_{\sigma=\uparrow,\downarrow}b^{\dagger}_{i;a\sigma}b_{i;a\sigma}=1$.
$b_{i;a\sigma}$ is a four-flavor boson labeled by layer index $a=t,b$ and spin
index $\sigma=\uparrow,\downarrow$.

In the Abrikosov fermion theory, the particle-hole operators are written as:

$\displaystyle\vec{S}_{i;a}$
$\displaystyle=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}}f^{\dagger}_{i;a\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}f_{i;a\sigma^{\prime}}$
$\displaystyle P^{\dagger}_{i}$
$\displaystyle=\sum_{\sigma}f^{\dagger}_{i;t\sigma}f_{i;b\sigma}$
$\displaystyle P^{-}_{i}$
$\displaystyle=\sum_{\sigma}f^{\dagger}_{i;b\sigma}f_{i;t\sigma}$
$\displaystyle P_{i;z}$
$\displaystyle=\frac{1}{2}\sum_{\sigma}(f^{\dagger}_{i;t\sigma}f_{i;t\sigma}-f^{\dagger}_{i;b\sigma}f_{i;b\sigma})$
$\displaystyle P^{\dagger}_{i}\vec{S}_{i}$
$\displaystyle=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}f^{\dagger}_{i;t\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}f_{i;b\sigma^{\prime}}$
$\displaystyle P^{-}_{i}\vec{S}_{i}$
$\displaystyle=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}f^{\dagger}_{i;b\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}f_{i;t\sigma^{\prime}}$
(10)

with the constraint
$\sum_{a=t,b}\sum_{\sigma=\uparrow,\downarrow}f^{\dagger}_{i;a\sigma}f_{i;a\sigma}=1$.
$f_{i;a\sigma}$ is a four-flavor boson labeled by layer index $a=t,b$ and spin
index $\sigma=\uparrow,\downarrow$.

In the above we only focus on the neutral spin and layer pseudospin sector. If
we are also interested in the dynamics of the charge, we can simply extend the
construction to a parton construction of the four-flavor electron: (I)
$c_{i;a\sigma}=h_{i}b_{i;a\sigma}$ or (II)
$c_{i;a\sigma}=\varphi_{i}f_{i;a\sigma}$. The statistics and physical quantum
numbers of the parton operators are listed in Table. 2. In either
construction, there is an emergent $U(1)$ gauge field $\mathbf{a}$
corresponding to gauge transformation: (I)$h_{i}\rightarrow
h_{i}e^{-i\alpha_{i}},b_{i;a\sigma}\rightarrow b_{i;a\sigma}e^{i\alpha_{I}}$
or
(II)$\varphi_{i}\rightarrow\varphi_{i}e^{-i\alpha_{i}},f_{i;a\sigma}\rightarrow
f_{i;a\sigma}e^{i\alpha_{I}}$.

Operator | Statistics | $Q$ | $P_{z}$ | $\vec{S}_{t}$ | $\vec{S}_{b}$ | $Q_{g}$
---|---|---|---|---|---|---
$c^{\dagger}_{t\sigma}$ | fermion | $1$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $0$
$c^{\dagger}_{b\sigma}$ | fermion | $1$ | $-\frac{1}{2}$ | $0$ | $\frac{1}{2}$ | $0$
$h$ | fermion | $-1$ | $0$ | $0$ | $0$ | $1$
$b^{\dagger}_{b\sigma}$ | boson | $0$ | $-\frac{1}{2}$ | $0$ | $\frac{1}{2}$ | $1$
$b^{\dagger}_{t\sigma}$ | boson | $0$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $1$
$c^{\dagger}_{t\sigma}h$ | boson | $0$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $1$
$\varphi$ | boson | $-1$ | $0$ | $0$ | $0$ | $1$
$f^{\dagger}_{b\sigma}$ | fermion | $0$ | $-\frac{1}{2}$ | $0$ | $\frac{1}{2}$ | $1$
$f^{\dagger}_{t\sigma}$ | fermion | $0$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $1$
$c^{\dagger}_{t\sigma}\varphi$ | fermion | $0$ | $\frac{1}{2}$ | $\frac{1}{2}$ | $0$ | $1$
Table 2: Table: quantum numbers of various operators. $c^{\dagger}_{a\sigma}$
is electron operator. There are two different parton constructions: (I)
$c_{i;a\sigma}=\varphi_{i}f_{i;a\sigma}$; (II)
$c_{i;a\sigma}=h_{i}b_{i;a\sigma}$. In both parton constructions, there is an
emergent U(1) gauge field. $Q$ is the physical charge. $P_{z}$ is the z
component of the layer pseudospin. $\vec{S}_{t}$ and $\vec{S}_{b}$ are spin of
the top and bottom layers separately. $Q_{g}$ is the gauge charge. $Q_{g}=0$
is the condition for gauge invariant operators. We note that the spinon
$b^{\dagger}_{t\sigma}$ or $f^{\dagger}_{t\sigma}$ can be identified as
electron-holon pair $c^{\dagger}_{t\sigma}h$ or
$c^{\dagger}_{t\sigma}\varphi$. We note that the assignment of the dipole
charge $P_{z}$ can be adjusted, as long as the total $P_{z}$ quantum number
for $f^{\dagger}_{t;\sigma}$ and $f_{b;\sigma}$ is $1$. Here we choose an
assignment of $1/2$+$1/2$, but one can also use $1$+$0$, which simply
corresponds to a redefinition of the internal gauge field $a_{\mu}$.

Ioffe-Larkin rule for the counterflow resistivity Before a detailed
application of either bosonic or fermionic parton theory, we first provide a
formula to calculate the counterflow resistivity. Because the two layers have
separate conserved charges, there are two U(1) probing gauge fields
$A^{t}_{\mu}$ and $A^{b}_{\mu}$, which couple to the currents $J^{t}_{\mu}$
and $J^{b}_{\mu}$ in the two layers separately in the form
$-J^{a}_{\mu}A^{a}_{\mu}$ with $a=t,b$. We can do a linear combination
$A^{c}_{\mu}=\frac{1}{2}(A^{t}_{\mu}+A^{b}_{\mu})$ and
$A^{s}_{\mu}=A^{t}_{\mu}-A^{b}_{\mu}$. Correspondingly we can define the
charge current $J^{c}_{\mu}=J^{t}_{\mu}+J^{b}_{\mu}$ and
$J^{s}_{\mu}=\frac{1}{2}(J^{t}_{\mu}-J^{b}_{\mu})$ so that the minimal
coupling is in the form $-J^{c}_{\mu}A^{c}_{\mu}-J^{s}_{\mu}A^{s}_{\mu}$.
Physically it is easier to see that the current $J^{c}_{\mu}$ is the current
associated with the total charge $Q=n^{t}+n^{b}$, while $J^{s}_{\mu}$ is the
current associated with the layer polarization
$P_{z}=\frac{1}{2}(n_{t}-n_{b})$. In our model, we assume that there is a
charge gap. Hence we always have $\vec{J}^{c}=0$. The only possible non-zero
current below the Mott gap is $\vec{J}^{s}$, which corresponds to the
configuration $\vec{J}^{t}=-\vec{J}^{b}$. Because the currents in the two
layers are opposite, this current is called counter-flow currentLiu  _et al._
(2017).

Let us use the fermionic parton as an example and derive a Ioffe-Larkin rule
for the counterflow resistivity defined as $\rho_{s;xx}=\sigma_{s;xx}^{-1}$,
with $\sigma_{s;xx}=\frac{J^{s}_{x}}{E^{s}_{x}}$.
$\vec{E}^{s}=-\vec{\partial}A^{s}_{0}+\frac{\partial\vec{A}^{s}}{\partial t}$
is the counter-flow electric field. Physically
$\vec{E}^{s}=\vec{E}^{t}-\vec{E}^{b}$, where $\vec{E}^{a}$ is the electric
field in the layer $a$. Suppose that we have fermionic partons $f_{t;\sigma}$
and $f_{b;\sigma}$ in both layers. There is a U(1) gauge field $\vec{a}$. From
Table. 2 we know $f^{\dagger}_{t;\sigma}$ couples to
$a_{\mu}+\frac{1}{2}A^{s}_{\mu}$, while $f^{\dagger}_{b;\sigma}$ couples to
$a_{\mu}-\frac{1}{2}A^{s}_{\mu}$. Then it is easy to write down a linear
response formula for each layer:

$\displaystyle J^{t}_{x}$
$\displaystyle=\sigma_{t;xx}(\frac{1}{2}E^{s}_{x}+e_{x})$ $\displaystyle
J^{b}_{x}$ $\displaystyle=\sigma_{b;xx}(-\frac{1}{2}E^{s}_{x}+e_{x})$ (11)

where $\vec{e}$ is the electric field associated with the internal gauge field
$a_{\mu}$. $\sigma_{a;xx}$ is the conductivity of the fermion $f_{a}$ in each
layer separately.

Because the charge current $\vec{J}^{c}=\vec{J}^{t}+\vec{J}^{b}=0$ below the
charge gap, we have $J^{t}_{x}=J^{s}_{x}$ and $J^{b}_{x}=-J^{s}_{x}$. Then we
get

$\displaystyle J^{s}_{x}$
$\displaystyle=\sigma_{t;xx}(\frac{1}{2}E^{s}_{x}+e_{x})$
$\displaystyle-J^{s}_{x}$
$\displaystyle=\sigma_{b;xx}(-\frac{1}{2}E^{s}_{x}+e_{x})$ (12)

From the above two equations, it is easy to obtain the Ioffe-Larkin rule:

$\rho_{s;xx}=\rho_{t;xx}+\rho_{b;xx}$ (13)

where $\rho_{a;xx}=\sigma_{a;xx}^{-1}$ is the resistivity of the parton at the
layer $a$. The above formula basically tells us that in a counter-flow
transport, the partons from the two layers form a series circuit.

The above formula applies to both fermionic and bosonic partons. One immediate
consequence of the formula is that the counter-flow transport is insulating as
long as parton in one of the layers is insulating. In the $x=0$ limit, we have
$\rho_{t}=\infty$ because $n_{t}=0$. In this limit there can not be any
counter-flow current. Once $x>0$, both layers may have finite or even zero
$\rho_{a}$, then the counter-flow resistivity can be finite or zero.

Counter-flow transport measures the response to $\vec{A}^{s}$. A similar
formula also exists for the compressibility
$\kappa_{s}=\frac{J^{s}_{0}}{A^{s}_{0}}$, which is essentially the inter-layer
polarizability $\frac{\partial P_{z}}{\partial D}$. Following the same
argument as above, the formula for $\kappa_{s}$ is:

$\kappa_{s}^{-1}=\kappa_{t}^{-1}+\kappa_{b}^{-1}$ (14)

where $\kappa_{a}$ is the compressibility of the parton at layer $a$. To have
a finite $\kappa_{s}$, the partons in both layers must be in a compressible
phase.

In the following we will apply the Schwinger boson theory or the Abrikosov
fermion theory to analyze the case where the layer polarized Mott insulator is
magnetically ordered or in a spin liquid phase with spinon Fermi surface
respectively. We need to emphasize that the usage of these fractionalized
parton operators does not necessarily leads to an exotic phase. Actually even
for a conventional symmetry breaking phase with magnetic order and simple
exciton condensation, it is convenient to start from the Schwinger boson
parton theory, as will be demonstrated in Sec. V.

##  V Inter-layer superfluids with bosonic excitons

Deep inside the Mott insulator with a large charge gap, we can restrict
ourselves to the spin-exciton model in Eq. 2 without ring-exchange terms. We
can also add next-nearest neighbor spin-spin coupling $J^{\prime}_{b}$ to tune
the spin phase in the layer polarized limit. In the limit
$J^{\prime}_{b}/J_{b}<<1$, we know the spin phase of the bottom layer is in
the $120^{\circ}$ ordered phase in the $x=0$ limit. We will see that doping
excitons can lead to various different inter-layer superfluid phases with
either magnetic order or spin gap.

To capture the magnetic order, it is convenient to use the four-flavor
Schwinger boson representation in Eq. 9. The spin model (Eq. 2) can be
rewritten as:

$\displaystyle H$ $\displaystyle=\frac{1}{4}J_{t}\sum_{\langle
ij\rangle}(b^{\dagger}_{i;t\sigma_{1}}\vec{\sigma}_{\sigma_{1}\sigma^{\prime}_{1}}b_{i;t\sigma^{\prime}_{1}})\cdot(b^{\dagger}_{j;t\sigma_{2}}\vec{\sigma}_{\sigma_{2}\sigma^{\prime}_{2}}b_{j;t\sigma^{\prime}_{2}})$
$\displaystyle+\frac{1}{4}J_{b}\sum_{\langle
ij\rangle}(b^{\dagger}_{i;b\sigma_{1}}\vec{\sigma}_{\sigma_{1}\sigma^{\prime}_{1}}b_{i;b\sigma^{\prime}_{1}})\cdot(b^{\dagger}_{j;b\sigma_{2}}\vec{\sigma}_{\sigma_{2}\sigma^{\prime}_{2}}b_{j;b\sigma^{\prime}_{2}})$
$\displaystyle+\frac{1}{2}J_{p}\sum_{\langle
ij\rangle}(b^{\dagger}_{i;t\sigma}b_{i;b\sigma^{\prime}}b^{\dagger}_{j;b\sigma^{\prime}}b_{j;t\sigma}+b^{\dagger}_{i;b\sigma}b_{i;t\sigma^{\prime}}b^{\dagger}_{j;t\sigma^{\prime}}b_{j;b\sigma})$
$\displaystyle-\frac{1}{4}J_{pz}\sum_{\langle
ij\rangle}(n_{i;t}n_{j;b}+n_{i;b}n_{j;t})$
$\displaystyle-D\sum_{i}(n_{i;t}-n_{i;b})$ (15)

where in the last term we added a constant term
$-\frac{1}{8}J_{pz}\sum_{\langle ij\rangle}(n_{i;t}+n_{i;b})(n_{j;t}+n_{j;b})$
by using the constraint $n_{i;t}+n_{i;b}=1$. The $D$ term is added if we use
grand canonical ensemble. Alternatively we can also fix the density to be
$n_{t}=x,n_{b}=1-x$ and then the $D$ term is not needed.

###  V.1 Classification of possible superfluid phases

| $\langle b_{b;\sigma}\rangle\neq 0$ |  $\langle b_{b;\sigma}\rangle=0$, $\langle\epsilon_{\sigma\sigma^{\prime}}b_{b;\sigma}b_{b;\sigma^{\prime}}\rangle\neq 0$
---|---|---
$\langle b_{t;\sigma}\rangle\neq 0$ | I: $\langle P^{\dagger}\rangle\neq 0$, $\Delta_{t}=0$, $\Delta_{b}=0$ | III: $\langle P^{\dagger}\rangle=0$, $\langle P^{\dagger}P^{\dagger}\rangle\neq 0$, $\Delta_{t}=0$, $\Delta_{b}>0$
$\langle b_{t;\sigma}\rangle=0$, $\langle\epsilon_{\sigma\sigma^{\prime}}b_{t;\sigma}b_{t;\sigma^{\prime}}\rangle\neq 0$ | II: $\langle P^{\dagger}\rangle=0$, $\langle P^{\dagger}P^{\dagger}\rangle\neq 0$, $\Delta_{t}>0$, $\Delta_{b}=0$ | IV: $\langle P^{\dagger}\rangle=0$, $\langle P^{\dagger}P^{\dagger}\rangle\neq 0$, $\Delta_{t}>0$, $\Delta_{b}>0$
Table 3: Four different types of phases described by the Schwinger boson mean
field theories. The Schwinger boson at each layer can be either single
condensed or pair condensed. If it is pair condensed, then there is a spin gap
$\Delta_{a}$ for the corresponding layer $a=t,b$. As long as one of
$\Delta_{t},\Delta_{b}$ is finite, the exciton condensation carries charge $2$
under $P_{z}$ with only $\langle P^{\dagger}P^{\dagger}\rangle\neq 0$. In this
case the condensed object is formed by two electrons and two holes. When both
$\Delta_{t}>0$ and $\Delta_{b}>0$, the paired exciton condensation coexists
with a $Z_{2}$ topological order.

We will show that the model in Eq. 15 hosts a rich phase diagram depending on
the parameters. Before going into details, we first list possible phases and
their physical properties. We have three quantum numbers: $P_{z}$, $S_{t}$,
$S_{b}$. From Eq. 15, we can obtain a mean field Hamiltonian in the form
$H_{M}=H_{t}+H_{b}$, where $H_{t},H_{b}$ are quadratic Hamiltonians for
$b_{t;\sigma}$ and $b_{b;\sigma}$ separately. Because of the conservation of
$P_{z}$, $\vec{S}_{t}$ and $\vec{S}_{b}$, terms like
$b^{\dagger}_{t;\sigma}b_{b;\sigma}$ and $b_{t;\sigma}b_{b;\sigma^{\prime}}$
are forbidden. Therefore $b_{t;\sigma}$ and $b_{b;\sigma}$ remain decoupled in
the mean field level.

The bosons at the two layers have densities $n_{t}=x$ and $n_{b}=1-x$. With
finite density, the bosons are known to form a condensation phase. But for
each layer $a$, there are two different possibilities: (1) Single boson
$b_{a;\sigma}$ is condensed, leading to magnetic ordering of the spin
$\vec{S}_{a}$. (2) Only a pair of boson is condensed:
$\langle\epsilon_{\sigma\sigma^{\prime}}b_{i;a\sigma}b_{j;a\sigma^{\prime}}\rangle\neq
0$. In this case, the spin $\vec{S}_{a}$ is gapped due to nearest neighbor
pairing of the spinons. Combining the two layers, we have four types of
superfluid phases as listed in Table. 3 and illustrated in Fig. 5. Note that
each type may contain subcategories specified by detailed symmetry realization
such as the momentum associated with the magnetic ordering. As long as spin in
one of the layer is gapped, the superfluid has a condensation of charge $2e$
under $P_{z}$, which supports paired exciton condensation with $\langle
P^{\dagger}P^{\dagger}\rangle\neq 0$ order instead of the simple exciton
condensation with $\langle P^{\dagger}\rangle\neq 0$.

Figure 5: Four possible exciton superfluid phases described by the Schwinger
boson mean field theories. Here the vertical red dashed square surrounds an
exciton. The horizontal blue solid square denotes a spin-singlet pairing term
within one layer, which causes a finite spin gap $\Delta_{a}$ in this layer
and also a gap for single exciton. (a) Single exciton condensation phase
coexisting with magnetic orders in both layers; (b) There is magnetic order in
the bottom layer, but the spin in the top layer is gapped. Only a pair of
exciton condenses. (c) Magnetic order in the top layer and spin gap in the
bottom layer. Only a pair of exciton condenses. (d) Both layers have a spin
gap. Paired exciton condensation coexists with $Z_{2}$ spin liquid. These four
phases in (a)(b)(c)(d) correspond to phases I, II, III, IV in Table. 3
respectively.

One particular case is the phase IV with both $\Delta_{t}>0$ and
$\Delta_{b}>0$. In the Schwinger boson theory, the spin gap is caused by the
paired condensation of the Schwinger boson. In this class, the emergent U(1)
gauge field is only higgsed down to $Z_{2}$. This means that there is still a
deconfined $Z_{2}$ gauge field coexisting with the inter-layer superfluid
order. The paired condensation of the Schwinger boson is the standard
description of the $Z_{2}$ spin liquid in the single layer model. Our bilayer
model provides an interesting realization of it. In contrast, if $b_{i;t}$ or
$b_{i;b}$ is condensed, the $Z_{2}$ gauge field would be completely higgsed
and we are left only with a conventional symmetry breaking order without
fractionalization.

In our case because $n_{b}=1-x$ is large, it is easy to have a magnetic order
($\langle b_{b;\sigma}\rangle\neq 0$) in the bottom layer. Then we can have
either the phase I or the phase II, which will be discussed in Sec. V.2. We
will find that the phase I is favored if $\frac{J_{t}}{J_{p}}$ is small, with
the spin in the top layer in a ferromagnetic phase due to kinetic energy. But
with a reasonably large value of $\frac{J_{t}}{J_{b}}$, the spin gap in the
top layer is finite and we get the phase II. In Sec.V.3 we will show that
increasing $x$ or adding a $J^{\prime}_{b}$ term can further destroy the
magnetic order in the bottom layer and favors the phase IV, which is a
fractional superfluid with deconfined Z2 topological order.

###  V.2 With magnetic order in the bottom layer: exciton condensation and
paired exciton condensation

We deal with the case with a robust $120^{\circ}$ order in the bottom layer
first. The magnetic order can be suppressed by $J^{\prime}_{b}/J_{b}$, which
will be discussed in the next subsection.

We can simplify our analysis by condensing the boson $b_{b}$ at momentum
$\mathbf{Q}=K,K^{\prime}$ to form the $120^{\circ}$ order. We represent the
Schwinger boson at the bottom layer as
$\chi_{i}=\begin{pmatrix}b_{i;b\uparrow}\\\ b_{i;b\downarrow}\end{pmatrix}$.
After the condensation of $\chi_{i}$ at momentum
$\mathbf{Q}=\mathbf{K},\mathbf{K^{\prime}}$, we get a three-sublattice
magnetic order. We can always rotate the spin on B and C sublattice so that
the ordered moment at each site is pointing along the $z$ direction. More
specifically, we can write the Schwinger boson as

$\displaystyle\chi_{A;i}$ $\displaystyle=\tilde{\chi}_{i}$
$\displaystyle\chi_{B;i}$
$\displaystyle=e^{-i\frac{2\pi}{3}\sigma_{y}}\tilde{\chi}_{i}$
$\displaystyle\chi_{C;i}$
$\displaystyle=e^{i\frac{2\pi}{3}\sigma_{y}}\tilde{\chi}_{i}$ (16)

with

$\tilde{\chi}_{i}=\left(\begin{array}[]{c}\sqrt{n_{b}-a^{\dagger}_{i}a_{i}}\\\
a_{i}\end{array}\right)$ (17)

where $a_{i}$ is the Holstein-Primakoff boson representing the spin wave on
top of the magnetic order. $n_{b}$ is the density of Schwinger boson at the
bottom layer.

We can do a large $n_{b}$ expansion to capture the dynamics of the goldstone
mode $a_{i}$. Substituting Eq. 17 into Eq. 15, we get

$\displaystyle H_{M}$ $\displaystyle=H_{b}+H_{t}+H_{int}$ (18)

with

$\displaystyle H_{t}$
$\displaystyle=-\frac{1}{4}J_{p}M\sum_{ij}b_{i;t\sigma}^{\dagger}b_{j;t\sigma}+J_{t}\sum_{\langle
ij\rangle}\vec{S}_{i;t}\cdot\vec{S}_{j;t}+\frac{1}{2}J_{pz}\sum_{\langle
ij\rangle}n_{i;t}n_{j;t}$ (19)

$\displaystyle H_{b}$ $\displaystyle=+\frac{1}{8}J_{b}M\sum_{\langle
ij\rangle}(a_{i}a_{j}+a_{i}^{\dagger}a_{j}^{\dagger})-\frac{3}{8}J_{b}M\sum_{\langle
ij\rangle}(a_{i}^{\dagger}a_{j}+a_{j}^{\dagger}a_{i})$
$\displaystyle+\frac{3}{2}J_{b}M\sum_{i}a_{i}^{\dagger}a_{i}+...$

and

$\displaystyle H_{int}$
$\displaystyle=-\frac{1}{4}J_{p}\sum_{ij}a_{j}^{\dagger}a_{i}b^{\dagger}_{i;t\sigma}b_{j;t\sigma}$
$\displaystyle+\frac{\sqrt{3}}{4}J_{p}\sqrt{M}\sum_{\langle
ij\rangle}\big{(}(a_{i}-a_{j}^{\dagger})b_{i;t\sigma}^{\dagger}b_{j;t\sigma}+h.c.\big{)}$
(21)

where the summation of $\langle ij\rangle$ is on $(i,j)=(A,B),(B,C),(C,A)$.
Note the order of $i,j$ matters for the second term of $H_{int}$. In $H_{t}$
we have rewritten the $J_{pz}$ term using the identity $n_{i;b}=1-n_{i;t}$ and
ignore the constant term and terms linear to $n_{i;t}$, which can be absorbed
into chemical potential. We have defined $M=n_{b}-\langle
a_{i}^{\dagger}a_{i}\rangle$, which is the magnetization of the magnetic order
in the bottom layer. For $H_{b}$ we have ignored higher order terms such as
$a_{i}a^{\dagger}_{j}a_{j}$.

In the following we will ignore the dynamics of the spin wave $a_{i}$ to
simplify our analysis, which represents the gapless spin waves on top of the
$120^{\circ}$ order in the bottom layer. The the main dynamical degree of
freedom is $b_{t\sigma}$. Because of the condensation of
$b_{i;b\sigma^{\prime}}$, now $b_{i;t\sigma}$ is a physical operator.
Physically $b^{\dagger}_{i;t\sigma}\sim b^{\dagger}_{i;t\sigma}\langle
b_{i;b\sigma^{\prime}}\rangle\sim
c^{\dagger}_{i;t\sigma}c_{i;b\sigma^{\prime}}$ now creates an exciton with
$S_{t}=\frac{1}{2}$ and $P_{z}=1$. Note that the exciton carries spin 1/2
under the SU(2) spin rotation in the top layer. The SU(2) spin rotation in the
bottom layer is already broken and $S_{b}$ is no longer well defined.

We can then just relabel $b_{i;\sigma}=b_{i;t\sigma}$, $n_{i}=n_{i;t}$ and
$\vec{S}_{i;t}=\vec{S}_{i}$. Then we basically have a gas of spin $1/2$ bosons
at density $n=x$, described by the Hamiltonian:

$\displaystyle H_{exciton}$
$\displaystyle=-t\sum_{ij}Pb^{\dagger}_{i;\sigma}b_{j;\sigma}P-\mu\sum_{i}b^{\dagger}_{i;\sigma}b_{i;\sigma}$
$\displaystyle+J\sum_{\langle
ij\rangle}\vec{S}_{i}\cdot\vec{S}_{j}+V\sum_{\langle ij\rangle}n_{i}n_{j}$
(22)

where $t=\frac{1}{4}(M_{0}-x)J_{p}$, $V=\frac{1}{2}J_{pz}$ and $J=J_{t}$. The
boson is at density $x$ per site. We have
$\vec{S}_{i}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}b^{\dagger}_{i;\sigma}\vec{\sigma}_{\sigma\sigma^{\prime}}b_{i;\sigma^{\prime}}$
and $n_{i}=\sum_{\sigma}b^{\dagger}_{i;\sigma}b_{i;\sigma}$. $M_{0}$ is the
magnetization of the $120^{\circ}$ order at the layer polarized limit with
$x=0$. With a finite x, the magnetization is further depleted and we have
$M=M_{0}-x$. For the model with $J_{b}^{\prime}=0$, we know $M_{0}\approx
1/2$. We also introduce the projection operator $P$ to implement the
constraint that $n_{i}=0,1$, inherited from the Schwinger boson formulation.

The above effective Hamiltonian works at the small x regime. We can see that
the dynamics of doped excitons on the background of $120^{\circ}$ magnetic
order is basically captured by a spin 1/2 hard core boson gas problem, an
observation already made in our previous paperZhang  _et al._ (2021a). Note
that Eq. 22 can be viewed as a bosonic t-J model. Next we analyze this
effective model. We will find a spin polarized BEC at small $J/t$ and a spin
gapped paired superfluid at large $J/t$. In the intermediate regime there is a
spiral phase where the boson $b$ condenses at non-zero momentum.

####  V.2.1 Spin polarized BEC at small J/t

If we set $J=0$, then it is natural to expect that the ground state of this
model is a spin polarized BEC, similar to what is usually found for BEC with
integer spinKawaguchi and Ueda (2012). For the spin polarized BEC, we take an
ansatz with $\langle b_{i;\uparrow}\rangle=\sqrt{n_{0}}$ and $\langle
b_{i;\downarrow}\rangle=0$. Next we analyze the goldstone modes above this
condensation using the Bogoliubov mean field theory and show that it is
unstable if $\frac{J}{t}$ is large enough.

In momentum space we assume
$b_{\uparrow}(\mathbf{k}=0)=\sqrt{N_{s}}\sqrt{n_{0}}$, where $N_{s}$ is the
number of sites. At mean field level we have energy:

$E_{M}=-(6t+\mu)n_{0}N_{s}+3(V+\frac{1}{4}J)n_{0}^{2}N_{s}$ (23)

Variation with respects to $n_{0}$ leads to:

$\mu=6(V+\frac{1}{4}J)n_{0}-6t$ (24)

Next we consider excitations by also including $b_{\sigma}(\mathbf{q}\neq 0)$.

the mean field Hamiltonian is

$\displaystyle H_{M}$
$\displaystyle=-t\sum_{ij}b_{i;\sigma}^{\dagger}b_{j;\sigma}-\mu\sum_{i}b^{\dagger}_{i;\sigma}b_{i;\sigma}-(6t+\mu)n_{0}N+3(V+\frac{1}{4}J)n_{0}^{2}N+6(V+\frac{1}{4}J)n_{0}\sum_{\mathbf{q}\neq
0}b^{\dagger}_{\uparrow}(\mathbf{q})b_{\uparrow}(\mathbf{q})$
$\displaystyle+\sum_{\mathbf{q}\neq
0}(6V-\frac{3}{2}J+\frac{J}{2}F(\mathbf{q}))n_{0}b^{\dagger}_{\downarrow}(\mathbf{q})b_{\downarrow}(\mathbf{q})+(V+\frac{1}{4}J)n_{0}\sum_{\mathbf{q}\neq
0}b^{\dagger}_{\uparrow}(\mathbf{q})b_{\uparrow}(\mathbf{q})F(\mathbf{q})$
$\displaystyle+\frac{Vn_{0}+\frac{1}{4}Jn_{0}}{2}\sum_{\mathbf{q}\neq
0}\big{(}F(\mathbf{q})b_{\uparrow}(\mathbf{q})b_{\uparrow}(-\mathbf{q})+F(\mathbf{q})b^{\dagger}_{\uparrow}(\mathbf{q})b^{\dagger}_{\uparrow}(-\mathbf{q})\big{)}$
(25)

where $F(\mathbf{q})=\sum_{j}\cos(\mathbf{q}\cdot\mathbf{r_{j}})$ and
$j=1,2,...,6$ lists the six nearest neighbors of site $i=(0,0)$.

It is easy to see that the $b_{\uparrow}$ and $b_{\downarrow}$ part decouple.
We can diagonalize the Hamiltonian using the Bogoliubov transformation and
express the Hamiltonian with new bosonic operator $\gamma(\mathbf{q})$.

In the end we get:

$H_{M}=E_{0}+\sum_{\mathbf{q}}\omega_{\uparrow}(\mathbf{q})\gamma^{\dagger}(\mathbf{q})\gamma(\mathbf{q})+\sum_{\mathbf{q}}\omega_{\downarrow}(\mathbf{q})b^{\dagger}_{\downarrow}(\mathbf{q})b_{\downarrow}(\mathbf{q})$
(26)

where,

$\displaystyle\omega_{\uparrow}(\mathbf{q})=$
$\displaystyle\sqrt{\big{(}t(6-F(q))+(V+\frac{1}{4}J)n_{0}F(\mathbf{q})\big{)}^{2}-\big{(}(V+\frac{1}{4}J)n_{0}F(\mathbf{q})\big{)}^{2}}$
(27)

and

$\omega_{\downarrow}(\mathbf{q})=(t-\frac{J}{2}n_{0})(6-F(\mathbf{q}))$ (28)

Using $F(\mathbf{q})=\sum_{j}\cos(\mathbf{q}\cdot\mathbf{r_{j}})\approx
6-\frac{3}{2}q^{2}$ at small $|\mathbf{q}|$, we find
$\omega_{\uparrow}(\mathbf{q})\approx 3\sqrt{2t(V+\frac{1}{4}J)}|\mathbf{q}|$
and
$\omega_{\downarrow}(\mathbf{q})=\frac{3}{2}(t-\frac{1}{2}Jn_{0})|\mathbf{q}|^{2}$
at small $|q|$. Physically the linear mode from $b_{\uparrow}$ corresponds to
the Goldstone mode breaking the U(1) symmetry generated by $P_{z}$. The
quadratic model from $b_{\downarrow}$ is the spin wave from breaking the
$SU(2)$ spin rotation symmetry of the top layer.

The stabilization of this spin polarized BEC requires that
$t-\frac{1}{2}Jn_{0}>0$. We expect the condensation $n_{0}\approx x$, and we
have $t=\frac{1}{4}(M_{0}-x)J_{p}$ with $M_{0}\approx\frac{1}{2}$ and
$J=J_{t}$. Then the stabilization condition becomes

$\frac{J_{t}}{J_{p}}<\frac{M_{0}}{2x}-1$ (29)

When $\frac{J_{t}}{J_{p}}>\frac{M_{0}}{2x}-1$, the spin polarized BEC is
unstable because of the antiferromagnetic spin-spin coupling. We will discuss
the possibility in the remaining of this section.

####  V.2.2 Paired superfluid at large J/t

When $J<<t$, we believe the ground state of Eq. 22 is a spin polarized BEC.
But when $J/t$ is large, it is obvious that the spin polarization is
disfavored. In this regime, a promising candidate of Eq. 22 is a paired
superfluid phase with
$\langle\epsilon_{\sigma\sigma^{\prime}}b_{i;\sigma}b_{j;\sigma^{\prime}}\rangle\neq
0$, but $\langle b_{i;\sigma}\rangle=0$. In the following we will describe
this paired condensation phase of the excitons and discuss its instability
when decreasing $J/t$.

We define operators
$\hat{T}_{ij}=\sum_{\sigma}b^{\dagger}_{i;\sigma}b_{j;\sigma}$ and
$\hat{\Delta}_{ij}=\sum_{\sigma\sigma^{\prime}}\epsilon_{\sigma\sigma^{\prime}}b_{i;\sigma}b_{j;\sigma^{\prime}}$.
A general symmetric mean field ansatz for the exciton from decoupling Eq. 22
is in the form

$\displaystyle H_{M}$ $\displaystyle=\sum_{ij}T_{ij}\hat{T}_{ij}-\sum_{\langle
ij\rangle}(\Delta_{ij}\hat{\Delta}_{ij}+h.c.)-\mu\sum_{i}n_{i}$ (30)

where $T_{ij}=T_{ji}^{*}$ and $\Delta_{ij}=-\Delta_{ji}$.

We have self consistent equations:

$\displaystyle T_{ij}$
$\displaystyle=\frac{3}{8}J\langle\hat{T}_{ji;a}\rangle-t$
$\displaystyle\Delta_{ij}$
$\displaystyle=\frac{3}{8}J\langle\hat{\Delta}^{\dagger}_{ij}\rangle$
$\displaystyle\langle n_{i}\rangle$ $\displaystyle=x$

where $J=J_{t}$ and $t=\frac{1}{4}J_{p}(M_{0}-x)$. $M_{0}$ is the
magnetization of the $120^{\circ}$ order in the bottom layer. $x$ is the
density of the excitons.

We will only consider uniform ansatz with $T_{ij}=T$. The pairing term needs
to be in the odd angular momentum channel. We consider ansatz with
$\Delta_{i,i\pm\hat{x}}=\pm\Delta$, while the other bonds can be generated by
$C_{3}$ rotation. This ansatz is in the $f$ wave channel and it still
preserves the time reversal symmetry. It is symmetric under $C_{6}$ rotation
if we apply a gauge transformation afterwards. The symmetry property of the
ansatz is the same as the zero-flux ansatz in the projection symmetry group
(PSG) of $Z_{2}$ spin liquids on triangular latticeWang and Vishwanath (2006).
This is expected as the sShwinger boson $b_{i}=b_{i;t}$ shares the same PSG
property as the Schwinger boson $b_{i;b}$ in the bottom layer, which is in the
$120^{\circ}$ ordered phase and is known to be proximate to the zero-flux
ansatzWang and Vishwanath (2006). The mean field Hamiltonian can be rewritten
in momentum space as:

$\displaystyle H_{M}$
$\displaystyle=\sum_{\mathbf{k}}(\xi(\mathbf{k})-\mu)b^{\dagger}_{\sigma}(\mathbf{k})b_{\sigma}(\mathbf{k})$
$\displaystyle+\Delta(\mathbf{k})(b_{\uparrow}(\mathbf{k})b_{\downarrow}(-\mathbf{k})-b_{\downarrow}(\mathbf{k})b_{\uparrow}(-\mathbf{k})+h.c.)$
(32)

with

$\displaystyle\xi(\mathbf{k})$ $\displaystyle=2T(\cos
k_{x}+2\cos\frac{1}{2}k_{x}\cos\frac{\sqrt{3}}{2}k_{y})$
$\displaystyle\Delta(\mathbf{k})$ $\displaystyle=2i\Delta(\sin
k_{x}-2\sin\frac{1}{2}k_{x}\cos\frac{\sqrt{3}}{2}k_{y})$ (33)

Using $\psi_{a}(\mathbf{k})=\begin{pmatrix}b_{\uparrow}(\mathbf{k})\\\
b^{\dagger}_{\downarrow}(-\mathbf{k})\end{pmatrix}$, the Hamiltonian can be
rewritten as:

$H_{a}=\sum_{\mathbf{k}}\psi^{\dagger}_{a}(\mathbf{k})\begin{pmatrix}\xi_{a}(\mathbf{k})-\mu_{a}&\Delta^{*}(\mathbf{k})\\\
\Delta(\mathbf{k})&\xi_{a}(-\mathbf{k})-\mu_{a}\end{pmatrix}\psi_{a}(\mathbf{k})+\sum_{\mathbf{k}}\mu$
(34)

Figure 6: The spin gap $\Delta_{t}$ for the paired superfluid phase from mean
field calculation using Eq. 37 and Eq.
LABEL:eq:paired_exciton_self_consistent_equations.

Using the Bogoliubov transformation, we can get:

$H_{a}=\sum_{\mathbf{k}}\omega_{a}(\mathbf{k})\sum_{\sigma}\gamma^{\dagger}_{\sigma}(\mathbf{k})\gamma_{\sigma}(\mathbf{k})+\sum_{\mathbf{k}}(\omega_{a}(\mathbf{k})+\mu)$
(35)

where $\gamma_{\sigma}(\mathbf{k})$ is bosonic operator which diagonalizes the
Hamiltonian.

We have

$\omega(\mathbf{k})=\sqrt{(TA_{a}(\mathbf{k})-\mu)^{2}-\Delta^{2}B(\mathbf{k})^{2}}$
(36)

where $A(\mathbf{k})=2(\cos
k_{x}+2\cos\frac{1}{2}k_{x}\cos\frac{\sqrt{3}}{2}k_{y})$ and
$B(\mathbf{k})=2(\sin
k_{x}-2\sin\frac{1}{2}k_{x}\cos\frac{\sqrt{3}}{2}k_{y})$.

$\langle\hat{T}_{ij}\rangle$ can be calculated as
$\sum_{ij}\langle\hat{T}_{ij}\rangle=\frac{\partial\langle
H_{M}\rangle}{\partial T}$. Similarly $\sum_{\langle
ij\rangle}\langle\hat{\Delta}_{ij}+\hat{\Delta}^{\dagger}_{ij}\rangle=-\frac{\partial\langle
H_{M}\rangle}{\partial\Delta}$ and $\sum_{i}n_{i}=-\frac{\partial\langle
H_{M}\rangle}{\partial_{\mu}}$. Eventually we get

$\displaystyle\langle\hat{T}_{ij}\rangle$
$\displaystyle=\frac{1}{6N_{s}}\sum_{\mathbf{k}}\frac{A(\mathbf{k})(TA(\mathbf{k})-\mu)}{\omega(\mathbf{k})}$
$\displaystyle\langle\hat{\Delta}_{i,i+\hat{x}}\rangle$
$\displaystyle=\frac{1}{6N_{s}}\sum_{\mathbf{k}}\frac{\Delta
B(\mathbf{k})^{2}}{\omega(\mathbf{k})}$ $\displaystyle\langle n\rangle$
$\displaystyle=\frac{1}{N_{s}}\sum_{\mathbf{k}}\frac{TA(\mathbf{k})-\mu}{\omega(\mathbf{k})}-1$
(37)

where $N_{s}$ is the number of sites in the mesh grid of the momentum space.

Combination of Eq. 37 and Eq.
LABEL:eq:paired_exciton_self_consistent_equations gives a set of self
consistent equations to obtain the mean field parameters $T,\Delta,\mu$ for
given values of $t,J,x$ in the Hamiltonian in Eq. 22. We show the obtained
spin gap $\Delta_{t}$ in Fig. 6. For each $x$, we find the spin gap closes
when $\frac{J_{t}}{J_{p}}<\frac{J_{t;c}}{J_{p}}$, indicating the instability
of the paired superfluid phase for smaller value of $J_{t}$. Intuitively the
single exciton condensation without spin gap is favored when
$\frac{J_{t}}{J_{p}}$ is small because the kinetic energy dominates.

####  V.2.3 Intermediate spiral phase

In the previous two parts, we show that the exciton is in a spin polarized BEC
phase when $\frac{J_{t}}{J_{p}}<\frac{M_{0}}{2x}-1$ and a paired superfluid
phase when $\frac{J_{t}}{J_{p}}>\frac{J_{t;c}}{J_{p}}$. Generically, we expect
$\frac{M_{0}}{2x}-1<\frac{J_{t;c}}{J_{p}}$. For example, when $x=0.1$, the
spin polarized BEC is in the regime $\frac{J_{t}}{J_{p}}<1.5$ and the paired
superfluid phase is in the regime $\frac{J_{t}}{J_{p}}>2$. In this case, we
are left with an intermediate regime. Note that the spin gap $\Delta_{t}$ must
be zero when $\frac{J_{t}}{J_{p}}<\frac{J_{t;c}}{J_{p}}$. So the single
exciton must be condensed with $\langle b_{i;t\sigma}\rangle\neq 0$. But it
can not be a uniform spin polarized BEC because the full spin polarization is
unstable due to the analysis in Sec. V.2.1. A natural possibility is then a
BEC phase at a non-zero momentum $\mathbf{Q}$, so the spin correlation has
anti-ferromagnetic component. The momentum $\mathbf{Q}$ may be generically
incommensurate and we can dub this phase as a spiral BEC phase. Note that both
the spin polarized BEC and the spiral BEC are labeled as phase I in Table. 3,
as we do not distinguish the momentum of the exciton condensation and magnetic
order there. We will show numerical evidence for such spiral phase in the next
subsection.

####  V.2.4 Numerical results

From the theoretical analysis, we know that the ground state goes through spin
polarized exciton condensation, spiral exciton condensation and paired exciton
condensation when we increasing $\frac{J_{t}}{J_{p}}$ while assuming the
bottom layer is in the $120^{\circ}$ ordered phase. Here we offer numerical
evidences for this picture based on infinite density matrix renormalization
group (DMRG) simulationMcCulloch (2008); White (1992, 1993). The infinite DMRG
is performed using the TeNPy Library (version 0.4.0)Hauschild and Pollmann
(2018). We define the triangular lattice unit vectors as:
$\mathbf{a}_{1}=(1,0)$ and $\mathbf{a}_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})$.
We will use a system with periodic boundary along $\mathbf{a}_{2}$ direction:
$O(\mathbf{r}+L_{y}\mathbf{a}_{2})=O(\mathbf{r})$, where $O$ is an arbitrary
operator and $L_{y}$ is the system size along the $\mathbf{a}_{2}$ direction.
We will use $L_{y}=6$. Along $\mathbf{a}_{1}$ the system size is infinite.
However we need to choose a unit cell with size $L_{x}\times L_{y}$. We will
use $L_{x}=6$ so that we can fix exciton density $x$ to be a small value such
as $\frac{1}{18}$. We use bond dimension up to $m=4000$. Typical truncation
error is at order $10^{-4}$.

We fix $J_{b}=J_{p}=J_{pz}=1$ and then vary $J_{t}$ in our DMRG simulation of
Eq. 2. The details of the implementation are shown in the Append. B. We label
the four states as
$\ket{1}=\ket{t\uparrow},\ket{2}=\ket{t\downarrow},\ket{3}=\ket{b\uparrow},\ket{4}=\ket{b\downarrow}$.
The operators can be labeled as $S_{ab}=\ket{a}\bra{b}$ with $a,b=1,2,3,4$. We
used three U(1) quantum numbers
$Q_{1}=S_{22}-S_{11},Q_{2}=S_{33}-S_{11},Q_{3}=S_{44}-S_{11}$. We can extract
correlation lengths within a sector fixed by quantum number using the transfer
matrix technique. We label $\xi_{S_{t}}$ for the largest correlation length in
the sector $(Q_{1},Q_{2},Q_{3})=(2,1,1)$, where the typical operator is
$S_{t}^{-}$. Similarly we label $\xi_{S_{b}}$ for the correlation length in
the sector $(Q_{1},Q_{2},Q_{3})=(0,-1,1)$, where the typical operator is
$S_{b}^{-}$. We label $\xi_{P}$ for the largest correlation length in the
sector $(Q_{1},Q_{2},Q_{3})=(1,2,1)$, where the typical operator is
$S_{31}=P^{-}(I+S_{z})$, which carries charge $-1$ under $P_{z}$. We label
$\xi_{PP}$ for the largest correlation length in the sector
$(Q_{1},Q_{2},Q_{3})=(0,2,2)$, where the typical operator is $S_{31}S_{42}$,
which carries charge $-2$ under $P_{z}$.

We plot the inverse of these various correlation lengths in Fig. 7 with
exciton density $x=\frac{1}{18}$ for a $L_{y}=6$ cylinder. $\xi^{-1}$ is never
zero at a finite bond dimension $m$ even if it is zero when $m=\infty$. The
extrapolation is hard given that we only have results for $m\leq 4000$ due to
computation cost. Therefore it is tricky to determine whether there is a spin
gap $\Delta_{t}$, $\Delta_{b}$ for the spin in the two layers. However, it is
useful to compare $\xi^{-1}_{P}$ and $\xi^{-1}_{PP}$. If the phase is in an
exciton condensation phase, single exciton is more elementary and we expect
$\xi^{-1}_{PP}=2\xi^{-1}_{P}$ because $\xi^{-1}_{PP}$ is roughly the gap of a
pair of exciton when the effective system size is finite (controlled by the
bond dimension $m$). This is exactly what we found for $J_{t}=1$ as shown in
Fig. 7(b). Thus we believe that the phase is in the single exciton
condensation phase for $J_{t}=1$. This would imply that
$\Delta_{t}=\Delta_{b}=0$ because single exciton
$c^{\dagger}_{t;\sigma}c_{b;\sigma^{\prime}}$ carries spin indexes of both
layers and its condensation must lead to gapless spin waves for both layers.
In contrast, when $J_{t}=4$, we can clearly see that
$\xi^{-1}_{PP}<\xi^{-1}_{P}$ and $\xi^{-1}_{PP}$ decreases significantly
faster than $\xi^{-1}_{P}$ when we increase the bond dimension $m$. This
strongly suggests that the phase is in a paired exciton condensation phase
with a gap for single exciton. This would imply that at lease one of
$\Delta_{t},\Delta_{b}$ is finite. We note that
$\xi^{-1}_{S_{b}}<<\xi^{-1}_{S_{t}}$ for the entire range of $J_{t}$, so we
conclude that $\Delta_{t}>0$, while $\Delta_{b}$ may still be zero, though it
really needs data at $m\rightarrow\infty$ to prove $\xi^{-1}_{S_{b}}=0$. A
finite $\Delta_{t}$ is further supported by the fact that $\xi^{-1}_{S_{t}}$
grows almost linearly with $J_{t}$ after $J_{t}\geq 2.2$ (see Fig. 7(a)),
which resembles the behavior of $\Delta_{t}$ from mean field calculation shown
in Fig. 6. Thus we conclude that the ground state is in the phase II (see
Table. 3) when $J_{t}\geq 2.2$ with $\Delta_{t}>0$. In the Schwinger boson
theory, this corresponds to the phase with $\langle b_{b;\sigma}\rangle\neq 0$
and $\langle b_{t;\sigma}\rangle=0$.

To further support this picture, we show the spin-spin structure factors
$\vec{S_{a}}(\mathbf{q})\cdot\vec{S_{a}}(-\mathbf{q})$ for the two layers
$a=t,b$ in Fig. 8. In the bottom layer, there are peaks at momentum
$K,K^{\prime}$ for both $J_{t}=1,4$, consistent with the $120^{\circ}$ order.
In the top layer, the spin structure factor has a peak at $\mathbf{q}=0$ for
$J_{t}=1$, consistent with the spin polarized state discussed in Sec. V.2.1.
When $J_{t}=4$, there is no peak in the spin structure factor in the top
layer, consistent with a spin gap $\Delta_{t}>0$. This agrees with the paired
superfluid phase discussed in Sec. V.2.2.

Figure 7: Correlation length from infinite DMRG. We fix
$J_{b}=J_{p}=J_{pz}=1$. We use system size $L_{y}=6$ and a unit cell size
$L_{x}=6$. The exciton density is fixed at $x=\frac{1}{18}$. (a) Inverse
correlation lengths with $J_{t}$ using bond dimension $m=2000$; (b) Results at
$J_{t}=1$. We used bond dimension $m=2000,3000,4000$. (c) Results at
$J_{t}=4$. We used bond dimension $m=2000,3000,4000$. Figure 8: Spin-spin
structure factors in the two layers. We fix $J_{b}=J_{p}=J_{pz}=1$. We use
system size $L_{y}=6$ and a unit cell size $L_{x}=6$. Bond dimension $m=4000$
is used. The exciton density is fixed at $x=\frac{1}{18}$.
(a)(b)$\langle\vec{S_{b}}(\mathbf{q})\cdot\vec{S_{b}}(-\mathbf{q})\rangle$ in
the bottom layer for $J_{t}=1,4$.
(c)(d)$\langle\vec{S_{t}}(\mathbf{q})\cdot\vec{S_{t}}(-\mathbf{q})\rangle$ in
the top layer for $J_{t}=1,4$.

We also show the single exciton correlation and paired exciton correlation in
the momentum space in Fig. 9. $\langle
P^{\dagger}(\mathbf{q})P^{-}(-\mathbf{q})\rangle$ is the structure factor for
single exciton operator $P^{\dagger}\otimes I$, where $I$ is the identity in
the spin sector. $\langle PP^{\dagger}(\mathbf{q})PP^{-}(-\mathbf{q})\rangle$
is the structure for the paired exciton operator
$PP^{\dagger}_{i}=(\epsilon_{\sigma\sigma^{\prime}}c^{\dagger}_{i;t\sigma}c^{\dagger}_{i+\hat{x};t\sigma^{\prime}})(\epsilon_{\alpha\beta}c_{i;b\alpha}c_{i+\hat{x};b\beta})$(see
Sec. B). For $J_{t}=1$, we can see that $\langle
P^{\dagger}(\mathbf{q})P^{-}(-\mathbf{q})\rangle$ has a peak at
$\mathbf{q}=K,K^{\prime}$. Because $P^{\dagger}\sim
b^{\dagger}_{t;\sigma}b_{b;\sigma^{\prime}}$ and we have $b_{t}$ condenses at
momentum $\mathbf{q}=0$ and $b_{b}$ condenses at momentum $K,K^{\prime}$ to
get the $120^{\circ}$ order, $P^{\dagger}$ needs to carry a momentum
$K,K^{\prime}$. On the other hand, $\langle
PP^{\dagger}(\mathbf{q})PP^{-}(-\mathbf{q})$ is featureless. In contrast, for
$J_{t}=4$, we find that $\langle
P^{\dagger}(\mathbf{q})P^{-}(-\mathbf{q})\rangle$ only has broad peaks at
momentum close to $K,K^{\prime}$ and $\langle
PP^{\dagger}(\mathbf{q})PP^{-}(-\mathbf{q})\rangle$ has a peak at momentum
$\mathbf{q}=0$. Combined with the correlation length
$\xi^{-1}_{P},\xi^{-1}_{PP}$ in Fig. 7, we conclude that the phase is in the
paired exciton condensation phase at momentum $\mathbf{q}=0$, while the single
exciton is gapped.

Figure 9: Exciton and paired exciton structure factors in the two layers. We
fix $J_{b}=J_{p}=J_{pz}=1$. We use system size $L_{y}=6$ and a unit cell size
$L_{x}=6$. Bond dimension $m=4000$ is used. The exciton density is fixed at
$x=\frac{1}{18}$. (a)(b)$\langle
P^{\dagger}(\mathbf{q})P^{-}(-\mathbf{q})\rangle$ in the bottom layer for
$J_{t}=1,4$. (c)(d)$\langle PP^{\dagger}(\mathbf{q})\cdot
PP^{-}(-\mathbf{q})\rangle$ in the top layer for $J_{t}=1,4$.

When $J_{t}$ is in the intermediate regime ($J_{t}=1.2$), we notice that the
spin structure factor in the top layer has peaks at momentum away from
$\mathbf{q}=0$ as shown in Fig. 10(a), consistent with the spiral phase
discussed in Fig. V.2.3 with $b_{t}$ condensed at a non-zero momentum. Because
$P^{\dagger}\sim b^{\dagger}_{t}b_{b}$, $\langle
P^{\dagger}(\mathbf{q})P^{-}(-\mathbf{q})\rangle$ will also have peaks at
momentum away from $K,K^{\prime}$. But this spiral phase seems to be gone
already when $J_{t}=1.5$ (see Fig. 10(b)) and then it is replaced by the
paired exciton condensation phase.

Figure 10: Spin-spin structure factors in the top layers. We fix
$J_{b}=J_{p}=J_{pz}=1$. We use system size $L_{y}=6$ and a unit cell size
$L_{x}=6$. (a)$J_{t}=1.2$ with bond dimension $m=2000$. One can see peaks at
momentum away from $\mathbf{q}=0$, consistent with the spiral phase discussed
in Sec. V.2.3 for intermediate value of $\frac{J_{t}}{J_{p}}$. (b)$J_{t}=1.5$
with bond dimension $m=4000$. No clear peak is seen, suggesting that the spin
in the top layer may already be gapped at $\frac{J_{t}}{J_{p}}=1.5$. Then it
enters the paired exciton condensation phase.

In the above we have assumed that the spin $\vec{S_{b}}$ in the bottom layer
is always in the $120^{\circ}$ ordered phase. This assumption should be valid
in the $x\rightarrow 0$ limit if there is a $120^{\circ}$ order with ordering
moment $M_{0}$ at $x=0$. However, a finite $x$ will deplete the ordering
moment $M\sim M_{0}-x$ to be zero when $x$ is large enough. We will discuss
this interesting case in the following part of this section.

###  V.3 Paired superfluid coexisting with $Z_{2}$ spin liquid

Previously we discussed the case with a magnetic order in the bottom layer. In
this limit the Schwinger boson $b_{b\sigma}$ in the bottom layer is condensed
at momentum $K,K^{\prime}$ with a condensation density $M$, then the low
energy dynamics is captured by a spin 1/2 boson gas represented by the
Schwinger boson $b_{t}$ in the top layer. We can get phase I and phase II in
Table. 3 depending on the competition between the $J_{t}$ and the $J_{p}$
term. Especially we find that there is a finite spin gap at the top layer when
$J_{t}$ is large enough. The next natural question is whether we can also get
rid of the single Schwinger boson condensation in the bottom layer and get a
spin gap also for $\vec{S}_{b}$.

We expect the condensation density for the Schwinger boson $b_{b\sigma}$ at
the bottom layer is $M=M_{0}-x$, where $M_{0}$ is the ordering moment at
$x=0$. Then we can imagine that $M$ is going to vanish when increasing $x$ to
be larger than a critical value $x_{c}=M_{0}$, especially in the case when
$M_{0}$ is small already in the layer polarized limit. This leads to the
possibility of the phase IV, which has spin gap for both layers and a $Z_{2}$
topological order.

We can obtain a mean field ansatz for such a phase in the Schwinger boson mean
field theory. We define operators
$\hat{T}_{ij;a}=\sum_{\sigma}b^{\dagger}_{i;a\sigma}b_{j;a\sigma}$ and
$\hat{\Delta}_{ij;a}=\sum_{\sigma\sigma^{\prime}}\epsilon_{\sigma\sigma^{\prime}}b_{i;a\sigma}b_{j;a\sigma^{\prime}}$.
A general spin rotation symmetric mean field ansatz is in the form

$\displaystyle H_{M}$
$\displaystyle=\sum_{a=t,b}\sum_{ij}T_{ij;a}\hat{T}_{ij;a}-\sum_{\langle
ij\rangle}(\Delta_{ij;a}\hat{\Delta}_{ij;a}+h.c.)$
$\displaystyle-(\mu+D)\sum_{i}n_{i;t}-(\mu-D)\sum_{i}n_{i;b}$ (38)

where $T_{ij}=T_{ji}^{*}$ and $\Delta_{ij}=-\Delta_{ji}$. $D$ is the
displacement field between the two layers. $\mu\pm D$ basically gives the
chemical potential for the two layers separately. There is no hybridization
between the two layers such as the term
$b^{\dagger}_{a\sigma}b_{b\sigma^{\prime}}$ or
$b_{t;\sigma}b_{b;\sigma^{\prime}}$, which are forbidden by the $SU(2)\times
SU(2)$ spin rotation symmetry.

We have self consistent equations:

$\displaystyle T_{ij;a}$
$\displaystyle=\frac{3}{8}J_{a}\langle\hat{T}_{ji;a}\rangle+\frac{1}{2}J_{p}\langle\hat{T}_{ji;\bar{a}}\rangle$
$\displaystyle\Delta_{ij;a}$
$\displaystyle=\frac{3}{8}J_{a}\langle\hat{\Delta}^{\dagger}_{ij;a}\rangle$
$\displaystyle\langle n_{i;t}\rangle$ $\displaystyle=x$ $\displaystyle\langle
n_{i;b}\rangle$ $\displaystyle=(1-x)$ (39)

where $\bar{\sigma}$ labels the opposite spin of $\sigma$. $\bar{a}$ is the
opposite layer index of $a$. The last equations are used to solve the chemical
potential $\mu$ and the displacement field $D$ at fixed exciton density $x$.

It is known that there are two different possible symmetric ansatz for the
Schwinger boson in the spin $1/2$ model: the zero flux ansatz and the $\pi$
flux ansatzWang and Vishwanath (2006). They are distinguished by the
projective symmetry group (PSG). In our case the Schwinger bosons in the two
layers share the same gauge transformation and it is easy to show that the PSG
classification is the same as the spin $1/2$ case and we still only have two
possible ansatz. The zero-flux ansatz is just a uniform ansatz for
$b_{a\sigma}$. The $\pi$ flux ansatz has a $2\times 1$ unit cell and a
projective translation symmetry $T_{x}T_{y}=-T_{y}T_{x}$. If the density of
the Schwinger boson $n_{a}$ is too large, it is known that the ansatz is
unstable to single boson condensation $\langle b_{a\sigma}\rangle\neq 0$,
resulting magnetic order in the layer $a$. Because $n_{t}=x<<1$, we do not
need to worry about the top layer, which should be in a paired condensation
phase provided $J_{t}/J_{p}$ is reasonably large according to the discussion
in the previous subsection. Then the main question is whether $b_{b\sigma}$
condenses or not.

In the $x=0$ limit, it is known that the single Schwinger boson condensation
leads to the $120^{\circ}$ ordered phase with $\mathbf{Q}=K$ and the stripe
ordered phase with $\mathbf{Q}=M$ respectively from the zero flux and the
$\pi$ flux ansatzWang and Vishwanath (2006). The zero flux ansatz and the
$\pi$ flux ansatz are separated by a first order transition by tuning the next
nearest neighbor coupling $\frac{J^{\prime}_{b}}{J_{b}}$ in the
$J_{b}-J^{\prime}_{b}$ modelWang and Vishwanath (2006). Indeed previous
numerical studies find the $120^{\circ}$ and the stripe ordered phases in the
small and large $J^{\prime}_{b}$ regimes. Our goal is to deplete the Schwinger
boson condensation $\langle b_{b\sigma}\rangle$ from $M_{0}$ to $M=M_{0}-x$ by
reducing the density $n_{b}$ from $n_{b}=1$ to $n_{b}=1-x$. When $x>x_{c}$
with $x_{c}=M_{0}$, we believe the symmetric ansatz such as in Eq. 38 is
stable. In principle we can get $x_{c}$ by solving self consistent equations
in Eq. 39. However, the Schwinger boson self consistent equations do not give
quantitatively precise results due to the lack of the quantum fluctuations.
Therefore we mainly use the Schwinger boson mean field ansatz as a guidance
and relies on DMRG simulation to decide the phase diagram.

Figure 11: Illustrated phase diagram with next nearest neighbor coupling
$J^{\prime}_{b}$ and the exciton density $x$. We assume $J_{b}=1$. In the
$x=0$ line, AF I and AF II stand for the 120∘ ordered phase and the stripe
phase respectivelyZhu and White (2015); Hu  _et al._ (2015). SL stands for the
spin liquid in the spin $1/2$ model studied previouslyZhu and White (2015); Hu
_et al._ (2015); Saadatmand and McCulloch (2016); Iqbal  _et al._ (2016). We
assume that $J_{t}$ is large enough so that there is always a spin gap in the
top layer and we can only obtain a paired superfluid (PSF) when doping
excitons. After the ordering moment in the bottom layer is depleted, we
naturally obtain two different Z2 spin liquid phases coexisting with paired
exciton superfluid (PSF). Z2 SL I and Z2 SL II are described by the zero flux
and $\pi$ flux ansatz of Schwinger bosons respectively and the condensation of
$b_{b\sigma}$ leads to the AF I and AF II phases nearby. The dashed line is a
first order transition separating the two different Z2 spin liquids+ PSF
phases.  Figure 12: Energy $E$ vs $J^{\prime}_{b}$ at fixed $J_{t}$. We use
$J_{b}=J_{p}=J_{pz}=1$ and $x=\frac{1}{18}$. The results are from infinite
DMRG with bond dimension $m=2000$ for system size $L_{y}=6$. (a) $J_{t}=1$.
The dashed line is at $J^{\prime}_{b}=0.18$. (b) $J_{t}=4$. The dashed line is
at $J^{\prime}_{b}=0.12$.

We assume that there is also a next nearest neighbor hopping $t^{\prime}$ in
the bottom layer, which results in a next nearest neighbor spin spin coupling
$J^{\prime}_{b}$ entering as $J^{\prime}_{b}\sum_{\langle\langle
ij\rangle\rangle}\vec{S}_{i;b}\cdot\vec{S}_{j;b}$ as included in Eq. 2. In the
$x=0$ limit, we can focus just on the $J_{b}-J^{\prime}_{b}$ model for the
spin 1/2 moment in the bottom layer. Such a model has been studied
beforeKaneko  _et al._ (2014); Li  _et al._ (2015); Zhu and White (2015); Hu
_et al._ (2015); Saadatmand and McCulloch (2016); Iqbal  _et al._ (2016);
Wietek and Läuchli (2017); Misumi  _et al._ (2017); Gong  _et al._ (2017);
Jiang and Jiang (2022). It was found that the ground state is in the
$120^{\circ}$ order ($\mathbf{Q}=K$) in $J^{\prime}_{b}/J_{b}\in[0,0.07]$, a
stripy ordered phase with momentum $\mathbf{Q}=M$ in
$J^{\prime}_{b}/J_{b}\in[0.15,\infty)$ and in a spin liquid phase when
$J^{\prime}_{b}/J_{b}\in[0.07,0.15]$. When approaching the intermediate spin
liquid phase from either magnetic ordered phase, we expect that the ordered
moment $M_{0}$ decreases and vanishes. Then around the intermediate
$J^{\prime}_{b}/J_{b}$, the ordered moment $M_{0}$ can be easily depleted by a
finite $x_{c}=M_{0}$, as illustrated in Fig. 11. The two magnetic ordered
phases are proximate to two different $Z_{2}$ spin liquids corresponding to
zero and $\pi$ flux phase of the Schwinger boson $b_{b}$ respectively at
$x=0$Wang and Vishwanath (2006). We expect $\langle
b_{b\sigma}\rangle=M_{0}-x$ and $\langle b_{t\sigma}\rangle=0$. When
$x>x_{c}=M_{0}$, we should have a phase in the class IV of Table. 3. Given
that both
$\langle\epsilon_{\sigma\sigma^{\prime}}b_{t\sigma}b_{t\sigma^{\prime}}\rangle\neq
0$ and
$\langle\epsilon_{\sigma\sigma^{\prime}}b_{b\sigma}b_{b\sigma^{\prime}}\rangle\neq
0$, we have the paired exciton condensation $\langle PP^{\dagger}\rangle\neq
0$, which coexists with the $Z_{2}$ spin liquids. There are two such Z2 spin
liquids in the small $J^{\prime}_{b}$ and large $J^{\prime}_{b}$, which are
connected by first order transition. With the above argument, we plot an
illustrated phase diagram in Fig. 11.

Figure 13: Spin structure factor at the bottom layer
$\langle\vec{S}_{b}(\mathbf{q})\cdot\vec{S}_{b}(-\mathbf{q})\rangle$ at
different values of $J^{\prime}_{b}$ at fixed $J_{t}=4$. We use
$J_{p}=J_{pz}=1$ and $x=\frac{1}{18}$. The results are from infinite DMRG with
bond dimension $m=2000$. Around $J^{\prime}_{b}=0.12$, there is a clear change
of the momentum of the peak from $\mathbf{Q}=K$ to $\mathbf{Q}=M$, consistent
with a first order phase transition indicated by energy in Fig. 12(b).  Figure
14: Spin structure factor at the bottom layer
$\langle\vec{S}_{b}(\mathbf{q})\cdot\vec{S}_{b}(-\mathbf{q})\rangle$ at
different values of $J^{\prime}_{b}$ at fixed $J_{t}=1$. We use
$J_{p}=J_{pz}=1$ and $x=\frac{1}{18}$. The results are from infinite DMRG with
bond dimension $m=2000$. Around $J^{\prime}_{b}=0.18$, there is a clear change
of the momentum of the peak from $\mathbf{Q}=K$ to $\mathbf{Q}=M$, consistent
with a first order phase transition indicated by energy in Fig. 12(a).

To support the phase diagram, we show DMRG results along the line of changing
$J^{\prime}_{b}$ at fixed $x=\frac{1}{18}$. First, in Fig. 12 we find that the
energy $E$ vs $J^{\prime}_{b}$ has a kink consistent with first order
transition for both $J_{t}=1$ and $J_{t}=4$. The difference between the two
sides across this first order transition line is most easily manifested in the
spin structure factor of the bottom layer as shown in Fig. 13. At $J_{t}=4$,
when changing $J^{\prime}_{b}$ from $0.1$ to $0.12$, momentum of the peak of
the structure factor
$\langle\vec{S}_{b}(\mathbf{q})\cdot\vec{S}_{b}(-\mathbf{q})\rangle$ suddenly
changes from $\mathbf{Q}=K$ to $\mathbf{Q}=M$. The same change also happens
for $J_{t}=1$ shown in Fig. 14. Therefore it is quite obvious that there are
two different phases separated by a first order transition. These two phases
origin from the $120^{\circ}$ order and the stripe order phase at $x=0$, which
are proximate to the zero flux and the $\pi$ flux ansatz of the Schwinger
bosonWang and Vishwanath (2006). The main question now is whether the single
boson $b_{a\sigma}$ is condensed or not, or equivalently whether the spin gap
$\Delta_{a}$ is finite for both layers.

Figure 15: Inverse of correlation lengths with the next-nearest neighbor spin-
spin coupling $J^{\prime}_{b}$ from infinite DMRG with bond dimension
$m=2000,3000,4000$. We fix $J_{b}=J_{p}=J_{pz}=1$ and $J_{t}=4$. We use system
size $L_{y}=6$. The exciton density is fixed at $x=\frac{1}{18}$. The
correlation length is obtained using the transfer matrix method (see Fig. 7
and texts around there). The value at $m=\infty$ is extrapolated using the
formula $\xi^{-1}(m)=\xi^{-1}(m=\infty)+a(\frac{1}{m})^{2}+b\frac{1}{m}$. (a)
$\xi_{S_{t}}$ is the correlation length of the spin in the top layer; (b)
$\xi_{P}$ is the correlation length for the single exciton; (c) $\xi_{PP}$ is
the correlation length of a pair of exciton; (d) $\xi_{S_{b}}$ is the
correlation length for the spin in the bottom layer.

We show the inverse of the correlation lengths from different operators in
Fig. 15, using $J_{b}=J_{p}=J_{pz}=1$ and $J_{t}=4$. In DMRG the correlation
length $\xi$ is always finite for a finite bond dimension $m$, but the scaling
of the inverse of the correlation length with $\frac{1}{m}$ can provide the
information on whether $\frac{1}{\xi}$ is finite for a given sector specified
by the quantum numbers. We find that the correlation lengths for both
$\vec{S}_{t}$ and the single exciton $P^{\dagger}$ are quite small with
$\frac{1}{\xi}>1$ as shown in Fig. 15(a)(b). Meanwhile, in Fig. 15(c) we show
that $\frac{1}{\xi}$ for the paired exciton decreases rapidly with increasing
the bond dimension $m$ and the extrapolated value $\frac{1}{\xi}$ at
$m=\infty$ is smaller than $0.2$ (almost $0$ at some value of
$J^{\prime}_{b}$). We believe the finite value of $\frac{1}{\xi_{PP}}$ is a
numerical artifact and $\frac{1}{\xi_{PP}}=0$ in the whole range if we really
take the bond dimension to infinite. This means that the phase has a paired
exciton condensation order, consistent with the expectation of the paired
superfluid (PSF) order in Fig. 11. The only remaining problem is whether it is
in the AF+PSF phase or in the Z2 SL+PSF phase illustrated in Fig. 11.

Figure 16:
$\langle\vec{S}_{b}(\mathbf{q})\cdot\vec{S}_{b}(-\mathbf{q})\rangle$ at
various values of $J^{\prime}_{b}$. The parameters are the same as in Fig. 15
and we use bond dimension $m=4000$.

In Fig. 15(d) we show the inverse of the correlation length of the spin in the
bottom layer. We find that $\frac{1}{\xi_{S_{b}}}$ increases with
$J^{\prime}_{b}$ and is maximized at the first order critical point at
$J^{\prime}_{b}\approx 0.12$. $\frac{1}{\xi_{S_{t}}}$ extrapolated to
$m=\infty$ at $J^{\prime}_{b}=0.10$ and $J^{\prime}_{b}=0.12$ are as large as
$0.6$ and $0.9$ respectively, suggesting short correlation length
($\xi_{S_{b}}<2$). In Fig. 16 we show the spin-spin structure factor in the
bottom layer for $J^{\prime}_{b}=0$, $0.1$, $0.12$, $0.4$. Obviously
$J^{\prime}_{b}=0.1$ and $J^{\prime}_{b}=0.12$ are separated by the first
order transition denoted as the dashed line in Fig. 11.
$\langle\vec{S}_{b}(\mathbf{q})\cdot\vec{S}_{b}(-\mathbf{q})\rangle$ shows
peaks at $\mathbf{Q}=K$ and $\mathbf{Q}=M$ for $J^{\prime}_{b}=0.10$ and
$J^{\prime}_{b}=0.12$ respectively, but the peaks are blurred compared to
$J^{\prime}_{b}=0.0$ and $J^{\prime}_{b}=0.4$, consistent with the rapid
increase of $\frac{1}{\xi_{S_{b}}}$ towards the intermediate regime. These
results strongly suggest that there is a spin gap $\Delta_{b}>0$ at least at
$J^{\prime}_{b}=0.10$ and $J^{\prime}_{b}=0.12$ and these two points have only
short ranged antiferromagnetic correlations and belong to the Z2 SL I + PSF
phase and the Z2 II + PSF phase in Fig. 11. This is also consistent with our
theoretical expectation that the magnetic order moment $M=M_{0}-x$ must vanish
when the doping density $x$ is finite given that $M_{0}$ should approach zero
towards the intermediate regime. Note that for $J^{\prime}_{b}=0.12$ the peaks
of the spin structure factor are at only two of the three M points, which
breaks the $C_{3}$ rotation symmetry. However, this may be because our system
in DMRG has a small $L_{y}=6$ and breaks the $C_{3}$ symmetry explicitly. A
larger size calculation is needed to determine whether this is a nematic phase
or a $C_{3}$ symmetric phase in the $L_{y}\rightarrow\infty$.

In our DMRG calculation it is not easy to decide the phase boundary between AF
+ PSF phase and Z2 SL + PSF phase because it is always a crossover at finite
bond dimension. For the particular parameter in Fig. 15,
$\frac{1}{\xi_{S_{b}}}(m\rightarrow\infty)$ becomes as small as $0.1$ at
$J^{\prime}_{b}=0.40$, but remains larger than $0.3$ at $J^{\prime}_{b}=0$.
Thus we conjecture that $\Delta_{b}=0$ at $J^{\prime}_{b}=0.4$, but
$\Delta_{b}>0$ even at $J^{\prime}_{b}=0$. This means that the AF I + PSF
phase is absent for this parameter and we are in the Z2 SL I + PSF phase even
with $J^{\prime}_{b}=0$. However, we note that it is tricky to distinguish
$\Delta_{b}>0$ and $\Delta_{b}=0$ when $\Delta_{b}$ is small. Therefore we
leave it to future to determine the precise boundary between Z2 SL and AF
phase using larger bond dimension and system size. In Appendix. D we show more
results to demonstrate that the fractional superfluid phase (Z2 SL + PSF) also
exist with smaller value of $J_{t}=1$ or $J_{t}=2$ and with other parameters
such as $J_{p}=2,J_{pz}=5$. Especially it survives a large $J_{pz}$ term,
which is the repulsion between the excitons and does not matter much at the
small $x$ regime.

Lastly we want to comment on the evolution from the spin liquid in the $x=0$
line to the Z2 SL I + PSF phase in the Fig. 11. Right now there are debates on
whether this spin liquid in the intermediate $\frac{J^{\prime}_{b}}{J_{b}}$ at
the $x=0$ limit is in a Z2 spin liquidZhu and White (2015); Hu  _et al._
(2015); Saadatmand and McCulloch (2016); Jiang and Jiang (2022) or in a U(1)
Dirac spin liquidIqbal  _et al._ (2016); Hu  _et al._ (2019). If it is a Z2
spin liquid, then it should be described by the zero flux ansatz of the
Schwinger boson theory because it has short ranged $120^{\circ}$ correlations.
In this case, doping the exciton induces a chemical potential tuned insulator
to superfluid phase transition for the paired excitons on top of the Z2 spin
liquid which leads to the Z2 SL I + PSF phase found in our numerical
calculation. If the layer polarized spin liquid is a U(1) Dirac spin liquid,
then we should use fermionic spinons such as in Eq. 10 to describe the spin
liquid phase. Then doping excitons will create spinon Fermi surfaces for both
layers. But there is no signature of spinon Fermi surface in our numerical
results at $x=\frac{1}{18}$. To reach the Z2 SL I + PSF phase, we need the
fermion spinons to further have pairing terms
$\epsilon_{\sigma\sigma^{\prime}}f_{a;\sigma}f_{a;\sigma^{\prime}}$ for both
layer $a=t,b$, which higgses the U(1) gauge field down to Z2. Because this is
a $\pi$ flux ansatz of fermionic spinon, it is equivalent to the zero flux
phase of the Schwinger boson222In this $Z_{2}$ spin liquid, the bosonic $e$
particle is a bound state of fermionic spinon and the vison. The vison feels a
$\pi$ flux and hence the $\pi$ flux of fermionic spinon corresponds to the
zero-flux of the bosonic spinon. and is in the same Z2 SL I phase. Thus we see
that there need to be two phase transitions starting from the U(1) Dirac spin
liquid phase to reach the Z2 SL I +PSF phase. Future calculations at
infinitesimal $x$ regime may be able to resolve the debates between the Z2
spin liquid and Dirac spin liquid in the $x=0$ limit.

Experimentally, in moiré + monolayer system, we expect $J_{t}$ to be larger
than $J_{b}$ and $J_{p}$ because the top layer is less correlated. Then
according to our analysis above it is easy to reach the paired superfluid
phase with a spin gap $\Delta_{t}>0$ in the top layer. As long as there is a
small $J^{\prime}_{b}\sim 0.1J_{b}$ in the bottom layer, we can realize either
the Z2 SL I + PSF phase or the Z2 SL II + PSF phases. These phases have spin
gapes in both layers and have preformed cooper pairs, the same as the
resonating-valence-bond (RVB) state first proposed for high Tc
cupratesAnderson (1987). Therefore it is very interesting to search for
possible superconductivity by doping the bilayer Mott insulator with the
$Z_{2}$ spin liquids orders.

##  VI Neutral Fermi surface formed by fermionic excitons

In Sec. IV we argue that the doped exciton can be either bosonic or fermionic
due to possible fractionalization. In Sec. V we did a survey of various
possible superfluid phases built with bosonic excitons, which are either local
excitations on top of a magnetic order or fractional excitations coupled to a
deconfined $Z_{2}$ gauge field. In this section, we turn to the more exotic
possibility of fermionic excitons. In this case, the fermionic exciton
obviously can not be local excitations and necessarily couples to a deconfined
gauge field.

As shown in Sec. IV, the spin operator $\vec{S}$ and exciton operator
$\vec{P}$ can be expressed either with bosonic partons $b_{a\sigma}$ or
fermionic parton $f_{a\sigma}$. When the layer polarized Mott insulator is in
the magnetically ordered phase, the natural parton description should be the
Schwinger boson theory shown in Sec. V. When $t/U$ is in the intermediate
regime, the layer polarized Mott insulator may be in a spinon Fermi surface
stateMotrunich (2005); Senthil (2008). The exact nature of the intermediate
weak Mott insulator of the spin 1/2 Hubbard model is still under debate and
may depend on details like whether further neighbor hoppings are included or
not. In the experimental study of the metal-insulator transition in the AA
stacked MoTe2/WSe2 moiré systemLi  _et al._ (2021), no magnetic order is found
down to the lowest temperature. The spin susceptibility in the weak Mott
regime seems to be just like the metallic phase, suggesting a possible spinon
Fermi surface ground state. However, such a neutral spinon Fermi surface is
hard to detect even if it already exists in the MoTe2/Wse2 moiré system. In
this section we propose to study a MoTe2-hBN-MoTe2/WSe2 system similar to the
moiré+monolayer setting up in Ref. Gu  _et al._ , 2021; Zhang  _et al._ ,
2021b. If the layer polarized Mott insulator indeed hosts a spinon Fermi
surface, then at finite exciton density $x>0$, the most natural phase is a
$U(1)$ spin liquid with spinon Fermi surfaces in both layers. Especially the
neutral fermion in the top layer can be viewed as a fermionic exciton formed
by electron $c_{t\sigma}^{\dagger}$ in the top layer bound to the holon
$\varphi$ in the bottom layer, as illustrated in Fig. 4. Possibility of
fermionic excitons has also been discussed previously in Kondo insulator
SmB6Chowdhury  _et al._ (2018) or in bilayer Landau levelsBarkeshli  _et al._
(2018); Zaletel  _et al._ (2018). In our realization the neutral fermion
carries a charge under a $U(1)$ symmetry generated by $P_{z}$ (the same is
true for the Fermionic exciton in the bilayer quantum Hall systemBarkeshli
_et al._ (2018); Zaletel  _et al._ (2018)). This makes it possible to detect
the movement of the neutral fermionic excitons through counter-flow transport.

To favor a U(1) spin liquid with spinon Fermi surface state, it is necessary
to include higher order ring exchange terms beyond our simple model in Eq. 2.
The task to establish the spinon Fermi surface state as a ground state of a
microscopic lattice model is very challenging and we leave it to future work.
Here our focus is to explore the fate under doping excitons while assuming the
$x=0$ limit indeed hosts a spinon Fermi surface. As said before, this
assumption is encouraged by the recent experimentLi  _et al._ (2021). But it
is hard to detect a neutral Fermi surface even if it already exists. The
purpose of this section is to demonstrate that doping excitons can lead to
smoking gun evidences of neutral Fermi surface from electric measurements.

Under the assumption that the $x=0$ limit has a spinon Fermi surface in the
bottom layer, we should use the Abrikosov fermion partons $f_{t\sigma}$ and
$f_{b\sigma}$ to represent the spin operators, as introduced in Eq. 10. A
typical mean field theory in terms of these fermionic partons is in the form:

$\displaystyle H_{M}$
$\displaystyle=-\sum_{ij}t_{b;ij}f^{\dagger}_{i;b\sigma}f_{j;b\sigma}-(\mu-D)\sum_{i}f^{\dagger}_{i;b\sigma}f_{i;b\sigma}$
$\displaystyle~{}~{}~{}-\sum_{ij}t_{t;ij}f^{\dagger}_{i;t\sigma}f_{j;t\sigma}-(\mu+D)\sum_{i}f^{\dagger}_{i;t\sigma}f_{i;t\sigma}$
(40)

where $\mu$ is the chemical potential to fix $n_{t}+n_{b}=1$ and $D$ is the
displacement field.

When $D\rightarrow-\infty$, the system is in a layer polarized Mott insulator
with spinon Fermi surface only in the bottom layer, as illustrated in Fig.
17(a). Then when increasing $D$, there is a Lifshitz transition at $D_{c}$
(see Fig. 17(b)), after which we have spinon Fermi surfaces in both layers,
illustrated in Fig. 17(c).

Figure 17: Illustration for the displacement field $D$ tuned transition. The
red line is the dispersion of the spinon in the bottom layer and the blue is
the dispersion of the spinon in the top layer. The dashed line is the chemical
potential as defined in Eq. 40. (a) A layer polarized Mott insulator with
spinon Fermi surface in the bottom layer. $\frac{\partial
P_{z}}{\partial_{D}}=0$ and counterflow transport is insulating. (b) At
critical $D_{c}$, the gap of the exciton is now zero. (c) A $U(1)$ spin liquid
phase with spinon Fermi surfaces in both layers. $\frac{\partial
P_{z}}{\partial D}$ is finite and the counterflow transport is now metallic.

A low energy theory can be written down to describe this spinon Fermi surface
phase which needs to couple to an internal gauge field:

$\displaystyle L$
$\displaystyle=\psi^{\dagger}_{t;\sigma}(\tau,x)\big{(}\partial_{\tau}-(\mu+D)-ia_{0}(\tau,x)-\frac{1}{2}iA^{s}_{0}(\tau,x)\big{)}\psi_{t;\sigma}(\tau,x)-\frac{\hbar^{2}}{2m_{t}}\psi^{\dagger}_{t;\sigma}(\tau,x)\big{(}-i\vec{\partial}-\vec{a}(\tau,x)-\frac{1}{2}\vec{A^{s}}(\tau,x)\big{)}^{2}\psi_{t;\sigma}(\tau,x)$
$\displaystyle~{}+\psi^{\dagger}_{b;\sigma}(\tau,x)\big{(}\partial_{\tau}-(\mu-D)-ia_{0}(\tau,x)+\frac{1}{2}iA^{s}_{0}(\tau_{x})\big{)}\psi_{b;\sigma}(\tau,x)-\frac{\hbar^{2}}{2m_{b}}\psi^{\dagger}_{b;\sigma}(\tau,x)\big{(}-i\vec{\partial}-\vec{a}(\tau,x)+\frac{1}{2}\vec{A^{s}}(\tau,x)\big{)}^{2}\psi_{b;\sigma}(\tau,x)$
(41)

where $\psi_{a;\sigma}(\tau,x)$ is the fermionic field in the layer $a=t,b$.
$a_{\mu}(\tau,x)$ is the internal U(1) gauge field associated with the gauge
symmetry: $f_{i;a\sigma}\rightarrow f_{i;a\sigma}e^{i\alpha_{i}}$ for $a=t,b$.
$A^{s}_{\mu}(\tau,x)$ is a probing field for the global U(1) symmetry
associated with $P_{z}$. Physically its electric field
$\vec{E}_{s}=-\vec{\partial}A^{s}_{0}+\frac{\partial\vec{A}^{s}}{\partial t}$
corresponds to a counter-flow electric field configuration:
$\vec{E}_{t}=-\vec{E}_{b}$ in the two layers. The probing field $A^{s}_{\mu}$
minimally couples to the counter-flow current: $-J^{\mu}_{s}A^{s}_{\mu}$. Here
$J^{s}_{0}=P_{z}$ is the layer polarization and $\vec{J}^{s}$ is its current.
We do not include the coupling to $A^{c}_{\mu}$ because it can always be
absorbed to the coupling to $a_{\mu}$ by redefinition $a_{\mu}\rightarrow
a_{\mu}-A^{c}_{\mu}$. This is another way to understand that the above theory
describes a Mott insulator without response to $A^{c}_{\mu}$.

In the following we provide two clear experimental signatures to detect this
exotic phase with neutral Fermi surfaces in both layers.

Counter-flow resistivity We first show that the counter-flow transport is
insulating when $D<D_{c}$ and metallic when $D>D_{c}$. Thus the displacement
field tunes an counter-flow insulator to metal transition. The counter-flow
resistivity is defined as $\rho_{s;xx}=\frac{E^{s}_{x}}{J^{s}_{x}}$, where
$\vec{E}=\vec{E}^{t}-\vec{E}^{b}$ and
$\vec{J}^{s}=\frac{1}{2}(\vec{J}^{t}-\vec{J}^{b})$. As shown in Eq. 13, the
counter-flow resistivity $\rho_{s}=\rho_{t}+\rho_{b}$, where
$\rho_{t},\rho_{b}$ are the resistivities of the fermionic partons
$f_{t\sigma}$ and $f_{b\sigma}$. When $D<D_{c}$, $\rho_{t}=\infty$ and hence
$\rho_{s}=\infty$. When $D>D_{c}$, we expect that $\rho_{t}=A_{t}T^{\alpha}$
and $\rho_{b}=A_{b}T^{\alpha}$. Therefore the counter-flow resistivity is
metallic $\rho_{s}(T)=T^{\alpha}$. The exponent $\alpha$ on the temperature
dependence is likely deviating from the Fermi liquid result $\alpha=2$ due to
the coupling to the gauge field $a_{\mu}$Lee and Nagaosa (1992); Maslov  _et
al._ (2011); Hartnoll  _et al._ (2014); Lee (2018, 2021), though the exact
value of $\alpha$ is still not well established theoretically. We will not try
to solve this problem in this paper. Instead, we propose to measure the
exponent $\alpha$ through counter-flow transport in experiments.

In addition to a metallic counter-flow transport, the spinon Fermi surface
state when $D>D_{c}$ also has a metallic spin susceptibility and inter-layer
polarizability. As shown in Eq. 14, the inter-layer polarizability
$\kappa_{s}=(\frac{1}{\kappa_{t}}+\frac{1}{\kappa_{b}})^{-1}$, where
$\kappa_{t}$ and $\kappa_{b}$ are the compressibilities of the fermionic
partons in the two layers. When $D>D_{c}$, $\kappa_{a}$ is a constant for
$a=t,b$. Therefore the inter-layer polarizability $\kappa_{s}$ is also a
constant just like a Fermi liquid. Note that the fermionic parton
$f_{a;\sigma}$ acquires a non-Fermi-liquid self energy
$\Sigma(i\omega)\propto|\omega|^{\frac{2}{3}}$Lee and Nagaosa (1992) from
coupling to the gauge field $a_{\mu}$, but it does not enter the spin
susceptibility because of the cancellation with the vertex correctionKim  _et
al._ (1994). Thus the spin susceptibility (or inter-layer polarizability)
behaves as a Fermi liquid. We need to emphasize that an exciton superfluid as
discussed in Sec. V also has constant inter-layer polarizability. Then a
constant inter-layer polarizability does not distinguish the superfluid and
neutral Fermi surface phase. To better characterize these two different kind
of phases, we really need counter-flow transport, which is either superfluid
or metallic in the exciton condensation and neutral Fermi surface phase.

Friedel oscillation in layer polarization It is known that the existence of a
Fermi surface can lead to Friedel oscillation because of the $2k_{F}$
scattering. In a Mott insulator with only spinon Fermi surface, the dominant
Friedel oscillation is in the magnetization and hence is not easy to detect
using the electric probe such as Scanning Tunneling Microscope (STM). An
electric signal may be found in the weak Mott regime where the holon density
is not exactly frozen to be one per siteMross and Senthil (2011), but the
amplitude is suppressed if the charge gap is large.

In our case at $D>D_{c}$, there are active neutral fermions $f_{a;\sigma}$ in
both layers and the fermionic parton carries not only spin $\vec{S}_{a}$ but
also the layer polarization $P_{z}$. Because
$P_{z}(\mathbf{r})=\frac{1}{2}\big{(}n_{t}(\mathbf{r})-n_{b}(\mathbf{r})\big{)}$,
its correlation function manifests $2k_{F}$ singularity of spinon Fermi
surfaces in both layers. Therefore, there will also be Friedel oscillations in
terms of the layer polarization:

$\langle
P_{z}(\mathbf{r})P_{z}(0)\rangle=\sum_{\mathbf{Q}=2\mathbf{k_{F}^{t}}}\frac{A_{t}}{|\mathbf{r}|^{\alpha_{t}(\mathbf{Q})}}e^{i\mathbf{Q}\cdot\mathbf{r}}+\sum_{\mathbf{Q}=2\mathbf{k_{F}^{b}}}\frac{A_{b}}{|\mathbf{r}|^{\alpha_{b}(\mathbf{Q})}}e^{i\mathbf{Q}\cdot\mathbf{r}}$
(42)

where $2\mathbf{k_{F}^{a}}$ is a vector connection two points in the Fermi
surface in layer $a=t,b$.

We can calculate the Fourier transformation of $\langle
P_{z}(\mathbf{r})P_{z}(0)\rangle$ in the mean field level:

$\displaystyle\chi_{zz}(\mathbf{q})=\langle
P_{z}(q_{0}=0,\mathbf{q})P_{z}(-q_{0}=0,-\mathbf{q})\rangle$
$\displaystyle=-\frac{1}{4}\sum_{a\sigma}\sum_{\omega}G_{a\sigma}(i\omega,\mathbf{k}+\mathbf{q})G_{a\sigma}(i\omega,\mathbf{k})$
$\displaystyle=\frac{1}{2}\sum_{a=t,b}\sum_{\mathbf{k}}\frac{f(\xi_{a}(\mathbf{k}))-f(\xi_{a}(\mathbf{k}+\mathbf{q}))}{\xi_{a}(\mathbf{k}+\mathbf{q})-\xi_{a}(\mathbf{k})}$
(43)

where $f(\xi(\mathbf{k}))=\theta(-\xi(\mathbf{k}))$ is the Fermi-Dirac
distribution at $T=0$.

We show the plot of $\chi_{zz}(\mathbf{q})$ in Fig. 18. As expected, there are
features along two circles centering at the $\Gamma$ point. They correspond to
$2k_{F;t}$ and $2k_{F;b}$. The one from $2k_{F;b}$ is outside the Brillouin
zone (BZ) and needs to be folded back to the first BZ. We have ignored the
gauge fluctuation in this calculation, which can further enhance the Friedel
oscillationsMross  _et al._ (2010).

Figure 18:  The Friedel oscillation in terms of the layer polarization. The
color plot shows $\chi_{zz}(\mathbf{q})=\langle
P_{z}(q_{0}=0,\mathbf{q})P_{z}(q_{0}=0,-\mathbf{q})\rangle$ at the mean field
level without gauge fluctuation. $q$ is in units of $2\pi$. The solid black
line is the boundary of the first Brillouin zone (BZ). We have used
$t_{t}=t_{b}=1$ in Eq. 40. There are two circles corresponding to $2k_{F;b}$
and $2k_{F;t}$, which are fixed to be $\frac{\pi
k_{F;t}^{2}}{4\pi^{2}}=\frac{x}{2}$ and $\frac{\pi
k_{F;b}^{2}}{4\pi^{2}}=\frac{1-x}{2}$.

The Friedel oscillation of the layer polarization $P_{z}$ may be detected in
the following way: one can apply a local electric field in z direction
$E_{z}(\mathbf{r}_{0})$, then we expect that $\langle
P_{z}(\mathbf{r})\rangle\sim\sum_{\mathbf{Q}=2\mathbf{k_{F}^{t}}}\frac{A_{t}}{|\mathbf{r}-\mathbf{r}_{0}|^{\alpha_{t}(\mathbf{Q})}}e^{i\mathbf{Q}\cdot(\mathbf{r}-\mathbf{r}_{0})}+\sum_{\mathbf{Q}=2\mathbf{k_{F}^{b}}}\frac{A_{b}}{|\mathbf{r}-\mathbf{r}_{0}|^{\alpha_{b}(\mathbf{Q})}}e^{i\mathbf{Q}\cdot(\mathbf{r}-\mathbf{r}_{0})}$.
The local profile of layer polarization $P_{z}(\mathbf{r})$ will further
induce electric field profile $E_{z}(\mathbf{r})\sim P_{z}(\mathbf{r})$, which
may be detected by electric probes. The Fourier transformation of
$E_{z}(\mathbf{r}-\mathbf{r}_{0})$ should be similar to the plot in Fig. 18.

In addition to the counter-flow transport and the Friedel oscillation, the
spinon Fermi surface state also has $C/T\sim T^{-\frac{1}{3}}$Lee and Nagaosa
(1992) where $C$ is the specific heat. Meanwhile the neutral Fermi surface may
show quantum oscillations under external magnetic field due to coupling
between the internal gauge flux $b=\nabla\times\vec{a}$ and the external
magnetic fieldMotrunich (2006); Chowdhury  _et al._ (2018); Sodemann  _et al._
(2018). But there is a large g factor in TMD material and whether spinon Fermi
surface state can survive under a finite field is not clear. In this sense the
electric probes such as counter-flow or Friedel oscillation at zero magnetic
field may be a better detection scheme.

In the above we start from the assumption that the layer polarized Mott
insulator at the weak Mott regime is a spinon Fermi surface state and then
argue that the small $x$ regime naturally hosts neutral Fermi surfaces in both
layers. Whether the assumption is correct is a hard problem both theoretically
and experimentally. Although there are suggestive signatures of constant spin
susceptibility at the $x=0$ limit in the MoTe2/WSe2 systemLi  _et al._ (2021),
it is hard to prove or rule out the possibility of spinon Fermi surface in the
$x=0$ sample because of the lack of the probe of neutral spin excitations. Our
proposal at finite $x$ in the MoTe2-hBN-MoTe2/WSe2 system is thus also a good
test of the nature of weak Mott insulator at the $x=0$ limit.

###  VI.1 Numerical results in one dimension

In the above we discussed the experimental signatures of the spinon Fermi
surface state, but a microscopic calculation of a 2D lattice model is avoided
due to its technical challenge. Here we provide numerical simulation for
analog of the spinon Fermi surface state in the one dimension chain. We will
show the 1D model has qualitatively the same physics as discussed in the
above. Note that the existence of exotic spinon Fermi surface-like phase in
one dimension is very special and can not be easily generalized to two
dimension. Indeed, we already showed in previous sections that the simple
model in Eq. 2 in two dimension hosts either single exciton or paired exciton
superfluid and we never found a neutral Fermi surface. To favor a spinon Fermi
surface in two dimension, more complicated higher order ring exchange terms
must be included and it is still not clear when a spinon Fermi surface is a
grond state. In this paper we avoid this challenging energetical problem. The
purpose of the 1D calculation is to demonstrate the behavior we expect when
doping excitons into a spinon Fermi surface Mott insulator at the x=0 limit.
We leave it to future work to understand when and why a spinon Fermi surface
state can be stabilized in a two dimension model. As said in the introduction
of this section, there is already encouraging experimental evidence for such a
state at the $x=0$ limit in the MoTe2/WSe2 systemLi  _et al._ (2021) and we
hope our calculation of the 1D model below can provide some insights on the
fate of doping excitons into the weak Mott insulator in the MoTe2/WSe2 system.

In one dimension, the $x=0$ limit of Eq. 2 is a spin $1/2$ chain and is in a
gapless phase without magnetic order. Such a phase described by SU(2)1
conformal field theory (CFT) is known to be an analog of the spinon Fermi
surface in higher dimension. Actually a good model wavefunction for the state
is a Gutzwiller projection of Fermi seaHaldane (1988), the same as the spinon
Fermi surface state in higher dimension. According to our analysis above, we
should expect the phase at finite x to host spinon Fermi surfaces in both
layers. Such a phase has a metallic counter-flow conductance and also show
$2k_{F}$ singularities in spin-spin correlation and dipole-dipole correlation.