Resonant superfluidity in the Rabi-coupled spin-dependent Fermi-Hubbard model

The RSFHM represented in terms of $\hat{c}_{\sigma,j}^{\dagger}$ and $\hat{c}_{\sigma,j}$ as in Eq. (1) is physically meaningful in the sense that the component index $\sigma$ indeed corresponds to two internal states of an atom. However, some of the essential properties of the model become clearer by considering a unitary rotation of the Hamiltonian with an angle which diagonalizes the on-site Hamiltonian at a given value of $\delta$ and $\Omega_{0}$. The unoccupied and doubly occupied states ($|0\rangle$ and $|2\rangle=\hat{c}^{{\dagger}}_{2}\hat{c}^{{\dagger}}_{1}|0\rangle$) are eigenstates of this Hamiltonian with eigenvalues 0 and $-U$ respectively, while the two remaining eigenstates are given by the creation operators $\hat{c}_{+}^{\dagger}=\alpha\hat{c}_{1}^{\dagger}+\beta\hat{c}_{2}^{\dagger}$ and $\hat{c}_{-}^{\dagger}=-\beta\hat{c}_{1}^{\dagger}+\alpha\hat{c}_{2}^{\dagger}$, where $\alpha=\cos(\theta),\beta=\sin(\theta)$. Here $\theta$ is an angle determined by the values of $\delta$ and $\Omega_{0}$ by the equations $\alpha\beta=\frac{\Omega_{0}}{2\sqrt{\Omega_{0}^{2}+\delta^{2}}}$ and $\beta^{2}-\alpha^{2}=\frac{\delta}{\sqrt{\Omega_{0}^{2}+\delta^{2}}}$. The corresponding on-site eigenvalues are given by

Rotating the full Hamiltonian at the angle $\theta$, it can be written in terms of the operators $\hat{c}_{+,j}$ and $\hat{c}_{-,j}$ as

with

The interaction and chemical potential terms are covariant under this transformation, while the $\delta$ and $\Omega_{0}$ fields combine as a Zeeman-type field for the new Hamiltonian with the direction determined by the angle $\theta$. For $t_{1}=t_{2}$ the rotated Hamiltonian is just the FHM with a Zeeman field that has a symmetric contribution from the Rabi and Zeeman terms of the original Hamiltonian. The SU(2) symmetry, which is present before the application of the on-site fields, remains such that the physics is essentially the same regardless of the orientation of the field. This model has been investigated in the literature, famously leading to the emergence of a Fulde–Ferrell–Larkin–Ovchinnikov (FFLO) phase between the BCS phase at small fields and the fully polarized phase at large fields [33, 34]. See also [35, 36] for exact numerical studies of the FFLO phase in a population-imbalanced one-dimensional FHM.

On the other hand, since $t_{1}\neq t_{2}$ breaks the initial SU(2) symmetry, the orientation of the on-site field becomes important. For $\Omega_{0}=0$ the rotated Hamiltonian corresponds to the original Hamiltonian, that is, a spin-dependent Hubbard model with a Zeeman field $\delta$ which has been studied in, e.g., [37]. This case has been investigated in order to understand how the FFLO phase is affected by asymmetric hopping. For $\Omega_{0}\neq 0$ the effective hoppings of the rotated components change but are still asymmetric and an extra symmetric hopping term which transforms between the components is added. This means that $[\hat{H}_{\theta},\hat{n}_{\pm,j}]\neq 0$ breaking the U(1) symmetry of these components, unlike the $t_{1}=t_{2}$ or $\Omega_{0}=0$ case. So even if the system is defined in the rotated frame, the Hilbert space cannot be divided into separate spin sectors. On the other hand, the system conserves the total number of particles which is therefore still a good quantum number.

For attractive interactions, quasi-long-range order (QLRO) of the pairing and/or CDW can be formed in the ground states of the RSFHM model. The properties of the pairs and density waves can be probed by the correlation functions

