Topological orbital superfluid with chiral d-wave order in a rotating optical lattice

Orbital degrees of freedom play a significant role in producing various exotic quantum states in complex condensed-matter systems, such as high temperature superconductors and quantum magnetic insulators. Recent experimental realizations of multi-orbital systems with ultra-cold atoms [1–4] have promoted the theoretical studies of high orbital physics in optical lattices, where a series of exotic quantum states have been proposed [5–11] . Among them, one of the remarkable characteristics is that the orbital hybridization can play the same role as spin–orbital couplings or artificial gauge fields which are the key ingredients needed to drive topologically insulating or superconducting states [12, 13]. Therefore, topologically nontrivial many-body states can be implemented in multi-orbital systems in the absence of spin–orbital coupling. Several methods exist to induce orbital hybridization in the context of cold-atom systems, including the many-body interaction effect [5], lattice shaking [14–17], and local rotation [18]. The relevant quantum states including topological semimetal [5] and topological band insulators [10, 15, 20, 21] have been proposed.

Recently, the superfluid of bosons with chiral odd-frequency orders, i.e., p + ip-wave and f + if-wave, have been experimentally realized in multi-orbital cold-atom systems [2, 22, 23]. For fermions, however, it is still a big challenge to realize superfluid states with chiral odd-frequency orders because the atom loss is strong near the Feshbach resonance in high-frequency channels [24]. Theoretically, thanks to the Rashba spin–orbital couplings, the topological superfluids of fermions with chiral odd-frequency orders have been proposed to emerge in the s-wave channel of the Feshbach resonance [25–28]. In comparison with the well-studied chiral odd-frequency superfluids of fermions, the superfluids of fermions with chiral even-frequency orders are rarely studied, and only some candidate materials are proposed to have the chiral even-frequency orders due to the unconventional superconducting pairing in condensed-matter systems [29–32]. More recently, a checkerboard lattice in a periodic Floquet driving field was proposed to support the chiral d-wave superfluid, where the sublattice degrees of freedom plays a key role and the periodic Floquet driving field induces the hybridization of two sublattices [33]. In this paper, we propose that a superfluid state of fermions with a chiral d-wave order can be implemented in a rotating multi-orbital optical lattice. In our proposal, the key ingredients needed to drive the underlying nontrivial topology of the multi-orbital superfluid state with a chiral d-wave order come from the two orbitals that are the counterparts of spin degrees of freedom in spin–orbital coupling, and the inter-orbital hybridization is induced by the local rotation with same frequency for every individual lattice site, which can be experimentally realized [18]. Interestingly, different from conventional chiral d-wave topological superfluid which supports two chiral Andreev edge modes at the boundary of the system, the topological orbital superfluid here supports four chiral Andreev edge modes due to the conservation of spin. More importantly, we find that the spatial barrier structure spontaneously formed by the intrinsic harmonic confining potential separates the trivial and nontrivial superfluid states, accumulates cold atoms and supports a mass current under injection of small angular momentum as the external driving force. These features can be experimentally adopted to verify the topologically nontrivial superfluid states. In comparison with the chiral p-wave and f-wave topological superconductor and superfluid [25, 26, 28, 34–40], where the spin–orbital couplings are essential, the chiral d-wave topological superfluid here only requires orbital hybridization. Therefore, our proposal provides a possible route to exploring topological superfluids with chiral d-wave order in multi-orbital cold-atom systems.

The paper is organized as follows. In section 2, we discuss the implementation of the multi-orbital system with a specific configuration of laser beams and construct the effective Hamiltonian to describe the multi-orbital system. In section 3, we study the homogeneous superfluid state with self-consistent mean-field approximation, and discuss the topological properties of the homogeneous superfluid state. In section 4, we discuss the inhomogeneous superfluid state modulated by the harmonic confining potential. In section 5, we discuss the experimental scheme and present a brief summary.

We consider a balanced mixture of fermion atoms with two internal states labeled by the spin index σ. The atoms are loaded in isotropic 2D square optical lattices. To introduce the couplings between different p orbital bands, one effective approach is to rotate the optical lattice with the same rotation frequency ${{\rm{\Omega }}}_{z}$ for every individual lattice site [18]. An alternative approach would be to directly couple the states with a drive laser [19]. Finally, the trapped atoms are tuned close to a Feshbach resonance to produce attractive s-wave interactions. The lattice potential takes the form,

$$\begin{eqnarray}&&V(x,y)={V}_{1}[\cos {k}_{L}x+\cos {k}_{L}y]+2{V}_{2}\cos {k}_{L}x\cos {k}_{L}y.\end{eqnarray} \tag{ 1 }$$

Here, V₁ and V₂ are the optical lattice potentials and k_L is the wave-vector of the laser fields. The realization of lattice potential $V(x,y)$ in equation (1) has been proposed for the case ${V}_{2}/{V}_{1}\gt 1/2$ [5]. Here, we consider the case ${V}_{2}/{V}_{1}\lt 1/2$, and the configuration of optical lattices under the condition that ${V}_{2}/{V}_{1}\lt 1/2$ can be implemented through four retro-reflected laser beams as shown in figure 1(a). The electric field generated by each laser beam is

$$\begin{eqnarray}&&{\vec{E}}_{j}(\vec{r},t)={E}_{j,0}{\vec{e}}_{j}\cos ({\vec{k}}_{j}\cdot {\vec{r}}_{j}){{e}}^{-{i}({\omega }_{j}t+{\varphi }_{j})},\end{eqnarray} \tag{ 2 }$$

where ${\vec{e}}_{j}$, ${\omega }_{j}$, and ${\varphi }_{j}$ are the polarization vector, the frequency, and the phase of the laser field, respectively. The parameters for each laser beams are summarized in table 1. The corresponding light-shift potential is

$$\begin{eqnarray}&&V(x,y)=-\chi {\left|\displaystyle \sum _{j}{\vec{E}}_{j}(\vec{r},t)\right|}^{2},\end{eqnarray} \tag{ 3 }$$

