Topological superfluid phases of attractive Fermi-Hubbard model in narrow-band cold-atom optical lattices

Topologically non-trivial superconductivity [1, 2] or superfluidity has been predicted in two-dimensional (2D) spin-half fermionic systems with Rashba spin-orbit coupling (SOC), Zeeman field, and on-site attractive interaction [3, 4]. The on-site attractive interaction aims to produce s-wave superconductivity or superfluidity in the Rashba system, which becomes topological when the applied Zeeman field exceeds a critical value. Such a spin-orbit coupled Fermi-Hubbard model with an applied Zeeman field can naturally occur in non-centrosymmetric superconductors in a magnetic field [5, 6] and can potentially be realized in cold atom systems in optical lattices [3, 4, 7, 8] with artificially engineered SOC [9, 10, 11, 12] and Zeeman field. Such a setup would leverage the significant experimental progress in realizing artificial SOC with ultracold fermions [13, 14, 15, 16, 17, 18, 19, 20, 21, 22].

A self-consistent mean field calculation for the pairing potential in the spin-orbit coupled Fermi-Hubbard model under a Zeeman field, with fermions following a quadratic dispersion, has shown that although a topological superconducting phase can emerge when the Zeeman field exceeds a critical value, the magnitude of the self-consistent pair potential in the TS phase is extremely small [23]. This makes it challenging to experimentally realize the TS phase in non-centrosymmetric superconductors where the fermion density is low enough to follow a quadratic dispersion. An alternative approach for experimentally achieving a topological superfluid state in spin-orbit coupled systems is to use two-component cold fermionic atoms in optical lattices with artificially engineered SOC and Zeeman field [4, 7, 8, 10, 11]. In this setup, an attractive Hubbard-type onsite interaction can be induced by an $s$-wave Feshbach resonance. The fermion density can be sufficiently high so that the (relevant) fermion dispersion relation corresponds to the lattice dispersion, which can be narrowed by adjusting the hopping parameter and the spin-orbit coupling, nearly reaching the flat-band limit [24, 25, 26, 27].

In flat bands, the critical temperature for Cooper pair formation is predicted to be linearly proportional to the attractive interaction between the Cooper pair constituents [28, 29]. This contrasts sharply with the Bardeen–Cooper–Schrieffer (BCS) theory of superconductivity, valid for quadratic fermion dispersion, where $T_{c}$ is proportional to the exponential of the inverse of the interaction strength. This indicates that the superconducting critical temperature, and thus the pair potential, is exponentially enhanced in flat-band systems, as compared to dispersive systems within the BCS formalism. The enhancement arises from the high density of states (DOS) and the dominance of interactions over kinetic energy. The band does not need to be perfectly flat to benefit from this; any band where the Hubbard interaction strength is larger than the bandwidth will suffice. In this spirit, we explore in this paper the possibility of realizing topological superfluid states in a narrow-band spin-orbit coupled Fermi-Hubbard system with a Zeeman field, where the relevant system parameters have values smaller than (but comparable to) the interaction strength, so that the pairing potential is enhanced by the narrowness of the band making the realization of robust TS phases more feasible experimentally.

We investigate the effects of attractive Hubbard interactions on two-component fermionic atoms in narrow two-dimensional energy bands with Rashba spin-orbit coupling (SOC) and applied Zeeman field. Using self-consistent mean field theory for evaluating the pairing potential, we uncover a complex phase diagram featuring various topological superfluid phases that depend on the chemical potential and the Zeeman field. We examine whether the self-consistent pairing potential characterizing the TS phases, which is notably small for a simple quadratic dispersion near the $\Gamma$-point, as found in earlier work [23], increases as the bandwidth narrows in the optical lattice, as predicted for systems near the flat-band limit [28, 29]. In addition, we determine the (topological) quasiparticle gap that protects different TS phases and identify the optimal system parameters that maximize this gap, ensuring the realization of a robust topological superfluid. We find that, in general, reducing the fermionic bandwidth enhances the pairing potential, consistent with the expected behavior for narrow- or flat-band systems. However, the gap characterizing the TS phases does not exceed maximum values evaluated within our mean-field approximation at about $12.5\%$ of the interaction strength. As discussed latter in the work, these maximum gap values are realized within specific parameter regimes characterized by finite values of the hopping amplitude, Rashba coupling, and Zeeman field. A topological gap $10-12\%$ of the attractive interaction strength $V_{0}$, which can be controlled by a Feshbach resonance, is a significant improvement over the topological gaps potentially achievable in naturally occurring systems, such as noncentrosymmetric superconductors [23]. In recent work [30], the phase diagram of this model was discussed as a function of Rashba coupling and Zeeman field but only near the half-filled limit. Thus, the TS phase found in this work is the Chern number $C=2$ phase near chemical potential equal to $4t$ (which corresponds to half-filling in our model) in Fig. 2 below. By contrast, in the present work, we discuss a variety of TS phases and discuss the phase diagram as a function of the control parameters chemical potential and Zeeman field.

