2D gapless topological superfluids generated by pairing phases

We consider a 2D bi-layer system with a square optical lattice within the x–y plane (see figure 1(a)), where a single-species Fermi gas is loaded with finite tunneling amplitudes between lattice sites in the same or opposite layers. As a result, the system could also be effectively understood as a 2D pseudo-spin $1/2$ system with an effective in-plane magnetic field via the inter-layer tunneling. In order to generate an effective interaction between these single-species fermions, an off-resonant two-photon transition is introduced to weakly couple the electronic ground state to a Rydberg excited state via an intermediate state. In the far detuning and weak coupling limit, it is known that the effective Rydberg-dressed interaction between these dressed-state atoms could be approximated by a finite ranged soft-core interaction (see figure 1(b)): $V_\mathrm{RD}(\mathbf{r})\,=\,\frac{U_{0}}{1\,+\,(r/R_{c})^{6}}$ within the standard perturbation and adiabatic approximation [28–33]. Here the effective interaction strength, U₀, and the interaction range, R_c , could be calculated explicitly from the Rabi coupling and detuning from the atomic excitations [28]. Many proposals discussed how to manipulate such a Rydberg-dressed interaction via various external fields [34, 35].

Zoom In Zoom Out Reset image size

We focus on the parameter regime to generate an effective attraction ($U_0\lt0$) and a finite-ranged Rydberg blockade radius ($R_c\,\sim\,a,d$, with a and d being the intra-layer and inter-layer lattice spacing respectively). As a result, these single-species fermions could be effectively interacting with each other with the following Hamiltonian in real space

$$\begin{align} \skew{2.4}\hat{H}& = -\,t\sum_{\mathbf{r},\sigma}\left[\hat{c}^\dagger_{\mathbf{r},\sigma}\hat{c}^{}_{\mathbf{r}+\hat{x},\sigma}\,+\,\hat{c}^\dagger_{\mathbf{r},\sigma}\hat{c}^{}_{\mathbf{r}+\hat{y},\sigma} \,+\,\textit{h.c.}\right]\nonumber\\ & \quad - t_z\sum_{\mathbf{r},\sigma} \hat{c}^\dagger_{\mathbf{r},\sigma}\hat{c}^{}_{\mathbf{r},-\sigma} \,-\,\mu\sum_{\mathbf{r},\sigma} \hat{c}^\dagger_{\mathbf{r},\sigma}\hat{c}^{}_{\mathbf{r},\sigma } \nonumber\\ & \quad - \frac{1}{2} \sum_{\mathbf{r}, \mathbf{r}^{^{\prime}},\sigma} V_{\parallel}(\mathbf{r}\,-\,\mathbf{r}^{^{\prime}}) \hat{c}^{\dagger}_{\mathbf{r},\sigma}\hat{c}^{\dagger}_{\mathbf{r}^{^{\prime}},\sigma}\hat{c}^{}_{r^{^{\prime}},\sigma} \hat{c}^{}_{r,\sigma}\nonumber\\ & \quad - \frac{1}{2} \sum_{\mathbf{r}, \mathbf{r}^{^{\prime}},\sigma} V_{\perp}(\mathbf{r}\,-\,\mathbf{r}^{^{\prime}}) \hat{c}^{\dagger}_{\mathbf{r},\sigma}\hat{c}^{\dagger}_{\mathbf{r}^{^{\prime}},-\sigma}\hat{c}^{}_{\mathbf{r}^{^{\prime}},-\sigma} \skew{2.4}\hat{c}^{}_{\mathbf{r},\sigma} \end{align} \tag{ 1 }$$

where $\hat{c}_{\mathbf{r},\sigma}^{\dagger}$($\hat{c}_{\mathbf{r},\sigma}^{}$) is the creation(annihilation) operator of fermions with the layer index $\sigma\,=\,(\uparrow ,\downarrow)$, µ is the chemical potential, and $\mathbf{r}\,=\,\left(i_{x}, i_{y}\right)$ is the in-plane coordinate. The intra-layer and inter-layer tunneling amplitude are t and t_z respectively. Both tunneling amplitudes can be tuned independently by the lattice strength and lattice spacing. Here $V_{\perp}(\mathbf{r}\,-\,\mathbf{r}^{^{\prime}})$ and $V_{||}(\mathbf{r}\,-\,\mathbf{r}^{^{\prime}})$ are the interlayer and intra-layer interaction provided by the Rydberg-dressed interaction, $V_\mathrm{RD}(\mathbf{r})$, and defined positively for the convenience of latter discussion.

To investigate the topological features, we further apply Fourier transform of the field operator, $\hat{c}_{{\mathbf{r}},\sigma}\,=\,\frac{1}{\Omega}\sum_{\mathbf{k}} \hat{c}_{\mathbf{k},\sigma} \mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}$, to transform the Hamiltonian into the momentum space and apply the BCS mean-field approximation between pairs of fermions. Here Ω is the system area. It is well-known that the resulting mean-field Hamiltonian ($\skew{2.4}\hat{H}_\mathrm{MF}$) can be written to be a $4\,\times\,4$ matrix form within the Nambu spinor representation, $\hat{\Psi}_\mathbf{k}\,\equiv\,\left[\hat{c}_{\mathbf{k} \uparrow}, \hat{c}_{\mathbf{k} \downarrow}, \hat{c}_{-\mathbf{k} \uparrow}^{\dagger}, \hat{c}_{-\mathbf{k} \downarrow}^{\dagger}\right]^\mathrm{T}$, and the final results can be divided in two parts: $\skew{2.4}\hat{H}_\mathrm{MF}\,=\,\skew{2.4}\hat{H}_\mathrm{BCS}\,+\,E_{\mathrm{c}}$, where