with χ denoting the real part of the polarizability. By adopting the parameters in table 1, we can get the lattice potential shown in equation (1) with an irrelevant constant shift ${V}_{0}=-\chi ({\epsilon }_{1}^{2}+{\epsilon }_{2}^{2})$. Here, ${V}_{1}=-\chi ({\epsilon }_{1}^{2}/2+{\epsilon }_{2}^{2})$, and ${V}_{2}=-\chi {\epsilon }_{2}^{2}/2$. The condition ${V}_{2}/{V}_{1}\lt 1/2$ can be achieved for arbitrary nonzero ${\epsilon }_{1}$ and ${\epsilon }_{2}$ and blue detuning with $\chi \lt 0$. Here, we set ${V}_{1}=1.2{E}_{R}$ and ${V}_{2}=0.4{E}_{R}$ . ${E}_{R}=\tfrac{{h}^{2}}{2{{ma}}^{2}}$ is the recoil energy and a is the lattice constant.

Zoom In Zoom Out Reset image size

Table 1.  The parameters of the electric fields of four laser beams shown in figure 1(a).

j	Ej0	${\vec{e}}_{j}$	${\vec{k}}_{j}$	${\omega }_{j}$	${\varphi }_{j}$
1	${\epsilon }_{1}$	$(1,0,0)$	$({k}_{L}/2,0)$	${\omega }_{0}$	0
2	${\epsilon }_{1}$	$(0,1,0)$	$(0,{k}_{L}/2)$	${\omega }_{0}$	0
3	${\epsilon }_{2}$	$(0,0,1)$	$({k}_{L}/2,{k}_{L}/2)$	${\omega }_{0}$	0
4	${\epsilon }_{2}$	$(0,0,1)$	$({k}_{L}/2,-{k}_{L}/2)$	${\omega }_{0}$	0

The contour of $V(x,y)$ is shown in figure 1(b). The lowest four band structures from the plane wave expansion approximation upon the potential $V(x,y)$ in equation (1) are shown in figure 1(d). It is straightforward to check that the splitting between two middle p_x and p_y bands off the high-symmetry point are induced by the coupling to the higher ${d}_{{x}^{2}-{y}^{2}}$ band [5]. Consider the three orbitals of p_x, p_y and ${d}_{{x}^{2}-{y}^{2}}$ shown in figure 1(e), a tight-binding (TB) Hamiltonian can be constructed to described the band structures of the fermionic square lattice, i.e.,

$$\begin{eqnarray}&&{H}_{{\rm{T}}{\rm{B}}}={H}_{d}+{H}_{p}+{H}_{{dp}},\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&{H}_{d}=\displaystyle \sum _{\langle i,j\rangle \sigma }[-{t}_{{dd}}+(\delta -{\mu }_{i}){\delta }_{{ij}}]{d}_{i,\sigma }^{\dagger }{d}_{j,\sigma },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{H}_{p} & = & -\displaystyle \sum _{i\sigma ,l=x,y}{\mu }_{i}{p}_{l,i,\sigma }^{\dagger }{p}_{l,i,\sigma }+{ih}{{\rm{\Omega }}}_{z}\displaystyle \sum _{i\sigma }{p}_{x,i,\sigma }^{\dagger }{p}_{y,i,\sigma }+H.c.\\ & & +{t}_{{pp}}\displaystyle \sum _{i\sigma ,l=x,y}{p}_{l,i,\sigma }^{\dagger }{p}_{l,i+{e}_{l},\sigma }+H.c.\\ & & -{t}_{{pp}}^{{\prime} }\displaystyle \sum _{i\sigma ,l=x,y,\bar{l}=-l}{p}_{l,i,\sigma }^{\dagger }{p}_{l,i+{e}_{\bar{l}},\sigma }+H.c.,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{H}_{{dp}}={t}_{{pd}}\displaystyle \sum _{i\sigma ,l=x,y}[{p}_{l,i+{e}_{l},\sigma }^{\dagger }{d}_{i,\sigma }-{p}_{l,i-{e}_{l},\sigma }^{\dagger }{d}_{i,\sigma }]+H.c.\end{eqnarray} \tag{ 7 }$$

Here, ${\mu }_{i}={\mu }_{0}+{V}_{{\rm{trap}}}({i}_{x},{i}_{y})$ with

$$\begin{eqnarray}&&{V}_{{\rm{trap}}}({i}_{x},{i}_{y})={V}_{t}\left[{\left({i}_{x}-\displaystyle \frac{{N}_{x}+1}{2}\right)}^{2}+{\left({i}_{y}-\displaystyle \frac{{N}_{y}+1}{2}\right)}^{2}\right]\end{eqnarray} \tag{ 8 }$$

being the weak harmonic confining potential to stabilize the optical lattice. ${p}_{x/y,i,\sigma }^{\dagger }$, and ${d}_{i,\sigma }^{\dagger }$ are the fermion creation operators for atoms in the relevant p_x, p_y and ${d}_{{x}^{2}-{y}^{2}}$ orbitals. We first set ${V}_{t}=0$ to simplify the discussions and recover it later. Note that all the energy scales are measured in the unit of t_pp as explained in the caption of figure 1 in the following parts of the paper if not specified. The energy spectra of TB Hamiltonian in equation (4) are shown in figures 1(e) and (f). It can be found that the TB Hamiltonian in equation (4) gives a good description of the band structures of lattice potential, and the on-site rotation in the second term in equation (6) induces the orbital hybridization to break the degeneracy of p_x and p_y bands around the Γ and M points shown in figure 1(c).

When the fermion atoms are loaded into the two p_x and p_y bands, the attractive s-wave interactions from the Feshbach resonance provide the two-orbital attractive Hubbard interactions as follows [41]:

$$\begin{eqnarray}{H}_{{\rm{int}}} & = & U\displaystyle \sum _{{il}}{n}_{{il}\uparrow }{n}_{{il}\downarrow }-\displaystyle \frac{J}{2}\displaystyle \sum _{i}\left[2{{\bf{S}}}_{{ix}}\cdot {{\bf{S}}}_{{iy}}+\displaystyle \frac{1}{2}{n}_{{ix}}{n}_{{iy}}\right]\\ & & +\displaystyle \frac{J}{2}\displaystyle \sum _{i}{n}_{{ix}}{n}_{{iy}}+{J}_{{\rm{\Delta }}}\displaystyle \sum _{i}{p}_{{ix}\uparrow }^{\dagger }{p}_{{ix}\downarrow }^{\dagger }{p}_{{iy}\downarrow }{p}_{{iy}\uparrow }+H.c.\end{eqnarray} \tag{ 9 }$$