Here $\bar{n}$ is the average density and is subtracted in order to remove the trivial lattice ordering. Note that both functions are invariant under the rotation of the Hamiltonian and can therefore be calculated in either frame, using the relevant creation and annihilation operators. The numerical calculations in this manuscript are performed in the original frame of the Hamiltonian given in Eq. (1). In accordance with general convention for one-dimensional systems (see, e.g., [39, 40]), we refer to the presence of the off-diagonal QLRO, signified by a power-law decay of $P(|j-k|)$, as a SF (quasi-)order throughout the manuscript. The density-density correlations $C(j,k)$ can similarly display off-diagonal QLRO corresponding to a CDW (quasi-)order. The model under consideration can also display true long-range CDW order depending on the parameters.

A simple way to gauge the degree of pairing is the average number of pairs $n_{\text{pair}}=2N_{\text{pair}}/N$ where $N$ is the number of particles and

while polarization along the $\theta$ direction

quantifies the imbalance between $+$ and $-$ particles.

The presence of QLRO in the real-space correlation functions can also be probed by the eigenvalues of $P(j,k),C(j,k)$ or the Fourier transforms of $P(j,k),C(j,k)$, see, e.g., [41, 39, 40]. The Fourier transforms of $P(j,k)$ and $C(j,k)$ correspond to the pair momentum distribution and static structure factor respectively

where $L$ is the number of sites, and these are of particular interest as both are experimentally measurable in cold atom setups [13] allowing for the experimental determination of QLRO. If the system is homogeneous and the boundary condition is periodic, the peaks in Fourier space are equivalent to the maximum eigenvalues of the real-space correlation functions, which can be diagonalized as $P(j,k)/N_{\text{pair}}=\sum_{n}\lambda_{P,n}\phi_{P,n}(j)\phi_{P,n}(k)$ and $C(j,k)/\sum_{j}C(j,j)=\sum_{n}\lambda_{C,n}\phi_{C,n}(j)\phi_{C,n}(k)$, where $\phi_{P(C),n}(j)$ are the eigenfunctions of a pair (density) correlation function. We denote these eigenvalues as $\lambda_{P,0}$ and $\lambda_{C,0}$ for the SF and CDW correlations respectively. Even in the case of an open boundary condition, which we adopt in our DMRG analyses shown below, the maximal values of $n(k)/N_{\text{pair}}$ and $S(\nu)/\sum_{j}C(j,j)$ are well approximated by $\lambda_{P,0}$ and $\lambda_{C,0}$.

The fact that maximum eigenvalues are proportional to $L^{a}$ (with $0<a<1$) while the remaining eigenvalues do not scale with system size indicates the presence of QLRO in the system [41]. If the maximum eigenvalues are constant with system size, however, it indicates the absence of order. In order to test this rigorously and extrapolate the eigenvalues to the thermodynamic limit, one must investigate $\lambda_{P,0}N_{\text{pair}}$ or $\lambda_{C,0}\sum_{j}C(j,j)$ as a function of the system size. For constant $L$, however, a large value of $\lambda_{P,0}$ or $\lambda_{C,0}$ compared to the other eigenvalues signifies QLRO in the real-space correlation function at that system size. With the above normalization, the sum of eigenvalues is unity such that $\lambda_{P,0}$ and $\lambda_{C,0}$ correspond to the relative size of the largest eigenvalues and can measure the QLRO. Comparing $\lambda_{P,0}$ and $\lambda_{C,0}$ also allows us to probe the dominant order at a constant $L$ as the same normalization is employed for both.

In this subsection, we show results obtained by DMRG simulations of Eq. 1 for $t_{2}=0$ and half filling. All results are calculated for $N=100$ and $L=100$. In Fig. 1(a,b) we plot $\lambda_{P,0}$ and $\lambda_{C,0}$ as functions of the interaction strength $U/t$ and the Rabi coupling $\Omega_{0}/U$ for $\delta=0$. While the CDW QLRO is always large for $\Omega_{0}<U$, a clear maximum in $\lambda_{P,0}$ emerges when $\Omega_{0}$ is slightly smaller than $U$ (although the CDW QLRO is still dominant), particularly for strong interactions. The area in which this enhancement is observed is very small, however, as a transition to an almost polarized band insulator (in which both the SF and CDW off-diagonal correlations decay exponentially) takes place at $\Omega_{0}\approx U$. The transition to the band insulator takes place at smaller values of $\Omega_{0}/U$ for smaller interactions $U$.

