#  Topological superfluids with time-reversal symmetry from $s$-wave
interaction in a bilayer system

Beibing Huang∗ Department of Physics and Center for quantum Coherence, The
Chinese University of Hong Kong Shatin, N.T., Hong Kong, China  Department of
Physics, Yancheng Institute of Technology, Yancheng, 224051, P. R. China  Pak
Hong Chui∗ Department of Physics and Center for quantum Coherence, The Chinese
University of Hong Kong Shatin, N.T., Hong Kong, China  Jia Liu  Department of
Physics and Center for quantum Coherence, The Chinese University of Hong Kong
Shatin, N.T., Hong Kong, China  Chuanwei Zhang  Department of Physics, The
University of Texas at Dallas, Richardson, TX 75080, USA  Ming Gong
Department of Physics and Center for quantum Coherence, The Chinese University
of Hong Kong Shatin, N.T., Hong Kong, China

###### Abstract

Topological superconducting phases with time-reversal (TR) symmetry have been
widely explored in recent years. However the involved unconventional pairings
are generally implausible in realistic materials. Here we demonstrate via
detailed self-consistent calculation that these topological phases with TR
symmetry in DIII and BDI classes can be realized in a spin-orbit coupled
bilayer system with only $s$-wave interaction. The bilayer freedom enables the
definition of TR symmetry between the layers even in the presence of local
Zeeman fields, which we propose to be realized using four laser beams. The
gapped phase in DIII class is characterized by $\mathbb{Z}_{2}$, while all the
gapless phases in these two classes are characterized by nontrivial winding
numbers and are also manifested from the Majorana flat bands. We also unveil
the intimate relation between TR symmetry and mirror symmetry due to phase
locking effect between the two layers, which harbors the mirror symmetry
protected topological phases. We finally demonstrate that these phases will
not be spoiled by interlayer pairings.

There is a major effort in realizing unpaired Majorana fermions (MFs) in both
condensed matter physics and cold atom physics, which have important
applications in topological quantum computationsNayak . These systems include
$p$-wave superconducting phases in Sr2RuO4SRO ; SRO2 , Helium-3HeB  and fully
spin polarized ultracold atomsColdP1 ; ColdP2 , and the $\nu=5/2$ fractional
quantum Hall stateNayak ; ReadGreen . Unfortunately, whether these systems can
be used to realize the elusive MFs is still a controversial issue due to the
lacking of convincing evidences. The conceptional breakthrough in recent years
is based on the idea that these exotic quantum phases can be realized using
the conventional $s$-wave interaction, spin-orbit coupling (SOC) and Zeeman
fieldFu ; DSarma1 ; DSarma2 ; Alicea . Very recently some pioneering
experiments both in one-dimensional (1D) nanowires and iron chain on the
realization of MFs has been carried outExp1 ; Exp2 ; Exp3 ; Exp4 ; Exp5 .
Remarkably, the localized wave functions at the two ends with zero eigenenergy
have been identifiedExp5 , which provide important evidence for the
realization of MFs. These systems belong to topological D class from the
Cartan classification due to absence of time-reversal (TR) symmetrySchnyder08
; Kitaev . In all these systems only the topological protected gapped phase
can be realized, while the topological phase transition is characterized
either by the Chern number of the occupied bands, or somewhat equivalently the
Pfaffian (Pf) of the Hamiltonian at the particle-hole (PH) invariant
point(s)Kitaev1dd .

Figure 1: Bilayer system for TR invariant superfluids. (b) The bilayer
structure can be realized by imposing a strong double-well potential along the
$z$ direction. The possibility of individual tunability of the local Zeeman
fields, spin-dependent tunnelings and the SOC coefficients can be used to
realize different topological models with TR symmetry. (b) Realization of
local Zeeman fields with opposite sign in the bilayer system for $T_{xy}$
model via a spin-dependent lattice; see more details in Methods.

These progress strongly stimulates the persuit of topological superconducting
phases with TR symmetry. However, these fantastical phases generally involve
some unconventional $p_{x}+ip_{y}$Qi09 , $d_{x^{2}-y^{2}}$FuL ; Sato10 ;
Nakosai12 ; Nakosai ; KTLaw  and $s_{\pm}$shusa ; zhangfan2 ; zhangfan1
pairings, which are implausible in realistic materials. Therefore we have the
intriguing question that whether there phases can exist with only $s$-wave
interaction? Here, instead of seeking for other new materials, we explore
these phases in a synthetic bilayer system where the bilayer freedom enables
the definition of nonlocal TR symmetry between the layers. While the single
layer physics is well-known in some sense, we show that the interlayer
tunnelings and TR symmetry can lead to fundamentally different superconducting
phases. We explore the topological protected gapped and gapless phases in both
DIII and BDI classes in this platform, in which the fully gapped phase in DIII
phase is characterized by $\mathbb{Z}_{2}$, while all the gapless phases are
characterized by nontrivial winding numbers yxzhao ; beri  around the
topological defects — the Dirac cones — in momentum space and are also
manifested from the Majorana flat bands (MFBs)satofb . We have proposed a
scheme to realize some of these phases in experiments. These results are
confirmed using detailed self-consistent calculation, in which the order
parameters and chemical potential are determined by the extrema of the free
energy. We also unveil the intimate relation between TR and mirror symmetries
in our models due to phase locking effect between the pairings in the two
layers, from which we realize some mirror symmetry protected gapped and
gapless phases. We finally demonstrate that these phases will not be spoiled
by the interlayer pairings. We expect these new phases to greatly enrich our
understanding of topological superconducting phases and related transitions.

Results

Physical Model. We first briefly discuss the topological superfluids in a
single layer system, which has been intensively studied in recent yearsMGong11
; MGong12 ; HH11 ; ZQY ; He12 . In this model the Zeeman field is used to
break the Kramers degeneracy at zero momentum, so that the chemical potential
can fill just one band, around which an effective $p$-wave pairing can be
induced. We found that the topological phase can be realized in a narrow
parameter regime even at finite temperature though the Zeeman field is (much)
larger than the pairing strengthMGong11 ; MGong12 , so the induced $p$-wave
pairing provides a new mechanism to go beyond the Chandrasekhar-Clogston
limitCC1 ; CC2 . The recent interests in this system also include the
topological inhomogeneous superfluidsFF1 ; FF2 ; huhuiff , the quench dynamics
of the topological superfluidsDongY ; Victor  and the topological superradiant
stateJian . Thus the exploring of these exotic phases in this model is still
an important topic in modern research.

To explore an experimental realization of the topological phases with TR
symmetry, we consider a two-dimensional (2D) bilayer system, see the schematic
structure in Fig. 1, which can be realized by imposing a strong double-well
potential along the $z$ directionDW1 ; DW2 . This synthetic structure has been
widely used in recent literatures for realizing some exotic phases in solid
systemsBL1 ; BL2 . We first consider the case that both layers have identical
SOC coefficients,
$V_{\text{so}}^{i}=(s_{i}k_{y}\sigma_{x}-\lambda_{i}k_{x}\sigma_{y})$, where
$\lambda_{i}=\alpha^{\prime}$ and $s_{i}=\alpha$, $i=1,2$, are the SOC
strength and $\sigma_{\alpha}$ are the Pauli matrices, and identical spin-
dependent tunnelings $t_{\uparrow,\downarrow}=t$, but with different Zeeman
fields $\Gamma_{i}$. The SOC in ultracold atoms both for fermions and bosons
has been extensively explored using two coupled Raman beamsLin11 ; Zhang12 ;
Wang12 ; Cheuk12 ; Qu13 ; Chris14 , and very recently, the required 2D SOC has
been realized in 40KZhang15 . The barrier width between the two layers can be
controlled by the laser intensity in experiments thus we can safely assume
that only the intralayer pairings, which can be defined as
$\Delta_{i}=\sum_{\bf k}g\langle c_{i{\bf k}\uparrow}c_{i-{\bf
k}\downarrow}\rangle$, are important. Here $g$ is the effective 2D interaction
strength, $c_{i{\bf k}s}$ is the annihilate operator for particle with
momentum k and spin $s$. The corresponding Bogoliubov-de Gennes (BdG) equation
can be written as

$\mathcal{H}=\begin{pmatrix}\mathcal{H}_{0}({\bf k})&\bar{\Delta}\\\
\bar{\Delta}^{\dagger}&-\mathcal{H}_{0}^{*}(-{\bf k})\end{pmatrix},$ (1)

with $\mathcal{H}_{0}({\bf k})$ being the single particle Hamiltonian

$\mathcal{H}_{0}({\bf k})=\begin{pmatrix}\epsilon_{{\bf
k}}+\Gamma_{1}&\rho_{1{\bf k}}&-t_{\uparrow}&0\\\ \rho_{1{\bf
k}}^{*}&\epsilon_{{\bf k}}-\Gamma_{1}&0&-t_{\downarrow}\\\
-t_{\uparrow}&0&\epsilon_{{\bf k}}+\Gamma_{2}&\rho_{2{\bf k}}\\\
0&-t_{\downarrow}&\rho_{2{\bf k}}^{*}&\epsilon_{{\bf
k}}-\Gamma_{2}\end{pmatrix}.$ (2)

Here $\epsilon_{{\bf k}}={\bf k}^{2}/2m-\mu$, and $\rho_{i{\bf
k}}=s_{i}k_{y}+i\lambda_{i}k_{x}$, and the concrete form of $\bar{\Delta}$ can
be found in Methods (Eq. 31). We construct the BdG Hamiltonian in the Nambu
spinor, $\Psi=(c_{1{\bf k}\uparrow},c_{1{\bf k}\downarrow},c_{2{\bf
k}\uparrow},c_{2{\bf
k}\downarrow},c_{1{\bf-k}\uparrow}^{\dagger},c_{1{\bf-k}\downarrow}^{\dagger},c_{2{\bf-k}\uparrow}^{\dagger},c_{2{\bf-k}\downarrow}^{\dagger})^{T}$.
Notice that in ultracold atoms the order parameters and chemical potential
should be solved self-consistently.

We define the TR operator $\mathcal{T}$ for the single particle Hamiltonian
$\mathcal{H}_{0}$ in Eq. 2 and $\mathbb{T}$ for the BdG equation in Eq. 1,
while both $\mathcal{T}^{2}=-1$ and $\mathbb{T}^{2}=-1$ are required for the
topological superfluids in DIII class. The realization of topological BDI
class superfluids with $\mathcal{T}^{2}=+1$ and $\mathbb{T}^{2}=+1$ will be
discussed later. Obviously, without the Zeeman fields, this model has a
nonlocal TR operator $\mathcal{T}=\rho_{x}\otimes\sigma_{y}K$, where
$\rho_{\alpha}$ are Pauli matrices acting on the layer space and $K$ is the
complex conjugate operator. In this TR operator, the spin up (down) in one
layer with momentum ${\bf k}$ is mapped to spin down (up) in the other layer
with opposite momentum. This symmetry is still respected by the Zeeman fields
when $\Gamma_{1}=-\Gamma_{2}=\Gamma$, which is the basic reason for the local
Zeeman fields discussed before. These local Zeeman fields may be challenging
to be realized in solid state system, however, in ultracold atoms — in the
ideal condition — it can be realized using the following spin-dependent
optical lattice,

$V(z)=V_{\text{e}}(z)+V_{\text{o}}(z)\sigma_{z},$ (3)

