Anisotropic pair superfluidity of trapped two-component Bose gases in an optical lattice

Ultracold gases in optical lattices are promising and flexible quantum simulators for the study of quantum phases that are not easily accessible in condensed-matter physics [1]. Experimental realization of the superfluid-to-Mott transition has paved the way for studies of strongly interacting Bose or Fermi gases in optical lattices [2]. Recently, Bose–Bose mixtures of ⁸⁷Rb and ⁴¹K have been produced and loaded into an optical lattice [3], which provide a new playground for investigating the interplay between kinetic energy, interaction and the spin degree of freedom. In a further line of investigation, a Bose–Bose mixture of two hyperfine states of ⁸⁷Rb has been used as a spin-gradient thermometer for measuring the temperature of ultracold gases in an optical lattice [4, 5]. Moreover, the effect of a second species on bosonic superfluidity in an optical lattice has been studied [6]. In these mixtures of different species or different hyperfine states, the most fundamental physics associated with quantum magnetism and the spin degree of freedom can be explored [7].

One interesting ground state of bosonic gases is the pair-superfluid phase (PSF) for attractive interactions, which has been predicted and studied theoretically both in free space [8, 9] and in optical lattices [10–19]. However, it is widely believed that the stability of the bosonic many-body system is questionable when interactions between atoms are attractive. Recently, it was found that a three-body hard-core constraint can stabilize a system of bosons in an optical lattice with attractive two-body interactions [20], and numerical simulations have been performed to study the pair superfluidity of three-body-constrained bosons in an optical lattice [21–23]. For two-component bosons, stable heteronuclear ⁸⁷Rb–⁴¹K mixtures, with negative inter-species interactions tuned via Feshbach resonances [24], have been realized recently. Now the question is how to theoretically understand the PSF of two-component bosonic gases in an optical lattice. Qualitatively, the PSF can be viewed as a condensate of bosonic pairs of different species or hyperfine states due to a second-order hopping of the pairs, but with a strongly suppressed first-order tunneling of single atoms. For hard-core bosons with total filling n = 1, the PSF for attractive interactions is equivalent to the XY -ferromagnetic phase (or super-counter-fluid state) for repulsive interactions, i.e. the PSF consists of particle–particle pairs of different species while the XY -ferromagnetic phase is composed of particle–hole pairs. It is apparent that the PSF can only exist at very low temperature compared to the critical temperature of quantum magnetic phases. The latter are also governed by second-order tunneling processes, leading to an effective spin-exchange coupling [25, 26], which has already been observed for a double-well system [27, 28]. It is expected that these phases can be detected via momentum-space correlations observed in time-of-flight measurements [16, 29] or by optical microscopy with single-site resolution [30, 31].

Previous studies of the PSF in a two- or three-dimensional optical lattice mostly focus on symmetric parameters for the two species (for an exception, see [10, 18]), since this is the most favorable condition for the PSF [16]. However, there is still a lack of detailed quantitative investigations of the PSF for the homogeneous system with asymmetric hopping amplitudes of the two species, for the trapped system and for imbalanced mixtures of the two species. Here we investigate the properties of the PSF of two-component ultracold bosons with attractive interspecies interaction, both in a homogeneous and a trapped optical lattice. For sufficiently low filling, this system can be well described by a single-band Bose–Hubbard model with pure on-site interaction. We investigate the homogeneous system by means of bosonic dynamical mean field theory (BDMFT), which is a non-perturbative approach toward strongly correlated bosonic systems [32], and the trapped system by real-space bosonic dynamical mean field theory (RBDMFT) [33], which includes arbitrary inhomogeneity such as a harmonic trap. For the homogeneous system, we focus on the phase diagram with filling number n = 1. In particular, we also present the phase diagram for the experimentally realized heteronuclear ⁸⁷Rb–⁴¹K mixtures, where double-species Bose–Einstein condensates with negative inter-species interactions have been observed [24]. For the trapped Bose–Bose mixture we study the coexistence of Mott insulator, superfluid and PSF.

The paper is organized as follows. In section 2 we give a detailed description of the Bose–Hubbard model and the BDMFT approach. Section 3 covers our results on the homogeneous Bose–Bose mixture in an optical lattice and the effect of the trap. We conclude in section 4.

