\section{Supplementary Materials}

\subsection{Momentum-space tight-binding characterization}

The Bloch states of the hexagonal lattice $V(\mathbf{r})=-V_{0}|%
\sum_{j=1}^{3}\mathbf{\epsilon }_{j}\exp [i\mathbf{k}_{L,j}\cdot (\mathbf{r}-%
\mathbf{r_{0}})]|^{2}$ can be written as $\phi _{l\mathbf{k}}(\mathbf{r}%
)=e^{i\mathbf{k}\cdot \mathbf{r}}u_{l\mathbf{k}}(\mathbf{r})$. The periodic
part can be expanded as $u_{l\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{p}}c_{l%
\mathbf{k}}^{\mathbf{p}}e^{i\mathbf{p}\cdot \mathbf{r}}$, where $\mathbf{p}%
=p_{1}\mathbf{e}_{1}+p_{2}\mathbf{e}_{2}$ with $\mathbf{e}_{i}$ the
primitive reciprocal lattice vectors and $p_{i}$ integers. By substituting
expansion of $\phi _{l\mathbf{k}}(\mathbf{r})$ into the Schr\"{o}dinger
equation $[-\frac{\mathbf{\nabla }^{2}}{2m}+V(\mathbf{r})]\phi _{l\mathbf{k}%
}=\mathcal{E}_{l\mathbf{k}}\phi _{l\mathbf{k}}$, the Bloch states (i.e., the
coefficients $c_{l\mathbf{k}}^{\mathbf{p}}$) and bands can be obtained. The
bare Bloch bands $\mathcal{E}_{sl\mathbf{k}_{s}}$ and states $\phi _{sl%
\mathbf{k}_{s}}$ of the two pseudospin states can be obtained similarly, as
shown in Fig.~\ref{fig:S1}. Here we keep the expansion coefficients up to $%
p_{i}=\pm 6$ in the calculation. We see that the two lowest bands (i.e., $%
\mathcal{E}_{sl\mathbf{k}_{s}}$ with $l=1,2$), which have a gap from higher
bands, form two Dirac cones at valleys $K_{s}$ and $K_{s}^{\prime }$~\cite%
{science.aad5812S}. For typical lattice depths, the two Dirac bands are
asymmetric with respect to Dirac-point energy $E=0$ due to the long-range
tunnelings beyond nearest neighbors that break the sublattice symmetry. Such
effect is more significant for shallower lattices.

\begin{figure}[b]
\includegraphics[width=0.85\linewidth]{FigS1.pdf}
\caption{Band structures of a hexagonal lattice with (a) $V_{0}=6E_{R}$, (b)
$V_{0}=4E_{R}$ and (c) $V_{0}=2E_{R}$. We have shifted the Dirac-point
energy to $E=0$. For a shallower lattice, the tunnelings (both short and
long range tunnelings) are stronger, therefore, the Dirac velocity is larger
and the nonlinearity of the dispersion around the Dirac points is stronger.
In addition, the two Dirac bands become more asymmetric with respect to the
Dirac-point energy $E=0$ due to the enhanced long-range tunnelings that
break the sublattice symmetry.}
\label{fig:S1}
\end{figure}

We are interested in the low-energy physics near the Dirac points, therefore
we only retain the two-lowest bare Dirac bands (i.e., $\mathcal{E}_{sl%
\mathbf{k}_{s}}$ with $l=1,2$) and drop all higher bands with $l>2$ safely.
That is, we keep only one Wannier orbital at each hexagonal lattice site,
and these Wannier orbitals form the full tight-binding basis set. We find
that it is more convenient to work in the Bloch basis $\{\phi _{sl\mathbf{k}%
_{s}}\}$, which is equivalent to the Wannier basis up to a Fourier
transformation. For commensurate twists $\cos (\theta )=\frac{n^{2}+m^{2}+4mn%
}{2(n^{2}+m^{2}+mn)}$, one can fold the bare Brillouin zone of each spin
states to the Moir\'{e} Brillouin zone. The inter-spin coupling coefficient
can be obtained as $\langle \phi _{\uparrow l\mathbf{k}_{\uparrow }}|\Omega
|\phi _{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}\rangle =\delta _{%
\mathbf{q},\mathbf{q}^{\prime }}J_{\mathbf{g}_{\uparrow }\mathbf{g}%
_{\downarrow }}^{ll^{\prime }}(\mathbf{q})$, where $\mathbf{k}_{\uparrow }=%
\mathbf{q}+\mathbf{g}_{\uparrow }$ and $\mathbf{k}_{\downarrow }=\mathbf{%
q^{\prime }}+\mathbf{g}_{\downarrow }$, with $\mathbf{g}_{s}$ the Moir\'{e}
reciprocal lattice vectors and $\mathbf{q}$ the Moir\'{e} Bloch momentum
that is a good quantum number. We obtain the Hamiltonian Eq.~(1) in the main
text, which incorporates all range real-space tunnelings with high accuracy.
Another advantage of this momentum-space approach is that we only need to
keep $\mathbf{g}_{s}$ around the Dirac cones to correctly characterize the
low-energy physics at small twist angles, leading to rather rapid
convergence of the basis set.