where the subscripts e and o represent an even and odd function of $z$,
respectively. The first term describes the spin-independent symmetric double-
well potential, which has been realized in experimentsDW1 ; DW2 , and the
second term means that the two layers have exactly opposite Zeeman fields.
Such a structure can be realized using four counterpropagating laser beams
(see details in Methods). An optical lattice that slightly deviates from the
above scenario can still be used to realize the model with TR symmetry by
fine-tuning the parameters.

The TR symmetry for the BdG equation is defined as
$\mathbb{T}=\text{diag}(\mathcal{T},e^{i\phi}\mathcal{T})$, where $\phi$ is an
arbitrary phase. The TR symmetry requires that
$\mathcal{T}\bar{\Delta}\mathcal{T}^{-1}e^{-i\phi}=\bar{\Delta}$, namely,
$\Delta_{1}=\Delta_{2}^{*}e^{-i\phi}$. The phase $\phi$ is used to compensate
the global phase of the order parameters such that only the relative phase
between the two order parameters $\Delta_{1,2}$ is important. After self-
consistent calculation (see details in Methods) we find that the two order
parameters can be treated as real numbers simultaneously. More precisely, we
find that the two layers can have zero relative phase and
$\Delta_{1}/\Delta_{2}>0$ when the two tunnelings have the same sign;
otherwise, the relative phase is $\pi/2$ and $\Delta_{1}/\Delta_{2}<0$. This
effect is rather robust, and this feature is independent of other parameters
such as binding energy, Zeeman fields etc. We further find that the relative
phase between the two order parameters can only be introduced to the
Hamiltonian by the complex tunneling terms. Hereafter, this effect is called
the phase locking effect. For these reasons, throughout this work, we treat
the two order parameters as real numbers in the numerical simulation. For this
particular model we have

$\mathbb{T}=\Lambda K=\tau_{0}\otimes\mathcal{T},$ (4)

where $\tau_{0}$ is a Paul matrix acting on PH space. The intrinsic PH
operator is defined as $\mathbb{C}=\tau_{x}K$, which ensures that
$\mathbb{C}\mathcal{H}({\bf k})\mathbb{C}^{-1}=-\mathcal{H}(-{\bf k})$.

Hereafter, this model is dubbed as $T_{xy}$ model to distinct itself from the
other models studied in the following.

The $\Lambda$ operator in Eq. 4 is essentially the mirror symmetry operator
since in our model only the $k_{x}$ component contains the imaginary number,
which can change sign upon complex conjugation $K$. Thus we can define the
mirror symmetry operator as $\mathcal{M}_{y}=\Lambda$, where

$\mathcal{M}_{y}\mathcal{H}(k_{x},k_{y})\mathcal{M}_{y}^{-1}=\mathcal{H}(k_{x},-k_{y}).$
(5)

We can also define another mirror operator about $k_{y}$-axis as
$\mathcal{M}_{x}=\mathbb{I}\cdot\mathcal{M}_{y}$, where
$\mathbb{I}=\tau_{z}\otimes\rho_{0}\otimes\sigma_{z}$ is the inversion
symmetry operator,
$\mathbb{I}\mathcal{H}(\textbf{k})\mathbb{I}^{-1}=\mathcal{H}(-\textbf{k})$.
In this case,

$\mathcal{M}_{x}\mathcal{H}(k_{x},k_{y})\mathcal{M}_{x}^{-1}=\mathcal{H}(-k_{x},k_{y}).$
(6)

Notice that both $\mathcal{M}_{x}$ and $\mathcal{M}_{y}$ are unitary hermitian
operators and their eigenvalues are either $-1$ or $+1$ since
$\mathcal{M}_{x,y}^{2}=1$.

These new symmetries can fundamentally change the topological invariant of the
Hamiltonian in the mirror symmetric invariant 1D subspaces along $k_{x}$ and
$k_{y}$ axes. To this end, we can diagonalize the mirror symmetry operators
via a unitary transformation
$U_{\alpha}\mathcal{M}_{\alpha}U_{\alpha}^{-1}=(\bf{I}_{4},-\bf{I}_{4})$,
where $\bf{I}_{4}$ is a $4\times 4$ unity matrix and $\alpha=x,y$. Under this
transformation we have

$U_{\alpha}\mathcal{H}_{\alpha}U_{\alpha}^{-1}=\text{diag}(M_{1}^{\alpha},M_{2}^{\alpha}),$
(7)

where $\mathcal{H}_{x}=\mathcal{H}(0,k_{y})$ and
$\mathcal{H}_{y}=\mathcal{H}(k_{x},0)$. In general $M_{1,2}^{x,y}$ can have
totally different symmetries than the original Hamiltonian $\mathcal{H}$. The
concrete form of these matrices and their topological invariants can be found
in Methods. For the $T_{xy}$ model considered here we find that both
$M_{1,2}^{x,y}$ belong to AIII class, while these two matrices are connected
by TR symmetry and PH symmetry to recover the symmetries of the original
Hamiltonian. In fact the consequence of mirror symmetry for any general
Hamiltonian has been intensively studied in literaturesMS1 ; MS2 ; Ueno ; Ken
by defining
$\mathcal{M}_{ss^{\prime}}^{\alpha}:=\mathcal{M}_{\alpha}\mathbb{T}\mathcal{M}_{\alpha}^{-1}=s\mathbb{T},\mathcal{M}_{\alpha}\mathbb{C}\mathcal{M}_{\alpha}^{-1}=s^{\prime}\Sigma$,
where $s,s^{\prime}=\pm 1$. This quantity,
$\mathcal{M}_{ss^{\prime}}^{\alpha}$, is crucial to determine the possible
symmetries of the mirror symmetric invariant subspaces, which are fully
consistent with our direct calculations above. This method will be adopted to
find all the possible mirror symmetric invariant Hamiltonians and their
topological classes for all the models presented in Table 1.

Figure 2: Topological superfluids for $T_{xy}$ model. (a) The phase diagram in
the parameter space by $\epsilon_{\text{b}}$ and $\Gamma$, in which the order
parameters and chemical potential are determined self-consistently. The number
in each colored regime marks the number of gapless Dirac cones in momentum
space, thus $0$ defines a fully gapped phase. The two solid line means that
the energy gap closing and reopening take place at $k_{x}$ and $k_{y}$ axes.
The dashed line means gap closing and reopening at nonzero nonzero momenta
($k_{x}k_{y}\neq 0$). PS denotes the regime for phase separation, which is
visible when $\epsilon_{\text{b}}>0.3E_{\text{F}}$. Parameters are:
$t=0.6E_{\text{F}}$, $\alpha=1.0E_{\text{F}}$,
$\alpha^{\prime}=0.5E_{\text{F}}$. (b) - (d) are band structures at the
critical boundaries when $\epsilon_{\text{b}}=0.1E_{\text{F}}$.

Topological phases in $T_{xy}$ model. We plot the phase diagram as a function
of binding energy $\epsilon_{\text{b}}$ and Zeeman field $\Gamma$ in Fig. 2a,
in which the number in each colored regime represents the number of robust
Dirac cones. We find that the two fully gapped phases are separated by two
gapless phases. The origin of these gapless phases can be understood using the
following way. Due to the presence of both PH and TR symmetries, we can define
a unitary sublattice symmetry $\mathbb{S}=\mathbb{T}\cdot\mathbb{C}$, which
satisfies $\\{\mathbb{S},\mathcal{H}\\}=0$. We then diagonalize this
sublattice operator via a unitary transformation
$V\mathbb{S}V^{-1}=\text{diag}(\bf{I}_{4},-\bf{I}_{4})$, under which we have

$V\mathcal{H}V^{\dagger}=\begin{pmatrix}0&Q\\\
Q^{\dagger}&0\end{pmatrix},\quad\text{Det}\mathcal{H}=|\text{Det}(Q)|^{2}.$
(8)

The closing of energy gap means that both the real and imaginary part of
Det$(Q)$ should be zero simultaneously, which give the following two
conditions,

$\displaystyle(\bm{\gamma}\cdot{\bf k})^{2}=t^{2},\quad\epsilon_{{\bf
k}}^{2}=\Gamma^{2}-\Delta^{2},$ (9)

where $\bm{\gamma}=(\alpha^{\prime},\alpha)$ and
$\Delta_{1}=\Delta_{2}=\Delta$. The first equation defines an ellipse and the
second equation defines one circle or two circles, depending strongly on the
values of the parameters. The robust gapless phases thus are determined by the
overlaps between the ellipse and the circle(s), which may always have
intersections in some appropriate parameter regimes. A topological explanation
for the robust accidental degeneracy at the Dirac cones and the physical
meaning of Eq. 9 will be presented shortly later. When $t\neq 0$, we find the
first equation prohibits the gap closing and reopening at zero momentum, thus
the topological phase transitions in this model will always take place at
nonzero momentum (${\bf k}\neq 0$). This is in sharp contrast to the results
in topological D class superconducting phasesExp1 ; Exp2 ; Exp3 ; Exp4 ; Exp5
; MGong11 ; MGong12  where the critical boundary is only determined by the gap
closing and reopening at ${\bf k}=0$. When $t=0$, the critical boundary is
reduced to,

${\bf k}=0,\quad\Gamma^{2}=\mu^{2}+\Delta^{2},$ (10)

which is well-known in previous literaturesExp1 ; Exp2 ; Exp3 ; Exp4 ; Exp5 ;
MGong11 ; MGong12 .

Figure 3: Phase separation (PS) in $T_{xy}$ model. (a) The evolution of order
parameters as a function of Zeeman field and binding energy. The regime for
sudden jump of order parameters is denoted by the two dashed lines, which
happens when $\epsilon_{\text{b}}>0.3E_{\text{F}}$. (b) The origin of the
sudden jump of order parameters from the perspective of free energy $F$ using
$\epsilon_{\text{b}}=0.8E_{\text{F}}$. (c) The mixed free energy
$F_{\text{m}}$ (defined in Eq. 18) can have slightly lower energy than $F$ by
mixing of two different superfluids. (d) The mixing coefficient $x$ as a
function of Zeeman field. For parameters see Fig. 2.

The phase diagram in Fig. 2a can then be understood from how the ellipse
intersects with the circle(s). Starting from the fully gapped trivial phase at
$\Gamma\sim 0$ and gradually increases the Zeeman fields, we find that the
ellipse will first intersect with two circles in the small binding energy
regime ($\epsilon_{\text{b}}<0.3E_{\text{F}}$), in which we observe a
topological gapless phase with 8 Dirac cones. This topological phase
transition will not take place at the $k_{x}$ and $k_{y}$ axes (see Fig. 2b).
With the further increasing of Zeeman fields, one of the circles will
disappear when $|\mu|<\sqrt{\Gamma^{2}-\Delta^{2}}$ and the system transits to
the topological gapless phase with 4 Dirac cones when $\epsilon_{{\bf
k}=(t/\alpha,0)}^{2}\epsilon_{{\bf
k}=(0,t/\alpha^{\prime})}^{2}<(\Gamma^{2}-\Delta^{2})^{2}$ with gap closing at
the $k_{y}$ axis for the particular parameters used here (see Fig. 2c).
Finally, this phase will evolve to the fully gapped phase with gap closing and
reopening along the $k_{x}$ axis (see Fig. 2d). In the large binding energy
regime, however, only one circle is allowed in the second equation of Eq. 9 so
the fully gapped trivial phase directly enters the phase with 4 Dirac cones.
However, this transition will be accompanied by a sudden jump of the order
parameter during the topological phase transition, see Fig. 3a. The
corresponding jump of the chemical potential is not shown. To understand this
effect we calculate the free energy $F$ as a function of order parameters in
Fig. 3b, in which two degenerate local minima in $F$ are clearly shown in the
critical boundary. The interplay of these two minima gives rise to the jump of
parameters in Fig. 3a. The jump of these parameters is an important signal for
phase separation (PS)caldas ; WYI11 . Fortunately using the method in Refs.
[caldas, ; WYI11, ] (see Eq. 18 in Methods) we find that the PS phase can only
be survived in a narrow regime near the boundary between trivial gapped phase
and the gapless phase with 4 Dirac cones.