$$\begin{align} \skew{2.4}\hat{H}_\mathrm{BCS}\,=\,\frac{1}{2}\sum_{\mathbf{k}} \hat{\Psi}_\mathbf{k}^\dagger \left[\begin{array}{cccc}{\varepsilon_{k}} & {-t_{z}} & {{2\Delta}_{p \uparrow }^{*}} & {\Delta_{s}^{*}} \\ {-t_{z}} & {\varepsilon_{k}} & {-\Delta_{s}^{*}} & {{2\Delta}_{p\downarrow}^{*}} \\ {{2\Delta}_{p \uparrow }} & {-\Delta_{s}} & {-\varepsilon_{k}} & {t_{z}} \\ {\Delta_{s}} & {{2\Delta}_{p \downarrow}} & {t_{z}} & {-\varepsilon_{k}}\end{array}\right]\hat{\Psi}_\mathbf{k} \end{align} \tag{ 2 }$$

and the constant energy term is given by

$$\begin{align} E_\mathrm{c } &\,=\,\frac{1}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{^{\prime}}} V_{\mathbf{k}, \mathbf{k}^{^{\prime}}}^{\perp}\langle \hat{c}_{\mathbf{k}^{^{\prime}}, \uparrow}^{\dagger} \hat{c}_{-\mathbf{k}^{^{\prime}}, \downarrow}^{\dagger}\rangle\left\langle \hat{c}_{-\mathbf{k}, \downarrow} \hat{c}_{\mathbf{k},\uparrow}\right\rangle \nonumber \\ & \quad + \frac{1}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{^{\prime}}} V_{\mathbf{k},\mathbf{k}^{^{\prime}}}^{\|}\langle \hat{c}_{\mathbf{k}^{^{\prime}}, \uparrow}^{\dagger} \hat{c}_{-\mathbf{k}^{^{\prime}}, \uparrow}^{\dagger}\rangle\left\langle \hat{c}_{-\mathbf{k} \uparrow} \hat{c}_{\mathbf{k},\uparrow}\right\rangle \nonumber \\ & \quad + \frac{1}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{^{\prime}}} V_{\mathbf{k}, \mathbf{k}^{^{\prime}}}^{\|}\langle \hat{c}_{\mathbf{k}^{^{\prime}}, \downarrow}^{\dagger} \hat{c}_{-\mathbf{k}^{^{\prime}}, \downarrow}^{\dagger}\rangle\left\langle \hat{c}_{-\mathbf{k}, \downarrow} \hat{c}_{\mathbf{k},\downarrow}\right\rangle \nonumber \\ & \quad + 2\sum_{\mathbf{k}} \varepsilon_{\mathbf{k}} \end{align} \tag{ 3 }$$

with $\langle \cdots\rangle$ being an expectation value and $\varepsilon_{\mathbf{k}}\,=\,-2 t( \cos k_{x}\,+\,\cos k_{y})\,-\,\mu$ being the bare kinetic energy. For simplicity, we use the layer index as the pseudo-spin index and define the inter-layer pairing order parameters to be the 's-wave' pairing ($\Delta_s$) and the intra-layer pairing order to be the 'p-wave' pairing ($\Delta_p$). These order parameters could be calculated through the inter- and intra-layer interaction strength as following:

$$\begin{align} \Delta_{s}(\mathbf{k}) &\,=\,\frac{1}{\Omega}\sum_{\mathbf{k}^{^{\prime}}} V_{\mathbf{k}, \mathbf{k}^{^{\prime}}}^{\perp}\left\langle \hat{c}_{-\mathbf{k}^{^{\prime}},\downarrow} \hat{c}_{\mathbf{k}^{^{\prime}}, \uparrow}\right\rangle \end{align} \tag{ 4 }$$

$$\begin{align} \Delta_{p,\sigma}(\mathbf{k}) &\,=\,\frac{1}{\Omega} \sum_{\mathbf{k}^{^{\prime}}} V_{\mathbf{k},\mathbf{k}^{^{\prime}}}^{\|}\left\langle \hat{c}_{-\mathbf{k}^{^{\prime}}, \sigma} \hat{c}_{\mathbf{k}^{^{\prime}},\sigma}\right\rangle \end{align} \tag{ 5 }$$

for $\sigma\,=\,\pm\,=\,\uparrow/\downarrow$. Here $V_{\mathbf{k},\mathbf{k}^{^{\prime}}}^{\|,\perp}$ is the Fourier transform of the intra-layer and the inter-layer interaction matrix elements $V_{\|,\perp}(\mathbf{r}\,-\,\mathbf{r}^{^{\prime}})$ respectively.

Note that, the inter-layer pairing could be in principle either s-wave (singlet) or p-wave (triplet). But the singlet pairing must be energetically more favorable, because its symmetric orbital wavefunction always lowered the attractive inter-layer interaction energy, similar to the systems of polar molecules with an electric field perpendicular to the layer plane [21, 22, 36]. On the other hand, the intra-layer pairing in the x–y plane must be triplet since the orbital wavefunction in the momentum space must be anti-symmetric as shown in equation (5).