Now we show how to evaluate the inter-spin coupling coefficients
\begin{eqnarray}
\langle \phi _{\uparrow l\mathbf{k}_{\uparrow }}|\Omega |\phi _{\downarrow
l^{\prime }\mathbf{k}_{\downarrow }}\rangle &=&\Omega \int d^{2}\mathbf{r}%
e^{-i\mathbf{k}_{\uparrow }\cdot \mathbf{r}}u_{\uparrow l\mathbf{k}%
_{\uparrow }}^{\ast }(\mathbf{r})e^{i\mathbf{k}_{\downarrow }\cdot \mathbf{r}%
}u_{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}(\mathbf{r})  \nonumber \\
&=&\Omega \sum_{\mathbf{M}}e^{i(\mathbf{q^{\prime }}-\mathbf{q})\cdot
\mathbf{M}}\langle u_{\uparrow l\mathbf{k}_{\uparrow }}|e^{i(\mathbf{k}%
_{\downarrow }-\mathbf{k}_{\uparrow })\cdot (\mathbf{r-M})}|u_{\downarrow
l^{\prime }\mathbf{k}_{\downarrow }}\rangle _{\mathbf{M}}  \nonumber \\
&=&\Omega \delta _{\mathbf{q},\mathbf{q}^{\prime }}\langle u_{\uparrow l%
\mathbf{k}_{\uparrow }}|e^{i(\mathbf{g}_{\downarrow }-\mathbf{g}_{\uparrow
})\cdot \mathbf{r}}|u_{\downarrow l^{\prime }\mathbf{k}_{\downarrow
}}\rangle .
\end{eqnarray}%
Here $\mathbf{M}$ denotes the Moir\'{e} lattice vectors, and the term $%
\langle \cdots \rangle _{\mathbf{M}}$ in the second line (with integral over
the $\mathbf{M}$-th Moir\'{e} unit cell) is independent of $\mathbf{M}$. The
coefficients $J_{\mathbf{g}_{\uparrow }\mathbf{g}_{\downarrow }}^{ll^{\prime
}}(\mathbf{q})$ is
\begin{eqnarray}
J_{\mathbf{g}_{\uparrow }\mathbf{g}_{\downarrow }}^{ll^{\prime }}(\mathbf{q}%
) &=&\Omega \langle u_{\uparrow l\mathbf{k}_{\uparrow }}|e^{i(\mathbf{g}%
_{\downarrow }-\mathbf{g}_{\uparrow })\cdot \mathbf{r}}|u_{\downarrow
l^{\prime }\mathbf{k}_{\downarrow }}\rangle  \nonumber \\
&=&\Omega \sum_{\mathbf{p}_{\uparrow },\mathbf{p}_{\downarrow }}c_{\uparrow l%
\mathbf{k}_{\uparrow }}^{\mathbf{p}_{\uparrow }\ast }c_{\downarrow l^{\prime
}\mathbf{k}_{\downarrow }}^{\mathbf{p}_{\downarrow }}\int d^{2}\mathbf{r}%
e^{i(\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}_{\downarrow
}-\mathbf{g}_{\uparrow })\cdot \mathbf{r}}  \nonumber \\
&=&\Omega \sum_{\mathbf{p}_{\uparrow },\mathbf{p}_{\downarrow },\mathbf{M}%
}c_{\uparrow l\mathbf{k}_{\uparrow }}^{\mathbf{p}_{\uparrow }\ast
}c_{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}^{\mathbf{p}_{\downarrow
}}\frac{e^{i(\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}%
_{\downarrow }-\mathbf{g}_{\uparrow })\cdot \mathbf{L}_{1}}-1}{(\mathbf{p}%
_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}_{\downarrow }-\mathbf{g}%
_{\uparrow })\cdot \mathbf{L}_{1}}\frac{1-e^{i(\mathbf{p}_{\downarrow }-%
\mathbf{p}_{\uparrow }+\mathbf{g}_{\downarrow }-\mathbf{g}_{\uparrow })\cdot
\mathbf{L}_{2}}}{(\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}%
_{\downarrow }-\mathbf{g}_{\uparrow })\cdot \mathbf{L}_{2}}  \label{eq:J} \\
&=&\Omega \sum_{\mathbf{p}_{\uparrow },\mathbf{p}_{\downarrow },\mathbf{M}%
}c_{\uparrow l\mathbf{k}_{\uparrow }}^{\mathbf{p}_{\uparrow }\ast
}c_{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}^{\mathbf{p}_{\downarrow
}}\delta _{\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow },\mathbf{g}%
_{\uparrow }-\mathbf{g}_{\downarrow }}.  \nonumber
\end{eqnarray}%
To obtain the last step, we have used $\mathbf{p}_{s}\cdot \mathbf{L}%
_{1}=2\pi (mp_{s1}+np_{s2})$ (similarly for $\mathbf{L}_{2}$), since $%
\mathbf{p}_{s}=p_{s1}\mathbf{e}_{s\mathbf{k}_{1}}+p_{s2}\mathbf{e}_{s\mathbf{%
k}_{2}}$ and $\mathbf{L}_{1}=m\mathbf{e}_{s1}+n\mathbf{e}_{s2}$, with $%
\mathbf{e}_{s1},\mathbf{e}_{s2}$ ($\mathbf{e}_{s\mathbf{k}_{1}},\mathbf{e}_{s%
\mathbf{k}_{2}}$) the primitive (reciprocal) lattice vectors of trap $V_{s}$.

\begin{figure}[b]
\includegraphics[width=1.0\linewidth]{FigS2.pdf}
\caption{(a) and (b) Typical density distributions of the Moire flat band
Bloch states $\protect\psi _{j\mathbf{q}}(\mathbf{r},\downarrow )$ in one
Moir\'{e} unit cell for $\mathbf{q}=K$ and $\mathbf{q}=\Gamma $,
respectively. Atoms are mainly distributed around the $AA$ region for $%
\mathbf{q}=K$ and become more uniformly distributed for $\mathbf{q}=\Gamma $%
. Other parameters are $\protect\theta =5.086^{\circ }$, $V_{0}=6E_{R}$ and $%
\Omega =0.116$. (c) Moir\'{e} bands at incommensurate twist angle $\protect%
\theta =4^{\circ }$ based on the generalized one-valley continuum model,
with $V_{0}=6E_{R}$ and $\Omega =0.0915$. The red and blue lines are bands
from $K_{s}$ and $K_{s}^{\prime }$ valleys, respectively. We have retained a
shell (containing 37 MBZs up to the third Moir\'{e} reciprocal lattice
vector shell) around the valley to construct the basis set. (d) Moir\'{e}
bands at commensurate twist $\protect\theta =5.086^{\circ }$ based on the
generalized one-valley continuum model, with $V_{0}=6E_{R}$ and $\Omega
=0.116$. The green thin dashed lines are Moir\'{e} bands obtained using the
full tight-binding basis.}
\label{fig:S2}
\end{figure}