Here, the first term is the intra-orbital attractive interaction, and the second term is the Hund's coupling with the spin operator ${{\bf{S}}}_{{il}}=\tfrac{1}{2}{p}_{{il},\alpha }^{\dagger }{\sigma }_{\alpha \beta }{p}_{{il},\beta }$ and $l=x,y$. U and J take the following forms:

$$\begin{eqnarray}&&U=4\pi {{\hslash }}^{2}{a}_{s}/m\int {\rm{d}}{r}| {\omega }_{x/y}(r){| }^{4},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&J=4\pi {{\hslash }}^{2}{a}_{s}/m\int {\rm{d}}{r}| {\omega }_{x}(r){| }^{2}| {\omega }_{y}(r){| }^{2}.\end{eqnarray} \tag{ 11 }$$

Here, a_s is the s-wave scattering length with negative value, i.e., ${a}_{s}\lt 0$. ${\omega }_{x/y}(r)$ are the Wannier functions of the ${p}_{x/y}$ orbitals. The third term in equation (9) is the inter-orbital attractive interaction with ${n}_{{il}}={n}_{{il},\uparrow }+{n}_{{il},\downarrow }$. The fourth term is the pair hopping term. Furthermore, we have $J=2U/3$ and ${J}_{{\rm{\Delta }}}=U/3$ [41]. Note that the Hund's coupling and inter-orbital interaction have same amplitudes, which are different from the electron system. The interaction terms shown in equations (9)–(11) are obtained under harmonic approximation. It is shown that the anharmonicity of the optical lattice can affect the properties of the multi-orbital system [42, 43]. In particular, the intra-orbital interaction U_xx is not equal to U_yy, and the inter-orbital interaction J is off $2{U}_{{xx}}/3$. Such imbalance can induce the modulations of superfluid order parameters. However, the topological superfluid is robust against such small modulations because nontrivial topology is the global feature of superfluid. For simplification, we neglect the irrelevant anharmonic effects in the present work.

Now, we turn to consider the homogeneous superfluid state with ${V}_{t}=0$ in equation (8) and the superfluid state is driven by the attractive interaction in equation (9). The spin–singlet superfluid pairing operators are defined as

$$\begin{eqnarray}&&{\hat{{\rm{\Delta }}}}_{s,{{ll}}^{{\prime} }}(k)=\displaystyle \sum _{\sigma {\sigma }^{{\prime} }}\displaystyle \frac{{[i{\sigma }_{y}]}_{\sigma {\sigma }^{{\prime} }}}{4}[{p}_{l,k\sigma }{p}_{{l}^{{\prime} },-k{\sigma }^{{\prime} }}+{p}_{{l}^{{\prime} },k\sigma }{p}_{l,-k{\sigma }^{{\prime} }}].\end{eqnarray} \tag{ 12 }$$

Then, we have

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=U\displaystyle \sum _{l}{\hat{{\rm{\Delta }}}}_{s,{ll}}^{\dagger }{\hat{{\rm{\Delta }}}}_{s,{ll}}+{J}_{{\rm{\Delta }}}\displaystyle \sum _{l\ne {l}^{{\prime} }}{\hat{{\rm{\Delta }}}}_{s,{ll}}^{\dagger }{\hat{{\rm{\Delta }}}}_{s,{l}^{{\prime} }{l}^{{\prime} }}+2J\displaystyle \sum _{l\gt {l}^{{\prime} }}{\hat{{\rm{\Delta }}}}_{s,{{ll}}^{{\prime} }}^{\dagger }{\hat{{\rm{\Delta }}}}_{s,{{ll}}^{{\prime} }}\end{eqnarray} \tag{ 13 }$$

with

$$\begin{eqnarray}&&{\hat{{\rm{\Delta }}}}_{s,{{ll}}^{{\prime} }}=\displaystyle \sum _{k}{\hat{{\rm{\Delta }}}}_{s,{l}^{{\prime} }{l}^{{\prime} }}(k).\end{eqnarray} \tag{ 14 }$$

Note that the spin–triplet pairing parts disappear because the Hund's coupling and inter-orbital interaction have the same amplitudes. Through the mean-field approximation, ${{\rm{\Delta }}}_{s,{{ll}}^{{\prime} }}=\langle {\hat{{\rm{\Delta }}}}_{s,{{ll}}^{{\prime} }}\rangle$, H_int can be decoupled to be

$$\begin{eqnarray}&&{H}_{{\rm{int}}}^{p}=\displaystyle \sum _{l,k}(U{{\rm{\Delta }}}_{s,{ll}}+{J}_{{\rm{\Delta }}}{{\rm{\Delta }}}_{s,\bar{l}\bar{l}}){\hat{{\rm{\Delta }}}}_{s,{ll}}^{\dagger }(k)+H.c.\,+\,\displaystyle \sum _{k}2J{{\rm{\Delta }}}_{s,{xy}}{\hat{{\rm{\Delta }}}}_{s,{xy}}^{\dagger }(k)+H.c.+{h}_{{\rm{con}}}\end{eqnarray} \tag{ 15 }$$

with

$$\begin{eqnarray}&&{h}_{{\rm{con}}}=-U\displaystyle \sum _{l}| {{\rm{\Delta }}}_{s,{ll}}{| }^{2}-2J| {{\rm{\Delta }}}_{s,{xy}}{| }^{2}-2{J}_{{\rm{\Delta }}}\mathrm{Re}({{\rm{\Delta }}}_{s,{xx}}{{\rm{\Delta }}}_{s,{yy}}^{* }).\end{eqnarray} \tag{ 16 }$$

The homogeneous superfluid state can be described by the mean-field Hamiltonian on the Nambu basis: ${\rm{\Psi }}(k)={[{d}_{k\uparrow },{p}_{x,k\uparrow },{p}_{y,k\uparrow },{d}_{k\downarrow }^{\dagger },{p}_{x,k\downarrow }^{\dagger },{p}_{y,k\downarrow }^{\dagger },{d}_{k\downarrow },{p}_{x,k\downarrow },{p}_{y,k\downarrow },{d}_{k\uparrow }^{\dagger },{p}_{x,k\uparrow }^{\dagger },{p}_{y,k\uparrow }^{\dagger }]}^{t}$,