$n_{\text{pair}}$ and $m=M/N$ (Fig. 1(c)) as well as the energy (Fig. 1(d)) as functions of $\sqrt{\delta^{2}+\Omega_{0}^{2}}$ for $\delta=0$ and $\delta=\frac{1}{2}U$ are also shown. Unlike the symmetric case ($t_{1}=t_{2}$) where the Rabi coupling is equivalent to a Zeeman field and $[\hat{H}_{\theta},\hat{n}_{\pm,j}]=0$ as $\gamma_{\theta}=0$, the magnetization is not a good quantum number for $t_{1}\neq t_{2}$ and $\Omega_{0}\neq 0$, i.e., $[\hat{H}_{\theta},\hat{n}_{\pm,j}]\neq 0$. In the former case, the system becomes fully polarized, while it only asymptotically approaches 1 after the transition in the latter case (see the inset of Fig. 1(c)) due to the presence of spin fluctuations. The behavior as a function of $\sqrt{\delta^{2}+\Omega_{0}^{2}}$ is very similar for any combination of $\Omega_{0},\delta$ (as long as $\Omega_{0}\neq 0$). Finally we plot the pair momentum distribution in Fig. 2. Before the transition a clear peak due to the presence of SF QLRO is visible. See the next two sections for a more detailed explanation of both the momentum distribution (e.g. why the peak is shifted from $k=0$ to $k=\pi$) and the physics in the paired regime in general.

In the limit of strong attractive interactions, an effective Hamiltonian in terms of the fermion-fermion composite bosons can be derived using second-order degenerate perturbation theory. The perturbation starts from the ground state of the Hamiltonian without hopping terms, given by $N/2$ pairs with an energy $E_{\rm pgs}=-NU/2$ considering the hopping terms as a perturbation. Note that, for $U<\sqrt{\delta^{2}+\Omega_{0}^{2}}$, the state consisting of $N$ spin-polarized fermions becomes the ground state, and the perturbation definitively breaks down. This approximately corresponds to the transition point identified in section III.B and the energies plotted in Fig. 1(d) correspond to $-NU/2$ (blue dotted line) before the transition, while the energy after the transition corresponds to $-\frac{N}{2}\sqrt{\delta^{2}+\Omega_{0}^{2}}$ (red dashed line). The alignment of the energies shown in Fig. 1(d) and the high degree of pairing displayed in Fig. 1(c) for $\sqrt{\delta^{2}+\Omega_{0}^{2}}<U$ indicates that the perturbation theory should describe the system well for large values of $U$.

Carrying out the perturbation (see the appendix for details), an effective Hamiltonian given in terms of the hard-core bosons $\hat{a}_{j}=\hat{c}_{j,+}\hat{c}_{j,-}$ can be found as

where

and

with $\tilde{\Omega}=\sqrt{\Omega_{0}^{2}+\delta^{2}}/U$.

The hardcore bosonic model can be equivalently described by an XXZ model, utilizing the mapping $\hat{a}^{\dagger}\rightarrow\hat{S}^{+}$, $\hat{a}\rightarrow\hat{S}^{-}$ and $\hat{n}\rightarrow\hat{S}^{z}+\frac{1}{2}$, which results in