The manuscript is organized as follows. In Sec. II we define the attractive Hubbard model and briefly describe our mean-field approach. In Sec. III we present numerical solutions to the mean field equations, discussing the topological phase diagram and the dependence of the pairing potential and quasiparticle gap on relevant parameters and identifying the optimal regimes that maximize the topological gaps characterizing different TS phases. We summarize in Sec. IV.

We consider a two-dimensional (2D) interacting system with Rashba-type spin-orbit coupling (SOC) and perpendicular Zeeman field and we describe it using a tight-binding model defined on a square lattice. In the basis corresponding to the eigenstates $|\phi_{{\bm{k}}\sigma}\rangle$ of the (non-interacting) system with no SOC and no Zeeman field, the corresponding (second quantized) Hamiltonian has the form $H=H_{0}+H_{int}$, with the non-interacting component having the form

$$H_{0}=\sum_{{\bm{k}},\sigma}(\xi_{{\bm{k}}}+\sigma\Gamma)\leavevmode\nobreak\ \!c_{{\bm{k}}\sigma}^{\dagger}c_{{\bm{k}}\sigma}+\sum_{\bm{k}}\left(\alpha_{\bm{k}}\leavevmode\nobreak\ \!c_{{\bm{k}}\uparrow}^{\dagger}c_{{\bm{k}}\downarrow}+h.c.\right),$$

(1)

where the operator $c_{{\bm{k}}\sigma}^{\dagger}$ creates a particle in the single-particle state $|\phi_{{\bm{k}}\sigma}\rangle$ with momentum ${\bm{k}}=(k_{x},k_{y})$ and spin projection $\sigma=\pm 1\equiv\uparrow\downarrow$ along the (perpendicular) $z$-axis, $\xi_{\bm{k}}=2t\leavevmode\nobreak\ \!(2-\cos k_{x}-\cos k_{y})-\mu$ is the energy spectrum relative to the chemical potential $\mu$, assuming nearest-neighbor hopping with amplitude $t$, $\Gamma$ is the Zeeman field, and $\alpha_{\bm{k}}=\alpha\leavevmode\nobreak\ \!(\sin k_{y}+i\sin k_{x})$ is the SOC contribution. For convenience, we have chosen the lattice constant $a$ as the unit for length, i.e., we have $a=1$. Considering purely local attractive interactions, the second term of the Hamiltonian becomes

$$H_{int}=-V_{0}\sum_{{\bm{k}},{\bm{k}}^{\prime},{\bm{q}}}c_{{\bm{k}}+{\bm{q}}\uparrow}^{\dagger}c_{{\bm{k}}\uparrow}\leavevmode\nobreak\ \!c_{{\bm{k}}^{\prime}-{\bm{q}}\downarrow}^{\dagger}c_{{\bm{k}}^{\prime}\downarrow}$$

(2)

with $V_{0}>0$ being the interaction strength. Note that, using a real space basis, $|\phi_{i\sigma}\rangle$, with $i$ labeling lattice sites, the Hamiltonian in Eq. (2) takes the familiar Hubbard form, $H_{int}=-V_{0}\sum_{i}n_{i\uparrow}n_{i\downarrow}$, with $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$.

In this study, we are interested in the narrow band regime characterized by system parameters ($t$, $\alpha$, and $V_{0}$) and control parameters ($\Gamma$ and $\mu$) having comparable values and, for convenience, we chose the interaction strength, $V_{0}$ as the energy unit. We focus on investigating the possible emergence of topological superfluid phases and finding the optimal regime(s) associated with these phases. To obtain a clear general understanding of the relevant phases and to efficiently explore the rather large parameter space, we use a mean-field approach. Note, however, that in the narrow band regime correlation effects could be significant, hence testing the validity of our results beyond mean-field remains an important future task.