$$\begin{eqnarray}{H}_{{\rm{mf}}}=\displaystyle \sum _{k}\displaystyle \frac{1}{2}{{\rm{\Psi }}}^{\dagger }(k)\left[\begin{array}{cccc}{H}_{{\rm{T}}{\rm{B}}}(k) & {\boldsymbol{\Delta }} & & \\ {{\boldsymbol{\Delta }}}^{\dagger } & -{H}_{{\rm{T}}{\rm{B}}}^{* }(-k) & & \\ & & {H}_{{\rm{T}}{\rm{B}}}(k) & -{\boldsymbol{\Delta }}\\ & & -{{\boldsymbol{\Delta }}}^{\dagger } & -{H}_{{\rm{T}}{\rm{B}}}^{* }(-k)\end{array}\right]{\rm{\Psi }}(k)+C.\end{eqnarray} \tag{ 17 }$$

Here, C is an operator-independent constant term. ${\boldsymbol{\Delta }}$ is a 3 × 3 matrix and takes the following form:

$$\begin{eqnarray}{\boldsymbol{\Delta }}=\left[\begin{array}{ccc}0 & 0 & 0\\ 0 & U{{\rm{\Delta }}}_{s,{xx}}+{J}_{{\rm{\Delta }}}{{\rm{\Delta }}}_{s,{yy}} & 2J{{\rm{\Delta }}}_{s,{xy}}\\ 0 & 2J{{\rm{\Delta }}}_{s,{xy}} & U{{\rm{\Delta }}}_{s,{yy}}+{J}_{{\rm{\Delta }}}{{\rm{\Delta }}}_{s,{xx}}\end{array}\right].\end{eqnarray} \tag{ 18 }$$

The mean-field Hamiltonian in equation (17) can be self-consistently solved with respect to the minimum of ground state energy, i.e.,

$$\begin{eqnarray}&&{E}_{g}={h}_{{\rm{con}}}-\displaystyle \frac{1}{4{\pi }^{2}}\displaystyle \sum _{n=1}^{3}\int {{\rm{d}}}^{2}{\bf{k}}[| {E}_{n}^{(s)}(k)| -| {E}_{n}^{(0)}(k)| ],\end{eqnarray} \tag{ 19 }$$

where, ${E}_{n}^{(s)}(k)$ and ${E}_{n}^{(0)}(k)$ are the eigen-energy spectra of the superfluid state and normal state. Here, we focus on the filling lying in the band splitting around the M point induced by the orbital hybridization as shown in figure 2(a). The typical Fermi surface is shown in figure 2(b). From equation (18), we can find that the superfluid order parameter in the intra-p_x orbital channel is ${{\rm{\Delta }}}_{22}=U{{\rm{\Delta }}}_{s,{xx}}+{J}_{{\rm{\Delta }}}{{\rm{\Delta }}}_{s,{yy}}$ while the superfluid order parameter in the intra-p_y orbital channel is ${{\rm{\Delta }}}_{33}=U{{\rm{\Delta }}}_{s,{yy}}+{J}_{{\rm{\Delta }}}{{\rm{\Delta }}}_{s,{xx}}$. To maximize the superfluid gap, one can find that ${{\rm{\Delta }}}_{s,{xx}}{{\rm{\Delta }}}_{s,{yy}}\gt 0$ is favorable to obtain the largest amplitudes of ${{\rm{\Delta }}}_{22}$ and ${{\rm{\Delta }}}_{33}$. The numerical results for the ground state energy and superfluid order parameters as functions of chemical potential ${\mu }_{0}$ and interaction amplitude $| U|$ are shown in figures 2(c) and (d), from which the intra-orbital ${{\rm{\Delta }}}_{22}$ and ${{\rm{\Delta }}}_{33}$ are degenerate in the whole parameter regime. This means that ${{\rm{\Delta }}}_{22}={{\rm{\Delta }}}_{33}$, and the only choice is ${{\rm{\Delta }}}_{s,{xx}}{{\rm{\Delta }}}_{s,{yy}}\gt 0$ thanks to ${{UJ}}_{{\rm{\Delta }}}\gt 0$. The aforementioned analyses are consistent, and one can achieve that the superfluid ground states favor ${{\rm{\Delta }}}_{s,{xx}}$ and ${{\rm{\Delta }}}_{s,{yy}}$ with the same sign to maximize the superfluid gap and to minimize the ground state energy. Furthermore, we can find that the inter-orbital ${{\rm{\Delta }}}_{23}$, which is also the matrix element in equation (18), is purely imaginary, and much smaller than ${{\rm{\Delta }}}_{22/33}$. The reason lies in that the inter-orbital ${{\rm{\Delta }}}_{23}$ is induced by the orbital hybridization and modulated by ${{\rm{\Omega }}}_{z}$. It is conceivable that the strength of inter-orbital ${{\rm{\Delta }}}_{23}$ could be comparable to intra-orbital ${{\rm{\Delta }}}_{22/33}$ when ${{\rm{\Omega }}}_{z}$ is large enough. However, the ${{\rm{\Delta }}}_{23}$ has no relation to the topological nature of the superfluid state, so we only focus on the case with ${{\rm{\Omega }}}_{z}$ set here.

Zoom In Zoom Out Reset image size

In order to reveal the underlying topological nature of the superfluid states, we first investigate the band characteristics of the normal states. As shown in figure 1(e), the full separation between the d band and p bands guarantees the feasibility of downfolding the Hamiltonian from the space spanned by the d and p orbitals to the space spanned by the two effective $\tilde{p}$ orbitals shown in figure 1(f). When ${V}_{t}=0$, the translation symmetry allows ones to write the TB Hamiltonian in momentum space under the effective basis ${\tilde{\psi }}_{\sigma }(k)={[{\tilde{p}}_{x,k,\sigma },{\tilde{p}}_{y,k,\sigma }]}^{t}$, i.e.,

$$\begin{eqnarray}&&{\tilde{H}}_{{\rm{T}}{\rm{B}}}=\displaystyle \sum _{k\sigma }{\tilde{\psi }}_{\sigma }^{\dagger }(k){\tilde{H}}_{{\rm{T}}{\rm{B}}}(k){\tilde{\psi }}_{\sigma }(k).\end{eqnarray} \tag{ 20 }$$

Here,