Figure 4: Pictures for topological gapped superfluids and MFB. (a) Topological
gapped superfluids with $\mathbb{Z}_{2}=-1$ can be realized when the chemical
potential $\mu$ fills the lowest two bands, thus only one TR invariant point
is enclosed by the Fermi surface. This picture holds in the weak pairing
limit. (b) The MFB is used to connect the two defects with opposite
chiralities ($\pm 1$), which can be computed using $\mathcal{W}({\bf
k}_{\parallel})$ in Eq. 27. The MFB is shrunk to a single point when the line
connecting the two defects perpendiculars to the strip direction.

Next we try to characterize the topology of each phase. We first focus on the
fully gapped topological phase at strong Zeeman field, which is characterized
by $\mathbb{Z}_{2}$ from the Cartan classificationSchnyder08 ; Kitaev . We
employ the method developed by Qi et alQi10 ,
$\mathcal{N}_{\text{2d}}=\prod_{j}\text{sgn}(\delta_{j{\bf k}})^{p_{j}}$,
where $\delta_{j{\bf k}}=\langle
j\text{k}|i\mathcal{U}\Delta|j\text{k}\rangle$ is the pairing gap at the $j$th
Fermi surface with $|j\textbf{k}\rangle$ being the eigenfunction of
$\mathcal{H}_{0}$, and $p_{j}=1$ if the TR invariant point at ${\bf k}=0$ is
enclosed by the Fermi surface and $0$ otherwise. In this formula,
$\mathcal{U}=\rho_{x}\otimes\sigma_{y}$ and
$i\mathcal{U}\bar{\Delta}=\rho_{x}\otimes{\bf I}_{2}$ is a Hermitian operator
thus $\delta_{j{\bf k}}\in\mathbb{R}$. Notice that for a relative large
pairing strength in this regime, we can always adiabatically deform this state
to the weak pairing state near zero binding energy or large Zeeman field
without closing the energy gap, thus this criteria can always be applied in
this work to characterize the topology of this fully gapped phase. A lengthy
but straightforward calculation shows that $\delta_{j{\bf
k}}=4\Delta(|t|\pm|\bm{\gamma}\cdot{\bf k}|)/(\Gamma\pm\Upsilon_{s})$, where
$\pm$ defines the upper and lower two bands and
$\Upsilon_{s}=\sqrt{\Gamma^{2}+(|t|+s|\bm{\gamma}\cdot{\bf k}|)^{2}}\neq 0$.
We immediately see from this result that $\delta_{j{\bf k}}$ can change sign
only when the energy gap is closed at $t\pm|\bm{\gamma}\cdot{\bf k}|=0$, which
is exactly the first equation in Eq. 9. The second equation therefore
determines the position of Fermi surface when $\Delta=0$. Notice that the sign
of $\delta_{j{\bf k}}$ is fixed in each Fermi surface. Then we find that in
this regime, $\mathcal{N}_{\text{2d}}=-1$. In this regime the chemical
potential just fills the lower two bands while the upper two bands are
unoccupied (see Fig. 4a). The same method yields $\mathcal{N}_{\text{2d}}=+1$
in the fully gapped trivial phase with small Zeeman field.

Figure 5: Edge states in $T_{xy}$ model along different directions. Edge
states for a strip with $L=200/k_{\text{F}}$ along $k_{x}$ (left), $[110]$
(middle) and $k_{y}$ (right) directions in the fully gapped topological phase
(a), gapless phase with 4 (b) and 8 (c) Dirac cones for the three points
marked by asterisk (*) in Fig. 2. For parameters see Fig. 2.

This topological gapped phase can support edge states for a strip along
different directions; see Fig. 5a. Note that this gapped phase along $k_{x}$
and $k_{y}$ axes may also be protected by mirror symmetry. For this reason, we
calculate the topological invariant — the winding number in 1D — in the mirror
symmetric invariant subspaces $M_{1,2}^{\alpha}$, which belong to AIII class
with only sublattice symmetry. This mirror winding number
$\mathcal{W}^{\alpha}$ is calculated using Eq. 27 (see Methods), which are
$\mathcal{W}^{x}=1$ and $\mathcal{W}^{y}=-1$. However, this mirror symmetry
can not protect the edge states along other directions.

The topological gapless phases emerged in our model can not be characterized
by the previous $\mathbb{Z}_{2}$ index. Instead, we should use the winding
numbers, or equivalently the chiralities, to characterize these topological
defects in momentum space. In the vicinity of the defect, the effective
Hamiltonian can be written as $H_{\text{eff}}=\sum_{i,j=x,y}v_{ij}({\bf
k})\delta k_{i}\sigma_{j}$, where ${\bf k}$ denotes the position of the
defect. The chirality thus is defined as $\nu({\bf
k})=\text{sgn}(\text{Det}[v])=\pm 1$, where $[v]$ is a matrix constructed by
matrix elements $v_{ij}$. We may also calculate the winding number
$\mathcal{W}({\bf k})$ using Eq. 21 for the BdG equation by considering a loop
$l:|{\bf k-q}|=r$ around the topological defect. We have a general relation
between these two invariants, $\mathcal{W}({\bf k})=-\nu({\bf k})$ (see Eq. 25
in Methods). The TR symmetry plays a crucial role here because the two
topological defects with opposite momenta ($\pm{\bf k}$) should have the same
chiralities and winding numbers. These defects can only be annihilated by
fusing two defects with opposite chiralities. The neutrality of topological
charge in the entire space means that some other two topological defects with
opposite chiralities should be presented at some other momenta. For this
reason, the system can only admit gapless phases with 4 or 8 Dirac cones. It
is worth to emphasize that these gapless phases may be a general feature in
all topological superconducting phaseszhangfan1 ; zhangfan2 ; shusa .

The edge states in these gapless phases will exhibit some salient anisotropic
features along different directions; see Fig. 5b. When the strip is along
$[110]$ direction, we observe two MFBs at finite momenta and one in-gap
gapless linear excitation near zero momentum. The two MFBs are a typical
feature in gapless phaseszigzag , which are used to connect the two Dirac
cones with opposite chiralities. These MFBs are protected by topology thus is
robust against perturbations. To this end, we can calculate the winding number
$\mathcal{W}(\textbf{k}_{\parallel})$ along an infinity straight line
perpendicular the strip direction using Eq. 27 (see Fig. 4b), where
$\textbf{k}_{\parallel}$ is the momentum along the strip direction. We find
that $\mathcal{W}(\textbf{k}_{\parallel})=-1$, which disappears only when a
topological defect is encountered at Det$Q({\bf k})=0$. When the strip is
along $k_{x}$ and $k_{y}$ directions, the Dirac cones with opposite
chiralities are projected to the same points, thus the MFBs are shrunk to some
single gapless points.

The second salient feature is about the topological invariant for a fully
gapped 1D line in momentum space, which is also characterized by
$\mathbb{Z}_{2}$. The definition of this topological index is identical to
that in $\mathcal{N}_{\text{2d}}$Qi09 . So in the gapless phase regime where
the $\mathbb{Z}_{2}$ index in 2D is ill-defined, we can still define the
$\mathbb{Z}_{2}$ index along some particular directions when the gapless
defects are carefully avoided. We find that $\mathcal{N}_{\text{1d}}=-1$ along
both $k_{x}$ and $[110]$ directions, thus for a strip along these directions
we can observe some in-gap linear excitations near zero momentum as
demonstrated in Fig. 5b. We also find $\mathcal{N}_{\text{1d}}=+1$ when along
$k_{y}$ direction since the gap closing and reopening in the 2D bulk takes
place at the $k_{x}$ axis ($k_{y}=0$). As a consequence the in-gap linear
excitation near zero momentum is absent. This result is consistent with the
mirror winding number analysis which are $\mathcal{W}^{x}=+1$ and
$\mathcal{W}^{y}=0$.

The same analysis can be applied to the regime with 8 Dirac cones, where the
major observations are quite similar to the other gapless phase. However, in
this regime it is possible to realize $\mathcal{W}({\bf k}_{\parallel})=2$
along some particular directions — for instance $[110]$ direction — with
doubly degenerate MFBs near zero momentum; see Fig. 5c. Moreover, two MFBs can
be found at nonzero momenta. We also found the $\mathbb{Z}_{2}$ invariant is
trivial ($\mathcal{N}_{\text{1d}}=+1$) for a strip along any direction, while
in the mirror invariant subspaces we find that the mirror winding number is
$\mathcal{W}^{x}=\mathcal{W}^{y}=0$, so we can not observe any gapless linear
excitations near zero momentum in this phase.

While some of these phases are protected by both TR and mirror symmetries, we
need to point out that the basic role of mirror symmetry played in this model
is quite subtle. In the topological gapped phase we indeed find that the
winding number in the mirror invariant subspaces are nonzero. Nevertheless if
we artificially break the mirror symmetry by introducing a relative phase to
the order parameters, we can still observe these edge states. This means that
the TR symmetry is a more important protection in our model. However, in the
BDI class phases below, we will show that the mirror symmetry will play the
primary role in protecting the edge states in some parameter regimes.

Figure 6: Topological superfluids for $T_{x0}$ model. (a) The phase diagram in
the parameter space by $\epsilon_{\text{b}}$ and $\Gamma$. The meaning of
number is the same as that in Fig. 2a. The two solid lines means that the
energy gap closing and reopening take place at ${\bf k}=0$, and the dashed
line means the gap closing and reopening at $k_{x}=0$ and $k_{y}\neq 0$.
Parameters are: $t_{\uparrow,\downarrow}=0.5E_{\text{F}}$,
$\alpha=\alpha^{\prime}=1.0E_{\text{F}}$. (b) - (d) show the band structures
at the critical boundaries when $\epsilon_{\text{b}}=0.5E_{\text{F}}$.  Figure
7: Edge states in $T_{x0}$ model along different directions. Edge states for a
strip with width $L=200/k_{\text{F}}$ along $k_{x}$ (left), $[110]$ (middle)
and $k_{y}$ (right) directions in the fully gapped topological phase (a),
gapless phase with 2 (b), 4 (c) and 6 (d) Dirac cones for the four points
marked by asterisk (*) in Fig. 6. For parameters see Fig. 6.