In order to calculate the ground state energy within the mean-field approximation, we first diagonalize the BCS mean-field Hamiltonian (equation (2)), and obtain the low energy excitations for the Bogoliubov quasi-particles:

$$\begin{align} \skew{2.4}\hat{H}_\mathrm{BCS}\,=\,\frac{1}{2}\sum_{\mathbf{k}}\hat{\Gamma}_\mathbf{k}^\dagger \left[\begin{array}{cccc}{E_{1, \mathbf{k}}} & {0} & {0} & {0} \\ {0} & {E_{2, \mathbf{k}}} & {0} & {0} \\ {0} & {0} & {E_{3, \mathbf{k}}} & {0} \\ {0} & {0} & {0} & {E_{4, \mathbf{k}}}\end{array}\right] \hat{\Gamma}_\mathbf{k} \end{align} \tag{ 6 }$$

where $\hat{\Gamma}_\mathbf{k}\,\equiv\,\left[\skew{2.4}\hat{\gamma}_{1,\mathbf{k}}, \skew{2.4}\hat{\gamma}_{2,\mathbf{k}}, \skew{2.4}\hat{\gamma}_{3,-\mathbf{k}}^\dagger, \skew{2.4}\hat{\gamma}_{4,-\mathbf{k}}^\dagger\right]^\mathrm{T}$ is the Nambu spinor in Bogoliubov eigenmode basis with $E_{j,\mathbf{k}}$ ($j\,=\,1,\ldots 4$) being their excitation energies. Their analytic forms may not be available for arbitrary $\Delta_s$ and $\Delta_p$, and therefore will be evaluated numerically in the rest of this paper.

Using the anti-commutation relationship between fermion operator, $\skew{2.4}\hat{\gamma}_{j,\mathbf{k}}$, and $\langle \skew{2.4}\hat{\gamma}_{j,\mathbf{k}}^\dagger \skew{2.4}\hat{\gamma}_{j,\mathbf{k}}\rangle\,=\,0$ for the ground state expectation value, we could easily derive the mean-field ground state energy to be (combined with the constant term, $E_\mathrm{c}$ in equation (3))

$$\begin{align} E_\mathrm{G}&\,=\,E_\mathrm{C}\,+\,\frac{1}{2} \sum_{j,\mathbf{k}} E_{j,\mathbf{k}}\,=\,\frac{1}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{^{\prime}}} V_{\mathbf{k}, \mathbf{k}^{^{\prime}}}^{\perp}\langle \hat{c}_{\mathbf{k}^{^{\prime}}, \uparrow}^{\dagger} \hat{c}_{-\mathbf{k}^{^{\prime}}, \downarrow}^{\dagger}\rangle\left\langle \hat{c}_{-\mathbf{k}, \downarrow} \hat{c}_{\mathbf{k},\uparrow}\right\rangle\nonumber\\ & \quad \,+\,\frac{1}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{^{\prime}}} V_{\mathbf{k},\mathbf{k}^{^{\prime}}}^{\|}\langle \hat{c}_{\mathbf{k}^{^{\prime}}, \uparrow}^{\dagger} \hat{c}_{-\mathbf{k}^{^{\prime}}, \uparrow}^{\dagger}\rangle\left\langle \hat{c}_{-\mathbf{k} \uparrow} \hat{c}_{\mathbf{k},\uparrow}\right\rangle \nonumber\\ & \quad \,+\,\frac{1}{\Omega}\sum_{\mathbf{k},\mathbf{k}^{^{\prime}}} V_{\mathbf{k}, \mathbf{k}^{^{\prime}}}^{\|}\langle \hat{c}_{\mathbf{k}^{^{\prime}}, \downarrow}^{\dagger} \hat{c}_{-\mathbf{k}^{^{\prime}}, \downarrow}^{\dagger}\rangle\left\langle \hat{c}_{-\mathbf{k}, \downarrow} \hat{c}_{\mathbf{k},\downarrow}\right\rangle\nonumber\\ & \quad \,+\,2\sum_{\mathbf{k}}\varepsilon_{\mathbf{k}} \,+\,\frac{1}{2}\sum_{j,\mathbf{k}} E_{j,\mathbf{k}}.\end{align} \tag{ 7 }$$

From the expression above, one could see that the ground state energy depends on the ground state expectation value of the pairing operators, $\langle \hat{c}_{-\mathbf{k},\sigma}\hat{c}_{\mathbf{k},\sigma^{^{\prime}}}\rangle$, which also appears in the definition of the order parameter, $\Delta_{s}(\mathbf{k})$ and $\Delta_{p,\sigma}(\mathbf{k})$, see equations (4) and (5). In this paper, we will use variational method to calculate the ground state energy by parametrizing these order parameters according to the lattice symmetry. More precisely, we could expand the intra- and inter-layer interaction as following:

$$\begin{align} V^{\|,\perp}_{\mathbf{k},\mathbf{k}^{^{\prime}}} \,=\,\sum_{mn\in Z}V_{mn}^{\|,\perp} \cos\left[m(k_x\,-\,k_x^{^{\prime}})a\,+\,n(k_y\,-\,k_y^{^{\prime}})a\right] \end{align} \tag{ 8 }$$