$$\begin{eqnarray}&&{\tilde{H}}_{{\rm{T}}{\rm{B}}}(k)=\displaystyle \frac{1}{2}{\xi }_{+}(k)-{\mu }_{0}+{\xi }_{{xy}}(k){\sigma }_{x}-h{{\rm{\Omega }}}_{z}{\sigma }_{y}+\displaystyle \frac{1}{2}{\xi }_{-}(k){\sigma }_{z},\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&{\xi }_{\pm }(k)=2({\tilde{t}}_{{pp}}\mp {\tilde{t}}_{{pp}}^{{\prime} })(\cos {k}_{x}\pm \cos {k}_{y}),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\xi }_{{xy}}(k)=4{\tilde{t}}_{{xy}}\sin {k}_{x}\sin {k}_{y}.\end{eqnarray} \tag{ 23 }$$

The Pauli matrices ${\sigma }_{i}$ with $i=x,y,z$ span the two effective ${\tilde{p}}_{x}$ and ${\tilde{p}}_{y}$ orbital space. The effective TB Hamiltonian ${\tilde{H}}_{{tb}}$ can be rewritten in the basis spanned by the orbital angular momentum eigenstate, i.e.,

$$\begin{eqnarray}&&{\bar{H}}_{{\rm{T}}{\rm{B}}}=\displaystyle \sum _{k\sigma }{\bar{\psi }}_{\sigma }^{\dagger }(k){\bar{H}}_{{\rm{T}}{\rm{B}}}(k){\bar{\psi }}_{\sigma }(k).\end{eqnarray} \tag{ 24 }$$

Here, ${\bar{\psi }}_{\sigma }^{\dagger }(k)=[{\bar{p}}_{+,k,\sigma }^{\dagger },{\bar{p}}_{-,k,\sigma }^{\dagger }]$ with ${\bar{p}}_{\pm ,k\sigma }^{\dagger }=\tfrac{1}{\sqrt{2}}[{\tilde{p}}_{x,k\sigma }^{\dagger }\pm {i}{\tilde{p}}_{y,k\sigma }^{\dagger }]$, and

$$\begin{eqnarray}&&{\bar{H}}_{{\rm{T}}{\rm{B}}}(k)=\displaystyle \frac{1}{2}{\xi }_{+}(k)-{\mu }_{0}+\displaystyle \frac{1}{2}{\xi }_{-}(k){s}_{x}+{\xi }_{{xy}}(k){s}_{y}-h{{\rm{\Omega }}}_{z}{s}_{z}.\end{eqnarray} \tag{ 25 }$$

The Pauli matrices s_i with $i=x,y,z$ span the two effective ${\bar{p}}_{+}$ and ${\bar{p}}_{-}$ orbital space. In the absence of ${{\rm{\Omega }}}_{z}$, $[{F}_{x}(k)=\tfrac{1}{2}{\xi }_{-}(k),{F}_{y}(k)={\xi }_{{xy}}(k)]$ forms a vector field in momentum space shown in figure 2(b). Then, the band degeneracy point at the M point can be mapped into a vortex in the momentum space with integer winding number [5], i.e.,

$$\begin{eqnarray}&&{W}_{\sigma }=\oint \displaystyle \frac{{\rm{d}}{\bf{k}}}{2\pi }\left[\displaystyle \frac{{F}_{x}(k)}{F(k)}{\boldsymbol{\nabla }}\displaystyle \frac{{F}_{y}(k)}{F(k)}-(x\,\longleftrightarrow \,y)\right],\end{eqnarray} \tag{ 26 }$$

with $F(k)=\sqrt{{F}_{x}^{2}(k)+{F}_{y}^{2}(k)}$. The direct calculation gives ${W}_{\sigma }=2$ in agreement with the pattern of the vector field $[{F}_{x}(k),{F}_{y}(k)]$ as shown in figure 2(b). Note that the total winding number $W={W}_{\uparrow }+{W}_{\downarrow }$ should be 4 when the spin degree of freedom is taken into account. In the presence of ${{\rm{\Omega }}}_{z}$, the induced orbital hybridization lifts the degeneracy at M point. Then, the above mapping does not work.

In the superfluid states, quasi-particle spectra are fully gapped and the nonzero ${{\rm{\Omega }}}_{z}$ breaks the pseudo-time-reversal symmetry. It is natural to introduce the Chern number to characterize the topological properties of the superfluid states. To show it, we consider the effective superfluid Hamiltonian spanned in the effective Nambu basis: $\bar{{\rm{\Psi }}}(k)={[{\bar{p}}_{+,k,\uparrow },{\bar{p}}_{-,k,\uparrow },{\bar{p}}_{+,-k,\downarrow }^{\dagger },{\bar{p}}_{-,-k,\downarrow }^{\dagger },{\bar{p}}_{+,k,\downarrow },{\bar{p}}_{-,k,\downarrow },{\bar{p}}_{+,-k,\uparrow }^{\dagger },{\bar{p}}_{-,-k,\uparrow }^{\dagger }]}^{t}$,

$$\begin{eqnarray}&&{\bar{H}}_{{\rm{mf}}}=\displaystyle \sum _{k}{\bar{{\rm{\Psi }}}}^{\dagger }(k)[{\bar{H}}_{{\rm{T}}{\rm{B}}}(k)+{\bar{H}}_{{\rm{int}}}^{p}(k)]\bar{{\rm{\Psi }}}(k),\end{eqnarray} \tag{ 27 }$$

with

$$\begin{eqnarray}{\bar{H}}_{{\rm{int}}}^{p}(k)={s}_{z}\otimes \left[\begin{array}{cc} & \bar{{\boldsymbol{\Delta }}}\\ {\bar{{\boldsymbol{\Delta }}}}^{\dagger } & \end{array}\right],\end{eqnarray} \tag{ 28 }$$

and

$$\begin{eqnarray*}\bar{{\boldsymbol{\Delta }}}=\left[\begin{array}{cc}| {{\rm{\Delta }}}_{{\rm{inter}}}| & {{\rm{\Delta }}}_{{\rm{intra}}}\\ {{\rm{\Delta }}}_{{\rm{intra}}} & -| {{\rm{\Delta }}}_{{\rm{inter}}}| \end{array}\right].\end{eqnarray*}$$