Topological superfluids with other TR symmetries. Now we try to extend this
idea to more general Hamiltonians by assuming that all coefficients — $s_{i}$,
$\lambda_{i}$, $\Gamma_{i}$ and $t_{\sigma}$ — can be independently tuned in
each layer. In this condition we can perform an exhaustive search strategy to
find all possible models with TR symmetry, which belong to either DIII class
with $\mathbb{T}^{2}=-1$ or BDI class with $\mathbb{T}^{2}=+1$. We first
assume that the single particle $\mathcal{H}_{0}$ has time-reversal symmetry
$\mathcal{T}=\mathcal{U}K$, where $\mathcal{U}=\rho_{i}\otimes\sigma_{j}$ and
$\mathcal{T}\mathcal{H}_{0}({\bf k})\mathcal{T}^{-1}=\mathcal{H}_{0}(-{\bf
k})$. With this operator, we can construct the TR operator for the BdG
Hamiltonian, which can be defined as either
$\mathbb{T}=\tau_{0}\otimes\mathcal{U}K$ or
$\mathbb{T}=\tau_{z}\otimes\mathcal{U}K$. For convenience this system is
dubbed as $T_{ij}$ model. All the possible results using this method are
summarized in Table 1, in which only the models with nonzero tunnelings and
nonzero SOC coefficients are presented. We find ten different models with TR
symmetries, six of which can support nontrivial topological phases. Notice
that some of these models require spin-dependent tunnelings, which may be
realized by spin selective laser-assisted hopping between the two layersLAH
and by modulating the spin-dependent optical latticeshaking1 ; shaking2 . The
two layers with different SOC coefficients may need some more complicated
laser configurations. There is an interesting relationship between the sign of
SOC coefficients and spin-dependent tunnelings, namely, a sign flip for the
SOC coefficients $s_{1}/s_{2}$ and $\lambda_{1}/\lambda_{2}$ is equivalent to
a sign flip in $t_{\uparrow}/t_{\downarrow}$, which can be mapped to each
other by a unitary transformation. Therefore the $T_{xy}$ and $T_{yx}$,
$T_{xx}$ and $T_{yy}$, and $T_{x0}$ and $T_{xz}$ are mathematically
equivalent, although their physical realizations are totally different. Here
we are interested in these models because all these phases may in principle be
realized based on $s$-wave interaction in this versatile platform, where
different phases are possible to be tuned to each other.

Table 1: Symmetry table for all the possible superfluids with TR symmetry. The single particle TR operator is defined as $\mathcal{T}=\mathcal{U}K$, where $\mathcal{U}$ is shown in the second column. The relation between the coefficients $\Gamma_{i}$, $s_{i}$, $\lambda_{i}$, $t_{i}$ and $\Delta_{i}$ are shown from column three to seven, in which the asterisk (*) means that their values can be arbitrary real numbers while $0/0$ means that both the two parameters should be zero simultaneously. $\mathbb{T}$ is the TR operator for the BdG Hamiltonian and $\mathcal{M}_{ss^{\prime}}^{\alpha}$ defines the symmetry in the mirror invariant subspace. The last column shows whether topological gapped and/or gapless phases can exist (Y) or not (N) in some appropriate parameter regimes. Model |  $\mathcal{U}$ |  $\Gamma_{1}/\Gamma_{2}$ |  $s_{1}/s_{2}$ |  $\lambda_{1}/\lambda_{2}$ |  $t_{\uparrow}/t_{\downarrow}$ |  $\Delta_{1}/\Delta_{2}$ |  $\mathbb{T}$ ($\mathbb{T}^{2}$) |  Mirror $\mathcal{M}_{ss^{\prime}}^{x}$ |  Mirror $\mathcal{M}_{ss^{\prime}}^{y}$ |  Topo.
---|---|---|---|---|---|---|---|---|---|---
$T_{xx}$ |  $\rho_{x}\otimes\sigma_{x}$ |  -1 |  -1 |  -1 |  +1 |  +1 |  $\tau_{z}\otimes\mathcal{T}$ (+1) |  $\mathcal{M}_{+-}^{x}$ (AI) |  $\mathcal{M}_{+-}^{y}$ (AI) |  N
$T_{yy}$ |  $\rho_{y}\otimes\sigma_{y}$ |  -1 |  +1 |  +1 |  -1 |  -1 |  $\tau_{z}\otimes\mathcal{T}$ (+1) |  $\mathcal{M}_{+-}^{x}$ (AI) |  $\mathcal{M}_{+-}^{y}$ (AI) |  N
$T_{xy}$ |  $\rho_{x}\otimes\sigma_{y}$ |  -1 |  +1 |  +1 |  +1 |  +1 |  $\tau_{0}\otimes\mathcal{T}$ (-1) |  $\mathcal{M}_{--}^{x}$ (AIII) |  $\mathcal{M}_{--}^{y}$ (AIII) |  Y
$T_{yx}$ |  $\rho_{y}\otimes\sigma_{x}$ |  -1 |  -1 |  -1 |  -1 |  -1 |  $\tau_{0}\otimes\mathcal{T}$ (-1) |  $\mathcal{M}_{--}^{x}$ (AIII) |  $\mathcal{M}_{--}^{y}$ (AIII) |  Y
$T_{xz}$ |  $\rho_{x}\otimes\sigma_{z}$ |  +1 |  +1 |  -1 |  * |  -1 |  $\tau_{0}\otimes\mathcal{T}$ (+1) |  $\mathcal{M}_{+-}^{x}$ (AI) |  $\mathcal{M}_{++}^{y}$ (BDI) |  Y
$T_{x0}$ |  $\rho_{x}\otimes\sigma_{0}$ |  +1 |  -1 |  +1 |  * |  +1 |  $\tau_{0}\otimes\mathcal{T}$ (+1) |  $\mathcal{M}_{+-}^{x}$ (AI) |  $\mathcal{M}_{++}^{y}$ (BDI) |  Y
$T_{zy}$ |  $\rho_{z}\otimes\sigma_{y}$ |  0/0 |  * |  * |  -1 |  * |  $\tau_{0}\otimes\mathcal{T}$ (-1) |  $\mathcal{M}_{--}^{x}$ (AIII) |  $\mathcal{M}_{--}^{y}$ (AIII) |  N
$T_{0y}$ |  $\rho_{0}\otimes\sigma_{y}$ |  0/0 |  * |  * |  +1 |  * |  $\tau_{0}\otimes\mathcal{T}$ (-1) |  $\mathcal{M}_{--}^{x}$ (AIII) |  $\mathcal{M}_{--}^{y}$ (AIII) |  N
$T_{0z}$ |  $\rho_{0}\otimes\sigma_{z}$ |  * |  * |  0/0 |  * |  * |  $\tau_{z}\otimes\mathcal{T}$ (+1) |  — |  $\mathcal{M}_{+-}^{y}$ (AI) |  Y
$T_{00}$ |  $\rho_{0}\otimes\sigma_{0}$ |  * |  0/0 |  * |  * |  * |  $\tau_{0}\otimes\mathcal{T}$ (+1) |  $\mathcal{M}_{+-}^{x}$ (AI) |  — |  Y

In all these models, our detailed self-consistent calculation demonstrates
that the intimate relation between TR and mirror symmetries is always
respected due to the phase locking effect. We determine all the topological
classes for the Hamiltonians in the mirror symmetric invariant subspaces using
the method in Eq. 7, which are summarized in Table 1. The analysis based on
the parameter $\mathcal{M}_{ss^{\prime}}^{\alpha}$ is completely consistent
with the direct analysis from the HamiltoniansMS1 ; MS2 ; Ueno ; Ken . In all
these models the subspaces can belong to AI, AIII and/or BDI classes, where
the later two classes can support nontrivial topological phases in
1DSchnyder08 ; Kitaev .

We utilize the $T_{x0}$ model as an example to highlight the unique features
of topological superfluids in BDI class. This is also a physical model that
may support topological protected MFs in 1DSumanta . Unfortunately, this
physical model may not have counterpart in solid state systems with 2D SOC, so
the realization of topological gapped and gapless superfluids in BDI class is
an important extension of our understanding of topological phases. For this
model the TR symmetry requires $s_{1}=-s_{2}=\alpha$,
$\lambda_{1,2}=\alpha^{\prime}$ and $\Gamma_{1,2}=\Gamma$ (see table I).
Hereafter we focus on the case $t_{\uparrow,\downarrow}=t$. The phase diagram
from the self-consistent calculation is presented in Fig. 6a. To understand
the origin of the gapless phases, we calculate the determinant of the
Hamiltonian using the method in Eq. 8, which can be written as

$\text{Det}(Q)=\mathcal{G}({\bf k})+8i\alpha^{\prime}k_{x}\epsilon_{{\bf
k}}t\Delta,$ (11)

where $\mathcal{G}({\bf k})=(\epsilon_{\bf
k}^{2}-t^{2}-k_{y}^{2}\alpha^{2}-k_{x}^{2}\alpha^{\prime
2}-\Gamma^{2}+\Delta^{2})^{2}-4t^{2}(k_{x}^{2}\alpha^{\prime
2}+\Gamma^{2})+4(t^{2}+k_{y}^{2}\alpha^{2}+k_{x}^{2}\alpha^{\prime
2})\Delta^{2}$, with $\Delta_{1}=\Delta_{2}=\Delta$. A direct calculation
shows that the gapless points can only happen along the $k_{y}$ axis
($k_{x}=0$). Furthermore we can verify easily that when $\Gamma=0$,
$\text{Det}(Q)$ is always greater than zero when $t\neq 0$, so these
topological superfluids are also driven by the Zeeman field. In fact we need
to emphasize that this is a quite general feature for all the topological
phases summarized in Table 1. Mathematically, the polynominal equation
$\text{Det}(Q)=0$ can always have real solutions in some parameter regimes. In
our numerical simulation, we find that the critical boundaries are also
determined by three lines as a function of binding energy and Zeeman field.
The two solid lines give topological transitions with gap closing and
reopening at ${\bf k}=0$ and $(t\pm\mu)^{2}+\Delta^{2}=\Gamma^{2}$; see Fig.
6b-c. Obviously these two critical boundaries are reduced to the well-known
result in Eq. 10 when $t=0$. The dashed line is the gap closing and reopening
at finite momentum ($k_{x}=0,k_{y}\neq 0$); see Fig. 6d. These critical
boundaries are totally different from the $T_{xy}$ model shown in Fig. 2,
which can have profound influence on the topological phases with slightly TR
symmetry breaking; see below.

We find that in this model the order parameters and chemical potential are
always a smooth function of Zeeman field and binding energy, so the PS phase
discussed in Fig. 3 is not shown up. However, this model can exhibit some
interesting features that not seen in the previous $T_{xy}$ model. We notice
that when $\mathbb{T}^{2}=+1$ the two Dirac cones with opposite momenta have
opposite chiralities (see an exact proof in Methods), so in principle, all the
topological protected gapless superfluids with even number ($\leq 8$) of Dirac
cones are admitted. Our crucial observations thus are in order. (1) We first
focus on the fully gapped phases which in principle can not support
topological phases in 2D from the standard tableSchnyder08 ; Kitaev . However,
mirror symmetry protected phases are still allowed, which can be characterized
by the winding number defined in Eq. 27 along the mirror symmetric invariant
subspaces. For the fully gapped phase at small Zeeman field ($\Gamma\sim 0$),
we find that the mirror winding numbers along $k_{y}$ axis are
$\mathcal{W}_{1}^{y}\oplus\mathcal{W}_{2}^{y}=0\oplus 0$. Notice that the
mirror symmetric subspaces along $k_{x}$ direction belong to AI class thus is
always trivial in 1DSchnyder08 ; Kitaev . For this reason we do not observe
gapless excitation for a strip along any directions (not shown in Fig. 7).
However, in the fully gapped phase at strong Zeeman field, we find that the
mirror winding numbers are
$\mathcal{W}_{1}^{y}\oplus\mathcal{W}_{2}^{y}=-1\oplus-1$ and
$\mathcal{W}({\bf k}_{\parallel}=k_{y})=+1$, so we can observe robust gapless
excitations for a strip along $k_{y}$ direction in Fig. 7a. (2) In the left
phase with 2 Dirac cones, the mirror winding numbers are
$\mathcal{W}_{1}^{y}\oplus\mathcal{W}_{2}^{y}=-1\oplus 0$. thus we observe
robust MFBs along both $k_{y}$ and $[110]$ directions in Fig. 7b. The MFB is
shrunk to a single point for a strip along $k_{x}$ direction since all the
Dirac cones are located at the $k_{y}$ axis. (3) In the regime with 4 Dirac
cones, the mirror winding numbers are
$\mathcal{W}_{1}^{y}\oplus\mathcal{W}_{2}^{y}=-1\oplus-1$, which are exactly
the same as the topological gapped phase at strong Zeeman field due to gap
closing and reopening at $k_{y}\neq 0$. This explains why the similar in-gap
excitations in Fig. 7a can also be observed in Fig. 7c. Meanwhile we can
observe two MFBs at finite momenta due to the same reason discussed in
$T_{xy}$ model. (4) In the regime with 6 Dirac cones, the mirror winding
numbers are $\mathcal{W}_{1}^{y}\oplus\mathcal{W}_{2}^{y}=-1\oplus 0$, which
are the same as the left gapless phase with 2 Dirac cones in Fig. 7b. In this
phase regime we can observe three MFBs for a strip along $k_{y}$ and $[110]$
directions in Fig. 7d, all of which will shrink to a single point for a strip
along $k_{x}$ direction. (5) These mirror symmetry protected phases will be
destroyed by mirror symmetry breaking terms, such as the relative phase
between the order parameters.