$$\begin{align} \Delta_{s}(\mathbf{k})\,=\,\sum_{m,n\in Z}\Delta^{s}_{mn} \left[\cos\left(mk_xa\,+\,nk_ya\right) \,+\,\cos\left(nk_xa\,-\,mk_ya\right)\right] \end{align} \tag{ 9 }$$

$$\begin{align} \Delta_{p,\sigma}(\mathbf{k})\,=\,\sum_{m,n\in Z}\Delta^{p,\sigma}_{mn}&\left[\sin\left(mk_xa\,+\,nk_ya\right) \right.\nonumber\\ & \;\;\left. \,-\,i\sin\left(nk_xa\,-\,mk_ya\right)\right] \end{align} \tag{ 10 }$$

$$\begin{align} \langle \hat{c}_{-\mathbf{k},\uparrow}\hat{c}_{\mathbf{k},\downarrow}\rangle \,=\,\sum_{m,n\in Z}c^s_{mn}& \left[\cos\left(mk_xa\,+\,nk_ya\right)\right.\nonumber\\ & \;\;\left.\,+\,\cos\left(nk_xa\,-\,mk_ya\right)\right] \end{align} \tag{ 11 }$$

$$\begin{align} \langle \hat{c}_{-\mathbf{k},\sigma}\hat{c}_{\mathbf{k},\sigma}\rangle\,=\,\sum_{m,n\in Z}c^{\,p,\sigma}_{mn}&\left[\sin\left(mk_xa\,+\,nk_ya\right)\right.\nonumber\\ & \;\;\left.\,-\,i\sin\left(nk_xa\,-\,mk_ya\right)\right].\end{align} \tag{ 12 }$$

Here $\Delta^s_{mn}$, $\Delta^{p,\sigma}_{mn}$, $c^s_{mn}$ and $c^{\,p,\sigma}_{mn}$ are all complex numbers, which could be determined later through a variational approach. m and n are integers for the index used in the reciprocal lattice. In the literature before, it is usually assumed that the pairing order parameters are proportional to the interaction strengths for simplicity, but this assumption cannot be justified when considering the competition between order parameters. In order to consider these order parameters correctly, we note that $\Delta^s_{mn}$ and $\Delta^{p,\sigma}_{mn}$ could be calculated from $c^s_{mn}$ and $c^{\,p,\sigma}_{mn}$ through the definition of gap function in equations (4) and (5). We will express their relationship in more details below for variational methods.

Although the expression of order parameters above could be applied to a general finite ranged-ranged interaction, it is still much more intuitive and transparent to start with a finite range interaction range, i.e. $R_c\,\sim\,a, d$, with a and d being the intra-layer and inter-layer lattice spacing respectively. For such a finite ranged attractive interaction, it is reasonable for us to first consider the nearest neighboring terms of pairing only, i.e. $|m|,|n|,|m\,\pm\,n|\,\leqslant\,1$, and neglect the longer-ranged pairing. Since our work is to emphasize the effects of these pairing phases on the ground state properties, it is much more instructive to concentrate on the coupling between the nearest neighboring intra-layer and inter-layer pairing only. Including longer-ranged pairing requires much more pairing phases to be calculated self-consistently, and makes it difficult to interpret the mechanism. We will investigate the pairing phase dependence for this long-ranged interaction in the future. As a result, we could simplify the general expression of the gap functions in the last section and obtain,

$$\begin{align} \Delta_{s}(\mathbf{k})\,=\,\Delta^s_{0}\,+\,\Delta^s_{1}\left[\cos (k_xa)\,+\,\cos (k_ya)\right] \end{align} \tag{ 13 }$$

$$\begin{align} \Delta_{p,\sigma}(\mathbf{k})\,=\,2i\mathrm{e}^{\mathrm{i}\phi_\sigma}\Delta^{p}_{1}\left[\sin\left({k}_xa\right)\,+\,i\mathrm{e}^{\mathrm{i}\alpha_\sigma}\sin\left({k}_ya\right)\right] \end{align} \tag{ 14 }$$

$$\begin{align} \langle \hat{c}_{-\mathbf{k},\uparrow}\hat{c}_{\mathbf{k},\downarrow}\rangle\,=\,c^s_{0}\,+\,c^s_{1}\left[\cos(k_xa)\,+\,\cos(k_ya)\right] \end{align} \tag{ 15 }$$

$$\begin{align} \langle \hat{c}_{-\mathbf{k},\sigma}\hat{c}_{\mathbf{k},\sigma}\rangle \,=\,2ic^{\,p}_{1}\left[\sin\left(k_xa\right)\,+\,i\sin\left(k_ya\right)\right]. \end{align} \tag{ 16 }$$

Here we have set the pairing phase of the s-wave pairing order parameter to be zero, and separate the phases of theses p-wave order parameters, so that $\Delta^s_{0/1}$, $\Delta^p_{1}$, $c_{0/1}^{s}$ and $c_{1}^{p}$ are all defined to be positive. α_σ and φ_σ are the four relative phases for the pairing order parameters as indicated above. Extension to a longer interaction range is straightforward but will not be considered in this paper.

Using the expression above, it is easy to directly connect the relationship between $\Delta^{s,p}_{0,1}$ and $c^{s,p}_{0,1}$ through the definition of order parameters in equations (4) and (5). We have

$$\begin{align} \Delta_{0}^s \,=\,V_{0}^{\perp}c_{0}^{s} \end{align} \tag{ 17 }$$