Within a mean-field approach, the interaction term can be approximated by a sum of pairing contributions. The corresponding effective Hamiltonian is $H_{eff}=H_{0}+H_{\Delta}$, with the pairing Hamiltonian having the form

$$H_{\Delta}=\sum_{\bm{k}}\left(\Delta\leavevmode\nobreak\ \!c_{{\bm{k}}\uparrow}^{\dagger}c_{{-\bm{k}}\downarrow}^{\dagger}+\Delta^{*}\leavevmode\nobreak\ \!c_{{-\bm{k}}\downarrow}c_{{\bm{k}}\uparrow}\right),$$

(3)

with the pairing potential $\Delta$ being determined self-consistently by solving a mean-field self-consistent gap equation. A derivation of the gap equation can be found, for example, in Ref. 23. At zero temperature, the pairing $\Delta$ is a solution of the gap equation $\theta(\Delta)=0$, with

	$\displaystyle\theta(\Delta)$	$\displaystyle=$	$\displaystyle-1-\frac{V_{0}}{4N_{k}}\sum_{\bm{k}}\left[\frac{1}{E_{1}}+\frac{1}{E_{2}}\right.$
		$\displaystyle-$	$\displaystyle\left.\frac{\Gamma^{2}}{\sqrt{\xi_{\bm{k}}^{2}|\alpha_{\bm{k}}|^{2}+\Gamma^{2}(\xi_{\bm{k}}^{2}+|\Delta|^{2})}}\left(\frac{1}{E_{1}}-\frac{1}{E_{2}}\right)\right],$

where the summation is done over the 2D Brillouin zone containing $N_{k}$ momentum values. The energies $E_{1}\geq 0$ and $E_{2}>0$, which depend on the system parameters and the pairing $\Delta$, are the (positive) eigenvalues of the effective Hamiltonian, $H_{eff}$. Explicitly, we have

$$E_{1(2)}^{2}=\xi_{\bm{k}}^{2}+|\alpha_{\bm{k}}|^{2}+\Gamma^{2}+|\Delta|^{2}\mp 2\sqrt{\xi_{\bm{k}}^{2}|\alpha_{\bm{k}}|^{2}\!+\!\Gamma^{2}(\xi_{\bm{k}}^{2}\!+\!|\Delta|^{2})}.$$

(5)

Note that the solutions of the pairing equation $\theta(\Delta)=0$ corresponding to $d\theta/d\Delta>0$ are unphysical [23]. For a given set of system and control parameters, ($t$, $\alpha$, $V_{0}$, $\Gamma$, and $\mu$), we determine the pairing $\Delta$ by solving numerically the equation $\theta(\Delta)=0$ and we characterize the low-energy physics of the system (at the mean-field level) by solving the eigenvalue problem associated with the effective Hamiltonian $H_{eff}$. In particular, we calculate the quasiparticle gap $\Delta_{qp}({\bm{k}})\leq\Delta$ characterising the system in different regimes. Note that $\Delta_{qp}$ vanishes at a topological quantum phase transition (TQPT), although the pairing $\Delta$ is finite, or in the absence of pairing, $\Delta=0$.

The main objective of this study is to determine to what degree narrowing the bandwidth, which enhances the density of states (DOS) and, therefore, is expected to strengthen the pairing potential, results in an enhancement of the topological gap. Two examples of narrow-band non-interacting energy spectra and the corresponding DOS are shown in Fig. 1. Note that the density of states is measured in units corresponding to the DOS at the bottom of the band for a square lattice model with nearest-neighbor hopping $t=V_{0}$, where $V_{0}$ can be viewed as an unspecified energy scale to be determined (in an interacting system) by the interaction strength. In the wide band limit ($t\gg V_{0}$), e.g., for a wide-band model with quadratic dispersion near the $\Gamma$-point, the corresponding DOS is $1/2\pi t\ll 1/2\pi V_{0}$. The DOS shown in Fig. 1 has values that are orders of magnitude higher that those corresponding to the large bandwidth regime. Below we investigate the impact of this large DOS on the stability of topological superfluid phases that emerge in the presence of onsite attractive interaction.

In this section we present the results of our numerical analysis, starting with a general discussion of the topological phase diagram for an effective Hamiltonian $H_{eff}$ with constant pairing (Sec. III.1). We continue in Sec. III.2 with self-consistent calculations of the topological phase diagram corresponding to two different parameter regimes. Finally, in Sec. III.3 we identify the optimal parameter regimes that maximize the topological gap, i.e., the minimum over the Brillouin zone of $\Delta_{qp}({\bm{k}})$.