where $J=2t_{\text{eff}}$ and $\Delta=\frac{V_{\text{eff}}}{2t_{\text{eff}}}$. Note that the model can be equivalently described by the antiferromagnetic model by applying the unitary transformation $-\hat{H}(J,-\Delta)=\hat{U}^{\dagger}\hat{H}(J,\Delta)\hat{U}$ and $(\Delta,J)\rightarrow(-\Delta,-J)$, where $\hat{U}=e^{-i\sum_{j=1}^{M}\hat{S}^{z}_{j}}$ [42]. For $t_{1}=t_{2}$ the model corresponds to the SU(2)-symmetric Heisenberg model for any values of $\delta$ and $\Omega_{0}$. The transverse and longitudinal spin correlations, respectively, correspond to the SF and CDW correlations in the full model.

For $\Omega_{0}=0$ the above Hamiltonian corresponds to earlier derivations [43, 27]. In this case $t_{2}=0$ (or equivalently $t_{1}=0$) maps to the pure Ising model. Finite asymmetry interpolates between the two limits, with the SU(2) symmetry being slowly restored as $t_{2}\rightarrow t_{1}$. A Kosterlitz-Thouless (KT) type transition takes place at $\Delta=1$ at half filling, with long-range CDW (antiferromagnetic order in the spin chain) present for $\Delta>1$ in the thermodynamic limit [44, 42, 43, 27]. Due to the exponentially slow opening of the gap, however, QLRO is still observed for relatively large system sizes even at half filling.

Similar observations can be made for the finite $\Omega_{0}$ case. The sign of $t_{\text{eff}}$ changes when $\sqrt{\Omega_{0}^{2}+\delta^{2}}$ is larger than some value, resulting in a change of sign for both $J$ and $\Delta$ which means that the model corresponds directly to the antiferromagnetic XXZ model. In the limit of $\sqrt{\Omega_{0}^{2}+\delta^{2}}=U$ ($\tilde{\Omega}=1$), $\Delta=-1$ so that the SU(2) symmetry is restored for any combination of $t_{1}$ and $t_{2}$. This corresponds to the maximum possible value of the ratio between the effective hopping $t_{\rm eff}$ and repulsive interaction $V_{\rm eff}$ and leads to the same degree of superfluidity as the $t_{1}=t_{2}$ case. The effective model breaks down in this limit, however, as this is the exact point at which the energy of the polarized state becomes smaller than that of the paired state. The exact point at which SU(2) symmetry is restored can therefore not be observed in the original FHM. The ratio $\frac{|t_{\text{eff}}|}{V_{\text{eff}}}$ for $t_{2}=0$ is plotted as a function of $\delta$ and $\Omega_{0}$ in Fig. 3. This explains why the results for finite $\delta$ and $\delta=0$ in section III.2 are so similar, in both cases the resonant enhancement is observed as long as $\Omega_{0}\neq 0$. From this, we conclude that there is a resonant enhancement of the superfluidity for $\sqrt{\Omega_{0}^{2}+\delta^{2}}\rightarrow U$ as observed in Fig. 1(a). $\tilde{\Omega}=\sqrt{\Omega_{0}^{2}+\delta^{2}}/U$ determines the importance of the Rabi contribution to the effective hopping with $\tilde{\Omega}=1$ corresponding to the resonant limit.

The change of sign for $t_{\rm eff}$ in the effective model shifts the energy bands and the lowest energy is therefore expected at $k=\pi$ rather than at $k=0$. From the pair momentum distribution plotted in Fig. 2 (based on the DMRG calculations of the full FHM), we see that the peak is indeed observed at $k=\pi$.

Owing to the presence of the U(1) symmetry breaking for the individual components $\zeta=+$ and $-$ in the Hamiltonian, the system will always have pair-breaking fluctuations. This cannot be captured by the effective model, which describes a definite, conserved number of fermion pairs which leads to an XXZ model with a well-defined and finite spin gap. However, the contribution beyond the effective theory is relatively less important in the regime of applicability of the perturbation theory, which explained the numerical results at half filling fairly well as discussed in the previous section. We therefore expect to be able to understand the important physical properties in the paired regime from the effective model, but will also confirm by numerical calculations of the full model.