Discussion and conclusion. Here we investigate the topological phases in
ultracold atom systems due to their flexibility of tuning all the parameters
in a wide range in experiments, which are challenging in solid state system.
In Methods we have shown a concrete scheme to realize the bilayer structure
with exactly opposite Zeeman fields. However, it is still very interesting to
ask the basic question that what will happen under weak TR symmetry breaking?
In this condition both the DIII and BDI class phases will be collapsed to D
class phasesMGong11 ; MGong12 , in which the topological phase transitions can
be characterized by either Chern number or Pfaffian aforementioned and all
gapless phases will be destroyed. In the $T_{xy}$ model we find that all the
phases under weak TR symmetry breaking have Chern number $\mathcal{C}=0$ and
Pfaffian $\mathcal{\chi}=\text{sgn Pf}(\mathcal{H}(0)\tau_{x})=+1$. This is
due to the fact that all the energy gap closings and reopenings in this model
take place at nonzero momenta, so the Pfaffian at zero momentum is unchanged.
However, the BDI model can exhibit totally different behaviors due to the
possibility of gap closing and reopening at zero momentum. So we find that the
trivial gapped phase will have $\mathcal{C}=0$ and $\mathcal{\chi}=+1$, and
the regime with 2 and 6 Dirac cones will have $\mathcal{C}=1$ and
$\mathcal{\chi}=-1$, and the fully gapped topological phase and the regime
with 4 Dirac cones will have $\mathcal{C}=2$ and $\mathcal{\chi}=+1$. In
general, $\mathcal{\chi}=(-1)^{\mathcal{C}}$HBB . If all the parameters in
experiments can be changed adiabatically, it is possible to observe the
topological transitions among DIII, D and BDI classes without gap
closingNagaosa .

We finally emphasize that these phases will not be spoiled by the presence of
interlayer pairings. We have included all the interlayer pairings in our model
(see Methods) and minimize the total free energy $F$ with respect to these
parameters. The detailed self-consistent simulation shows that the TR symmetry
is always respected. This is due to the phase locking effect between the two
layers with real tunnelings. These pairings only slightly modify the real and
imaginary part of Det$Q$ discussed in Eq. 9 and Eq. 11, thus Det$Q=0$ can
always have real solutions in some parameter regimes. Instead, these pairings
just slightly modify the topological boundaries without qualitatively
influence all the conclusions.

To conclude, the bilayer structure defined in this work provides an
experimentally controllable platform to realize new topological protected
phases in ultracold atom systems with TR symmetry that belonging to DIII and
BDI classes. These phases can be realized using the conventional $s$-wave
interactions, thus is in stark contrast to previous physical proposals based
on unconventional pairings. In these models, the fully gapped phase in DIII
class is characterized by $\mathbb{Z}_{2}$, while all the gapless phases in
these two classes are characterized by winding numbers around the topological
defects. The phase locking effect between the two layers ensures the
realization of mirror symmetry protected gapped and gapless phases. These new
phases are expected to greatly enrich our understanding of topological
superconducting phases and related transitions. This novel idea can also be
applied to explore the MFs in 1D and the Weyl superfluids in 3D, and their
crossover from the Bose-Einstein condensation in the strong coupling limit to
the Bardeen-Cooper-Schrieffer in the weak coupling limitMGong11 ; MGong12 ;
HH11 ; ZQY ; He12 , which will be published elsewhere.

Methods

Generation of local Zeeman fields in $T_{xy}$ model. The spin-dependent
double-well potential in Eq. 3 is essential for the first four models in Table
1, which can be realized using the state-of-the-art techniques. Here we follow
the basic idea in Ref. [Ivan, ]. For the cold atom optically pumped from the
$|J=1/2\rangle$ manifold to the $|J^{\prime}=3/2\rangle$ manifold and in the
large detuning limit the rank 2 polarizability tensor takes the following
formIvan ,

$\alpha_{ij}(J\rightarrow J^{\prime})=\tilde{\alpha}({2\over
3}\delta_{ij}-{i\over 3}\epsilon_{ijk}\sigma_{k}),$ (12)

where $\tilde{\alpha}$ is a constant that inversely proportional to the
detuning. For transparency we first consider two counterpropagating laser
beams with momenta $\pm k$ and frequency $\omega$ along $z$ direction. Then
the electric field can be written as

${\bf E}=(E_{1}{\bf e}_{1}e^{i(kz+\phi)}+E_{2}{\bf e}_{2}e^{-ikz})e^{i\omega
t},$ (13)

where ${\bf e}_{1,2}$ represent the polarization direction for the two laser
beams with field strength $E_{1,2}$, and $\phi$ is the relative phase between
the two laser beams. Let ${\bf e}_{1}=\cos({\theta\over 2}){\bf
e}_{x}+\sin({\theta\over 2}){\bf e}_{y}$ and ${\bf e}_{2}=\cos({\theta\over
2}){\bf e}_{x}-\sin({\theta\over 2}){\bf e}_{y}$, where $\theta$ is the
relative polarization angle between the two laser beams, then the spin-
dependent optical lattice can be written asIvan

$U(z)={\bf E}^{*}\cdot\alpha\cdot{\bf E}=\tilde{\alpha}[{2\over 3}|{\bf
E}|^{2}+{i\over 3}({\bf E}^{*}\times{\bf E})\cdot\bm{\sigma}],$ (14)

where in the second term only the $\sigma_{z}$ component is nonzero since
${\bf E}$ and ${\bf E}^{*}$ lie in the $x$-$y$ plane. Plugging the expression
of ${\bf E}$ into Eq. 14 we find

$\displaystyle
U(z)=U[2\cos(\theta)\cos(2kz+\phi)+\sin(\theta)\sin(2kz+\phi)\sigma_{z}],$
(15)

where $U=-{2\over 3}\tilde{\alpha}E_{1}E_{2}$. Thus these two laser beams
create a spin-dependent periodic optical lattice along $z$ direction.

To realize an isolated double-well structure in Fig. 1a, we can apply two
extra laser beams along the same direction with momenta $\pm 2k$DW1 ; DW2 ,
polarizability $\tilde{\alpha}^{\prime}$, electric field strength
$E_{1,2}^{\prime}$ (hence $U^{\prime}=-{2\over
3}\tilde{\alpha}^{\prime}E_{1}^{\prime}E_{2}^{\prime}$), polarization angle
$\theta^{\prime}$ and relative phase $\phi^{\prime}$. These four colinear
laser beams can create a spin-dependent potential in the following form,

$\displaystyle U(z)$
$\displaystyle=2U\cos(\theta)\cos(2kz)+2U^{\prime}\cos(\theta^{\prime})\cos(4kz+\delta\phi)+$
(16)
$\displaystyle(U\sin(\theta)\sin(2kz)+U^{\prime}\sin(\theta^{\prime})\sin(4kz+\delta\phi))\sigma_{z},$

where a shift of $z$ by $-\phi/2k$ has been made and
$\delta\phi=\phi^{\prime}-2\phi$. We plot the band structure in Fig. 1b using
the above equation with $\delta\phi=0$, which is exactly the ideal potential
presented in Eq. 3. In this potential the first part —an even function of $z$
—gives the symmetric double-well potential while the second part — an odd
function of $z$ — gives the spin-dependent Zeeman fields. The two layers
therefore have exactly opposite Zeeman fields and this feature is independent
of the choice of other parameters. We also need to emphasize that a bilayer
system with TR symmetry can still be achieved even deviates from this idea
condition because five independent parameters — $U$, $U^{\prime}$, $\theta$,
$\theta^{\prime}$, $\delta\phi$ — can be tuned to ensure just two constraints
— $t_{\uparrow}=t_{\downarrow}$ and $\Gamma_{1}=-\Gamma_{2}$ — in $T_{xy}$
model, which can always be fulfilled in a vast range of parameters.
Furthermore the external Zeeman field along $z$ direction, $B_{z}\sigma_{z}$,
can also contribute to the recovery of TR symmetry.

Self-consistent calculation. The pairings in ultracold atoms is induced by the
$s$-wave scattering between the particles thus should be determined in a self-
consistent manner. In the major numerical simulation, we only consider the
intralayer pairings by assuming that the interlayer pairings are sufficiently
suppressed by the strong barrier between the two layers (see Fig. 1). The
thermodynamic potential $\Omega$ is defined as
$Z=e^{-\Omega/T}=\text{Tr}e^{-H/T}$. When $T\rightarrow 0$ we have

$\Omega={1\over 2}\sum_{\eta<0,{\bf k}}(E_{\eta{\bf k}}+2\epsilon_{1{\bf
k}}+2\epsilon_{2{\bf k}})+{|\Delta_{1}|^{2}\over g}+{|\Delta_{2}|^{2}\over
g},$ (17)

where $\eta<0$ means summation over all the occupied bands with $E_{\eta{\bf
k}}\leq 0$. The condensation energy should be regularized using the
prescription: $g=\sum_{{\bf k}}1/({\bf k}^{2}/m+\epsilon_{\text{b}})$, where
$\epsilon_{\text{b}}$ defines the binding energy which can be controlled by
the Feshbach resonance. The corresponding free energy used in Fig. 3b is
defined as $F=\Omega+n\mu$. The order parameters and chemical potential are
determined by the extrema of $F$: $\partial F/\partial\Delta_{i}=0$, $\partial
F/\partial\Delta_{i}^{*}=0$ and $\partial F/\partial\mu=0$.

In our simulation, the two order parameters are assumed to
$\Delta_{1}e^{i\delta\theta}$, and $\Delta_{2}e^{-i\delta\theta}$
respectively, where $\Delta_{1,2}$, $\delta\theta$ $\in\mathbb{R}$. The global
phase has been gauged out, so $\delta\theta$ determines the relative phase
between the order parameters in the two layers. There three parameters
together with the chemical potential are used to minimize the total free
energy $F$. In our numerical simulation we have assumed
$k_{\text{F}}=\sqrt{n\pi}$, where $n$ is the fixed total particle density, and
$E_{\text{F}}=k_{\text{F}}^{2}/2m$. In all our simulation the energy and
momentum are rescaled by $E_{\text{F}}$ and $k_{\text{F}}$, respectively. We
find that in all our numerical simulation $\delta\theta=0$, which means that
the order parameters between the two layers always have fixed phase. This
basic feature, dubbed as phase locking effect throughout this work, is still
respected by the presence of interlayer pairings (see below).