$$\begin{align} \Delta_{1}^s \,=\,\frac{1}{2}V_{1}^{\perp}c_{1}^{s} \end{align} \tag{ 18 }$$

$$\begin{align} \Delta_{1}^p \,=\,\frac{1}{2}V_{1}^{\|}c_{1}^{\,p} \end{align} \tag{ 19 }$$

and therefore the mean-field ground state energy in equation (7) can be calculated directly to be

$$\begin{align*} E_\mathrm{G}\,=\,\sum_{\mathbf{k}}\left( 2\xi_{\mathbf{k}}\,-\,\frac{1}{2} \sum_{j} E_{j,\mathbf{k}} \right)\,+\,\frac{|\Delta_{0}^{s}|^2}{V_{0}^{\perp}}\,+\,\frac{2|\Delta_{1}^{s}|^2}{V_{1}^{\perp}}\,+\,\frac{2 |\Delta_{1}^{p}|^{2}}{V_{1}^{\|}}. \end{align*}$$

Here the three order parameters, $\Delta^s_{0,1}$ and $\Delta^p_1$, are then treated as independent variational parameters.

However, besides the magnitude of order parameters, the ground state energy also depends on the relative values between these pairing phases, $(\alpha_{\uparrow},\alpha_{\downarrow},\phi_{\uparrow},\phi_{\downarrow})$. They are embedded inside the expression of the Bogoliubov eigenstate energy, $E_{j,\mathbf{k}}$. We have examined all the possible combinations of these pairing phases, but present results only using their relative values, $\alpha\,\equiv\,\alpha_{\uparrow}\,-\,\alpha_{\downarrow}$ and $\phi\,\equiv\,\phi_{\uparrow}\,-\,\phi_{\downarrow}$ by setting $\alpha_{\downarrow}\,=\,\phi_{\downarrow}\,=\,0$ in the rest of this paper. It is because this could give the most representative results without missing other information. As an example, we could analytically calculate the eigenvalues of Bogoliubov excitations for $(\alpha,\phi)\,=\,(0,0)$ to be

$$\begin{align}\begin{array}{*{20}{c}} {{E_{j,{\boldsymbol{k}}}}}&{ = \pm \left[ {\varepsilon _{\boldsymbol{k}}^2 + t_z^2 + 4{{\left| {{\Delta }_1^p} \right|}^2}\left( {{{\sin }^2}\left( {{k_x}a} \right) + {{\sin }^2}\left( {{k_y}a} \right)} \right) + {{\left| {{\Delta}_0^s} \right|}^2}} \right.}\\ {}&{{{\left. { \pm 2\sqrt {\varepsilon _{\boldsymbol{k}}^2t_z^2 + t_z^2{{\left| {{\Delta}_0^s} \right|}^2} + 4{{\left| {{\Delta }_1^p} \right|}^2}{{\left| {{\Delta}_0^s} \right|}^2}{{\sin }^2}\left( {{k_x}a} \right)} } \right]}^{1/2}}.} \end{array} \end{align} \tag{ 20 }$$

For $(\alpha,\phi)\,=\,(\pi,\pi)$, we have

$$\begin{align}\begin{array}{*{20}{c}} {{E_{j,{\boldsymbol{k}}}}}&{ = \pm \left[ {\varepsilon _{\boldsymbol{k}}^2 + t_z^2 + 4{{\left| {{\Delta}_1^p} \right|}^2}\left( {{{\sin }^2}\left( {{k_x}a} \right) + {{\sin }^2}\left( {{k_y}a} \right)} \right) + {{\left| {{\Delta}_0^s} \right|}^2}} \right.}\\ {}&{{{\left. { \pm 2\left| {{t_z}} \right|\sqrt {\varepsilon _{\boldsymbol{k}}^2 + {{\left| {{\Delta}_0^s} \right|}^2} + 4{{\left| {{\Delta}_1^p} \right|}^2}{{\sin }^2}\left( {{k_x}a} \right)} } \right]}^{1/2}}.} \end{array} \end{align} \tag{ 21 }$$

