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

Topological superfluid is an exotic state of quantum matter that possesses a nodeless superfluid gap in the bulk and Andreev edge modes at the boundary of a finite system. Here, we study a multi-orbital superfluid driven by attractive s-wave interaction in a rotating optical lattice. Interestingly, we find that the rotation induces the inter- orbital hybridization and drives the system into topological orbital superfluid in accordance with intrinsically chiral d-wave pairing characteristics. Thanks to the conservation of spin, the topological orbital superfluid supports four rather than two chiral Andreev edge modes at the boundary of the lattice. Moreover, we find that the intrinsic harmonic confining potential forms a circular spatial barrier which accumulates atoms and supports a mass current under injection of small angular momentum as external driving force. This feature provides an experimentally detectable phenomenon to verify the topological orbital superfluid with chiral d-wave order in a rotating optical lattice.


I. INTRODUCTION
Orbital degrees of freedom play a significant role to produce various exotic quantum states in complex condensed-matter systems, such as high temperature superconductors and quantum magnetic insulators. Recent experimental realizations of multiorbital systems with ultra-cold atoms [1][2][3][4] have promoted the theoretical studies of high orbital physics in optical lattices, where a series of exotic quantum states have been proposed [5][6][7][8][9][10][11] . Among them, one of remarkable characteristics is that the orbital hybridization can play the same role as spin-orbital coupling or artificial gauge fields which are the key ingredient 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 couplings. There exist several methods to induce the orbital hybridization in the context of cold atom systems, including many-body interaction effect 5 , lattice shaking [14][15][16][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 the fermions, however, it is still a big challenge to realize the superfluid states with chiral odd-frequency orders, because the atom loss is strong near the Feshbach resonance in highfrequency 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 s-wave channel of the Feshbach resonance. [25][26][27][28] . In comparison with 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][30][31][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 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 non-trivial superfluid states. In comparison with the chiral p-wave and f -wave topological superconductor and superfluid 25,26,28,[34][35][36][37][38][39][40] , where the spin-orbital couplings are essential, the chiral d-wave topological superfluid here only requires the orbital hybridization. Therefore, our proposal provides a possible route to explore topological superfluids with chiral d-wave order in multi-orbital cold-atom systems.
The paper is organized as follows. In section II, 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 III, we study the homogeneous superfluid state with self-consistent mean-field approximation, and discuss the topological properties of the homogeneous superfluid state. In section IV, we discuss the inhomogeneous superfluid state modulated by the harmonic  Fig. 1(a).

II. OPTICAL LATTICE AND MODEL
We consider a balanced mixture of fermion atoms with two internal states labeled by the spin index σ. The atoms are loaded in an isotropic 2D square optical lattices. To introduce the couplings between different p orbital bands, one effective approach is to rotate the optical lattice with same rotation frequency Ω 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 turned close to a Feshbach resonance to produce attractive s-wave interactions. The lattice potential takes the form, Here, V 1 and V 2 are the optical lattice potentials and k L is the wave-vector of laser fields. The realization of lattice potential V (x, y) in Eq. (1) has been proposed for the case V 2 /V 1 > 1/2 5 . Here, we consider the case V 2 /V 1 < 1/2, and the configuration of optical lattices under the condition V 2 /V 1 < 1/2 can be implemented through four retro-reflected laser beams as shown in Fig. 1(a). The electric field generated by each laser beam is where e j , ω j , and ϕ 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 I. The corresponding light-shift potential is with χ denoting the real part of the polarizability. By adopting the parameters in Table I, we can get the lattice potential shown in Eq.(1) with an irrelevant constant shift V 0 = −χ(ǫ 2 1 + ǫ 2 2 ). Here, V 1 = −χ(ǫ 2 1 /2 + ǫ 2 2 ), and V 2 = −χǫ 2 2 /2. The condition V 2 /V 1 < 1/2 can be achieved for arbitrary nonzero ǫ 1 and ǫ 2 and blue detuning with χ < 0. Here, we set V 1 = 1.2E R and V 2 = 0.4E R . E R = h 2 2ma 2 is the recoil energy and a is the lattice constant. The contour of V (x, y) is shown in Fig. 1(b). The lowest four band structures from the plane-wave expansion approximation upon the potential V (x, y) in Eq. (1) are shown in Fig. 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 Fig. 1(e), a tight-binding (TB) Hamiltonian can be constructed to described the band structures of the fermionic square lattice, i.e., with Here, being the weak harmonic confining potential to stabilize the optical lattice. p † x/y,i,σ , and d † i,σ 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 Fig. 1 in the following parts of the paper if not special specified. The energy spectra of TB Hamiltonian in Eq. (4) are shown in Figs. 1(e) and 1(f). It can be found that the TB Hamiltonian in Eq. (4) gives a good description of the band structures of lattice potential, and the on-site rotation in the second term in Eq. (6) induces the orbital hybridization to break the degeneracy of p x and p y bands around the Γ and M points shown in Fig. 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 give the two-orbital attractive Hubbard interactions as follows 41 , Here, the first term is the intra-orbital attractive interaction, and the second term is the Hund's coupling with the spin operator S il = 1 2 p † il,α σ αβ p il,β and l = x, y. U and J take the following forms, Here, a s is the s-wave scattering length with negative value, i.e., a s < 0. ω x/y (r) are the Wannier functions of p x/y orbitals. The third term in Eq. (9) is the inter-orbital attractive interaction with n il = n il,↑ + n il,↓ . The fourth term is the pair hopping term. Furthermore, we have J = 2U/3 and J ∆ = 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 Eqs.(9)-(11) are obtained under the harmonic approximation. It is shown that the an-harmonicity of the optical lattice can affect the properties of the multiorbital 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 2U 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 an-harmonic effects in the present work.

III. HOMOGENEOUS SUPERFLUID STATES WITH CHIRAL D-WAVE ORDER
Now, we turn to consider the homogeneous superfluid state with V t = 0 in Eq. (8) and the superfluid state is driven by the attractive interaction in Eq. (9). The spin-singlet superfluid pairing operators are defined aŝ Then, we have 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, ∆ s,ll ′ = ∆ s,ll ′ , H int can be decoupled to be with The homogeneous superfluid state can be described by the mean-field Hamiltonian in the Nambu basis: Here, C is an operator-independent constant term. ∆ is a 3 × 3 matrix and takes the following form, The mean-field Hamiltonian in Eq. (17) can be self-consistently solved with respect to the minimum of ground state energy, i.e., where, E n (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 Fig. 2(a). The typical Fermi surface is shown in Fig. 2(b). From Eq. (18), we can find that the superfluid order parameter in the intra-p x orbital channel is To maximize the superfluid gap, one can find that ∆ s,xx ∆ s,yy > 0 is favorable to obtain the largest amplitudes of ∆ 22 and ∆ 33 . The numerical results for the ground state energy and superfluid order parameters as functions of chemical potential µ 0 and interaction amplitude |U | are shown in Figs. 2(c) and 2(d), from which the intra-orbital ∆ 22 and ∆ 33 are degenerate in the whole parameter regime. It means that ∆ 22 = ∆ 33 , and the only choice is ∆ s,xx ∆ s,yy > 0 thanks to U J ∆ > 0. The aforementioned analyses are consistent, and one can achieve that the superfluid ground states favor ∆ s,xx and ∆ s,yy with same sign to maximize the superfluid gap and to minimize the ground state energy. Furthermore, we can find that the inter-orbital ∆ 23 , which is also the matrix element in Eq. (18), is purely imaginary, and much smaller than ∆ 22/33 . The reason lies in that the inter-orbital ∆ 23 is induced by the orbital hybridization and modulated by Ω z . It is conceivable that the strength of inter-orbital ∆ 23 could be comparable to intra-orbital ∆ 22/33 when Ω z is large enough. However, the ∆ 23 has no relation with the topological nature of the superfluid state, we only focus on the case with Ω z set here.
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 Fig. 1(e), the full separation between the d band and p bands guarantees the feasibility to downfold the Hamiltonian from the space spanned by d and p orbitals to the space spanned by two effectivep orbitals shown in Fig. 1(f). When V t = 0, the translation symmetry allows ones to write the TB Hamiltonian in momentum space under the effective basis ψ σ (k) = [p x,k,σ ,p y,k,σ ] t , i.e.,H Here,H tb (k) = and ξ ± (k) = 2(t pp ∓t ′ pp )(cos k x ± cos k y ), ξ xy (k) = 4t xy sin k x sin k y .
The Pauli matrices σ i with i = x, y, z span the two effectivep x andp y orbital space. The effective TB HamiltonianH tb can be rewritten in the basis spanned by the orbital angular momentum eigen-state, i.e., Here, the parameters are same as these in Fig. 1(f). (c) The zero-temperature ground-state energy of superfluid state as change as chemical potential µ0 and interaction amplitude |U |. (d) The intra-and inter-orbital superfluid order parameters as change as chemical potential µ0 and interaction amplitude |U |. Here, ∆intra = ∆22 with ∆33 = ∆22, and ∆inter = ∆23. The explicit expressions of ∆22, ∆23 and ∆33 are shown in Eq. (18), which are the relevant matrix elements. The mesh of kx × ky = 51 × 51.
The Pauli matrices s i with i = x, y, z span the two effectivep + andp − orbital space. In the absence of Ω z , [F x (k) = 1 2 ξ − (k), F y (k) = ξ xy (k)] forms a vector field in momentum space shown in Fig. 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., with F (k) = F 2 x (k) + F 2 y (k). The direct calculation gives W σ = 2 in agreement with the pattern of the vector field [F x (k), F y (k)] as shown in Fig. 2(b). Note that the total winding number W = W ↑ + W ↓ should be 4 when the spin degree of freedom is taken into account. In the presence of Ω z , the induced orbital hybridization lifts the degeneracy at M point. Then, the above mapping does not work. Here, the y direction has periodic boundary condition while the lattice number along x direction is set to be Nx = 41. (c) The amplitudes of wave-function of the in-gap states labeled 1-8 in (a). Note that each point are double degeneracy by taking into account the spin degree of freedom. Here, the red and blue "o" marks label particle-like |u px ,↑/↓ (ky, ix)| 2 and hole-like |v px,↑/↓ (ky, ix)| 2 while the red and blue "⊳" marks label particle-like |u py ,↑/↓ (ky, ix)| 2 and hole-like |v py ,↑/↓ (ky, ix)| 2 . (d) The phase diagram as change as chemical potential µ0 and interaction amplitude |U |.
In the superfluid states, quasi-particle spectra are fully gapped and the nonzero Ω 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: Here ∆ intra = ∆ 22 and |∆ inter | = |∆ 23 |. Upon an unitary rotation 26 , we can obtain a dual form of the Hamiltonian, i.e., where In the dual HamiltonianH D± mf (k) shown in Eq. (31), [ ξ−(k) 2 , ξ xy (k)] resembles two components of pairing order parameters of the chiral d-wave superfluid and ∆ inter corresponds to the mixed s-wave component. "∆ intra ± hΩ z " is the pseudo-kinetic energy with k-independent, and resembles kinetic energy term " 2m " of the chiral d-wave superfluid when µ 0 is set to satisfy the condition µ 0 = 1 2 ξ + (π, π). Then, the dual HamiltonianH D± mf (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, ∆ inter by itself cannot drive the gap-closing condition, because it is much smaller than ∆ intra and Fermi energy. Therefore, the small ∆ inter can be absorbed and set to zero. The topological nontrivial superfluid states can be achieved under the condition 44 ∆ intra < hΩ z when µ 0 = 1 2 ξ + (π, π), which naturally corresponds to the weak-coupling condition For the general case with arbitrary µ 0 , one can obtain nontrivial superfluid states if f (Ω z , µ 0 , ∆ intra ) > 0 with f (Ω z , µ 0 , ∆ intra ) shown in Eq. (33), and trivial superfluid states if f (Ω z , µ 0 , ∆ intra ) < 0. The topological phase transition condition coincides with the gap-closing condition with f (Ω z , µ 0 , ∆ intra ) = 0. The phase diagram separating the topological trivial and non-trivial superfluid phases is plotted in Fig. 3(d) according to phase transition condition f (Ω z , µ 0 , ∆ intra ) = 0.
The nontrivial topological nature of the superfluid states can be characterized by the Chern number, with u s,n (k) the Bloch functions of occupied quasi-particle states with s = up and down to label the the up-block and downblock parts of Hamiltonian in Eq. (17). The straightforward calculations give C up = C down = 2 for hΩ z > 0 and C up = C down = −2 for hΩ z < 0 under the condition f (Ω z , µ 0 , ∆ intra ) > 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 Fig. 3(a) and no gapless edge states emerge in trivial superfluid state shown in Fig.3(b). The local feature of the edge states in the Fig. 3(a) are explicitly demonstrated through the amplitude distributions of the wave-functions shown in Fig. 3(c).

IV. MASS DENSITY MODULATION FROM THE HARMONIC CONFINING POTENTIAL
Now, we consider the realistic case with nonzero harmonic confining potential in Eq. (8), and the pattern of V trap (i x , i y ) is shown in Fig. 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 tb + H p int in Eqs. (4) and (15) in lattice space. The quasi-particle spectra and the distribution of superfluid order parameters are shown in Figs. 4(b), 4(d), and 4(e) for two different hΩ z = 0.2 and hΩ z = 0.05 under the periodic boundary condition. We find that the amplitudes of superfluid order parameters in both cases are similar from Fig. 4(d) and 4(e), but the quasi-particle spectra are quite different from Fig. 4(b) with in-gap fermion modes for hΩ z = 0.2 and without in-gap fermion modes for hΩ z = 0.05. The reason lies in that V trap (i x , i y ) forms a spatial barrier structure [The position is marked with red-dashed circle in Fig. 4(a)] separating the nontrivial superfluid state with f (Ω z , µ i , ∆ intra ) > 0 and trivial superfluid state with f (Ω z , µ i , ∆ intra ) < 0 for Ω z = 0.2. Note that µ i = µ 0 + V trap (i x , i y ), thus the position of spatial barrier coincides with the gap-closing condition with f (Ω z , µ i , ∆ intra ) = 0 . For fixed µ 0 and V t , one can find that f (Ω z , µ i , ∆ intra ) is always smaller than zero when Ω z = 0.05. The superfluid is always trivial, because Ω z = 0.05 is too small to overcome the gap-closing condition f (Ω z , µ i , ∆ 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 . 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 , where ω a = c a /R is the angular velocity of the rotation along the spatial barrier with R the radius of spatial barrier of V trap (i x , i y ) 48 , a labels the ath branch, and Q a = 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) = −ω a Q a . In the square lattice space, the circular rotation symmetry SO(2) for Eq. (35) is broken down to C 4 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 4 symmetry. Correspondingly, the energy levels of the in-gap fermion modes trapped by the spatial barrier are discrete [see Fig. 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 where u n i,lσ and v n i,lσ 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 Fig. 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 Fig. 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 Fig. 4(c) for details].
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 where we have assumed the thickness along z direction to be unity. The extra 1/2 in denominator is added to compensate the double count 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 cross the zero energy, such as the modulations of V 1 and V 2 in Eq. (1) to deform the ξ−(k) 2 and ξ xy (k) and introducing additional laser beam to modulation the trapping potential. Here, we consider a more convenient method. From Eq. (35), it is straightforward to inject non-zero Q into the superfluid state through slight modulation of local rotating frequency Ω z . As a consequence, the in-gap fermion modes can be driven to cross the zero energy by the non-zero δΩ z . If we further assume that all the in-gap fermion modes trapped in the spatial barrier have the relation ∼ hδΩ z , we can obtain that the response of change of mass current to the modulation of the rotating frequency δI M ∼ mδΩz 2 a sgn(c a ) with the summation involving all the in-gap fermion modes cross zero energy. However, in the square lattices, we can find that the k a is different for different a-th branch from Fig. 5. As a good approximation, we can define an effective k to remove the difference of different k a , and k can be replaced with the averge Fermi momentum k F . Then, we can obtain that the modulation of mass current density is proportional to the change of the LDOS, i.e., with The pattern of δρ(i x , i y ) for δΩ z = 0.02/h is shown in Fig. 5(c), from which we can find that the mass current is trapped around the spatial barrier.

V. DISCUSSIONS AND CONCLUSIONS
In terms of experiment, the fermion atoms can be selected as lithium 6 Li, two internal states can be selected as 2 S 1/2 with M=± 1 2 . The principal fluorescence line from 2 S 1/2 to 2 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 ∼ h × 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 Ω z can be turned with large flexibility 18 . From the energy bands in Fig. 1, we can estimate that it is enough for Ω z ∼ h × 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 p-orbital bands, and increasing the interactions and optical lattice potentials can suppress anisotropy of condensate density 49 . 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 one dimensional limit [50][51][52] . Thus, the boson with infinitely repulsive interactions is roughly equivalent to free fermion even in two dimensional system. Such effective "repulsive interactions"can suppress the anisotropy of condensate density, and guarantee the validity of LDA in 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 onsite 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 [53][54][55] . Besides the radio-frequency spectroscopy, the recently developed matter-wave interference technique 56 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 realize 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 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.