We consider a two-component bosonic mixture loaded into a two-dimensional (2D) or three-dimensional (3D) optical lattice. In experiments this mixture could consist of two different species, e.g. ⁸⁷Rb and ⁴¹K as in [3] or two different hyperfine states of a single isotope, e.g. ⁸⁷Rb as in [4]. Besides the optical lattice, we also impose an external harmonic trapping potential which introduces inhomogeneity in the system. For sufficiently low filling, the whole system can be effectively described by a two-component inhomogeneous Bose–Hubbard model within the single-band approximation

$$\begin{equation}\fl \mathcal{H}=- \sum_{\stackrel{\langle i,j\rangle}{\nu=b,d}} t_\nu (b^\dagger_{i\nu}b_{j\nu}+{\rm h.c}.)+\frac{1}{2}\sum_{i,\lambda\nu} U_{\lambda\nu} \hat{n}_{i,\lambda} (\hat{n}_{i\lambda}-\delta_{\lambda\nu}) +\sum_{i,\nu=b,d} (V_i-\mu_\nu)\hat{n}_{i\nu}, \end{equation} \tag{ 1 }$$

where 〈i,j〉 denotes summation over nearest-neighbor sites i, j, and the two bosonic species are labeled by the index λ(ν) = b, d. $\hat {b}^{\dagger }_{i\nu }$ ( $\hat {b}_{i\nu }$) denotes the bosonic creation (annihilation) operator for species ν at site i and $\hat {n}_{i,\nu }=\hat {b}^\dagger _{i\nu }\hat {b}_{i\nu }$ the corresponding local density. Due to possibly different masses or a spin-dependent optical lattice, these two species generally hop with non-equal amplitudes t_b and t_d. U_λν denotes the inter- and intra-species interactions, which can be tuned via Feshbach resonances or spin-dependent lattices [6, 34]. μ_ν is the global chemical potential for the two bosonic species, and V_i ≡ V^ν₀ r²_i denotes the external harmonic trap, where V^ν₀ is the strength of the harmonic trap for the ν component and r_i is the distance from the trap center.

To obtain the ground state of this system, we apply BDMFT for a homogeneous lattice [32], and its real-space generalization (RBDMFT) for the trapped system [33]. Within RBDMFT, the Hamiltonian (1) is mapped onto a set of individual single-site problems where the physics of lattice site i is described by a local effective action