Analytical forms of these eigenstate energies could be also obtained in other specific values of $(\alpha,\phi)$. However, for a general value of $(\alpha,\phi)$, we will evaluate it numerically for the variational calculation below. In figures 2(a1)–(a3) and (b1)–(b3), we show how the numerically calculated eigenenergies changes as a function of these pairing phases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the Anderson pseudo-spin point of view [37], the ${SU}(2)$ and Z₂ symmetries are broken explicitly in the BCS mean-field Hamiltonian (see equation (2)). Nevertheless, we find that, by selecting appropriate parameters, say $t_{z}\,=\,0$, $\phi_{\uparrow}\,=\,-\phi_{\downarrow}$, and $\alpha_{\uparrow}\,=\,\pi+\alpha_{\downarrow}$, the bilayer system we study here becomes time-reversal symmetric. More precisely speaking, the chiral directions of the two superfluids in the upper and lower layers are opposite, making the time-reversal symmetry restored if reversing the parity in the z axis. However, if the inter-layer tunneling is finite, such time-reversal symmetry disappears as expected, both in the original Hamiltonian and the mean-field one. In terms of the mean-field Hamiltonian, the pairing term in low energy limit can be expressed as $\Delta(k)\,=\,i \sigma_2 (\Delta_s \sigma_0\,+\,\Delta_p (k_x \sigma_1\,+\,i \mathrm{e}^{\mathrm{i}\alpha} k_y \sigma_2 ))$. The value of α determines the time-reversal symmetry of the system. We find that in the self-consistent mean-field treatment, α can be π. The time-reversal operation $T:\,=\,i \sigma_2 K$ maps $\Delta(k)$ to $\Delta(-k)$ indicates the preserving of the time-reversal symmetry. In our self-consistent mean-field treatment, we do not add other constraints besides the breaking of U(1) symmetry. However, since the relative values of these phases may be pinned by the interaction effects, some symmetry (say the C₄ rotational symmetry of original Hamiltonian in the x–y plane) may be also broken in the exotic phases we discussed in the context of figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 3(a), we show the calculated quantum phase diagram in terms of the two pairing phases, α and φ. The result is obtained by a given set of system parameters, $V_0^\perp/t\,=\,11$, $V_1^\perp/t\,=\,5.8$, $V_1^\|/t\,=\,6.2$, $t_z/t\,=\,0.8$, and $\mu\,=\,0.2$. The ground state is obtained by minimizing the ground state energy in the parameter space, $(\Delta_0^s,\Delta_1^s,\Delta_1^p)$. One can see that when the relative pairing phase between the two layers (i.e. α) is small, the ground state is mainly s-wave for $\phi\,\sim\,\pi/2$ and becomes p-wave otherwise. The phase boundary between these two phases depend on α weakly. The transition between these two phases are found to be first order, as shown in figure 3(b). Such a result is reasonable because, when only the onsite inter-layer s-wave pairing is present (i.e. $\Delta_0^s\,\neq\,0$, $\Delta_1^s\,=\,\Delta_1^p\,=\,0$), all the pairing phase dependence of the ground state energy disappears. On the other hand, when the p-wave order parameters dominant, their relative phase respect to the s-wave also becomes unimportant, while the ground state energy strongly depends on the relative phase, φ, between the upper and lower layers. The lowest energy for small α regime stays at $\phi\,=\,0,\pi$ as expected.

Note that, from energetic point of view, one could easily expect that the s-wave pairing will dominate as the two layers are very close (i.e. $d\,\ll \,a$), while the p-wave become dominant as the two layers are far away from each other (i.e. $d\,\gg\,a$). The phase diagram here is obtained by changing the pairing phases for a fixed value of $d/a$ (here given by $t_z/t\,=\,0.8$). Such pairing phase dependence in ground state energy could lead to a certain spontaneous phase transition, since these pairing phases are automatically determined by the interaction effects.

Besides the two traditional pairing mechanism shown above, we find that a hybrid ground state, where both $\Delta_0^s$ and $\Delta_1^p$ are finite (but within the parameter regime, $\Delta_1^s\,=\,0$). It appears in a regime near $(\alpha, \phi)\,=\,(\pi,\pi)$, which was not investigated before. From the self-consistent mean-field calculation demonstrated here, the phase transition between this sp-coexisting ground sate appears through a second order phase transition with a lower energy than either s-wave or p-wave states alone. The presence of such a sp-coexisting state is from the fact that near $(\alpha, \phi)\,=\,(\pi,\pi)$, $\Delta_0^s$ and $\Delta_1^p$ are coupled in a higher order terms in Landau's free energy expression through the complicated Bogoliubov excitation spectrum.

In order to investigate the phase diagram in terms of the realistic system parameters, we further express the quantum phase diagram in terms of interaction strengths, $V_0^\perp$ and $V_1^\|$. The ground states are then determined by minimizing the full variational energy ($E_\mathrm{G}$) in terms of three variational order parameter, ($\Delta_0^s$, $\Delta_1^s$, $\Delta_1^p$), and two pairing phases, $(\alpha,\phi)$. Note that, both $V_0^\perp$ and $V_1^\|$ can be directly calculated from the full effective Rydberg-dressed interaction between dressed states. More generally, we analytically express their values in terms of a general inter-layer and intra-layer interaction, i.e. $V_0^\perp\,=\,V_{0,0}^\perp$, $V_1^\perp\,=\,V_{\pm 1,0}^\perp\,=\,V_{0,\pm 1}^\perp$, and $V_{1}^\|\,=\,V_{\pm 1,0}^\|\,=\,V_{0,\pm 1}^\|$, where

$$\begin{align} V_{m,n}^{\perp}\,\equiv\,\frac{U_0}{\left[d^{2}\,+\,(n a)^{2}\,+\,(ma)^{2}\right]^{3}\,+\,R_{c}^{6}} \end{align} \tag{ 22 }$$

$$\begin{align} V_{m,n}^{\|} \,\equiv\,\frac{U_0}{[(na)^2\,+\,(ma)^2]^3\,+\,R_{c}^6} \end{align} \tag{ 23 }$$

with $n,m\in Z$ being the lattice site index. As a result, both $V_{1}^{\|}$ and $V_{0}^{\perp}$ can be tuned independently by changing U₀ and d, keeping the lattice constant a and the interaction range R_c the same. The value of inter-layer and inter-site interaction, $V_1^\perp$ is then changed accordingly.