The one-dimensional XXZ model can be solved exactly utilizing the Bethe ansatz approach [45, 46]. As the magnetization in the XXZ model is connected to the number of particles per site in the Hubbard model, we can utilize the results and formulas presented in [46] to understand the phase diagram of the effective model. Based on these prior results and numerical evaluations of the full model, the $n$ vs $\tilde{\Omega}=\sqrt{\Omega_{0}^{2}+\delta^{2}}/U$ schematic phase diagram of the RSFHM with $t_{2}=0$ at constant (relatively large) $U$ is given in Fig. 4, where $n=N/L$ is the number of particles per site.

For $n=1$ (half filling), the system displays long-range CDW order for any $\tilde{\Omega}<\tilde{\Omega}_{c}$ according to the effective model. Away from half filling, however, the effective model corresponds to a Luttinger liquid. In the Luttinger liquid regime, the decay of the two correlation functions can be expressed as functions of $r=|i-j|$ [46, 47]

where $\kappa_{\rho}$ and $\kappa_{\sigma}$ are, respectively, the Luttinger parameters for the charge and spin degrees of freedom and $k_{F}$ corresponds to the Fermi momentum. Within the effective theory corresponding to the paired regime, the XXZ model is spin-gapped, and we can consider $\kappa_{\sigma}=0$. The dominant QLRO can then be determined by calculating the Luttinger charge parameter $\kappa_{\rho}$. The shape of the $\kappa_{\rho}=1$ line in Fig. 4, which separates the two regimes of dominant SF or CDW QLRO in the system, is based on numerically solving the Bethe ansatz equations. The Luttinger parameter can also be approximately obtained from the numerical DMRG results in the full model as [48] $\kappa_{\rho}=C(k=2\pi/L)L/2$ and is expected to converge to a constant with increasing system size in the regime of applicability of the Luttinger theory. The SF QLRO for the Luttinger liquid corresponds to a standard SF (i.e. a single superfluid peak in the pair momentum distribution), and the resonantly enhanced superfluid region described by the effective model is therefore not of the FFLO-type.

At $\tilde{\Omega}\gtrsim 1$, the effective model breaks down as the system becomes mostly polarized, requiring analyses of the full model for understanding the physics. In order to understand this regime and confirm the predictions of the effective model in the regime of its applicability, we consider the decay of the correlation functions with increasing $r=|i-j|$ in the full system. In Fig. 5(a), we show the Luttinger parameter $\kappa_{\rho}$, which is closely related to the algebraic decay of the correlation functions through Eqs. (14) and (15), as a function of the system size (with a maximum of $L=200$) for the half filling ($n=1$) and quarter filling ($n=1/2$) cases at $U=30t_{1}$ and $\Omega_{0}=28t_{1}$ ($\tilde{\Omega}=0.93$), which is in the resonant region. Figure 5(b-d) show three representative parameter sets of the correlation functions as functions of $r=|i-j|$, which respectively correspond to the cases at half filling ($n=1$) and at quarter filling ($n=1/2$) in the resonant regime $\Omega_{0}=28t_{1}$ ($\tilde{\Omega}=0.93$), and the one at quarter filling ($n=1/2$) in the polarized regime $\Omega_{0}=34t_{1}$ ($\tilde{\Omega}=1.13$). In order to minimize the boundary effects, we consider $P(r)=P(\frac{L+r}{2},\frac{L-r}{2})$ and $C(r)=C(\frac{L+r}{2},\frac{L-r}{2})$. In general the boundary effects start to come into play as $r>L/2$ (see Fig. 5(b-d)) after which the behavior changes drastically. Indeed the calculations at larger system size clearly shows that this change is a boundary effect.