Here ${{\rm{\Delta }}}_{{\rm{intra}}}={{\rm{\Delta }}}_{22}$ and $| {{\rm{\Delta }}}_{{\rm{inter}}}| =| {{\rm{\Delta }}}_{23}|$. Upon an unitary rotation [26], we can obtain a dual form of the Hamiltonian, i.e.,

$$\begin{eqnarray}&&{\bar{H}}_{{\rm{mf}}}^{D}=S{\bar{H}}_{{\rm{mf}}}{S}^{\dagger },\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}S=\displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{cccc}1 & {s}_{x} & & \\ {{is}}_{y} & -{s}_{z} & & \\ & & 1 & -{s}_{x}\\ & & -{{is}}_{y} & -{s}_{z}\end{array}\right],\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray*}{\bar{H}}_{{\rm{mf}}}^{D}(k)=\left[\begin{array}{cc}{\bar{H}}_{{\rm{mf}}}^{D+}(k) & \\ & {\bar{H}}_{{\rm{mf}}}^{D-}(k)\end{array}\right]\end{eqnarray*}$$

$$\begin{eqnarray}{\bar{H}}_{{\rm{mf}}}^{D\pm }(k)=\left[\begin{array}{cc}{{\rm{\Delta }}}_{{\rm{intra}}}-h{{\rm{\Omega }}}_{z}{s}_{z} & \pm h(k)\\ \pm {h}^{\dagger }(k) & -{{\rm{\Delta }}}_{{\rm{intra}}}+h{{\rm{\Omega }}}_{z}{s}_{z}\end{array}\right],\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&h(k)={{is}}_{y}\left[-\displaystyle \frac{{\xi }_{+}(k)}{2}+{\mu }_{0}+\displaystyle \frac{{\xi }_{-}(k)}{2}{s}_{x}-{\xi }_{{xy}}(k){s}_{y}+{{is}}_{y}| {{\rm{\Delta }}}_{{\rm{inter}}}| \right].\end{eqnarray} \tag{ 32 }$$

In the dual Hamiltonian ${\bar{H}}_{{\rm{mf}}}^{D\pm }(k)$ shown in equation (31), $\left[\tfrac{{\xi }_{-}(k)}{2},{\xi }_{{xy}}(k)\right]$ resembles two components of pairing order parameters of the chiral d-wave superfluid and ${{\rm{\Delta }}}_{{\rm{inter}}}$ corresponds to the mixed s-wave component. '${{\rm{\Delta }}}_{{\rm{intra}}}\pm h{{\rm{\Omega }}}_{z}$' is the pseudo-kinetic energy with k-independent, and resembles kinetic energy term '$\tfrac{{k}^{2}\pm {k}_{F}^{2}}{2m}$' of the chiral d-wave superfluid when ${\mu }_{0}$ is set to satisfy the condition ${\mu }_{0}=\tfrac{1}{2}{\xi }_{+}(\pi ,\pi )$. Then, the dual Hamiltonian ${\bar{H}}_{{\rm{mf}}}^{D\pm }(k)$ resembles the standard Hamiltonian describing the chiral d-wave superconductors [30, 44], and belongs to class C according to the classification by Schnyder et al [45]. Here, ${{\rm{\Delta }}}_{{\rm{inter}}}$ by itself cannot drive the gap-closing condition because it is much smaller than ${{\rm{\Delta }}}_{{\rm{intra}}}$ and Fermi energy. Therefore, the small ${{\rm{\Delta }}}_{{\rm{inter}}}$ can be absorbed and set to zero. The topological nontrivial superfluid states can be achieved under the condition [44] ${{\rm{\Delta }}}_{{\rm{intra}}}\,\lt \,$ $h{{\rm{\Omega }}}_{z}$ when ${\mu }_{0}=\tfrac{1}{2}{\xi }_{+}(\pi ,\pi )$, which naturally corresponds to the weak-coupling condition$\tfrac{{k}^{2}-{k}_{F}^{2}}{2m}\lt 0$ [34]. For the general case with arbitrary ${\mu }_{0}$, one can obtain nontrivial superfluid states if $f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})\gt 0$ with $f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})$ shown in equation (33), and trivial superfluid states if $f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})\lt 0$. The topological phase transition condition coincides with the gap-closing condition with $f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})=0$. The phase diagram separating the topological trivial and nontrivial superfluid phases is plotted in figure 3(d) according to phase transition condition $f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})=0$.

$$\begin{eqnarray}&&f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})=| h{{\rm{\Omega }}}_{z}| -\sqrt{{{\rm{\Delta }}}_{{\rm{intra}}}^{2}+{\left[\displaystyle \frac{{\xi }_{+}(\pi ,\pi )}{2}-{\mu }_{0}\right]}^{2}}.\end{eqnarray} \tag{ 33 }$$

The nontrivial topological nature of the superfluid states can be characterized by the Chern number,

$$\begin{eqnarray}&&{{ \mathcal C }}_{s}=\displaystyle \frac{i}{2\pi }\displaystyle \sum _{{E}_{n}\lt 0}{\int }_{{BZ}}{\rm{d}}{\bf{k}}\langle {{\rm{\nabla }}}_{k}{u}_{s,n}({\bf{k}})| \,\times \,| {{\rm{\nabla }}}_{k}{u}_{s,n}({\bf{k}})\rangle ,\end{eqnarray} \tag{ 34 }$$

with ${u}_{s,n}({\bf{k}})$ the Bloch functions of occupied quasi-particle states with s = up and down to label the the up-block and down-block parts of the Hamiltonian in equation (17). The straightforward calculations give ${{ \mathcal C }}_{{\rm{up}}}={{ \mathcal C }}_{{\rm{down}}}=2$ for $h{{\rm{\Omega }}}_{z}\gt 0$ and ${{ \mathcal C }}_{{\rm{up}}}={{ \mathcal C }}_{{\rm{down}}}=-2$ for $h{{\rm{\Omega }}}_{z}\lt 0$ under the condition $f({{\rm{\Omega }}}_{z},{\mu }_{0},{{\rm{\Delta }}}_{{\rm{intra}}})\gt 0$, which means the inverse local rotation corresponds to reverse chirality. From the bulk-edge correspondence, the quasi-particle spectra have two chiral gapless edge states at the open boundary shown in figure 3(a) and no gapless edge states emerge in the trivial superfluid state shown in figure 3(b). The local feature of the edge states in figure 3(a) are explicitly demonstrated through the amplitude distributions of the wave-functions shown in figure 3(c).