In figure 4, we show the calculated quantum phase diagram in terms of $V_{1}^{\|}$ and $V_{0}^{\perp}$ for $t_z/t\,=\,0.8$ and $\mu/t\,=\,0.2$. Several interesting and important properties could be observed: first, in the limit of zero interaction strength (white regime), there is no superfluid order parameter, because the finite inter-layer tunneling plays as a magnetic field to open a gap between different pseudo-spin components, suppressing the formation of Cooper pairs. In principle, the ground state could be an FFLO state due to the Fermi wavevector mis-matching, while we did not include this order parameters in our mean-field approach to simplify the calculation. This simplification is justified here because the critical temperature of the FFLO state is known much lower than regular superfluids due to its non-zero condensate momentum. (For example, the estimated $T_\mathrm{c}$ of such an FFLO state of Na–Li polar molecules in a bilayer system is just about 10 nK [36], well below the temperature current experiments.) Therefore, we believe our calculated quantum phase diagram should not be affected much even considering the FFLO state at a finite temperature.

Secondly, in the inter-mediate and stronger interaction regime, we observe the s- (p-wave) superfluid in the regime when $V_0^\perp$($V_1^\|)$ becomes dominant in the blue (pink) regime. This reflects the fact that this bilayer system has a great flexibility to investigate the quantum phase transition between these two different superfluids via changing Rydberg-coupling strength or the inter-layer distance directly. As described above, the phase boundary between these two phases are first order without inter-mediate phase if the interaction is not strong.

However, in a stronger interaction regime, we do find the possibility to have a co-existing s-wave and p-wave superfluid, where the pairing phase becomes $(\alpha, \phi)\,=\,(\pi,\pi)$, as shown in the upper right corner of figure 3(a). We have to emphasize that, if one fixes the relative phases between s- and p-wave order parameters to be zero ($\alpha\,=\,\phi\,=\,0$) as one usually did in the literature, we will not be able to find such a co-existing multi-order superfluidity from the variational approach. Although the regime of such an sp-coexisting phase seems not very large, it still provides an important route to investigate a new topological superfluidity, i.e. the traditional chiral p-wave superfluids could be coupled and paired together, making a possible new topological superfluid, which we will study in more details later.

Since we are more interested in the phase with the co-existence of both s- and p-wave superfluids, it will be more instructive to investigate its band structure and edge state properties first. In figure 5(a) we show the calculated band structure for $L_x\,=\,L_y\,=\,100$ with a periodic boundary condition in both x and y directions. We have set $(\alpha, \phi)\,=\,(\pi,\pi)$ and find several gapless points in the band structure. Note that such gapless structure does not exist if the relative phases between these order parameters are the same (i.e. $\alpha\,=\,\phi\,=\,0$). In other words, the pairing phase modulation makes two gapped superfluids (s- and p-waves) co-exist and close their gaps, leading to interesting topological superfluid similar to nodal superconductors [38–42] as shown below.

To investigate if there could be any edge states due to the topological properties in the bulk, in figures 5(b) and (c), we show the calculated band structures with an open boundary condition in the y and x directions (i.e. k_x and k_y are good quantum numbers respectively). We could find that their band structures are very different due to the phase locking between the inter-layer pairing and the intra-layer pairing order parameters. More importantly, we could clearly find the localized state band energy in both sides. This reflects the bulk-edge corresponding principle of the topological matter. When the boundary is open in the y direction, the localized edge mode has a flat band, connecting different gapless points projected on the edge Brillouin zone, as shown in figure 5(b). As we will show below, this indicates the quantized topological charge for these gapless points.

In order to investigate the topological properties of this new phase in more details, here we calculate the possible topological charge of these gapless points, whose positions, $(k^0_x,k^0_y)$, could be calculated easily from the eigenmode spectrum, see equation (21). The gap is closed at $k^0_{y}\,=\,0,\pi$, and $\sqrt{\varepsilon_{\mathbf{k}^0}^{2}\,+\,\left|\Delta_{0}^{s}\right|^{2}\,+\,4\left|\Delta_{1}^{p}\right|^{2} \sin ^{2} (k_{x}^0a)}\,\pm\,t_{z}\,=\,0$, leading to four solutions distributed in the momentum space as shown in figure 5(d). The red/blue colors indicate positive/negative topological charges as we will show below.

It is well-known [43–45] that in three dimensions, the topological charge is defined to be the surface integral of Berry curvature around a topological monopole, i.e. the Chern number. In 2D, a loop that encircles the gapless point defines the topological charge, which is then equivalent to the winding number. There are several methods to calculate the topological charges around the gapless points. Here we note that the band disperses linearly in the momentum space through the gapless point, so that we could obtain the effective Hamiltonian around $\left(k_{x}^{0}, k_{y}^{0}\right)$ up to the linear terms, i.e. $\mathbf{k}\,=\,\mathbf{k}^{0}\,+\,\mathbf{q}\,+\,{\cal O}(q^2)$, and hence $\skew{2.4}\hat{H} (\mathbf{k} )\,=\,\skew{2.4}\hat{H}_{0}(\mathbf{k}^{0})\,+\, \skew{2.4}\hat{H}_{1}(\mathbf{q} )\,+\,{\cal O}(\mathbf{q}^2)$. Here zeroth order term of the Hamiltonian is

$$\begin{align} H_{0}\left( \mathbf{k} ^{0}\right)&\,=\,\varepsilon_{k }\left(\sigma_{3} \otimes \sigma_{0}\right)\,-\,\Delta_{s}\left(\sigma_{2} \otimes \sigma_{2}\right)\,-\,t_{z}\left(\sigma_{3} \otimes \sigma_{1}\right) \nonumber \\ &\quad\,-2\,\Delta_{p}\left[\sin \left(k_{x}^{0}a\right)\right]\left(\sigma_{2} \otimes \sigma_{0}\right)\nonumber \\ &\quad\,-\,2\Delta_{p}\left[\sin \left(k_{y}^{0}a\right)\right]\left(\sigma_{1} \otimes \sigma_{3}\right) \end{align} \tag{ 24 }$$