Phase separation. The instability of the uniform phases towards the PS phase
can be understood from the mixing of two different phases, which is captured
by the following mixed free energycaldas ; WYI11 ,

$F_{\text{m}}=xF(\mu,\Delta_{1},\Delta_{2})+(1-x)F(\mu,\Delta_{1}^{\prime},\Delta_{2}^{\prime}),$
(18)

where $0\leq x\leq 1$ defines the mixing coefficient (see numerical result in
Fig. 3d). The two phases should have the same chemical potential but different
pairing strengthes $\Delta_{1,2}$ and $\Delta_{1,2}^{\prime}$. We minimize the
total mixed free energy $F_{\text{m}}$ with respect to all these parameters to
map out the PS phase in Fig. 2a. We find that this phase is not favorable in
the $T_{x0}$ model.

Topological invariants in the gapless phases. Now we investigate the effect of
TR, PH and mirror symmetries on the winding numbers and chiralities for the
gapless phases.

(1) We first notice that the off-diagonal matrix $Q$ in Eq. 8 and Eq. 11 have
the following propertiesSchnyder08 ; Kitaev

$\displaystyle\text{Det}[Q_{\text{DIII}}^{T}({\bf k})]$
$\displaystyle=\text{Det}[-Q_{\text{DIII}}(-{\bf k})],\quad\text{ and }$ (19)
$\displaystyle\text{Det}[Q_{\text{BDI}}^{*}({\bf k})]$
$\displaystyle=\text{Det}[Q_{\text{BDI}}(-{\bf k})].$ (20)

In these models the winding number around a single topological defect is
defined as

$\displaystyle\mathcal{W}({\bf k})$ $\displaystyle={i\over 4\pi}\oint_{l}d{\bf
q}\text{Tr}\mathbb{S}\mathcal{H}^{-1}\partial_{\bf q}\mathcal{H}$ (21)
$\displaystyle={1\over 2\pi}\Im\oint_{l}d{\bf q}\partial_{\bf
q}\text{Det}Q({\bf q})\in\mathbb{Z},$

where ${\bf k}$ defines the position of the defect and the loop $l:=|{\bf
q-k}|=r$ should be small enough to enclose only one defect. We have

$\displaystyle\mathcal{W}_{\text{DIII}}({\bf
k})=\mathcal{W}_{\text{DIII}}(-{\bf k}),\mathcal{W}_{\text{BDI}}({\bf
k})=-\mathcal{W}_{\text{BDI}}(-{\bf k}).$ (22)

We see that the TR symmetry plays a fundamental role in deriving the above
results. In fact when TR is broken the $Q$ matrix is no longer well-defined
and the topological defects are then destroyed. We do not observe the robust
Dirac cones in topological D class phases in all our models in Table 1 by
introducing some TR symmetry breaking terms to Eq. 1.

(2) Next we study the chirality of the topological defect. We assume the
defect is located at ${\bf k}$ with degenerate zero energy eigenvectors
$\psi_{\pm}$, where $\psi_{+}=\mathbb{S}\psi_{-}$ due to chiral symmetry.
Based on these functions we can define two chiral basis
$\phi_{\pm}=\psi_{+}\pm\psi_{-}$, which satisfy
$\mathbb{S}\phi_{\pm}=\pm\phi_{\pm}$. We assume that $H({\bf k}+\delta{\bf
k})=H({\bf k})+M(\delta{\bf k})$ in the vicinity of the defect, where
$M(\delta{\bf k})=\delta k_{x}M_{x}+\delta k_{y}M_{y}$. The chiral symmetry
ensures that $\mathbb{S}M\mathbb{S}^{\dagger}=-M$. With these results we
readily have $\langle\phi_{\pm}|M|\phi_{\pm}\rangle=0$, and the off-diagonal
term,
$\langle\phi_{+}|M|\phi_{-}\rangle=\langle\psi_{+}|M|\psi_{+}\rangle-\langle\psi_{-}|M|\psi_{-}\rangle-\langle\psi_{+}|M|\psi_{-}\rangle+\langle\psi_{-}|M|\psi_{+}\rangle$.
Notice that
$\langle\psi_{+}|M|\psi_{+}\rangle=-\langle\psi_{-}|M|\psi_{-}\rangle$ and
$\langle\psi_{+}|M|\psi_{-}\rangle=-\langle\psi_{-}|M|\psi_{+}\rangle$, so the
last two terms should be exactly imaginary numbers. Therefore we can assume
$\langle\phi_{+}|M|\phi_{-}\rangle=a_{1}\delta k_{x}+a_{2}\delta
k_{y}+i(b_{1}\delta k_{x}+b_{2}\delta k_{y})$, where $a_{1,2}$, $b_{1,2}$
$\in\mathbb{R}$. The effective Hamiltonian can be written as

$H_{\text{eff}}=\begin{pmatrix}0&c_{1}\delta k_{x}+c_{2}\delta k_{y}\\\
*&0\end{pmatrix}=\sum_{ij}v_{ij}\delta k_{i}\sigma_{j},$ (23)

where $c_{i}=a_{i}+ib_{i}$. The effective Hamiltonian constructed in this way
can fully embody the importance of TR symmetry in this model. The chirality
for Eq. 23 is defined as

$\nu=\text{sgn}(\text{Det}[v])=\text{sgn}(b_{1}a_{2}-b_{2}a_{1}).$ (24)

We also find a general relation between these invariants,

$\mathcal{W}({\bf k})=-\nu({\bf k}).$ (25)

We turn to study the topological charge at $-{\bf k}$, which also have two
degenerate eigenvectors $\psi_{\pm}^{\prime}$ with
$\psi_{+}^{\prime}=\mathbb{S}\psi_{-}^{\prime}$. The PH and TR symmetries
ensure that $\psi_{-}^{\prime}=\mathbb{T}\psi_{-}=\eta\mathbb{C}\psi_{+}$, and
$\psi_{+}^{\prime}=\mathbb{C}\psi_{-}=\eta\mathbb{T}\psi_{+}$, where
$\eta=\mathbb{T}^{2}$. The perturbation term is assumed to be $M^{\prime}$.
Let $\phi_{\pm}^{\prime}=\eta\mathbb{T}\psi_{+}\pm\mathbb{T}\psi_{-}$, then

$\displaystyle\langle\phi_{+}^{\prime}|M^{\prime}|\phi_{-}^{\prime}\rangle=\langle\eta\psi_{+}-\psi_{-}|M|\eta\psi_{+}+\psi_{-}\rangle.$
(26)

For DIII model, $\eta=-1$, and
$\langle\phi_{+}^{\prime}|M^{\prime}|\phi_{-}^{\prime}\rangle=\langle\phi_{+}|M|\phi_{-}\rangle$,
so the two defects with opposite momenta have the same chirality. They will
have opposite chiralities when $\eta=1$ in BDI model since
$\langle\phi_{+}^{\prime}|M^{\prime}|\phi_{-}^{\prime}\rangle=(\langle\phi_{+}|M|\phi_{-}\rangle)^{\dagger}$.
This analysis is completely consistent with the results in Eqs. 22 and 25.

(3) To understand the origin of the MFBs localized at the boundary of the
strip, we consider another integer topological invariant defined as (see Fig.
4b)satofb ,

$\mathcal{W}(\textbf{k}_{\parallel})={1\over 2\pi}\Im\int
d\textbf{k}_{\perp}\partial\textbf{k}_{\perp}\ln\text{Det}Q(\textbf{k})\in\mathbb{Z},$
(27)

where $\textbf{k}_{\perp}$ is the direction perpendicular to the strip. If
this number is nonzero, it means topological protected edge states with zero
energy can be found at the two boundaries. The number of edge states is
defined by $|\mathcal{W}(\textbf{k}_{\parallel})|$.

Mirror symmetries and mirror winding numbers For the $T_{xy}$ model,
$M_{1,2}^{x}=-\epsilon_{k}^{x}\sigma_{z}\otimes\rho_{0}-\Gamma\sigma_{z}\otimes\rho_{z}\
-(k_{x}\alpha^{\prime}\pm
t)\sigma_{z}\otimes\rho_{y}-\Delta\sigma_{y}\otimes\rho_{y}$, and
$M_{1,2}^{y}=-\epsilon_{k}^{y}\sigma_{z}\otimes\rho_{0}-\Gamma\sigma_{z}\otimes\rho_{z}+(k_{y}\alpha\pm
t)\sigma_{0}\otimes\rho_{x}-\Delta\sigma_{y}\otimes\rho_{y}$. Here
$\epsilon_{k}^{x,y}=k_{x,y}^{2}/2m-\mu$, and $\rho$ and $\sigma$ are Pauli
matrices which do not have the meaning defined in the main text. These two
models belong to AIII model with only chiral symmetry and the topological
invariant is characterized by $\mathbb{Z}$ in 1D.

For the $T_{x0}$ model only the mirror spaces generated by $\mathcal{M}_{y}$
is important. For the case studied in the main text,
$t_{\uparrow}=t_{\downarrow}=t$, $\alpha=\alpha^{\prime}$ and
$M_{1,2}^{y}=-(\epsilon_{k}^{y}\pm
t)\sigma_{z}\otimes\rho_{0}+\Gamma\sigma_{z}\otimes\rho_{z}-k_{x}\alpha^{\prime}\sigma_{z}\otimes\rho_{y}-\Delta\sigma_{y}\otimes\rho_{y}$.
These two Hamiltonians belong to topological BDI class with TR, PH and chiral
symmetries, and the invariant is characterized by $\mathbb{Z}$ in 1D.

Since all the above mirror subspaces have chiral symmetry (AIII for $T_{xy}$
and BDI for $T_{x0}$), we assume $s_{\alpha}$ is the chiral symmetry operator
for $M_{\alpha}$ (see Eq. 5 and 6), where

$u_{\alpha}s_{\alpha}u_{\alpha}^{-1}=\begin{pmatrix}\bf{I}_{2}&0\\\
0&-\bf{I}_{2}\end{pmatrix},u_{\alpha}M^{\alpha}_{i}u_{\alpha}^{-1}=\begin{pmatrix}0&q^{\alpha}_{i}\\\
q_{i}^{\alpha\dagger}&0\end{pmatrix}.$ (28)

Then the mirror winding number in these mirror symmetric invariant subspaces
can be calculated using Eq. 27, in which $\mathcal{H}$ and $\mathbb{S}$ are
replaced by $M_{i}^{\alpha}$ and $s_{\alpha}$. For $T_{xy}$ model,
$s_{x}=\sigma_{y}\otimes\sigma_{z}$ and $s_{y}=\sigma_{x}\otimes\sigma_{0}$;
and for $T_{x0}$ model $s_{y}=\sigma_{x}\otimes\sigma_{0}$.

There is an important difference between these two models ($T_{xy}$ and
$T_{x0}$) thought the definition of mirror winding numbers is exactly the
same. In $T_{xy}$ model, the two mirror subspaces $M_{1}^{\alpha}$ and
$M_{2}^{\alpha}$, which belongs to AIII class in both $k_{x}$ and $k_{y}$
directions, are not independent, but are related with each other by TR and PH
symmetries, thus there is just only one independent winding number
$\mathcal{W}^{\alpha}$ along each direction. In the main text this number is
calculated using $M_{1}^{\alpha}$ since
$\mathcal{W}^{\alpha}=\mathcal{W}_{1}^{\alpha}\equiv\mathcal{W}_{2}^{\alpha}$.
However in $T_{x0}$ model the two mirror symmetric invariant subspaces belong
to BDI class when along $k_{y}$ direction, which have independent TR and PH
symmetries, thus the mirror winding numbers are independent. For this reason
the mirror topological properties in BDI model are labeled by two integer
numbers $\mathcal{W}^{y}_{1}\oplus\mathcal{W}^{y}_{2}$, where
$\mathcal{W}^{y}_{1,2}$ are calculated for $M_{1,2}^{y}$ respectively. We do
not have to consider the invariant along $k_{x}$ direction, which is trivial.