Zoom In Zoom Out Reset image size

Now, we consider the realistic case with nonzero harmonic confining potential in equation (8), and the pattern of ${V}_{{\rm{trap}}}({i}_{x},{i}_{y})$ is shown in figure 4(a) with ${V}_{t}=1.2/{N}_{x}{N}_{y}$. We perform the self-consistent calculations about the Bogoliubov–de Gennes (BdG) Hamiltonian ${H}_{{\rm{T}}{\rm{B}}}+{H}_{{\rm{int}}}^{p}$ in equations (4) and (15) in lattice space. The quasi-particle spectra and the distribution of superfluid order parameters are shown in figures 4(b), (d), and (e) for two different $h{{\rm{\Omega }}}_{z}=0.2$ and $h{{\rm{\Omega }}}_{z}=0.05$ under the periodic boundary condition. We find that the amplitudes of superfluid order parameters in both cases are similar to figures 4(d) and (e), but the quasi-particle spectra are quite different from figure 4(b) with in-gap fermion modes for $h{{\rm{\Omega }}}_{z}=0.2$ and without in-gap fermion modes for $h{{\rm{\Omega }}}_{z}=0.05$. The reason lies in that ${V}_{{\rm{trap}}}({i}_{x},{i}_{y})$ forms a spatial barrier structure (The position is marked with a red-dashed circle in figure 4(a)) separating the nontrivial superfluid state with $f({{\rm{\Omega }}}_{z},{\mu }_{i},{{\rm{\Delta }}}_{{\rm{intra}}})\gt 0$ and trivial superfluid state with $f({{\rm{\Omega }}}_{z},{\mu }_{i},{{\rm{\Delta }}}_{{\rm{intra}}})\lt 0$ for ${{\rm{\Omega }}}_{z}=0.2$. Note that ${\mu }_{i}={\mu }_{0}+{V}_{{\rm{trap}}}({i}_{x},{i}_{y})$, thus the position of spatial barrier coincides with the gap-closing condition with $f({{\rm{\Omega }}}_{z},{\mu }_{i},{{\rm{\Delta }}}_{{\rm{intra}}})=0$ . For fixed ${\mu }_{0}$ and V_t, one can find that $f({{\rm{\Omega }}}_{z},{\mu }_{i},{{\rm{\Delta }}}_{{\rm{intra}}})$ is always smaller than zero when ${{\rm{\Omega }}}_{z}=0.05$. The superfluid is always trivial because ${{\rm{\Omega }}}_{z}=0.05$ is too small to overcome the gap-closing condition $f({{\rm{\Omega }}}_{z},{\mu }_{i},{{\rm{\Delta }}}_{{\rm{intra}}})=0$. The spatial barrier traps in-gap fermion modes and accumulates atoms when the negative energy states are occupied [44, 46]. The in-gap fermion modes trapped by the spatial barrier have the same origin as the fermion modes in spectrum of the Caroli–de Gennes–Matricon bound states in the vortex core [47].

Zoom In Zoom Out Reset image size

In the low-energy limit, the spectrum of in-gap fermion modes in terms of the angular momentum Q takes the following form under the axisymmetric condition [44, 46]:

$$\begin{eqnarray}&&{E}_{a}(Q)={\omega }_{a}(Q-{Q}_{a}),\end{eqnarray} \tag{ 35 }$$

where ${\omega }_{a}={c}_{a}/R$ is the angular velocity of the rotation along the spatial barrier with R the radius of spatial barrier of ${V}_{{\rm{trap}}}({i}_{x},{i}_{y})$ ⁶ , a labels the ath branch, and ${Q}_{a}={\hslash }{k}_{a}R$. The total number of the branches is four according to the index theorem [46] when the spin degree of freedom is taken into account. In the absence of external driving, the energy of in-gap fermion modes is ${E}_{a}(0)=-{\omega }_{a}{Q}_{a}$. In the square lattice space, the circular rotation symmetry SO(2) for equation (35) is broken down to C₄ symmetry, and the Fermi velocity is strongly anisotropic and the superfluid order parameters are highly inhomogeneous. Q_a can only take the discrete values under the constraint of C₄ symmetry. Correspondingly, the energy levels of the in-gap fermion modes trapped by the spatial barrier are discrete (see figure 4(b) for details), and several energy levels close to zero usually correspond to in-gap fermion modes trapped by the spatial barrier.

The localization feature of the in-gap fermion modes trapped by the spatial barrier can be reflected by the local density of states (LDOS), which is calculated by

$$\begin{eqnarray}&&{\rho }_{i}(\omega )=\displaystyle \sum _{n,l,\sigma }[| {u}_{i,l\sigma }^{n}{| }^{2}\delta ({E}_{n}-\omega )+| {v}_{i,l\sigma }^{n}{| }^{2}\delta ({E}_{n}+\omega )],\end{eqnarray} \tag{ 36 }$$

where ${u}_{i,l\sigma }^{n}$ and ${v}_{i,l\sigma }^{n}$ are the particle-like and hole-like components of eigenstate with quasi-particle energy E_n at site i and orbital l. The LDOS of the five in-gap fermion modes with the highest negative energy are shown in figures 5(a1)–(a5), from which we can find that four levels with energy $-0.0017,-0.0145,-0.0208,-0.0282$ are the fermion modes which are trapped by the spatial barrier. To make a comparison, the level with energy −0.0322 is the extended state. We also plot the LDOS of the the five levels with the lowest positive energy in figures 5(b1)–(b5) for comparison. Furthermore, we find that the highest negative energy level and the lowest positive energy level approach zero energy with increasing the lattice size N × N (see figure 4(c) for details).

Zoom In Zoom Out Reset image size

In the presence of external driving, the spectrum of the in-gap fermion modes is a function of the angular momentum Q from the external driving, and the in-gap fermion modes could cross the zero energy and form the variation of the mass current. The change of the mass current trapped in the spatial barrier is [44]

$$\begin{eqnarray}&&\delta {I}_{M}=\displaystyle \frac{\hslash }{8\pi }\displaystyle \sum _{a}\delta ({k}_{a}^{2}),\end{eqnarray} \tag{ 37 }$$