$$\begin{eqnarray}&&\fl S_{{\rm eff}}^{(i)}\!=\!\int_0^\beta \mathrm{d}\tau\,{\rm d}\tau'\!\!\!\sum_{\lambda,\nu=\{b,d\}} \mathbf{b}_{\lambda}^{(i)}(\tau)^\dagger{\mathcal{G}}_{0,\lambda\nu}^{(i)}(\tau\!-\!\tau')^{-1} \mathbf{b}_{\nu}^{(i)}(\tau')+\int_0^\beta {\rm d}\tau\!\left\{\sum_{\lambda,\nu} \frac{1}{2}U_{\lambda,\nu}n^{(i)}_{\lambda}(\tau)\Big(n^{(i)} _{\nu}(\tau)\!-\!\delta_{\lambda\nu}\Big)\right. \nonumber\\ &&-\left.\sum_{\langle 0j\rangle , \nu} t_\nu \Big(b_{\nu}^{(i)}(\tau)^*\phi_{j,\nu}(\tau)+b_{\nu}^{(i)}(\tau)\phi_{j,\nu}^*(\tau)\Big)\right\}. \end{eqnarray} \tag{ 2 }$$

In this action τ denotes imaginary time and the function $\mathcal {{G}}_{0,\lambda \nu }^{(i)}(\tau -\tau ')$ is a local non-interacting propagator interpreted as a dynamical Weiss mean field which is determined in a self-consistent manner. We use the Nambu notation b⁽ⁱ⁾_0,ν(τ) ≡ (b⁽ⁱ⁾_0,ν(τ),b⁽ⁱ⁾_0,ν(τ)*). Moreover, the superfluid order parameters are defined as

$$\begin{eqnarray*} &&\phi_{j,\nu}(\tau)=\langle b_{j,\nu}\rangle_0, \end{eqnarray*}$$

where the index 0 indicates that the average is taken for the cavity system, i.e. excluding the impurity site.

Now, each of these local actions can be treated as an impurity in the presence of a bath (i.e. the rest of the lattice) and is obtained in a self-consistent manner. In practice, we start with an initial set of the local Weiss Green functions and local bosonic superfluid order parameters ϕ_i,ν(τ). After solving the local effective action (2), we obtain a set of local self-energies for each species Σ⁽ⁱ⁾_λν(iω_n) with ω_n being Matsubara frequencies. Then we employ the lattice Dyson equation in real-space representation in order to obtain the interacting lattice Green function

$$\begin{equation} \mathbf{G}({\rm i}\omega_n)^{-1}=\mathbf{G}_0({\rm i}\omega_n)^{-1}-\boldsymbol{\Sigma}({\rm i}\omega_n). \end{equation} \tag{ 3 }$$

The site dependence of the Green functions is indicated by boldface quantities which denote a matrix form with site-indexed elements. Here G₀(iω_n)⁻¹ is the inverse non-interacting Green function

$$\begin{equation} \mathbf{G}_0({\rm i}\omega_n)^{-1}=(\mu+\mathrm{i}\omega_n)\mathbf{1}-\mathbf{t}-\mathbf{V}. \end{equation} \tag{ 4 }$$

In this expression, the hopping amplitudes are given by matrix elements of t and the external potential is included via V_ij = δ_ijV_i. The diagonal elements of the lattice Green function are then identified with the interacting local Green functions, i.e. G⁽ⁱ⁾(iω_n) = (G(iω_n))_ii where the spin dependence of the local Green functions is shown by boldface quantities which denote a matrix form. Eventually, the self-consistency loop is closed by specifying the Weiss Green function via the local Dyson equation

$$\begin{equation} {\mathcal{G}}^{(i)}_{0}({\rm i}\omega_n)^{-1}= {G}^{(i)}({\rm i}\omega_n)^{-1}+{\Sigma}^{(i)}({\rm i}\omega_n). \end{equation} \tag{ 5 }$$

This self-consistency loop is repeated until the desired precision for the superfluid order parameter and the Weiss Green function is obtained.

In order to solve the local effective action, we map it onto a Hamiltonian described by an Anderson impurity model. This step is analogous to the solution of the local action in a homogeneous system [32]. The main difference to a homogeneous BDMFT calculation is that in the present case, the lattice Green function G⁽ⁱ⁾_λν(iω_n) is obtained via the real-space Dyson equation (3), which incorporates the effect of spatial inhomogeneity. We can now use standard techniques to solve the impurity problem. Here we adopt exact diagonalization as a solver [32, 35].

To include the effect of spatial inhomogeneity, we also employ an LDA approximation combined with single-site BDMFT. The advantage of this approach is the larger system size accessible. Within LDA + BDMFT, the local chemical potential for each species is set to μ_ν(r) = μ_ν − V (r). In this work, we apply both RBDMFT and LDA + BDMFT to the 2D trapped square lattice, and only LDA + BDMFT to the 3D trapped cubic lattice.

In this section, we will investigate properties of two-component bosonic mixtures with negative inter-species interactions in a homogeneous 3D optical lattice at zero and finite temperatures, and also in harmonically trapped, inhomogeneous systems, both in 2D and 3D. For the homogeneous system, we will study the stability of the PSF against asymmetric hopping and finite temperature. In the presence of an external harmonic trap, we will investigate finite-size effects of the PSF, both in 2D and 3D, within RBDMFT and LDA + BDMFT. We choose the absolute value of the inter-species interaction U_bd as our unit of energy in the following. For simplicity, we denote U_bd as U in the following.

3.1. Bose–Bose mixture in a three-dimensional (3D) cubic lattice

We first investigate Bose–Bose mixtures with negative inter-species interaction U < 0 in a 3D cubic optical lattice. The system is unstable for |U| > U_b,d, since the strongly attractive inter-species interaction cannot be compensated by repulsive intra-species interactions, leading to a collapse of the system. In the following discussion, we choose |U| < U_b,d. We explore the zero-temperature phase diagram with asymmetric hopping parameters at total filling n = 1 (with n_b = n_d = 0.5), as shown in figure 1. We observe three different phases: the PSF with 〈b〉 = 〈d〉 = 0 but ϕ_bd ≡ 〈bd〉 ≠ 0, the superfluid phase (SF) with 〈b〉,〈d〉 > 0 and a charge density wave (CDW) with Δ_CDW = |n_α − n_α'| > 0 and 〈b〉 = 〈d〉 = 0, where n denotes the filling and α the sublattice (α = −α'). In particular, we confirm the existence of a PSF for asymmetric hopping amplitudes t_b ≠ t_d. In the regime of comparably weak hopping for both species, first-order tunneling is suppressed by the strong interactions. Formation of bosonic pairs between different species can thus be energetically favored and will typically compete with single-species condensation, since the bosonic pairs can hop via second-order tunneling. As a result, the PSF will have a non-vanishing order parameter ϕ_bd but vanishing superfluid order 〈b〉 and 〈d〉. On the other hand, when both species acquire comparably large hopping, a superfluid phase will appear with 〈b〉 > 0 and 〈d〉 > 0. For asymmetric hopping, we observe that CDW order develops, since the small hopping amplitude for one component localizes bosonic pairs of the two species and CDW order is favored by the system. We remark here that the transition from a PSF to a superfluid phase for symmetric hopping amplitudes occurs at the same value of t_ν/U_bd as that from the XY -ferromagnetic to the superfluid phase in the corresponding system with repulsive U_bd > 0 of equal magnitude. In contrast to the XY -ferromagnetic phase, however, the PSF also exists for non-integer filling. Note that we have also investigated the influence of finite-size effects on the phase boundary due to a finite number of bath orbitals (4, 5 and 6 bath orbitals). We find that their effect is smaller than the symbol size.

Zoom In Zoom Out Reset image size

The effect of finite temperature is shown in figure 2. Generally, the PSF is sensitive to temperature, since the pairs are formed in the weak hopping regime and their coherence can be easily destroyed by thermal fluctuations, due to their small effective pair tunneling of order O(t²/U). At finite temperature, the PSF regime shrinks in the weak hopping regime in favor of developing a new unordered phase (UN) with vanishing values for ϕ_bd, 〈b〉 and 〈d〉. To further understand this unordered phase, we calculate the dependence of the local density–density correlator g⁽²⁾ ≡ 〈n_bn_d〉 − 〈n_b〉〈n_d〉 on temperature, as shown in the lower panel of figure 2. We observe that g⁽²⁾ starts to decrease noticeably only above temperatures of the order of 10⁻¹|U|, which indicates that local pairs still exist below this temperature. We therefore conclude that the unordered phase, shown in the upper panel of figure 2, consists of non-coherent pairs of different species. In sufficiently deep optical lattices, these pairs are localized, while for larger hopping the pairs delocalize over the whole lattice. As a result, the local density–density correlator decreases as a function of t_b,d, as shown in the lower panel of figure 2. Another interesting feature of the temperature dependence of g⁽²⁾ is the increasing (non-monotonic) trend at low temperatures, since thermal fluctuations first localize and then break the pairs. We remark here that the temperature regime of non-condensed pairs (≈ 0.1|U|) is experimentally accessible [36], and could be detected via radio frequency spectroscopy [37]. We also observe that the (single-particle) superfluid phase remains almost unchanged for the temperature considered here.

Zoom In Zoom Out Reset image size

One crucial question regarding the observation of the PSF is how fragile it is against finite-temperature effects. To address this issue, figure 3 shows T_c as a function of the hopping amplitudes t_b = t_d at different interactions. We notice that T_c rises as the hopping amplitudes increase, due to the growing second-order tunneling which stabilizes long-range order. The inset of figure 3 shows the temperature dependence of ϕ_bd. It indicates a second-order phase transition from the PSF to the unordered phase. We also observe that the critical temperatures for the PSF shown here are comparable with those of the XY -ferromagnetic phase [33] and notably smaller than the coldest temperatures which have been measured in most experiments until now, with the exception of the MIT group where temperatures as low as 350 pK (≈ 0.01U_bd with t_b/U_bd ≈ 0.029) have been achieved [5].

Zoom In Zoom Out Reset image size

To verify the validity of the BDMFT results, comparison has been made with a hard-core boson model solved by a tensor-product-state approximation [17], as shown in figure 4. We find excellent agreement between the two methods. We also plot the phase diagram for soft-core bosons (U_b = U_d = 2|U|), and observe that in this case, the phase boundary, between the PSF and the superfluid phase is shifted to lower hopping values.

Zoom In Zoom Out Reset image size

3.2. Rubidium–potassium mixture

Our investigations have so far focused on symmetric interactions U_b = U_d, which is a good approximation for mixtures of hyperfine states of ⁸⁷Rb [4]. However, this symmetry is not present for mixtures of ⁸⁷Rb and ⁴⁰K, where a negative inter-species interaction has been achieved via a Feshbach resonance [24]. Here we consider a mixture of ⁸⁷Rb and ⁴¹K loaded into a 3D cubic lattice with wavelength λ = 757 nm, which yields equal dimensionless lattice strength s for the two species. Due to different masses, the ratio of the intra-species interaction strengths is then fixed to U_Rb/U_K = m_Ka_Rb/m_Rba_K ≈ 0.72 and the ratio of the hopping amplitudes to t_Rb/t_K ≈ 0.47, where E_R is the recoil energy.

Now we explore the phase diagram of ⁸⁷Rb and ⁴¹K mixtures in a 3D cubic lattice and make predictions for ongoing experiments. Since the depth s of the optical lattice and the inter-species scattering length a_RbK are tunable with high accuracy, we show in figure 5 the phase diagram in the a_RbK–s plane for total filling n = 1 (n_b = n_d) at zero temperatures. At zero temperature, three phases appear: superfluid, PSF and CDW. When the scattering length a_RbK is small, the system is in a superfluid phase for a shallow lattice. When the depth of the lattice is increased, the ratio of t_Rb/U_Rb decreases, resulting in a strong suppression of first-order tunneling. The dominant process will then be hopping of composite pairs, which energetically favors PSF. We also observe CDW phase in the intermediate regime, where it is surrounded by PSF. It is expected that PSF and CDW shrink to a small region in parameter space at finite temperature due to thermal fluctuations, which is similar to results in [32].

Zoom In Zoom Out Reset image size

3.3. Trapped Bose–Bose mixtures in two-dimensional and 3D lattices

In this section, we simulate the two-component bosonic system in both 2D and 3D in the presence of a harmonic trap, as is relevant for most experiments. In particular, we investigate the stability of the PSF in the trapped system. Here we choose a 41 × 41 square lattice for the 2D case and a 41 × 41 × 41 cubic lattice for the 3D case. In 2D, in our simulations we apply both RBDMFT and BDMFT + LDA, while in 3D we only use BDMFT + LDA due to its lower computational effort. For simplicity, we investigate rubidium–rubidium mixtures with V^ν₀ ≡ V₀.

3.3.1. Balanced mixture

Figure 6 shows the density distributions n_b, order parameter ϕ_b and correlator ϕ_bd for the PSF versus radius r at different hopping amplitudes in a trapped 2D optical lattice. Since the PSF is stabilized only within a narrow region for the symmetric parameters (see figure 4), the harmonic trap should be very shallow and the hopping amplitudes need to be fine-tuned. Otherwise, the system will go through a phase transition directly from a Mott-insulating to a (single-component) superfluid phase. Here we choose completely symmetric parameters: t = t_b = t_d and U_b = U_d = 2|U| with balanced filling for the two components. Therefore only one value for n_b,d and ϕ_b,d is shown in figure 6. We observe that a wedding-cake structure appears in the trapped system, and the coexistence of different phases sensitively depends on the hopping amplitudes. At lower hopping t = 0.07|U|, only two phases appear, and the corresponding phase transition is from a Mott insulator with total filling n = 2 to a PSF with total filling 0 < n < 2 indicated by the non-vanishing value of ϕ_bd while ϕ_b = 0. If we increase the hopping amplitudes, the first-order tunneling of single atoms will increase, which induces large density fluctuations in the system, leading to a phase transition from the PSF to the superfluid phase. We clearly observe this effect from the middle panel of figure 6 where the superfluid phase starts to appear in the middle of the PSF. With further increase of the tunneling amplitudes the superfluid dominates, as shown in the lower panel of figure 6, where the PSF completely disappears at t = 0.7|U|.

Zoom In Zoom Out Reset image size

Figure 7 shows a comparison between the results of RBDMFT and those of BDMFT + LDA for a 2D square lattice. We observe good agreement between the two methods except close to the phase transition. In spite of the artificially sharp phase transition feature of LDA, the results of BDMFT + LDA are still reliable with sufficient accuracy in the regime away from the transition. We will therefore apply BDMFT + LDA to tackle the 3D case due to its higher computational efficiency compared to RBDMFT.

Zoom In Zoom Out Reset image size

Let us now investigate the stability of PSF in the 3D case in the presence of a harmonic trap. Results obtained within LDA are shown in figure 8. Here we choose completely symmetric parameters: t = t_b = t_d and U_b = U_d = 12|U| with balanced filling for the two components. Only one value for n_b,d and ϕ_b,d is shown in figure 8, respectively. Compared to 2D, we observe a similar scenario of phase coexistence in the trapped 3D cubic lattice: at lower hopping amplitudes the Mott insulator and the PSF are coexisting, at intermediate hopping amplitudes, the superfluid phase appears due to increased density fluctuations, and at even larger hopping amplitudes, the PSF will disappear.

Zoom In Zoom Out Reset image size

We are also interested in the effects of temperature on the PSF. Figure 9 shows the correlator ϕ_bd for different temperatures. We observe that the PSF is sensitive to thermal fluctuations. At finite T, the PSF is reduced in favor of developing an unordered phase characterized by ϕ_bd = 0. On the other hand, the density distribution is rather insensitive to small finite T.

Zoom In Zoom Out Reset image size

3.3.2. Imbalanced mixture

As pointed out above, asymmetry in the hopping of the two species does not destroy the PSF. However, filling imbalance hinders the formation of the pairs [16, 38]. We will now study this effect in more detail. The imbalance, N_b ≠ N_d, will be controlled by a non-zero chemical potential difference Δμ = μ_b − μ_d which can be viewed as an effective magnetic field. Results in 2D are shown in figure 10. Upon increasing Δμ, the PSF will cease to exist and will be replaced by a superfluid phase, since the chemical potential difference eventually exceeds the pairing gap for the PSF, allowing unpaired excess atoms to enter the PSF region. We find that the density distribution is almost unchanged for small Δμ, as shown in figure 10. When increasing the imbalance parameter Δμ further, the PSF will disappear already for a small population imbalance, since the particles form a conventional superfluid. Here we do not find any phase separation.

Zoom In Zoom Out Reset image size

Finally, we discuss the influence of population imbalance on the PSF in a trapped 3D cubic lattice using BDMFT + LDA. From our results shown in figure 11 we conclude that here the physics is qualitatively similar to the 2D case.

Zoom In Zoom Out Reset image size

We have investigated low-temperature properties of Bose–Bose mixtures with attractive inter-species interaction both in 2D and 3D optical lattices by means of BDMFT/RBDMFT. In particular, we found that the pair superfluid and charge density wave are stable for asymmetric hopping of the two species. We obtained the critical temperature of the PSF, which we found to be of the same order as that of the XY -ferromagnet in the corresponding system with repulsive interactions of equal magnitude. We have confirmed the stability of the PSF in a balanced Bose–Bose mixture in the presence of the harmonic trap both in 2D and 3D. On the other hand, we found that even a small population imbalance can destroy the PSF. This novel PSF quantum phase can be detected in future experiments via the momentum distribution of pairs, which shows signatures of the pair condensate [16, 38].

We acknowledge useful discussions with M Reza Bakhtiari. This work was supported by the China Scholarship Fund, the National Natural Science Foundation of China under grant numbers 11304386 (YL) and the Deutsche Forschungsgemeinschaft via SFB-TR/49, the DIP project HO 2407/5-1 and FOR 801. WH acknowledges the hospitality of KITP Santa Barbara, where this research was supported in part by the National Science Foundation under grant no. PHY05-25915.