Below we summarize the results for mirror winding numbers for all phases in
the phase diagrams in the main text. For $T_{xy}$ model, the topological
gapped phase has $\mathcal{W}^{x}_{1}=1$ and $\mathcal{W}^{y}_{1}=-1$; the
phase with 4 Dirac cones has $\mathcal{W}^{x}_{1}=1$ and
$\mathcal{W}^{y}_{1}=0$; and the regime with 8 Dirac cones has
$\mathcal{W}^{x}_{1}=0$ and $\mathcal{W}^{y}_{1}=0$. For $T_{x0}$ model, the
right phase with 2 Dirac cones has
$\mathcal{W}^{y}_{1}\oplus\mathcal{W}^{y}_{2}=0\oplus-1$; the left phase with
2 Dirac cones has $\mathcal{W}^{y}_{1}\oplus\mathcal{W}^{y}_{2}=-1\oplus 0$;
the regime with 4 Dirac cones has
$\mathcal{W}^{y}_{1}\oplus\mathcal{W}^{y}_{2}=-1\oplus-1$; the regime with 6
Dirac cones has $\mathcal{W}^{y}_{1}\oplus\mathcal{W}^{y}_{2}=-1\oplus 0$ and
the fully gapped topological phase has
$\mathcal{W}^{y}_{1}\oplus\mathcal{W}^{y}_{2}=-1\oplus-1$.

Effect of interlayer pairings. For convenience we use the following notations
to define all the possible scatterings between the two layers with
$g_{\alpha\beta\gamma\delta}$ being the scattering strength and the Greek
alphabet being the layer index,

$\displaystyle F_{\alpha\beta\gamma\delta}:$
$\displaystyle=g_{\alpha\beta\gamma\delta}\sum_{{\bf
k,q}}(c_{\alpha\textbf{k}\uparrow}^{\dagger}c_{\beta-\textbf{k}\downarrow}^{\dagger})(c_{\gamma\textbf{q}\downarrow}c_{\delta-\textbf{q}\uparrow})$
$\displaystyle\simeq{g_{\alpha\beta\gamma\delta}\over
g}\cdot[\sum_{\textbf{k}}\Delta_{\alpha\beta}^{*}c_{\delta\textbf{k}\uparrow}c_{\gamma-\textbf{k}\downarrow}+\text{h.c}+{\Delta_{\alpha\beta}^{*}\Delta_{\delta\gamma}\over
g}].$

We define the intralayer and interlayer pairings using the following way:
$\Delta_{\alpha\beta}=g\langle\sum_{\textbf{k}}c_{\alpha\textbf{k}\uparrow}c_{\beta-\textbf{k}\downarrow}\rangle$,
and
$\Delta_{\alpha\beta}^{*}=g\langle\sum_{\textbf{k}}c_{\beta-\textbf{k}\downarrow}c_{\alpha\textbf{k}\uparrow}\rangle$.
The leading pairings in the form $\Delta_{\alpha\alpha}=\Delta_{\alpha}$ have
been used in the main text. We have to sum up all possible Feynman diagrams to
construct the mean-field BdG equation. To this end, we define
$g_{1111}=g_{2222}=g$, $g_{1112}=g_{1121}=\cdots=g^{\prime}$, and
$g_{1122}=g_{1221}=\cdots=g^{\prime\prime}$, with
$g^{\prime\prime}<g^{\prime}\ll g$. Then we can write

$U=U^{\dagger}+U^{-}+U_{0},$ (30)

where
$U^{\dagger}=\sum_{\textbf{k}}\Delta_{1}c_{1\textbf{k}\uparrow}^{\dagger}c_{1-\textbf{k}\downarrow}^{\dagger}+\Delta_{2}c_{2\textbf{k}\uparrow}^{\dagger}c_{2-\textbf{k}\downarrow}^{\dagger}+\Delta_{3}c_{1\textbf{k}\uparrow}^{\dagger}c_{2-\textbf{k}\downarrow}^{\dagger}+\Delta_{4}c_{2\textbf{k}\uparrow}^{\dagger}c_{1-\textbf{k}\downarrow}^{\dagger}$,
$U^{-}=(U^{\dagger})^{\dagger}$, and
$U_{0}=-\sum_{\alpha,\beta,\gamma,\delta=1,2}{g_{\alpha\beta\gamma\delta}\over
g^{2}}\Delta_{\alpha\beta}^{*}\Delta_{\delta\gamma}$. Here we have defined
four order parameters due to all possible scatterings. In this condition, for
the BdG equation in Eq. 1, the off-diagonal block matrix $\bar{\Delta}$ have
to be replaced by

$\bar{\Delta}=\begin{pmatrix}0&\Delta_{1}&0&\Delta_{3}\\\
-\Delta_{1}&0&\Delta_{4}&0\\\ 0&-\Delta_{4}&0&\Delta_{2}\\\
-\Delta_{3}&0&-\Delta_{2}&0\end{pmatrix},$ (31)

where $\Delta_{1}=-\Delta_{11}-{g^{\prime}\over
g}(\Delta_{12}+\Delta_{21})-{g^{\prime\prime}\over g}\Delta_{22}$,
$\Delta_{2}=-\Delta_{22}-{g^{\prime}\over
g}(\Delta_{12}+\Delta_{21})-{g^{\prime\prime}\over g}\Delta_{11}$,
$\Delta_{3}=-{g^{\prime}\over
g}(\Delta_{11}+\Delta_{22})-{g^{\prime\prime}\over
g}(\Delta_{12}+\Delta_{21})$, $\Delta_{4}=-{g^{\prime}\over
g}(\Delta_{11}+\Delta_{22})-{g^{\prime\prime}\over
g}(\Delta_{21}+\Delta_{12})$. The TR symmetry for this new order parameter is
still determined by the operator $\mathbb{T}$ used in the main text, but now,
$\Delta_{1,2}$ and $\Delta_{3,4}$ should respect some restrict relations. If
we define $\mathcal{O}=\Delta_{12}+\Delta_{21}$, then we can simplify the
condensation energy

$\displaystyle-U_{0}=$ $\displaystyle{1\over
g}(|\Delta_{11}|^{2}+|\Delta_{22}|^{2})+{g^{\prime}\over
g^{2}}\mathcal{O}(\Delta_{11}+\Delta_{22}+\Delta_{11}^{*}+\Delta_{22}^{*})$
(32) $\displaystyle+{g^{\prime\prime}\over
g^{2}}(\mathcal{O}^{2}+\Delta_{11}\Delta_{22}^{*}+\Delta_{11}^{*}\Delta_{22}).$

This energy should be used to replace the last two terms in Eq. 17. The self-
consistent numerical sumulation demonstrate that all the conclusions obtained
in the main text are still valid in the presence of these interlayer pairings.

Acknowledgements We thank Y. X. Zhao, F. Zhang and Z. D. Wang for valuable
discussions. C.C., J.L. and M.G. are supported by Hong Kong RGC/GRF Projects
(No. 401011, 401113) and The Chinese University of Hong Kong Focused
Investments Scheme. B.H. is supported by Natural Science Foundation of Jiangsu
Province under Grant No. BK20130424. C.Z. is supported by ARO
(W911NF-12-1-0334) AFOSR (FA9550-13-1-0045) and NSF (PHY-1505496).

Author Contributions B.H. and C.C contributes equally to this work. M.G.
conceived the idea and supervised the project. B.H. and C.C. performed the
self-consistent calculation and symmetry analysis. M.G. and C.Z. wrote this
manuscript. All authors participated in discussions about this work.
Correspondence and requests for materials should be addressed to M.G.
<EMAIL_ADDRESS>and C.Z. (chuanweizhang@utdallas.com).

Competing Interests The authors declare no competing financial interests.