and the leading order Hamiltonian is

$$\begin{align} \skew{2.4}\hat{H}_{1}(\mathbf{q})&\,=\,2t [q_{x} \sin \left(k_{x}^{0}a\right)+q_{y} \sin \left(k_{y}^{0}a\right)]\left(\sigma_{3} \otimes \sigma_{0}\right)\nonumber\\ &\quad-2\,\Delta_{1}^{p} q_{x} \cos \left(k_{x}^{0}a\right)\left(\sigma_{2} \otimes \sigma_{0}\right). \end{align} \tag{ 25 }$$

Note that we use the direct product of two Pauli matrix for a general expression of the original $4\,\times\,4$ matrix.

In order to get the effective two-band Hamiltonian, we first numerically solve the $4\,\times\,4$ unperturbed Hamiltonian, $\skew{2.4}\hat{H}_{0}(\mathbf{k}^{0})$, and express the leading order term, $\skew{2.4}\hat{H}_{1}(q)$, in this eigenstate basis. We then obtain the $2\,\times\,2$ effective Hamiltonian by projecting the original Hamiltonian to two bands near zero energy

$$\begin{align} H _\mathrm{eff }(\mathbf{q} )&\,=\,\left[\begin{array}{cc} \left\langle v_{1}\left|H_{1}(\mathbf{q})\right| v_{1}\right\rangle & \left\langle v_{1}\left|H_{1}(\mathbf{q})\right| v_{2}\right\rangle \nonumber \\ \left\langle v_{2}\left|H_{1}(\mathbf{q})\right| v_{1}\right\rangle & \left\langle v_{2}\left|H_{1}(\mathbf{q})\right| v_{2}\right\rangle \end{array}\right]\nonumber\\ &\,=\,v_{x}\left( \mathbf{k} ^{0}\right) q_{x} \sigma_{x}^{^{\prime}}\,+\,v_{y}\left( \mathbf{k} ^{0}\right) q_{y} \sigma_{y}^{^{\prime}}. \end{align} \tag{ 26 }$$

In the last line, we have transformed the $2\,\times\,2$ matrix into a new basis, expressed by $\{\sigma_{x}^{^{\prime}}$ and $\sigma_{y}^{^{\prime}}\}$ in order to match the general form of gapless points with coefficients $v_{x}\left( \mathbf{k} ^{0}\right)$ and $v_{y}\left( \mathbf{k} ^{0}\right)$ respectively. In order to calculate the winding number around these gapless points more conveniently, we expresses the final effective Hamiltonian in polar coordinates by $v_{x}\left( \mathbf{k} ^{0}\right) q_{x}\,=\,R(\mathbf{q}) \sin \theta_\mathbf{q} \cos \phi_\mathbf{q}$, $v_{y}\left( \mathbf{k} ^{0}\right) q_{y}\,=\,R(\mathbf{q}) \sin \theta_\mathbf{q} \sin \phi_\mathbf{q}$, so that

$$\begin{align} H_\mathrm{eff }( \mathbf{q} )\,=\,R(\mathbf{q} )\left(\begin{array}{cc} 0 & \mathrm{e}^{\mathrm{i} \theta_{\mathbf{q} }} \\ \mathrm{e}^{-\mathrm{i} \theta_{\mathbf{q} }} & 0 \end{array}\right) \end{align} \tag{ 27 }$$

with $\scriptstyle{R(\mathbf{q} )\,\equiv\,\sqrt{v_{x}\left( \mathbf{k} ^{0}\right)^{2} q_{x}^{2}\,+\,v_{y}\left( \mathbf{k} ^{0}\right)^{2} q_{y}^{2}}}$ and $\scriptstyle{\theta_{ \mathbf{q} }\,\equiv\ \tan ^{-1}\left[\frac{v_{y}\left( \mathbf{k} ^{0}\right) q_{y}}{v_{x}\left( \mathbf{k} ^{0}\right) q_{x}}\right]}$. The corresponding winding number is then given by

$$\begin{align} W(\mathbf{k}^0)\,=\,\frac{1}{2\pi i} \oint \mathrm{d}\theta_{\tilde{\mathbf{q}}} \partial_{\theta_{\tilde{\mathbf{q}}}} \mathrm{e}^{\mathrm{i} \theta_{\tilde{\mathbf{q}}}}\,=\,\pm 1, \end{align} \tag{ 28 }$$

depending on the positions of Weyl points. Note that due to the fermion-doubling theorem [46–49], the sum of all topological charges needs to be zero since topological point nodes with positive and negative values come in pairs.

In figure 5(d), we show the calculated results of topological charges for each gapless point in the momentum space. The bulk energy spectrum is shown in figure 5(a). It is clear that the points of positive charge and negative charge appears alternatively in space, while the breakdown of rotational C₄ symmetry makes the projection of total charge to be zero on the k_y axis, but non-zero on the k_x axis (see figures 5(b) and (c)). It explains why a zero mode appears in the edge state when the boundary is open in the y direction (hence k_x is a conserved quantum number), as shown in figure 5(b).