We study a model of fermions with attractive Hubbard interaction in a square 2D optical lattice in the presence of SOC and a time-reversal symmetry breaking Zeeman splitting. We are interested in the large Zeeman energy regime, known to drive band flattening [24, 25, 26, 27] while simultaneously enabling topological superfluid phases. It has been suggested before that the superconducting/superfluid pairing potential $\Delta$ increases with decreasing bandwidth and in the limit of vanishing bandwidth (flat-band limit) $\Delta$ becomes linearly dependent on the attractive interaction $V_{0}$. This is in contrast to the case of the wide-band regime with a quadratic dispersion (near the $\Gamma$-point), where the BCS theory applies. According to the BCS theory, the superconducting pair potential $\Delta$ is proportional to the exponential of the inverse of $V_{0}$, and thus the pair potential should increase exponentially as the bandwidth narrows. This enhancement of the pair potential with narrowing of the bandwidth is ultimately connected to the enhancement of the density of states at the Fermi surface [28, 29]. In this paper, we explore whether the enhancement of the narrow/flat band density of states could enhance the pairing potential in the TS phase(s) and strengthen the topological superfluid phase by enhancing the corresponding topological gap.

Using the zeros of the energy eigenvalues (at ${\bm{k}}=0$), as given by Eq. (5), and the Chern number in Eq. (6), we first show that in the narrow band regime a model with constant pairing potential supports three distinct TS phases that occur in a parameter region significantly expanded as compared to the topological region associated with a wide-band system with quadratic dispersion [23]. Since the former applies to optical lattices and the latter applies to naturally-occurring noncentrosymmetric superconductors [23], we predict that the TS phases should be more easily accessible in optical lattice systems. However, greater accessibility, i.e., an expanded range of control parameters that support TS phases, does not automatically imply more robust topological phases. In general, the robustness of a TS phase is controlled by the topological (quasiparticle) gap that protects it. Using a self-consistent mean field theory to determine the pairing potential, we calculate the dependence of the topological gap on the relevant system parameters, including the hopping amplitude, Rashba coupling, Zeeman field, and chemical potential. As expected, we find that the pairing potential increases with reducing the bandwidth and, generally, decreases with the applied Zeeman field. However, accessing the topological regime requires a large-enough Zeeman field. In particular, $\Gamma$ should exceed the value of the (mean-field) pairing potential, $\Delta$. As a result, in the narrow lattice dispersion regime the pairing potential is significantly enhanced as compared to the wide-band, quadratic dispersion regime, but, ultimately, its values in the topological phase are constrained by the Zeeman field $\Gamma$ satisfying the topological condition. In addition, the pairing potential represents an upper bound for the topological gap. To determine the maximum possible values of this gap (within our mean-field approach), we systematically investigate its dependence on the system parameters and identify the optimal regimes. Our best-case scenarios correspond to topological superfluid phases with gaps on the order of $10-12.5\%$ of the attractive interaction strength, $V_{0}$. This is a significant improvement over the topological gaps achievable in naturally occurring systems, such as noncentrosymmetric superconductors.

Our work sets the stage to address additional important open questions regarding correlation effects and the thermal stability of topological superfluidity with broken time-reversal symmetry. The effect of correlations on the pairing potential can be investigated using non-perturbative approaches, such as, for example, the dynamical cluster approximation [34]. On the other hand, the Berenski-Kosterlitz-Thouless transition temperature, $T_{\text{BKT}}$, establishes the temperature below which we can expect to observe topological superfluidity. $T_{\text{BKT}}$ is proportional to the superfluid stiffness [8, 35]. In contrast to strictly flat-band systems, the superfluid stiffness is not zero in the present case where the band is narrow but not completely flat. Nonetheless, there should be a contribution to the superfluid stiffness from quantum geometry. Future work will examine the geometric contribution [36, 37, 38, 39] to superfluid stiffness in modeling the topological system considered in this paper with broken time-reversal symmetry.

ST and TDS acknowledge support from ONR-N000142312061. VS acknowledges support from AFOSR (FA2386-21-1-4081, FA9550-23-1-0034, FA9550-19-1-0272). ST and VS acknowledge support from ARO W911NF2210247.