## References

  * (1) Nayak, C., Simon, S. H., Stern, A., Freedman, M. & Das Sarma, S. Non-Abelian anyons and topological quantum computation. Rev. Mod. Phys. 80, 1083 (2008).
  * (2) Matzdorf, R., Fang, Z., Ismail, Zhang, J., Kimura, T., Tokura, Y. Terakura, K., & Plummer, E. W. Ferromagnetism stabilized by lattice distortion at the surface of the $p$-wave superconductor Sr2RuO4. Science, 289, 746-748 (2000).
  * (3) Hicks, C. W., Brodsky, D. O., Yelland, E. A., Gibbs, A. S., Bruin, J. A. N., Barber, M. E., Edkins, S. D., Nishimura, K., Yonezawa, S., Maeno, Y., & Mackenzie, A. P. Strong increase of $T_{c}$ of Sr2RuO4 under both tensile and compressive strain. Science, 344, 283-285 (2014).
  * (4) Volovik, G. E. The Universe in a Helium Droplet (Oxford Univ. Press, 2003).
  * (5) Gurarie, V., Radzihovsky, L., & Andreev, A.V. Quantum phase transitions across a $p$-wave Feshbach resonance. Phys. Rev. Lett. 94, 230403 (2005).
  * (6) Levinsen, J., Cooper, N. R., & Shlyapnikov, G. V. Topological $p_{x}+ip_{y}$ superfluid phase of fermionic polar molecules, Phys. Rev. A 84, 013603 (2011).
  * (7) Read, N., & Green, D. Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum hall effect. Phys. Rev. B 61, 10267 (2000).
  * (8) Fu, L., & Kane, C. L. Superconducting proximity effect and Majorana fermions at the surface of a topological insulator. Phys. Rev. Lett. 100, 096407 (2008).
  * (9) Lutchyn, R. M., Sau, J. D., & Das Sarma, S. Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures. Phys. Rev. Lett. 105, 077001 (2010).
  * (10) Lutchyn, R. M., Stanescu, T. D., & Das Sarma, S. Search for Majorana fermions in multiband semiconducting nanowires. Phys. Rev. Lett. 106, 127001 (2011).
  * (11) Alicea, J., Oreg, Y., Refael, G., Oppen, F. von, & Fisher, M. P. A. Non-Abelian statistics and topological quantum information processing in 1D wire networks. Nat. Phys. 7, 412 (2011).
  * (12) Mourik, V., Zuo, K., Frolov, S. M., Plissard, S. R., Bakkers, E. P. A. M., & Kouwenhoven, L. P. Signatures of Majorana fermions in hybrid superconductor-semiconductor nanowire devices. Science 336, 1003 (2012).
  * (13) Das, A., Ronen, Y., Most, Y., Oreg, Y., Heiblum, M., & Shtrikman, H. Zero-bias peaks and splitting in an Al–InAs nanowire topological superconductor as a signature of Majorana fermions. Nat. Phys. 8, 887 (2012).
  * (14) Finck, A. D. K., Harlingen, D. J. V., Mohseni, P. K., Jung, K., & Li, X. Anomalous modulation of a zero-bias peak in a hybrid nanowire-superconductor device. Phys. Rev. Lett. 110, 126406 (2013).
  * (15) Deng, M. T., Yu, C. L., Huang, G. Y., Larsson, M., Caroff, P., & Xu, H. Q. Anomalous Zero-Bias Conductance Peak in a Nb–InSb Nanowire–Nb Hybrid device. Nano. Lett. 12, 6414 (2012).
  * (16) Nadj-Perge, S., Drozdov, I. K., Li, J., Chen, H., Jeon, S., Seo, J., MacDonald, A. H., Bernevig, B. A., & Yazdani, A. Observation of Majorana fermions in ferromagnetic atomic chains on a superconductor. Science 346, 602 (2014).
  * (17) Schnyder, A. P., Ryu, S., Furusaki, A., & Ludwig, A. W. W. Classification of topological insulators and superconductors in three spatial dimensions. Phys. Rev. B 78, 195125 (2008).
  * (18) Kitaev, A. Periodic table for topological insulators and superconductors. AIP Conf. Proc. 1134, 22-30 (2009).
  * (19) Kitaev, A. Y. Unpaired Majorana fermions in quantum wires. Phys. Usp. 44, 131 (2001).
  * (20) Qi, X., Hughes, T. L., Raghu, S., & Zhang, S. Time-reversal-invariant topological superconductors and superfluids in two and three dimensions. Phys. Rev. Lett. 102, 187001 (2009).
  * (21) Fu, L., & Berg, E. Odd-parity topological superconductors: theory and application to CuxBi2Se3. Phys. Rev. Lett. 105, 097001 (2010).
  * (22) Sato, M. Topological odd-parity superconductors. Phys. Rev. B 81, 220504(R) (2010).
  * (23) Nakosai, S., Budich, J. C., Tanaka, Y., Trauzettel, B., & Nagaosa, N. Majorana bound states and non-local spin correlations in a quantum wire on an unconventional superconductor. Phys. Rev. Lett. 110, 117002(2013).
  * (24) Nakosai, S., Tanaka, Y. & Nagaosa, N. Topological superconductivity in bilayer Rashba system. Phys. Rev. Lett. 108, 147003 (2012).
  * (25) Wong, C. L. M., & Law, K. T. Majorana Kramers doublets in $d_{x^{2}-y^{2}}$-wave superconductors with Rashba spin-orbit coupling. Phys. Rev. B 86, 184516 (2012).
  * (26) Zhang, F., Kane, C. L., & Mele, E. J. Time-reversal-invariant topological superconductivity and Majorana Kramers pairs. Phys. Rev. Lett. 111, 056402(2013).
  * (27) Zhang, F., Kane, C. L., & Mele, E. J. Topological mirror superconductivity. Phys. Rev. Lett. 111, 056403 (2013).
  * (28) Deng, S., Viola, L., & Ortiz, G. Majorana modes in time-reversal invariant $s$-wave topological superconductors. Phys. Rev. Lett. 108, 036803 (2012)
  * (29) Zhao Y. X., & Wang, Z. D. Topological classification and stability of Fermi surfaces. Phys. Rev. Lett. 110, 240404 (2013).
  * (30) Béri, B. Topologically stable gapless phases of time-reversal-invariant superconductors. Phys. Rev. B 81, 134515 (2010).
  * (31) Sato, M., Tanaka, Y., Yada, K., & Yokoyama, T. Topology of Andreev bound states with flat dispersion. Phys. Rev. B 83, 224511 (2011).
  * (32) Gong, M., Tewari, S., & Zhang, C. BCS-BEC crossover and topological phase transition in 3D spin-orbit coupled degenerate Fermi gases. Phys. Rev. Lett. 107, 195303 (2011).
  * (33) Gong, M., Chen, G., Jia, S., & Zhang, C. Searching for Majorana fermions in 2D spin-orbit coupled fermi superfluids at finite temperature. Phys. Rev. Lett. 109, 105302 (2012).
  * (34) Hu, H., Jiang, L., Liu, X. J., & Pu, H. Probing anisotropic superfluidity in atomic Fermi gases with Rashba spin-orbit coupling. Phys. Rev. Lett. 107, 195304 (2011).
  * (35) Yu, Z. Q. & Zhai, H. Spin-orbit coupled Fermi gases across a Feshbach fesonance. Phys. Rev. Lett. 107, 195305 (2011).
  * (36) He, L. & Huang, X. G. BCS-BEC crossover in 2D Fermi gases with Rashba spin-orbit coupling. Phys. Rev. Lett. 108, 145302 (2012).
  * (37) Qu, C., Zheng, Z., Gong, M., Xu, Y., Mao, L., Zou, X., Guo, G., & Zhang, C. Topological superfluids with finite-momentum pairing and Majorana fermions,. Nat. Commun. 4, 2710 (2013).
  * (38) Zhang, W., & Yi, W. Topological Fulde-Ferrell-Larkin-Ovchinnikov states in spin-orbit-coupled Fermi gases. Nat. Commun. 4, 2711 (2013).
  * (39) Cao, Y., Zou, S. H., Liu, X.-J., Yi, S., Long, G. L., & Hu, H. Gapless topological Fulde-Ferrell superfluidity in spin-orbit coupled Fermi gases. Phys. Rev. Lett. 113, 115302 (2014).
  * (40) Dong, Y., Dong, L., Gong, M., & Pu, H. Dynamical phases in quenched spin-orbit-coupled degenerate Fermi gas. Nat. Comms. 6, 6103 (2015).
  * (41) Radic, J., Natu, S. S., & Galitski, V. Interaction-tuned dynamical transitions in a Rashba spin-orbit-coupled Fermi gas. Phys. Rev. Lett. 112, 095302 (2014).
  * (42) Pan, J. S., Liu, X. J., Zhang, W., Yi, W., & Guo, G. C. Topological superradiant states in a degenerate Fermi gas. Phys. Rev. Lett. 115, 045303 (2015).
  * (43) Clogston, A. M. Upper limit for the critical field in hard superconductors. Phys. Rev. Lett. 9, 266 (1962).
  * (44) Chandrasekhar, B. S. A note on the maximum ciritical field of high-field superconductor. Appl. Phys. Lett. 1, 7 (1962).
  * (45) Bloch, I. Ultracold quantum gases in optical lattices. Nat. Phys. 1, 23-30 (2005).
  * (46) Cheinet, P., Trotzky, S., Feld, M., Schnorrberger, U., Moreno-Cardoner, M., Fölling, S., & Bloch, I. Counting atoms using interaction blockade in an optical superlattice. Phys. Rev. Lett. 101, 090404 (2008).
  * (47) Eisenstein, J. P. & MacDonald, A. H. Bose-Einstein condensation of excitons in bilayer electron systems. Nature 432, 691 - 694 (2004).
  * (48) Qiao, Z., Tse, W.-K., Jiang, H., Yao, Y., & Niu, Q. Two-dimensional topological insulator state and topological phase transition in bilayer graphene. Phys. Rev. Lett. 107, 256801 (2011).
  * (49) Lin, Y., Jimenez-Garca, K., & Spielman, I. Spin-orbit-coupled Bose-Einstein condensates. Nature 471, 83 (2011).
  * (50) Zhang, J., Ji, S., Chen, Z., Zhang, L., Du, Z., Yan, B., Pan, G., Zhao, B., Deng, Y., Zhai, H., Chen, S. & Pan, J. Collective dipole oscillations of a spin-orbit coupled Bose-Einstein condensate. Phys. Rev. Lett. 109, 115301 (2012).
  * (51) Wang, P., Yu, Z., Fu, Z., Miao, J., Huang, L., Chai, S., Zhai, H., & Zhang, J. Spin-orbit coupled degenerate Fermi gases. Phys. Rev. Lett. 109, 095301 (2012).
  * (52) Cheuk, L., Sommer, A., Hadzibabic, Z., Yefsah, T., Bakr, W., & Zwierlein, M. Raman-induced interactions in a single-component Fermi gas near an $s$-wave Feshbach resonance. Phys. Rev. Lett. 111, 095301 (2013).
  * (53) Qu, C., Hamner, C., Gong, M., Zhang, C., & Engels, P. Observation of Zitterbewegung in a spin-orbit coupled Bose-Einstein condensate. Phys. Rev. A 88, 021604 (2013).
  * (54) Hamner, C., Qu, C., Zhang, Y., Chang, J., Gong, M., Zhang, C., & Engels, P. Dicke-type phase transition in a spin-orbit-coupled Bose-Einstein condensate. Nat. Commun. 5, 4023 (2014).
  * (55) Huang, L., Meng, Z., Wang, P., Peng, P., Zhang, S., Chen, L., Li, D., Zhou, Q., & Zhang, J. Experimental realization of a two-dimensional synthetic spin-orbit coupling in ultracold Fermi gases. arXiv:1506.02861.
  * (56) Chiu, C., Yao, H., & Ryu, S. Classification of topological insulators and superconductors in the presence of reflection symmetry. Phys. Rev. B 88, 075142(2013).
  * (57) Yao, H., & Ryu, S. Interaction effect on topological classification of superconductors in two dimensions. Phys. Rev. B 88, 064507 (2013).
  * (58) Ueno, Y., Yamakage, A., Tanaka, Y., & Sato, M. Symmetry-protected Majorana fermions in topological crystalline superconductors: theory and application to Sr2RuO4. Phys. Rev. Lett. 111, 087002 (2013).
  * (59) Shiozaki, K., & Sato, M. Topology of crystalline insulators and superconductors. Phys. Rev. B 90, 165114 (2014) .
  * (60) Tewari, S., & Sau, J. D. Topological invariants for spin-orbit coupled superconductor nanowires. Phys. Rev. Lett. 109, 150408 (2012).
  * (61) Bedaque, P. F., Caldas, H., & Rupak, G. Phase separation in asymmetrical fermion superfluids. Phys. Rev. Lett. 91, 247002 (2003).
  * (62) Yi, W., & Guo, G. C. Phase separation in a polarized Fermi gas with spin-orbit coupling. Phys. Rev. A 84, 031608(R) (2011).
  * (63) Qi, X., Hughes, T. L., & Zhang, S. Topological invariants for the Fermi surface of a time-reversal-invariant superconductor. Phys. Rev. B 81, 134508 (2010).
  * (64) Nakada, K., Fujita, M., Dresselhaus, G., & Dresselhaus, M.S. Edge state in graphene ribbons: nanometer size effect and edge shape dependence. Phys. Rev. B 54, 17954 (1996).
  * (65) Miyake, H., Siviloglou, G. A., Kennedy, C. J., Burton, W. C., & Ketterle, W. Realizing the Harper Hamiltonian with laser-assisted tunneling in optical lattices. Phys. Rev. Lett. 111, 185302 (2013).
  * (66) Jiménez-García, K., LeBlanc, L. J., Williams, R. A., Beeler, M. C., Qu, C., Gong, M., Zhang, C., & Spielman, I. B. Tunable spin-orbit coupling via strong driving in ultracold-atom systems. Phys. Rev. Lett. 114, 125301 (2015).
  * (67) Valle, G. D., Ornigotti, M., Cianci, E., Foglietti, V., Laporta, P., & Longhi, S. Visualization of coherent destruction of tunneling in an optical double well system. Phys. Rev. Lett. 98, 263601 (2007).
  * (68) Ezawa, M., Tanaka, Y., Nagaosa, N. Topological phase transition without gap closing. Scientific Reports 3, 2790 (2013).
  * (69) Huang, B. B., Chan, C. F., & Gong, M. Large Chern-number topological superfluids in a coupled-layer system. Phys. Rev. B 91, 134512 (2015).
  * (70) Deutsch, I. H. & Jessen, P. S. Quantum-state control in optical lattices. Phys. Rev. A 57, 1972 (1998).