With these coupling coefficients, we can diagonalize the single particle
Hamiltonian $H_{0}=H_{\uparrow }+H_{\downarrow }+H_{\uparrow \downarrow }$.
In the Moir\'{e} Bloch eigenbasis $\{\psi _{j\mathbf{q}}\}$, it reads
\begin{equation}
H_{0}=\sum_{j,\mathbf{q}}E_{j\mathbf{q}}\beta _{j\mathbf{q}}^{\dag }\beta _{j%
\mathbf{q}}.
\end{equation}%
The typical distributions of $\psi _{j\mathbf{q}}(\mathbf{r},s)$ are shown
in Figs.~\ref{fig:S2}a and \ref{fig:S2}b. Atoms are mainly distributed
around the $AA$ region for all $\mathbf{q}$ except a small area near $%
\mathbf{q}=\Gamma $, where atoms become more uniformly distributed. As a
result, the interaction and thereby the pairing is weak at $\Gamma $ point.

The above results and the Hamiltonian Eq.~(1) in the main text apply for any
commensurate twist angles. For incommensurate twist angles, there are no
well-defined Moir\'{e} patterns and Moir\'{e} bands if the twist angle is
too large. However, if the twist angle is small enough, Moir\'{e} patterns
can form even for incommensurate twist. In this case, the low-energy physics
is mainly determined by the states around the Dirac valleys, and the
inter-valley coupling is also negligible. Therefore, we can adopt a similar
approach as the continuum model by only keeping $\mathbf{g}_{s}$ around one
valley in the Hamiltonian, and the inter-spin couplings can be obtained
using Eq.~\ref{eq:J}. In Fig.~\ref{fig:S2}c, we plot the Moir\'{e} bands at
a small incommensurate twist angle using the generalized one-valley
continuum model mentioned above. We also plot the Moir\'{e} bands at a small
commensurate twist angle using the one-valley continuum model approach and
compare it with the results based on the full tight-binding basis (see Fig.~%
\ref{fig:S2}d). Their agreement confirms the validation of the one-valley
continuum model at small twist angles.