For quarter filling in the resonant regime, we clearly see in Fig. 5(a) that the Luttinger parameter is well converged to $\kappa_{\rho}\approx 1.3$ at $L>50$. This is consistent with the fact that both the SF and CDW correlations display a power-law decay with $P(r)\propto r^{-0.84}$ and $C(r)\propto r^{-1.2}$ as shown in Fig. 5(b). Indeed, these exponents are close to those predicted by the Luttinger theory, which is respectively $-0.76$ and $-1.3$. The actual decay exponents are slightly biased towards CDW QLRO compared to those expected from the Luttinger liquid, which is likely due to the presence of finite polarization in the full model. In the polarized regime, both $P(r)$ and $C(r)$ also display a power-law decay (Fig. 5(c)), implying that there is no phase transition between the two regimes. While the pairing fraction in the polarized regime never becomes zero, it asymptotically approaches zero as $\tilde{\Omega}$ increases and the power law decay for the correlations is very fast $P(r)\propto r^{-\eta}~{}(\eta\approx 4)$. This is in stark contrast to the resonant regime where the SF correlations are dominant at quarter filling and we can therefore clearly distinguish the SF properties of these two regimes. For half filling, both correlation functions display a power-law decay (Fig. 5(d)), although a slight flattening of the $C(r)$ curve is visible for large $r$ at $L=200$. The Luttinger parameter decreases with the system size (Fig. 5(a)), which suggests that the system is in fact not a Luttinger liquid. This is consistent with the analysis of the effective theory where a KT transition to a long-range ordered CDW phase takes place at the resonant point and long-range CDW is expected for $\tilde{\Omega}<1$ at half filling. We therefore expect that $P(r)$ would start to show exponential decay for larger system sizes and that $C(r)$ would become constant. Note that this means the resonant regime will display SF-like behavior even at half filling for system sizes relevant in most cold atom lattice experiments. At half filling, there is a phase transition to a band insulator as mentioned in section III.2.

From the effective model, we expect the values of $\tilde{\Omega}$ where superfluidity is dominant to be broadened as $n$ is lowered. However, observing this in the DMRG calculations is difficult because the low-energy states are energetically closer the further one goes from the resonant limit and half filling. Instead, we focus on the resonantly enhanced region, as the DMRG method gives well-converged results in this regime, and investigate the CDW and SF QLRO as a function of the particle number per site $n$, corresponding to vertical lines within the non-polarized regime of the diagram given by Fig. 4.

To investigate the dominant QLRO as a function of $n$, we compute the Luttinger parameter $\kappa_{\rho}$ as a function of $n$ for $t_{2}=0$ in the strongly-interacting region $U=30t_{1}$, $\Omega_{0}=25t_{1}$, $\delta=0$ ($\tilde{\Omega}=0.83$) where the effective model is expected to work well and for more moderate interactions $U=5t_{1}$, $\Omega_{0}=2.8t_{1}$, $\delta=0$ ($\tilde{\Omega}=0.56$). A smaller value of $\tilde{\Omega}$ is chosen in the latter case as the polarization becomes dominant at smaller $\tilde{\Omega}$ for smaller values of $U$. From Fig. 6(a), in which $S(\nu)/L$ is shown, it is clear that the single $\pi$-peak at half filling for the structure factor splits into two symmetric peaks around $\nu=\pi$ and that their absolute values decrease as $n$ is lowered. In the upper panels of Fig. 6(b,c), $\lambda_{P,0}$ and $\lambda_{C,0}$ are plotted. $\lambda_{P,0}$ remains bigger than $\lambda_{C,0}$ for a relatively large $n$ region; $0<n\lesssim 0.8$ in the case of $U=30t_{1}$ and $\Omega_{0}=25t_{1}$ ($0<n\lesssim 0.5$ in the case of $U=5,\Omega_{0}=2.8$). In Fig. 6(b,c), we also plot the Luttinger parameter $\kappa_{\rho}$ in the lower panels and show that these crossing points of $\lambda_{P,0}$ and $\lambda_{C,0}$ correspond to the point where $\kappa_{\rho}$ crosses 1 signifying the change in dominant QLRO. Overall this implies that for lower $n$ the SF QLRO is dominant. These results correspond to a vertical slice of the schematic phase diagram in Fig. 4 with $\tilde{\Omega}=0.83$ and $0.56$, respectively, where both are below the value at which polarization becomes dominant, which is smaller for the smaller value of $U$.