Superfluidity of fermions with repulsive on-site interaction in an anisotropic optical lattice near a Feshbach resonance

We present a numerical study on ground state properties of a one-dimensional (1D) general Hubbard model (GHM) with particle-assisted tunnelling rates and repulsive on-site interaction (positive-U), which describes fermionic atoms in an anisotropic optical lattice near a wide Feshbach resonance. For our calculation, we utilize the time evolving block decimation (TEBD) algorithm, which is an extension of the density matrix renormalization group and provides a well-controlled method for 1D systems. We show that the positive-U GHM, when hole-doped from half-filling, exhibits a phase with coexistence of quasi-long-range superfluid and charge-density-wave orders. This feature is different from the property of the conventional Hubbard model with positive-U, indicating the particle-assisted tunnelling mechanism in GHM brings in qualitatively new physics.

We present a numerical study on ground state properties of a one-dimensional general Hubbard model (GHM) with particle assisted tunnelling rates and repulsive on-site interaction (positive-U), which describes fermionic atoms in an anisotropic optical lattice near a wide Feshbach resonance. For our calculation, We utilize the time evolving block decimation (TEBD) algorithm, which is an extension of the density matrix renormalization group and provides a well controlled method for one-dimensional systems. We show that the positive-U GHM, when hole doped from half-filling, exhibits a phase with coexistence of quasi-long-range superfluid and charge-density-wave orders. This feature is different from the property of the conventional Hubbard model with positive-U, indicating the particle assisted tunnelling mechanism in GHM brings in qualitatively new physics. The combination of Feshbach resonance and optical lattice techniques has opened up possibilities to investigate strongly interacting ultracold atoms under tunable configurations [1]. Ability to control such strongly interacting systems provides an unprecedented opportunity to explore interesting states of matter. Many interesting physics has been predicted for ultracold atom systems with fundamental Hubbard model hamiltonians. For example, with Bose-Hubbard model and its derivations, people have studied superfluid to Mott-insulator transition [2], existence of supersolid order [3], etc, for ultracold bosons while with Fermi-Hubbard model, Luther-Emery [4] and FFLO [5] phases are predicted to be observable for ultracold Fermions with attractive interaction. In particular, it is well known that for repulsive (positive-U ) conventional (Fermi-)Hubbard model the susceptibility for superfluid and charge density wave (CDW) orders are suppressed at low temperature and the leading quasilong-range order is given by a spin density wave (SDW) at any filling fraction [6].
However, in this work, we show that coexistence of quasi-long-range superfluid and CDW orders can be observed for fermionic atoms with repulsive on-site interaction in an anisotropic optical lattice near a wide Feshbach resonance. The interactions in this strongly interacting system is described by a one-dimensional positive-U general Hubbard model (GHM) with particle assisted tunnelling rates [7]. The GHM is an effective one-band Hamiltonian that takes into account the multi-band populations and the off-site atom-molecule couplings in an optical lattice near a wide Feshbach resonance (see the detailed derivation in Ref. [7]). It is interesting to note that the GHM with similar particle assisted tunnelling also arises in different physical contexts, as proposed in Ref. [8]. In contrast with the case of conventional positive-U Hubbard model, we show that the superfluid and CDW emerge as dominant quasi-long range orders over spin orders for the positive-U GHM when the system is significantly hole-doped below half-filling, although at or very close to half-filling, the dominant correlation in GHM is still anti-ferromagnetic. This feature indicates that the particle assisted tunnelling in GHM brings in qualitatively new physics. It makes the effective interaction in GHM doping dependent, showing different behaviors with a possible phase transition in between. We get our results through numerical calculation based on the time evolving block decimation (TEBD) algorithm [9,10], which, as an extension of the density matrix renormalization group method [11], is a well controlled approach to deal with one-dimensional systems. We compare our numerical results with some known exact results for the conventional Hubbard model, and the remarkably precise agreement shows that the calculation here can make quantitatively reliable predictions.
As shown in Ref. [7], a generic Hamiltonian to describe strongly interacting two-component fermions in an optical lattice (or superlattice) is given by the following general Hubbard model: where n iσ ≡ a † iσ a iσ , n i ≡ n i↑ + n i↓ , µ is the chemical potential, i, j denotes the neighboring sites, and a † iσ is the creation operator to generate a fermion on the site i with the spin index σ. The symbol σ stands for (↓, ↑) for σ = (↑, ↓). The δg and δt terms in the Hamiltonian represent particle assisted tunnelling, for which the inter-site tunnelling rate depends on whether there is another atom with opposite spin on these two sites. The particle assisted tunnelling comes from the multi-band population and the off-site atom-molecule coupling for this strongly interacting system [7]. For atoms near a wide Feshbach resonance with the average filling number n i ≤ 2, each lattice site could have four different states, either empty (with state |0 ), or a spin ↑ or ↓ atom (a † iσ |0 ), or a dressed molecule (d † i |0 ) which is composed by two atoms with opposite spins. The two atoms in a dressed molecule can distribute over a number of lattice bands due to the strong on-site interaction, with the distribution coefficient fixed by solving the singlesite problem. One then can mathematically map the dressed molecule state d † i |0 to a double occupation state a † i↓ a † i↑ |0 by using the atomic operators a † iσ [7] . After this mapping, the effective Hamiltonian is transformed to the form of Eq. (1). The GHM in Eq. (1) reduces to the conventional Hubbard model when the particle assisted tunnelling coefficients δg and δt approaching zero, as one moves far away from the Feshbach resonance. Near the resonance, δg and δt can be significant compared with the atomic tunnelling rate t due to the renormalization from the multi-band populations and the direct neighboring coupling [7].
We consider in this work an anisotropic optical lattice for which the potential barriers along the x, y directions are tuned up to completely suppress tunnelling along those directions. The system becomes a set of independent one-dimensional chains. We thus solve the GHM in one dimension through numerical analysis. For this purpose, first we transfer all the fermion operators to the hard core boson operators through the Jordan-Wigner transformation [12]. In the one-dimensional case, we can get rid of the non-local sign factor, and after the transformation the hard core boson operators satisfy the same Hamiltonian as Eq. (1). On each site we then have two hard core boson modes which are equivalent to a spin-3/2 system with the local Hilbert space dimension d = 4. We can therefore use the TEBD algorithm to solve this pseudo-spin system [9]. Similar to the density matrix renormalization group (DMRG) method [11], the TEBD algorithm is based on the assumption that in the one-dimensional case the ground state |Ψ = d i1=1 · · · d in=1 c i1...in |i 1 · · · i n of the Hamiltonian with short-range interactions can be written into the following matrix product form: where Γ [s]is denotes the matrix associated with site-s with the matrix dimension χ. When χ = 1, the assumption reduces to the mean-field approximation, and for a larger χ, the matrix product state well approximates the ground state as it catches the right entanglement structure for 1D systems [9,10]. To use the TEBD algorithm, we just start with an arbitrary matrix product state in the form of Eq. lation functions. This calculation has a well controlled precision since at each time step to update the matrix product state, the Hilbert space truncation error can be suppressed by choosing an appropriate matrix dimension χ [9]. In this calculation, we use the infinite lattice algorithm by assuming that the lattice is bipartite and the ground state has a translational symmetry for each sublattice [9]. This allows us to directly calculate the system in the thermodynamic limit.
To show that our calculation is capable of making reliable predictions, we first test our results by comparing them with some known exact results of the Hubbard model in certain cases. For the Hubbard model at half-filling n i↑ = n i↓ = 0.5, the ground state energy per site is known to have the analytic expression E = −4 ∞ 0 J0(ω)J1(ω)dω ω[1+exp(ωU/2)] in the thermodynamic limit from the exact Bethe ansatz solution [13], where J 0 and J 1 are Bessel functions and we have chosen the tunnelling rate t as the energy unit. In Fig. 1 (a), we show our numerical results for the ground state energy of the Hamiltonian (1) with δg = δt = 0, and one can see that it agrees very well with the exact energy of the Hubbard model in particular when U > t. The error is in general smaller than 10 −3 as shown in Fig. 1(b). In this and the following calculations, we choose the matrix dimension χ = 40. We have tried larger χ which gives better precision, but we choose χ = 40 to have a faster speed and its precision is enough for our purpose.
We have also tested the final state from our calculation by comparing its correlation functions with some known results. It is difficult to get correlations analytically from the Bethe ansatz solution, but from the bosonization ap- proach to the one-dimensional Hubbard model, we know its correlation functions take certain asymptotic forms. For instance, one can look at the spin-spin correlation, defined as S r ≡ s i · s i+r , where the spin operator for the site i is given by s i ≡ a † iα σ αβ a iβ /2 with α and β =↓, ↑ and σ standing for the Pauli matrices. The correlation S r is independent of i because of the translational symmetry. The Hubbard spin correlation function has the following asymptotic form [14] where K ρ is the Luttinger parameter whose value has been determined before from the exact Bethe ansatz solution [13,14], k F is the fermi momentum related to the filling number n i through k F = n i π/2, and A is a non-universal model dependent constant. In Fig.1 (c) and (d), we compare our calculation results for S r with this asymptotic form for filling number n i = 0.5 and 0.75, and the agreement is again remarkable as long as r is not too small (the expression of S r in Eq. (3) is not accurate for small r).
With the confidence in numerics built from the above comparison, we now present our main calculation results for the repulsive GHM in Eq. (1) with U > 0. Apart from the spin correlation S r defined before, we also calculate the charge-density-wave (CDW) correlation, defined as D r ≡ n i n i+r − n i n i+r , and the pair (superfluid) correlation, defined as P r ≡ a i↑ a i↓ a † i+r↓ a † i+r↑ . The results are shown in Fig. 2 for different filling fraction n i and for models with different particle assisted tunnelling rates δg and δt. First at half filling with n i = 1, the correlation functions S r , D r , and P r for the GHM with different δg and δt all look qualitatively similar to the corresponding results for the conventional Hubbard model, although with increase of the coefficient δg the spin correlation reduces a bit while the CDW and superfluid correlations increase slightly. Clearly, the dominant correlation in this case is in spin which suggests a quasilong range anti-ferromagnetic order. In this and following calculations, we take U = 8t for all the cases, which corresponds to a significant on-site repulsion.
Qualitatively different results show up when the system is doped with holes. At the filling fraction n i = 0.75, although for the Hubbard model the spin correlation is still the dominant one (the spin density wave order has been pinned to the corresponding 2k F = 3π/4), for the GHM with a noticeable δg, the superfluid and the CDW emerge as the leading quasi-long-range orders, and their correlations increase significantly and decay much slower in space compared with the spin correlation when δg grows. These features become more evident when we further increase the doping. For instance, at the right column of Fig. (2), we show the correlations for the filling fraction n i = 0.5. The qualitative behavior is similar to the case with n i = 0.75, but the CDW and superfluid correlations for the GHM get significantly larger at long distance, and the contrast with the Hubbard model becomes sharper. One also note that for all these calculations, change of the coefficient δt in the GHM makes little difference to the result. This is understandable as a significant positive U suppresses the possibility of double occupation in the lattice, and the δt term in the GHM has no effect without double occupation. The δg term in the GHM, however, is critically important, which favors superfluidity in general and brings in the qualitatively different features mentioned above.
To show the spatial structure of these quasi-long range (QLR) orders, we plot in Fig. 3 the spin, the CDW, and superfluid correlations in the momentum space for the GHM with δg = 3t at different filling fractions n i . The momentum space correlations are defined by the Fourier transform X k = 1/ √ N N r=0 X r cos(kr), where X stands for the correlations S, D, or P . From these momentum space curves, one can clearly see that this GHM at half filling has a QLR anti-ferromagnetic order (characterized by the peak at k = π), and away from half filling a QLR superfluid order (peak at k = 0) and CDW order (peaks at k = 2k F and 2π − 2k F , where 2k F = 3π/4 (π/2) for the filling fractions n i = 0.75 (0.5), respectively). The peaks in Fig. 3 have finite widths because these orders in 1D are only quasi-long-range with algebraic decay. Note that if we turn on small tunnelling interaction between different 1D chains, a leading QLR order, such as superfluid order, could be stabilized to a true long range order [6]. The GHM thus provides an example of a microscopic Hamiltonian that with hole doping from half filling, an anti-ferromagnetic phase could be transferred to a superfluid phase (or a CDW phase in some case depending on which order becomes more dominant with the interchain coupling). The correlations that characterize these QLR orders can be detected for the cold atomic gas, for instance, through the method described in Ref. [16].
In summary, we have investigated the ground state properties of the general Hubbard model with repulsive on-site interaction in one dimension through well controlled numerical analysis. For the system with significant particle assisted tunneling rates δg and δt, we have found coexistence of quasi-long range superfluid and charge-density-wave orders when the system is holedoped from half filling. This feature is in sharp contrast with convention Hubbard model, in which case for positive-U the charge and superfluid orders are always suppressed regardless of the filling fraction. With a combination of the Bosonlization approach and the numerical method here, it may be possible to determine the compete phase diagram for the GHM. The model here describes strongly interacting fermionic atoms in an anisotropic optical lattice. The possibility of a transition from an anti-ferromagnetic phase to a superfluid phase for the GHM with hole doping may also have interesting indications for other areas.
This work is supported under the MURI program, and under ARO Award W911NF0710576 with funds from the DARPA OLE Program.