As we discussed in the main text, the energy scales of the magic physics is
small and at the order of $10^{-3}E_{R}$. Therefore, in the Wannier basis,
long-range (intra- and inter-spin) tunnelings beyond nearest neighbors that
are $\gtrsim 10^{-3}E_{R}$ should be taken into account to obtain correct
magic flat bands~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S,PNAS.108.12233S,PhysRevB.86.155449S}%
. In particular, the inter-spin site separations take various values and are
nearly continuously distributed for small twists, therefore the on-site, NN,
next-NN, etc., inter-spin tunnelings should be considered even when the
tunneling cutoff distance is small. The inter-spin tunnelings are also
highly anisotropic, which depend on both the relative orientation and
distance between the corresponding sites due to the three-fold rotational
symmetry of the Wannier orbitals. A small deviation in the tunnelling
coefficients may result in significant change in the band structures near
the `magic angle' due to the narrow bandwidths and approximate degeneracy of
the flat bands. The simplified real-space tight-binding model in~\cite%
{PhysRevA.100.053604S}, which assumes a simply isotropic Gaussian Wannier
function and includes only the nearest-neighbor intra-spin tunneling and
on-site inter-spin couplings, has significant deviation in determining the
bands at small twist angles.

To elaborate more on the effects of long-range real-space tunnelings, we
first solve for the maximally localized Wannier functions $w_{s,j}(\mathbf{r}%
)$ (with $j$ the site index) of the two spin-dependent lattices, then
calculate the band structures using the real-space tight-binding Hamiltonian
\[
H_{s}=\sum_{j,j^{\prime }}t_{jj^{\prime }}\alpha _{s,j}^{\dag }\alpha
_{s,j^{\prime }}
\]%
and
\[
H_{\uparrow \downarrow }=\sum_{j,j^{\prime }}J_{jj^{\prime }}\alpha
_{\uparrow ,j}^{\dag }\alpha _{\downarrow ,j^{\prime }}+h.c.,
\]%
where $\alpha _{s,j}^{\dag }$ is the creation operator of the Wannier
orbitals, and the tunneling coefficients are $t_{jj^{\prime }}=\int d^{2}%
\mathbf{r}w_{s,j}^{\ast }(\mathbf{r})[-\nabla _{\mathbf{r}}^{2}+V_{s}(%
\mathbf{r})]w_{s,j^{\prime }}(\mathbf{r})$ and $J_{jj^{\prime }}=\int d^{2}%
\mathbf{r}w_{\uparrow ,j}^{\ast }(\mathbf{r})\Omega w_{\downarrow ,j^{\prime
}}(\mathbf{r})$. For the lattice depth considered in this paper, we have $%
(t_{0};t_{1};t_{2};t_{3};t_{4};t_{5})=(-4.3458;-0.1967;0.0166;-0.0047;-0.0025;0.0009)E_{R}
$ for $V_{0}=6E_{R}$ and $%
(t_{0};t_{1};t_{2};t_{3};t_{4};t_{5})=(-2.2498;-0.2351;0.0245;-0.0089;-0.0045;0.002)E_{R}
$ for $V_{0}=4E_{R}$, with $t_{jj^{\prime }}=t_{\tilde{\imath}}$ if site $j$
is the $\tilde{\imath}$-th neighbor of site $j^{\prime }$ (where $t_{0}$ is
the on-site energy). Therefore, the cutoff distance for tunneling should be
up to the 5-th nearest neighbors of the bare hexagonal lattice. The
inter-spin coupling $J_{jj^{\prime }}$ is nearly continuously distributed
which decays exponentially with site separation (the on-site tunneling rate
is close to $\Omega $). In Figs.~\ref{figS:real_space_TB}a and ~\ref%
{figS:real_space_TB}b, we show the band structures for different cutoff
distance. We see that the flat band energy is largely shifted (by $\sim
0.25t_{1}$) for a cutoff only to the 1st nearest neighbors of the bare
hexagonal lattice. Nevertheless, magic flat bands are still observed with a
larger width (which can be reduced further by tuning $\Omega $), indicating
the robustness of the magic-angle physics.

\begin{figure}[b]
\includegraphics[width=0.95\linewidth]{FigS_real_space_TB_r.pdf}
\caption{(a) and (b) Moir\'{e} bands based on real-space tight-binding
Hamiltonian for $V_0=4E_R$ and $V_0=6E_R$, respectively. For cutoff distance
$I_c$, we keep the (intra- and inter-spin) tunnelings with site separation
not greater than the $I_c$-th neighbor separation of the bare hexagonal
lattice (e.g., we keep $t_{\tilde{\imath}}$ with $\tilde{\imath}\leq I_c$). The blue,
green and red solid lines correspond to the cutoff distance up to the $1$st,
$2$nd and $5$-th neighbors of the bare hexagonal lattice. The thin dashed
lines are the results of the momentum-space tight-binding Hamiltonian. Other
parameters are $\protect\theta =5.086^{\circ }$ and $\Omega =\Omega_f$. (c)
Schematic illustration of the pairing (indicated by dashed lines) in the
Brillouin zones, between $\uparrow$ states (green dots at $+\mathbf{k}$)
around valley $K^{\prime }_\uparrow$ and $\downarrow$ states (red dots at $-%
\mathbf{k}$) around valley $K_\downarrow$. The two pairing states are
located at the same side of the corresponding valleys (as indicated by the
black arrows), and thus the relative phase between $A$ and $B$ sublattice
sites (which is related to the chirality of the valley) leads to the
staggered pairing order.}
\label{figS:real_space_TB}
\end{figure}

\subsection{Superfluid orders}

The mean-filed interaction reads
\begin{eqnarray}
H_{\text{int}} &=&U_{0}\int d^{2}\mathbf{r}\hat{\Psi}_{\uparrow }^{\dag }(%
\mathbf{r})\hat{\Psi}_{\downarrow }^{\dag }(\mathbf{r})\hat{\Psi}%
_{\downarrow }(\mathbf{r})\hat{\Psi}_{\uparrow }(\mathbf{r})  \nonumber \\
&=&\int d^{2}\mathbf{r}[\hat{\Psi}_{\uparrow }^{\dag }(\mathbf{r})\hat{\Psi}%
_{\downarrow }^{\dag }(\mathbf{r})\Delta (\mathbf{r})+\hat{\Psi}_{\downarrow
}(\mathbf{r})\hat{\Psi}_{\uparrow }(\mathbf{r})\Delta ^{\ast }(\mathbf{r})-%
\frac{|\Delta (\mathbf{r})|^{2}}{U_{0}}]
\end{eqnarray}%
with local pairing amplitude $\Delta (\mathbf{r})=U_{0}\langle \hat{\Psi}%
_{\downarrow }(\mathbf{r})\hat{\Psi}_{\uparrow }(\mathbf{r})\rangle =\sum_{%
\mathbf{g}}\Delta _{\mathbf{g}}e^{i\mathbf{g}\cdot \mathbf{r}}$. We expand
the field operator in the Moir\'{e} Bloch basis $\hat{\Psi}_{s}(\mathbf{r}%
)=\sum_{j,\mathbf{q}}\beta _{j\mathbf{q}}\psi _{j\mathbf{q}}(\mathbf{r},s)$
and obtain the gap equation
\begin{equation}
\Delta _{\mathbf{g}}=U_{0}\sum_{j^{\prime },j,\mathbf{q}}\chi _{j^{\prime
}j}^{\mathbf{q}}(\mathbf{g})C_{\mathbf{q}}^{j^{\prime }j}[\Delta (\mathbf{r}%
)],  \label{eq:GapS}
\end{equation}%
where $\chi _{j^{\prime }j}^{\mathbf{q}}(\mathbf{g})=\frac{1}{L}\int d^{2}%
\mathbf{r}e^{-i\mathbf{g}\cdot \mathbf{r}}\psi _{j^{\prime }-\mathbf{q}}(%
\mathbf{r},\downarrow )\psi _{j\mathbf{q}}(\mathbf{r},\uparrow )$. The
correlation $C_{\mathbf{q}}^{j^{\prime }j}=\langle \beta _{j^{\prime }-%
\mathbf{q}}\beta _{j\mathbf{q}}\rangle $ can be obtained by solving the
Bogoliubov-de Gennes (BdG) Hamiltonian%
\begin{equation}
H_{\text{BdG}}=\sum_{j,\mathbf{q}}(E_{j\mathbf{q}}-\mu)\beta _{j\mathbf{q}}^{\dag
}\beta _{j\mathbf{q}}+\sum_{j,j^{\prime },\mathbf{q}}[\bar{\Delta}%
_{jj^{\prime }}(\mathbf{q})\beta _{j\mathbf{q}}^{\dag }\beta _{j^{\prime }-%
\mathbf{q}}^{\dag }+h.c.]
\end{equation}%
with $\bar{\Delta}_{jj^{\prime }}(\mathbf{q})=L\sum_{\mathbf{g}}\Delta _{%
\mathbf{g}}\chi _{j^{\prime }j}^{\mathbf{q}\ast }(\mathbf{g})$. We solve the
gap equation self-consistently by retaining only the four flat bands which
have much larger DOS than nearby bands.

As we discussed in the main text, at finite temperature, the relevant
physics is the Berezinskii-Kosterlitz-Thouless (BKT) transition since there
is no long-range superfluid order in two dimensions due to phase
fluctuations, and free vortex-antivortex pairs are formed spontaneously
below a critical temperature $T_{\text{BKT}}$~\cite{PhysRevLett.109.105302S,PhysRevLett.114.110401S}. The critical temperature $%
T_{c}$ by mean-field theory is often overestimated. Nevertheless, the
mean-field critical temperature $T_{c}$ provides an upper bound of $T_{\text{%
BKT}}$. Furthermore, the true BKT transition temperature could be obtained
from the mean-field superfluid weight $\rho _{s}$ which can be calculated
based on current response~\cite{PhysRevLett.112.086401S,julku2018superfluidS}%
. On the Moir\'{e} superlattice scale, our system corresponds to an $s$-wave
superfluid, we introduce a slowly varying superfluid phase $\Delta (\mathbf{r%
})\rightarrow \Delta (\mathbf{r})e^{i\Theta (\mathbf{r},\tau)}=\sum_{\mathbf{g}%
}\Delta _{\mathbf{g}}e^{i\mathbf{g}\cdot \mathbf{r}}e^{i\Theta (\mathbf{r}%
,\tau)}$, and the phase $\Theta (\mathbf{r},\tau)=\mathbf{Q}\cdot \mathbf{r}$
(with $\mathbf{Q}\rightarrow 0$) leads to a supercurrent $\mathbf{j}%
_{s}=\rho _{s}\mathbf{Q}$. The supercurrent can also be calculated using the
variation of the free energy $\mathbf{j}_{s}=\frac{1}{L}\langle \frac{\delta
H_{\text{BdG}}}{\delta \mathbf{Q}}\rangle $, the superfluid weight is
obtained by $\rho _{s}=\mathbf{j}_{s}/\mathbf{Q}$, and the BKT transition
temperature is given by $T_{\text{BKT}}=\frac{\pi }{2}\rho _{s}(T_{\text{BKT}%
})$~\cite{PhysRevLett.112.086401S,julku2018superfluidS}. Here the BdG
Hamiltonian reads
\begin{equation}
H_{\text{BdG}}(\mu ,\Delta ,\mathbf{Q})=\sum_{j,\mathbf{q}}(E_{j\mathbf{q}%
}-\mu)\beta _{j\mathbf{q}}^{\dag }\beta _{j\mathbf{q}}+\sum_{j,j^{\prime },%
\mathbf{q}}[\bar{\Delta}_{jj^{\prime }}(\mathbf{q,Q})\beta _{j\mathbf{q+Q/2}%
}^{\dag }\beta _{j^{\prime }\mathbf{-q+Q/2}}^{\dag }+h.c.]
\end{equation}%
where $\bar{\Delta}_{jj^{\prime }}(\mathbf{q,Q})=L\sum_{\mathbf{g}}\Delta _{%
\mathbf{g}}\chi _{j^{\prime }j}^{\mathbf{q,Q}\ast }(\mathbf{g})$ and $\chi
_{j^{\prime }j}^{\mathbf{q,Q}}(\mathbf{g})=\frac{1}{L}\int d^{2}\mathbf{r}%
e^{-i\mathbf{g}\cdot \mathbf{r}}\psi _{j^{\prime }\mathbf{-q+Q/2}}(\mathbf{r}%
,\downarrow )\psi _{j\mathbf{q+Q/2}}(\mathbf{r},\uparrow )$. The gap
equation is $\Delta _{\mathbf{g}}=U_{0}\sum_{j^{\prime },j,\mathbf{q}}\chi
_{j^{\prime }j}^{\mathbf{q,Q}}(\mathbf{g})C_{\mathbf{q,Q}}^{j^{\prime }j}$
with $C_{\mathbf{q,Q}}^{j^{\prime }j}=\langle \beta _{j^{\prime }\mathbf{%
-q+Q/2}}\beta _{j\mathbf{q+Q/2}}\rangle $.

Typically, $\rho _{s}(T)$ is
nearly a constant at low temperature $T$, and therefore we use the following
approximation $T_{\text{BKT}}\simeq \frac{\pi }{2}\rho _{s}(T=0)$. We
numerically calculate $T_{\text{BKT}}$ (the results are shown in Fig.~3a in
the main text) and find that $T_{\text{BKT}}\simeq 0.4T_{c}$. In particular,
we first set $\Theta (\mathbf{r},\tau)=0$ and solve for the order parameter $%
\Delta _{0}(\mathbf{r})$ with fixed chemical potential $\mu $, the
corresponding zero-temperature free energy is $\mathcal{F}(\mathbf{Q}%
=0)\equiv \frac{1}{L}\langle H_{\text{BdG}}(\mu ,\Delta _{0},\mathbf{Q}%
=0)\rangle $. Then we introduce the superfluid phase $\Theta (\mathbf{r},\tau)=%
\mathbf{Q}\cdot \mathbf{r}$ by setting the order parameter as $\Delta _{0}(%
\mathbf{r})e^{i\mathbf{Q}\cdot \mathbf{r}}$, and we calculate the free
energy $\mathcal{F}(\mathbf{Q})\equiv \frac{1}{L}\langle H_{\text{BdG}}(\mu
,\Delta _{0},\mathbf{Q})\rangle $. Since $\mathbf{j}_{s}=\rho _{s}\mathbf{Q}=%
\frac{\delta F(\mathbf{Q})}{\delta \mathbf{Q}}$ for $\mathbf{Q}\rightarrow 0$
at $T=0$, we have $F(\mathbf{Q})-F(0)=\frac{1}{2}\rho _{s}\mathbf{Q}^{2}=%
\frac{1}{2}\rho _{s}(\nabla \Theta )^{2}$ and the superfluid weight is
calculated as $\rho _{s}(T=0)=2\times \frac{F(\mathbf{Q})-F(0)}{\mathbf{Q}%
^{2}}$. The BKT critical temperature is obtained as $T_{\text{BKT}}\simeq
\pi \times \frac{F(\mathbf{Q})-F(0)}{\mathbf{Q}^{2}}$.

It can be verified that, the results based on current response are
consistent with those obtained through the effective action based on quantum
field theory. The partition function can be written as $Z=\text{Tr}%
(e^{-H/T})=\int \mathcal{D}[\bar{\Psi},\Psi ]e^{-S_{\text{eff}}[\bar{\Psi}%
,\Psi ]}$. The effective action is
\begin{equation}
S_{\text{eff}}[\bar{\Psi},\Psi ]=\int_{0}^{1/T}d\tau \left[ \int d^{2}%
\mathbf{r}\sum_{s}\bar{\Psi}_{s}\partial _{\tau }\Psi _{s}+\mathcal{H}(\bar{%
\Psi},\Psi )\right] ,
\end{equation}%
with $\mathcal{H}(\bar{\Psi},\Psi )=\int d^{2}\mathbf{r}\left\{ \sum_{s}[%
\bar{\Psi}_{s}(-\nabla ^{2}+V_{s}-\mu)\Psi _{s}+\bar{\Psi}_{s}\Omega \Psi _{\bar{%
s}}]+U_{0}\bar{\Psi}_{\uparrow }\bar{\Psi}_{\downarrow }\Psi _{\downarrow
}\Psi _{\uparrow }\right\} $ obtained by replacing field operator $\hat{\Psi}%
_{s}^{\dag }(\mathbf{r},t)$ and $\hat{\Psi}_{s}(\mathbf{r},t)$ with Grassman
field number $\bar{\Psi}_{s}$ and $\Psi _{s}$ in the Hamiltonian, and $\tau $
is the imaginary time. We can integrate out the quartic interaction term by
Hubbard-Stratonovich transformation, where the order parameter $\Delta $ is
defined. By further integrating out fermion fields, the partition function
reads $Z=\int \mathcal{D}[\Delta ^{\ast },\Delta ]e^{-S_{\text{eff}}[\Delta
^{\ast },\Delta ]}$ with
\begin{equation}
S_{\text{eff}}[\Delta ^{\ast },\Delta ]=\int_{0}^{1/T}d\tau \int d^{2}%
\mathbf{r}\frac{-|\Delta |^{2}}{U_{0}}-\frac{1}{2}\text{Tr ln}\left( i\omega
_{m}+M_{\text{BdG}}\right) +\frac{1}{T}\sum(\mathcal{E}_{sl\mathbf{k}_{s}}-\mu ).
\end{equation}%
Where $M_{\text{BdG}}$ is the BdG matrix in the Nambu basis, and $\omega _{m}
$ are Matsubara frequencies. For a static phase $\Theta (\mathbf{r},\tau)=%
\mathbf{Q}\cdot \mathbf{r}$, the effective action induced by the phase
fluctuation is
\begin{equation}
S_{\text{fluc}}=\frac{1}{2}\text{Tr ln}[i\omega _{m}+M_{\text{BdG}%
}(0)]-\frac{1}{2}\text{Tr ln}[i\omega _{m}+M_{\text{BdG}}(\mathbf{Q})].
\end{equation}
We can
further sum over Matsubara frequencies and obtain
\begin{equation}
S_{\text{fluc}}=\frac{1}{2%
}\sum_{j,\mathbf{q}}\text{ln cosh}\frac{E_{\text{BdG}}(j,\mathbf{q,0})}{2T}-\frac{%
1}{2}\sum_{j,\mathbf{q}}\text{ln cosh}\frac{E_{\text{BdG}}(j,\mathbf{q,Q})}{2T},
\end{equation}
with $E_{\text{BdG}}(j,\mathbf{q,Q})$ the BdG eigenfrequencies, which can be calculated numerically.
It is hard to proceed further to obtain an analytical
expression of $S_{\text{fluc}} (\mathbf{Q})$ due to the complicated setup.
Instead, we
assume  $S_{\text{fluc}}$ takes the form
\begin{equation}
S_{\text{fluc}%
}=\int_{0}^{1/T}d\tau \int d^{2}\mathbf{r}\frac{1}{2}\rho _{s}(\nabla\Theta)^{2}=\int_{0}^{1/T}d\tau \int d^{2}\mathbf{r}\frac{1}{2}\rho _{s}\mathbf{Q}^{2}
\end{equation}
for $\mathbf{Q}\rightarrow 0$ and a static phase $\Theta (\mathbf{r},\tau)=%
\mathbf{Q}\cdot \mathbf{r}$~\cite{PhysRevLett.109.105302S,PhysRevLett.114.110401S},
leading to $\rho _{s}=\frac{2T}{
\mathbf{Q}^{2}L}S_{\text{fluc}}$.
In the low temperature limit $T\rightarrow 0
$, we have
\begin{equation}\lim_{T\rightarrow 0}\frac{T}{L}S_{\text{fluc}}=\frac{1}{2L}\sum_{j,\mathbf{q}%
,E<0}[E_{\text{BdG}}(j,\mathbf{q,Q})-E_{\text{BdG}}(j,\mathbf{q,0})]=F(%
\mathbf{Q})-F(0).\end{equation}
Therefore, the results based on current response are
consistent with those obtained through the effective action. We have
numerically verified that, $\rho _{s}(T)=\frac{2T}{L\mathbf{Q}^{2}}S_{\text{fluc%
}}$ is nearly a constant at low temperature with $\rho _{s}(T_{\text{BKT}})\simeq\rho _{s}(T=0)$,
and $T_{\text{BKT}}=%
\frac{\pi }{2}\rho _{s}(T_{\text{BKT}})\simeq \frac{\pi }{2}\rho _{s}(T=0)$
is a very good approximation.


The attractive $s$-wave interaction pairs atoms from opposite valleys of the
bare hexagonal lattice and the intra-sublattice pairing is dominant. Note
that the pairing is between Moir\'{e} momentum $\pm \mathbf{q}$, which is
mainly determined by the bare Bloch states $\phi _{\uparrow l\mathbf{k}}$
and $\phi _{\downarrow l-\mathbf{k}}$ nearest to the valleys (which
contribute most to the flat bands). Due to the relative twist, $\pm \mathbf{k%
}$ are located at the same side of the corresponding valleys (as illustrated
in Fig.~\ref{figS:real_space_TB}c). Therefore, we have $\phi _{\uparrow l%
\mathbf{k}}\propto \lbrack 1,e^{i\gamma _{\uparrow \mathbf{k}}}]^{T}$ and $%
\phi _{\downarrow l-\mathbf{k}}\propto \lbrack 1,e^{i\gamma _{\downarrow -%
\mathbf{k}}}]^{T}$ in the $A$ and $B$ sublattice basis, with the relative
phases $\gamma _{\uparrow \mathbf{k}}\simeq -\gamma _{\downarrow -\mathbf{k}%
}+\pi $ (which are related to the chirality of the valleys, i.e., the Berry
phase on loops surrounding the valley). The relative phase $\gamma _{s%
\mathbf{k}}$ between $A$ and $B$ sublattice sites are responsible for the
Larkin-Ovchinnikov phase with staggered pairing order $\Delta (\mathbf{r}%
)\propto \langle \phi _{\uparrow l\mathbf{k}}\phi _{\downarrow l-\mathbf{k}%
}\rangle \propto \lbrack 1,-1]^{T}$.

\begin{figure}[tb]
\includegraphics[width=1.0\linewidth]{FigS3_r.pdf}
\caption{(a)-(c) The correlations $C_{\mathbf{q}}^{j^{\prime }j}$ for
different $j^{\prime },j$ with $U_{0}=-0.5$. (d) The superfluid band
structures (two middle particle/hole bands are not shown). (e) Superfluid
order $\Delta (\mathbf{r}=0)$ as a function of $U_{0}$ obtained by retaining
12 bands (circles) and 4 bands (plus signs) with $\protect\mu =0.00485$. (f)
Superfluid order $\Delta (\mathbf{r})$ at $U_{0}=-0.5$. Common parameters: $%
\protect\theta =5.086^{\circ }$, $V_{0}=6E_{R}$ and $\Omega =0.116$.}
\label{fig:S3}
\end{figure}

\begin{figure}[b]
\includegraphics[width=1.0\linewidth]{FigS4.pdf}
\caption{(a) and (b) Moir\'{e} bands for $AB$ stacking with twist axis at
one coinciding site (where $A$ site of $V_{\uparrow }$ coincides with $B$
site of $V_{\downarrow }$). All other parameters in (a) and (b) are the same
as Figs.~2a and 2b in the main text, respectively. The bands are almost
identical with that in Figs.~2a and 2b in the main text, except that the two
valence (conduction) flat bands become degenerate at $K$ and $K^{\prime }$
points and a tiny gap develops between the valence and conduction flat
bands. (c) The normalized pairing order $\Delta (\mathbf{r})$ (maximum is
normalized to 1) starting from $AA$ stacking with twist axis at one hexagon
center. Other parameters are $\protect\theta =5.086^{\circ }$, $V_{0}=6E_{R}$%
, $\Omega =\Omega _{f}=0.116$ and $U_{0}=-0.5$.}
\label{fig:S4}
\end{figure}

The correlations $C_{\mathbf{q}}^{j^{\prime }j}=\langle \beta _{j^{\prime }-%
\mathbf{q}}\beta _{j\mathbf{q}}\rangle $ for different $j,j^{\prime }$ are
shown in Figs.~\ref{fig:S3}a-c. The pairing is mainly between bands from
different valleys, and the conduction bands from different valleys become
degenerate along the high symmetric $\Gamma $-$K$ lines with avoided
crossing (a tiny gap) due to inter-valley couplings. As we discussed in the
main text, $C_{\mathbf{q}}^{11}$ changes from characterizing $K_{s}$-$K_{%
\bar{s}}^{\prime }$ to $K_{s}^{\prime }$-$K_{\bar{s}}$ correlations;
therefore, $C_{\mathbf{q}}^{11}$ changes the sign across the $\Gamma $-$K$
lines where $C_{\mathbf{q}}^{12}$ is mainly distributed.
This means that though $C_{\mathbf{q}}^{j^{\prime }j}$ show $f$-wave
structures, their combined effects lead to the nearly uniform superfluid gap
(see Fig.~\ref{fig:S3}d) and the pairing is $s$-wave. The small superfluid
gap (weak pairing) at $\Gamma$ is due to the nearly uniform Moir\'{e}
wavefunction there (see Fig.~\ref{fig:S2}b). Moreover, we also calculate the
results by including 8 nearby Moir\'{e} bands (12 bands in total with 4 flat
bands, 4 higher and 4 lower bands) for comparison, and find that these
nearby bands have very little effects on the phase diagram and superfluid
order (as shown in Figs.~\ref{fig:S3}e and \ref{fig:S3}f).

In the above mean-field approach, we have assumed a real-space pairing
order. Alternatively, we can assume the momentum-space pairing order. In
particular, we first write the interaction Hamiltonian as
\[
H_{\text{int}}=\sum U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{Q;q,q^{\prime }}%
}\beta _{j_{1}\mathbf{Q+q}}^{\dag }\beta _{j_{2}\mathbf{Q-q}}^{\dag }\beta
_{j_{3}\mathbf{Q-q^{\prime }}}\beta _{j_{4}\mathbf{Q+q^{\prime }}}
\]%
with $U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{Q;q,q^{\prime }}}=U_{0}\int d^{2}%
\mathbf{r}\psi _{j_{1}\mathbf{Q+q}}^{\ast }(\mathbf{r},\uparrow )\psi _{j_{2}%
\mathbf{Q-q}}^{\ast }(\mathbf{r},\downarrow )\psi _{j_{3}\mathbf{Q-q^{\prime
}}}(\mathbf{r},\downarrow )\psi _{j_{4}\mathbf{Q+q^{\prime }}}(\mathbf{r}%
,\uparrow )$. Here $\mathbf{Q}$, $\mathbf{q}$ and $\mathbf{q^{\prime }}$ are
superlattice momenta in the Moir\'{e} BZ. We restrict the interaction to the
$\mathbf{Q}=0$ BCS channel, and assume the momentum-space order $\bar{\Delta}%
_{j_{1}j_{2}}(\mathbf{q})=\sum_{j_{3},j_{4},\mathbf{q}^{\prime
}}U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{0;q,q^{\prime }}}C_{\mathbf{q}^{\prime
}}^{j_{3}j_{4}}$. The correlation $C_{\mathbf{q}^{\prime }}^{j_{3}j_{4}}$
can be obtained by solving the BdG Hamiltonian%
\begin{equation}
H_{\text{BdG}}=\sum_{j,\mathbf{q}}(E_{j\mathbf{q}}-\mu)\beta _{j\mathbf{q}}^{\dag
}\beta _{j\mathbf{q}}+\sum_{j_{1},j_{2},\mathbf{q}}[\bar{\Delta}%
_{j_{1}j_{2}}(\mathbf{q})\beta _{j_{1}\mathbf{q}}^{\dag }\beta _{j_{2}-%
\mathbf{q}}^{\dag }+h.c.],
\end{equation}%
which allows us to obtain the superfluid order $\bar{\Delta}_{j_{1}j_{2}}(%
\mathbf{q})$ self-consistently. Using this approach, we calculate superfluid
order by keeping only the four flat bands and find that the orders $\bar{%
\Delta}_{j_{1}j_{2}}(\mathbf{q})$ are the same (up to tiny numeric errors)
as those obtained by assuming a real-space pairing order [i.e., $%
\sum_{j_{3},j_{4},\mathbf{q}^{\prime }}U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{%
0;q,q^{\prime }}}C_{\mathbf{q}^{\prime }}^{j_{3}j_{4}}\simeq L\sum_{\mathbf{g%
}}\Delta _{\mathbf{g}}\chi _{j_{1}j_{2}}^{\mathbf{q}\ast }(\mathbf{g})$].
The two approaches lead to the same superfluid phase, correlation $C_{%
\mathbf{q}}^{j_{3}j_{4}}$, as well as $\Delta (\mathbf{r})$ (which is
determined by $C_{\mathbf{q}}^{j_{3}j_{4}}$).

\subsection{Effects of different stackings and twist axes}

We have focused on the twists starting from $AA$ stacking with the twist
axis at one sublattice site. Like the magic behaviors in TBG~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S,PNAS.108.12233S,PhysRevB.86.155449S}%
, here the twist axis or stacking position do not affect the appearance of
magic flat bands at small twist angles, as shown in Figs.~\ref{fig:S4}a and %
\ref{fig:S4}b. When $\theta $ is small, the Moir\'{e} bands for different
stackings and twist axes are almost identical. For $AB$ stacking with twist
axis at one coinciding site (where $A$ site of $V_{\uparrow }$ coincides
with $B$ site of $V_{\downarrow }$), the two valence (conduction) flat bands
become degenerate at $K$ and $K^{\prime }$ points and a tiny gap develops
between the valence and conduction flat bands. We find that the Moir\'{e}
bands for $AB$ ($AA$) stacking with twist axis at one hexagon center is the
same as that for $AA$ ($AB$) stacking with twist axis at one coinciding
site. For different stackings and twist axes, the superfluid orders and the
phase diagrams are similar, where the pairing order is staggered and
distributed mainly around the $AA$ region (see Fig.~\ref{fig:S4}c).

\subsection{Difference with twisted bilayer graphene}

In this section, we reemphasize some of the differences between our system
and the TBG system.

(1) The two twisted lattice potentials are state dependent, and one
potential does not affect atoms trapped by the other. This is different from
the electrons in TBG, where electrons in one layer can feel the potential of
atoms in the other layer.

(2) Our system is physically a single-layer system and we twist the lattice
for atomic (pseudo-)spin states (i.e., atomic internal energy levels). The $%
z $-direction is tightly confined by an additional state-independent
potential using the so-called magic-wavelength lasers. Therefore, the two
spin states have identical Wannier orbital along the $z$ direction. The
inter-spin tunnelings, realized by additional lasers, are different from the
inter-layer tunnelings between $p_{z}$ orbitals~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S} in TBG where a large
inter-layer distance exist. The existence of magic behaviors in our system
is not a straightforward derivative of TBG.

(3) The optical lattice potential here takes a cosine form which is much
simpler comparing to the atomic potential in graphene. Therefore the bare
bands and inter-spin couplings can be obtained accurately by directly
solving for the Bloch states $\phi _{sl\mathbf{k}_{s}}$ in our system. While
for the TBG, real-space tight-binding approximation based on Slater-Koster
parameters is usually adopted~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S}.

(4) Long-range tunnelings are more significant in our system because the
optical lattices considered here are relatively shallow. A shallow lattice
not only improves the atomic lifetime (through reducing the atomic decay
rate), but also increases the bare Dirac velocity $v_{D}$ (a larger $v_{D}$
leads to larger gaps and better isolation of the flat bands).

(5) The interactions are dominated by the $s$-wave scattering between atoms
in different spin states that are coupled and relatively twisted. In TBG,
the electronic interactions are more complex and include both Coulomb
repulsive interaction and/or phonon-mediated attractive interactions, which
mainly involve electrons in the same layer with no relative twist~\cite%
{PhysRevB.98.220504S,PhysRevLett.122.257002S,PhysRevLett.121.257001S,
PhysRevX.8.041041S,PhysRevLett.121.217001S,PhysRevLett.122.026801S}. The
unique interaction in our system can lead to interesting Larkin-Ovchinnikov
superfluid orders that do not exist in TBG.

(6) Finally, the advantage of cold atom system is that the parameters (e.g.,
inter-spin tunnelings, lattice depth, lattice constant, interactions, etc.)
are highly tunable. This not only leads to magic behaviors in a wide range
of parameter space, but also opens various possibilities for exploring novel
twistronics in cold atom systems.