where we have assumed the thickness along the z direction to be unity. The extra 1/2 in the denominator is added to compensate for the double counting due to the particle-hole symmetry. Generally, there are several external perturbations which can be introduced to be the driving force to move the in-gap fermion modes across the zero energy, such as the modulations of V₁ and V₂ in equation (1) to deform the $\tfrac{{\xi }_{-}(k)}{2}$ and ${\xi }_{{xy}}(k)$ and introducing an additional laser beam to modulate the trapping potential. Here, we consider a more convenient method. From equation (35), it is straightforward to inject nonzero Q into the superfluid state through the slight modulation of local rotating frequency ${{\rm{\Omega }}}_{z}$. As a consequence, the in-gap fermion modes can be driven to cross zero energy by the nonzero $\delta {{\rm{\Omega }}}_{z}$. If we further assume that all the in-gap fermion modes trapped in the spatial barrier have the relation $\tfrac{{\hslash }^{2}\delta ({k}_{a}^{2})}{2m}\sim h\delta {{\rm{\Omega }}}_{z}$, we can obtain that the response of the change in mass current to the modulation of the rotating frequency $\delta {I}_{M}$ $\sim \tfrac{m\delta {{\rm{\Omega }}}_{z}}{2}{\sum }_{a}{sgn}({c}_{a})$ with the summation involving all the in-gap fermion modes across the zero energy. However, in the square lattices, we can find that the k_a is different for different a-th branches from figure 5. As a good approximation, we can define an effective $\langle k\rangle$ to remove the difference of different k_a, and $\langle k\rangle$ can be replaced with the average Fermi momentum $\langle {k}_{F}\rangle$. Then, we can obtain that the modulation of mass current density is proportional to the change in the LDOS, i.e.,

$$\begin{eqnarray}&&\delta {j}_{M}({i}_{x},{i}_{y})\propto \delta \rho ({i}_{x},{i}_{y}),\end{eqnarray} \tag{ 38 }$$

with

$$\begin{eqnarray}&&\delta \rho ({i}_{x},{i}_{y})=\rho ({i}_{x},{i}_{y}){| }_{{{\rm{\Omega }}}_{z}+\delta {{\rm{\Omega }}}_{z}}-\rho ({i}_{x},{i}_{y}){| }_{{{\rm{\Omega }}}_{z}},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&\rho ({i}_{x},{i}_{y}){| }_{{{\rm{\Omega }}}_{z}}=\displaystyle \sum _{n,l,\sigma }| {v}_{l\sigma }^{n}({i}_{x},{i}_{y}){| }^{2}\theta (-{E}_{n}){| }_{{{\rm{\Omega }}}_{z}}.\end{eqnarray} \tag{ 40 }$$

The pattern of $\delta \rho ({i}_{x},{i}_{y})$ for $\delta {{\rm{\Omega }}}_{z}=0.02/h$ is shown in figure 5(c), from which we can find that the mass current is trapped around the spatial barrier.

In terms of experiment, the fermion atoms can be selected as lithium ⁶Li, two internal states can be selected as ${}^{2}{S}_{1/2}$ with M = $\pm \tfrac{1}{2}$. The principal fluorescence line from ²S${}_{1/2}$ to ²P is at 670.8 nm. Therefore, a Nd:YAG-laser with 532 nm could be selected to be the light source to realize the optical potential with the lattice constant $a=532$ nm. The recoil energy ${E}_{R}\sim h\times 100$ kHz. The local rotation around each potential minimum has been experimentally realized through inserting electrooptic phase modulators into the beams forming the 2D lattice potential, and the relevant rotating frequency ${{\rm{\Omega }}}_{z}$ can be turned with large flexibility [18]. From the energy bands in figure 1, we can estimate that it is enough for ${{\rm{\Omega }}}_{z}\sim h\times 2$ KHz to satisfy the topological superfluid condition.

In the presence of the harmonic trap, it has been shown that the local density approximation(LDA) breaks down for trapped non-interacting bosons in the p-orbital bands, and increasing the interactions and optical lattice potentials can suppress anisotropy of condensate density [48] . However, the picture is different for trapped non-interacting fermions in p-orbital bands due to the different statistics. It is shown that the hard-core boson known as Tonks–Girardeau boson with infinitely repulsive interactions can be mapped into non-interacting free fermion in the 1D limit [49–51]. Thus, the boson with infinitely repulsive interactions is roughly equivalent to free fermion even in a 2D system. Such effective 'repulsive interactions' can suppress the anisotropy of condensate density, and guarantee the validity of LDA in a system with trapped fermions in p-orbital bands. Furthermore, the tunability of the optical lattice potential and quite small trap potential can further reduce the anisotropy of condensate density. Though the breaking down of LDA can be suppressed, the particle density per site will inevitably vary and the s-orbital atoms will thereby shift the on-site energies for p-orbital atoms in the presence of the trap. Thanks to the small trap potential, one can expect that the density fluctuations of the both trapped s-orbital and p-orbital atoms should be small, and the main results throughout the paper are not changed qualitatively.

The change of the mass current and the accumulation of the atoms around the spatial barrier can be spatially resolved with the radio-frequency spectroscopy [52–54]. Besides the radio-frequency spectroscopy, the recently developed matter-wave interference technique [55] is a more powerful tool, which can directly represent the phase properties of the superfluid order parameter. More remarkably, one can reconstruct the spatial geometry of certain low-energy in-gap fermion modes and verify the formation of the spatial barrier structure, both of which are the key signatures in our proposal.

In summary, we propose that the superfluid states of fermions with a chiral d-wave order can be implemented in a rotating optical lattice where the orbital degrees of freedom play a key role. Our proposal presents an alternative route to realizing the topological superfluids with chiral even-frequency order in the absence of the spin–orbital coupling. Furthermore, we show that the intrinsic harmonic confining potential can form a circular spatial barrier structure which accumulates atoms and support a mass current under the injection of small angular momentum as the driving force. The mass current associated with the accumulated atoms can be experimentally detected, and provides a signature to verify the emergence of topological superfluid state with chiral d-wave order in a rotating optical lattice.

We thank S-Q Shen, S Z Zhang, D W Zhang, and G C Liu for useful discussions. This work was supported by the NSFC under Grants No. 11674331, No. 11274032, No.11625415, the Ministry of Science and Technology of China (Grant number: 2017YFA0303200), and by the 100 Talents Programme of the CAS.