Resonant Superfluidity in an Optical Lattice

We study a system of ultracold fermionic Potassium (40K) atoms in a three-dimensional optical lattice in the vicinity of an s-wave Feshbach resonance. Close to resonance, the system is described by a multi-band Bose-Fermi Hubbard Hamiltonian. We derive an effective lowest-band Hamiltonian in which the effect of the higher bands is incorporated by a self-consistent mean-field approximation. The resulting model is solved by means of Generalized Dynamical Mean-Field Theory. In addition to the BEC/BCS crossover we find a phase transition to a fermionic Mott insulator at half filling, induced by the repulsive fermionic background scattering length. We also calculate the critical temperature of the BEC/BCS-state and find it to be minimal at resonance.


Introduction
The first experimental realizations of Bose-Einstein condensates in dilute atomic gases of rubidium [1], lithium [2] and sodium [3] atoms initiated a new field of condensedmatter research, providing an ideal laboratory for comparing theoretical models and experimental results with high accuracy. In particular the important consequences of Bose-Einstein condensation could be investigated, which up to 1995 had remained an elusive and inaccessible phenomenon in experiments.
Not long after the first realization of BEC, ultracold fermionic gases were studied experimentally as well. The first important results of quantum degeneracy in trapped Fermi gases were obtained in 1999 at JILA [4] and later on by other groups [5,6]. A break-through experiment in this field was the investigation of fermionic superfluidity at the crossover between the BEC state and the Bardeen-Cooper-Schrieffer (BCS) state [7,8,9,10,11]. This is made possible by the use of Feshbach resonances, which have become an indispensable experimental tool for ultracold atom experiments. Feshbach resonances not only allow one to tune the interatomic interaction with high precision, but also make it possible to increase it to a level where the critical temperature becomes high enough for the investigation of attraction-induced superfluidity. On the contrary, away from resonance the critical temperature is usually exponentially suppressed and experimentally inaccessible. The regime of strong, even diverging interactions, the so-called unitarity region, on the other hand defines a new field of research where standard mean-field methods break down and the physics has to be described in a nonperturbative way. Moreover, this system allows for the experimental investigation of the BEC-BCS crossover: for negative scattering lengths the system is a BCS superfluid, whereas for positive scattering length fermionic atoms with opposite spin pair up to form a bosonic molecular bound state. More recent experimental work has focused on studying the effect of spin imbalance on the BCS state, i.e. the case when an unequal number of fermions occupies the two different spin states [12,13,14,15], as well as mixtures of fermions with unequal masses, such as 6 Li and 40 K [16,17].
Also the effect of periodic potentials on trapped Fermi gases has been studied experimentally [18,19,20,21]. Recently, evidence for a fermionic Mott insulator was obtained in a system of repulsively interacting 40 K fermions in an optical lattice [22,23]. Optical lattices and Feshbach resonances have been combined experimentally as well: the group at ETH Zürich reported the production of 40 K molecules in threedimensional cubic optical lattices using s-wave Feshbach resonances in early 2006 [19], but no evidence of a superfluid state in the lattice was found until later that year, when superfluid 6 Li was loaded in an optical lattice at MIT where both a condensate and an insulating state were observed [20]. The results where interpreted in terms of a superfluid-to-Mott insulator transition, for which a detailed theoretical description is still lacking.
In this work we study an ultracold mixture of fermionic atoms in two different hyperfine states in a three-dimensional optical lattice close to a Feshbach resonance. Schematic picture of the BEC-BCS crossover in an optical lattice. By tuning the interaction strength between the two fermionic spin states, one obtains a smooth crossover from a BEC regime of tightly bound bosonic molecules (a) to a BCS regime of large Cooper pairs (c). In between these two extremes, one encounters an intermediate crossover regime where the pair size is comparable to the interparticle spacing (b). For total fermionic filling one, the system can undergo a phase transition to the Mott insulator phase (d).
This system has all the characteristics of the continuum BEC/BCS crossover described above: For magnetic field values below the resonance, fermions with different spin form bosonic molecules (see Fig. 1). By varying the magnetic field the bosonic level is detuned relative to the fermionic one, which changes the ratio of the densities of fermions and molecular bosons as well as the the effective interaction between the fermions. On top of this BEC-BCS crossover physics, which is familiar from the system without lattice, new features emerge when an optical lattice is applied. The most prominent one is the occurrence of a fermionic Mott insulator for half-filled fermions deep in the BCS regime, which is stabilized by the repulsive fermionic background scattering length. As described below, we find a first order transition between the BEC/BCS state and the Mott insulator. For a total filling of two fermions per site the BCS state competes against a band insulating state [24,25,26].
The presence of an optical lattice allows to utilize powerful non-perturbative methods that are available for lattice systems. Here we apply Generalized Dynamical Mean-Field Theory (GDMFT) [27,28] to the system described above. However, GDMFT is a single-band approach, whereas Feshbach-resonant interactions in an optical lattice lead to a multi-band model [34,35,36,37,38,39]. We therefore first perform a mean-field decoupling of the higher bands, thereby deriving an effective single-band Hamiltonian which is self-consistently coupled to the higher bands.
The paper is organized as follows: in the next section we introduce the microscopic model and in Sec. 3 we introduce the Generalized Dynamical Mean-Field approach we use to solve this model. In Sec. 4 we present the result of our numerical calculations and in Sec. 5 we end up with concluding remarks. In the appendix, we describe in detail how the self-energy for the resonantly interacting Bose-Fermi mixture studied here can be calculated in the Dynamical Mean-Field framework.

Microscopic Model
Studying ultracold fermions close to a Feshbach resonance is a challenging problem. Due to the fact that exactly on resonance the scattering length is infinite, the standard fermionic Hubbard Hamiltonian cannot be defined. To solve this problem, it is necessary to formulate a two-channel Hamiltonian by separating out the resonance state and treating it explicitly [40]. The nonresonant contributions give rise to a background scattering length. As the Feshbach resonance occurs due to a coupling with the bosonic molecular state, the additional degrees of freedom introduced in the two-channel theory are bosons [40].
An ultracold atomic gas of fermionic atoms and molecular bosons close to a Feshbach resonance in the presence of an optical lattice is thus well described by a Bose-Fermi Hubbard model [35,41]. In our calculation we assume the molecular bosons to be in the lowest band. For the fermions, on the other hand, we have to take into account also the higher bands, in order to properly incorporate the two-body physics associated with the Feshbach resonance [35,36,42]. Since the bandwidth is much smaller than the band gap, we approximate the higher bands to be flat and only take into account the full band-structure for the lowest band. Moreover, we neglect the interaction between fermions in higher bands with each other and with the bosons. This is justified because the filling in the higher bands is very small, so that interaction effects are also small. The Hamiltonian thus has the following form: whereĉ † iσ,l is the creation operator of a fermion with spin σ for the l-th band on lattice site i.b † i is the creation operator of a boson at site i.n f iσ,l =ĉ † iσ,lĉiσ,l ,n f i,l =n i↑,l +n i↓,l are the fermionic number operators, andn b i =b † ibi is the bosonic number operator. t f and t b are the tunneling amplitude for fermions and bosons, respectively. U f , U b and U f b are Fermi-, Bose-, and Bose-Fermi-Hubbard interactions, respectively. These interactions arise due to the background scattering lengths. Furthermore µ is the chemical potential, δ is the detuning of the bosonic level, and ω is the frequency of the harmonic oscillator associated with an optical lattice well. Finally, g l = g 0 L (1/2) l (0) is the Feshbach coupling to the l-th band of the lattice, where g 0 is the Feshbach coupling for the lowest Hubbard band and L (1/2) l (0) is the generalized Laguerre polynomial. Choosing the Feshbach couplings in this way guarantees that the two-body physics associated with the Feshbach resonance is incorporated exactly [35,42].
The parameters of the generalized Hubbard Hamiltonian (1) are given by: Here is the recoil energy, V 0 is the amplitude of the optical lattice potential and λ is the laser wavelength. a f , a b , and a f b are the fermion-fermion, boson-boson, and fermion-boson background scattering lengths. In our calculation we approximate the background boson-boson and the Bose-Fermi scattering lengths by a b = 0.6a f [43] and a f b = 1.2a f [44]. Furthermore, B is the magnetic field, and B 0 and ∆B are the position of the Feshbach resonance and its width, respectively. ∆µ mag is the difference in the magnetic moment between the closed and open channel of the Feshbach resonance. Finally, m f and m b are the respective masses of the fermions and bosons.
The Hamiltonian (1) can be simplified by the following rescaling: such that the factor 3 ω 2 disappears.

Derivation of the effective single-band Hamiltonian
The multi-band Hamiltonian derived so far is very complicated, since it both involves strong correlations and many bands. Simply neglecting the higher bands would lead to an incorrect description close to the Feshbach resonance, since the Feshbach parameter g is very large and even exceeds the band gap [35]. However, the filling of fermions in the higher bands is strongly suppressed by the band gap. This allows us to perform a mean-field decoupling in the higher bands [35]. The lowest band is left untouched in this procedure, since the fermionic filling can be large there. We thus perform the following decoupling for l > 0 on each site: This step implies that the lowest band and the higher bands are only coupled in a meanfield way. They can thus be diagonalized separately, but are coupled by the mean-field self-consistency relations. The full Hamiltonian is now given bŷ where the following terms are added to the bosonic part of the lowest band Hamiltonian: where the mean-field ∆ has been defined as ∆ = − ∞ l=1 g l ĉ i↑,lĉi↓,l . For each of the higher bands l > 0 we obtain the following Hamiltonian (here we suppress the site index The system described by Eqs. (15), (16) and (17) needs to be solved self-consistently with respect to the mean fields b and ∆.
To diagonalize the Hamiltonian (17), one has to perform the following Bogoliubov transformation: where This leads to the following expectation values: where f (ω l ) is the Fermi function and T is the temperature. We use absolute values in the equation for ĉ l↑ĉl↓ because of the ambiguity of the sign, arising from the fact that still a divergence has to be subtracted (see below). The total number of fermions is equal to: This is a converging sum, which can be evaluated numerically. From Eq. (16) follows that we have to evaluate the sum which is divergent. This divergence always arises in the gap equation of the BCS model when the T-matrix is approximated by a delta-potential [35,42]. One way to solve this problem is by using a pseudo-potential instead of the delta-potential [42]. Here, however, we follow Ref. [35] and explicitly isolate the diverging contribution from the sum. First, we notice that for large l, ω l can be approximated by ω l = 2l ω −μ and tanh ω l 2kT 1. Therefore Here N is a large integer number (in our calculation we took N = 500). The first two terms of Eq. (27) are finite sums, but the last term diverges. This sum is known from the literature [42,45]. To separate the diverging part, we have to take the following limit: Since the diverging part √ π/r is independent of the model parameters, we can cure the divergence by neglecting this term [35,42]. Doing so, we obtain We now fix the sign by requiring ∆ > 0, since this solution minimizes the (free) energy.
Summarizing, we have reduced the multi-band problem to an effective single-band Hamiltonian:Ĥ whereĤ 0 f ,Ĥ b ,Ĥ 0 f b andĤ b are given by Eqs. (2), (3), (4) and (16), respectively. The chemical potential µ has to be adjusted such that the total filling is equal to the desired value n tot : This leads to the following self-consistency loop: we start from an initial guess of the superfluid order parameter b and calculate ∆ using Eq. (29). As a result we know all parameters in the Hamiltonian (30), and can find its eigenvalues and eigenvectors, and correspondingly calculate new correlation functions, including the superfluid order parameter b . With this step the self-consistency loop is closed.
It is worth noting that the effective single-band Hamiltonian we have derived here is different from the effective single-band model in terms of dressed particles derived in other approaches [36,39]: the bosons and fermions in our Hamiltonian correspond to the bare particles in the lowest band.

Generalized Dynamical Mean-Field Theory
To analyze the Hamiltonian (30) we use Generalize Dynamical Mean Field Theory (GDMFT) which is explained in detail in Ref. [27,28]. Here we only mention that within GDMFT one considers a single site which is self-consistently coupled to a dynamical fermionic bath corresponding to DMFT [29,30] and a static bosonic meanfield corresponding to bosonic Gutzwiller [31,32,33]. These are the leading order contributions in a 1/z-expansion of the effective action (z being the lattice coordination number). Hence GMDFT is exact in infinite dimensions and expected to be a good approximation for the cubic lattice considered here, for which z = 6. The typical accuracy for low-temperature expectation values is around 20 percent.
In the specific case considered here the system is described by a generalized single impurity Anderson model (GSIAM) with the following Anderson Hamiltonian: where z is the lattice coordination number and ϕ = b is the superfluid order parameter. k labels the noninteracting orbitals of the effective bath, V k are the corresponding fermionic hybridization matrix elements, W k describes the superfluid properties of the bath and a † kσ is the creation operator of a fermion in the k th orbital of the bath with spin σ.n f σ =ĉ † σ c σ is the fermionic number operator andn =n f ↑ +n f ↓ . To solve the Anderson Hamiltonian we use exact diagonalization as impurity solver [46,47,48,49]. In this algorithm the infinite number of orbitals in the Hamiltonian (32) is truncated and only a finite number of n s orbitals is considered. The resulting finitesize problem is fundamentally different from the finite-size problem of a finite number of lattice sites of the original Hubbard model, and the truncation procedure can be viewed as using a finite number of parameters (energy scales) to describe the local dynamics encoded in the Weiss Green's function: where β is the inverse temperature and ω n = (2n+1)π/β are the Matsubara frequencies.
To close the self-consistency loop by using the lattice Dyson equation we calculate the normal and superfluid Green's functions which can be written as follows [46] where ζ = iω n +μ − Σ(iω n ) and D(ε) is the non-interacting density of the states of the cubic lattice. Σ(iω n ) and Σ SC (iω n ) are the normal and superfluid self-energies, which as shown in Appendix A can be expressed via a set of higher order Green's functions: Here G(iω n ) = c σ,0 , c † σ,0 ω and F (iω n ) = c ↑,0 , c ↓,0 ω are the normal and superfluid single particle Green's functions. In addition we have also defined the following additional interacting Green's functions: and Z = n e −βEn (40) is the partition function. The relation between the Weiss field and the Green's function is given by the local Dyson equation: whereĜ is the matrix of interacting Green's functions, is the self-energy matrix and is the matrix of Weiss Green's functions.
To determine new parameters for the Anderson Hamiltonian, we fit the Weiss functions calculated from (33) and (34) to the ones calculated from the eigenstates of the Anderson Hamiltonian via the local Dyson equation (41). We use a steepest decent method with the following norm: where N max is the number of Matsubara frequencies taken into account. The minimization procedure works as follows: we start from an initial guess of the GSIAM parameters ( lσ , V lσ and W l ), and then knowing the local Green's functions calculated from Eq. (39) and the self-energies calculated from (37) and (38) we calculate the lattice Green's function according to Eqs. (35) and (36). Subsequently, using the Dyson equation (41) we can calculate the Weiss Green's functions G −1 σ,ex (iω n ) and F −1 σ,ex (iω n ). The next step is to fit this result by the parameterization in Eqs. (33) and (34) and thus to find a new set of parameters for the GSIAM. These new parameters serve as input for the next iteration. This procedure is repeated until convergence is reached.

Calculation of the critical temperature
The combination of the mean-field approximation in the higher bands and GDMFT explained so far, however, leads to a problem. Both approximations involve the superfluid order parameter b i . The mean-field approximation for the higher bands implies that the local correlator b † iĉi↑,lĉi↓,l is approximated by b † i ĉ i↑,lĉi↓,l . The GDMFT scheme, on the other hand, involves the approximation to replace the nonlocal correlator b † ib j by b † i b j . This means that b both measures the local phase coherence between bosons and fermions and the non-local bosonic long range order. However, these are two very different quantities which generally cannot be described by a single mean-field order parameter. At zero temperature, this problem is not too severe, because in this case one expects both long-range order and on-site boson-fermion coherence, so that b is large for both reasons. At finite temperature, however, this becomes a real problem, because the bosonic long range order is expected to vanish at temperatures of the order of the bosonic hopping t b . The local boson-fermion coherence, on the other hand, persists for much higher temperatures, since the coupling g is orders of magnitudes larger than t b . Indeed, we find that in the approximation outlined above the full GDMFT calculations are in good agreement with a single site approximation.
In the single-site approach the impurity site couples neither to the fermionic nor to the bosonic bath and long range order cannot be inferred. The critical temperature obtained from this calculation can be identified with the pair breaking temperature T pair , which is much higher than the relevant temperature in experiments. The scheme explained in the previous subsection can therefore not be used to infer the critical temperature for superfluid long range order. In order to do so, we have to modify the approximation and remove the ambiguous nature of the order parameter b . This is made possible by the observation that the term ∆b † in the Hamiltonian merely renormalizes the self-energy of the bosons: in the BEC regime the bosons are in a coherent state and this term is equivalent to a shift of the bosonic chemical potential. This is also clear from the treatment in [35], where terms from higher bands enter the bosonic self-energy. To make this more explicit, we write We can therefore replace the term in the Hamiltonian, by such that terms from the higher bands only renormalize the chemical potential. We remark here, that this might look like an additional approximation. However, one has to keep in mind that the term connecting ∆ to the bosonic creation operator originated from the mean-field approximation in the higher bands. By treating the contribution of higher bands within the bosonic self-energy we are therefore restoring part of the mean-field approximation made in the previous step. Indeed, by using second order perturbation theory in the couplings to higher bands (which is justified if the band energy exceeds the Feshbach coupling), we can also obtain this correction directly as part of the bosonic self-energy, without invoking a mean-field decoupling. This improved approximation for treating higher bands gives for T = 0 similar results as before; in particular the position of the transition to the Mott insulator is in good approximation the same. The superfluid order parameter is smaller, as expected. However, for nonzero temperatures this improved approximation scheme allows for a calculation of the critical temperature for superfluid long range order, which was not possible in the previous approximation.

Results
We study a mixture of potassium atoms ( 40 K) and Feshbach molecules in a threedimensional optical lattice. The on-site harmonic oscillator frequency is chosen to be ω = 2π × 58275Hz, which corresponds to a lattice with wavelength λ = 806nm and Rabi frequency of Ω R = 2π × 1.43GHz. The Feshbach resonance considered here is at B = 202.1G and the width of the resonance is 7.8 G [50]. The difference between the magnetic moments of the closed and open channels of the Feshbach resonance is ∆µ = 16/9µ B , where µ B is the Bohr magneton. The total filling per lattice site in our calculation is n tot = 1.

Zero temperature
First we consider the case of zero temperature. Our calculations for the ground state are summarized in Fig. 2. Deep in the BEC regime only bosonic molecules are present. When the magnetic field is increased, close to resonance the number of fermions is increasing and the number of bosons decreasing. Above the resonance we mainly have fermions and the number of bosons is small. Both fermions and bosons are superfluid. We remark again that here we describe the physics in terms of bare bosons and fermions: in terms of dressed particles as in Ref. [35], these are still molecular bosons and the BEC/BCS crossover takes place when the bosonic self-energy crosses twice the Fermi energy [35]. However, in the case of half-filled fermions, this crossover is intercepted by a first order phase transition to a fermionic Mott insulator state, which happens at a critical value of the magnetic field of B = 249G. Calculations which include only the lowest band of the Bose-Fermi-Hubbard model (as well as with one and two exited bands), yield this transition into the Mott insulator phase already close to the Feshbach resonance at B 205G (results not shown). This implies that to capture the superfluid region 205T B < 249T higher bands which renormalize the bosonic self-energy are crucial.
Another point worth noting is the first-order nature of the transition to the Mottinsulating state. In contrast, if one integrates out the bosonic degree of freedom and describes the Feshbach resonance in terms of an effective, attractive interaction between the fermions within a single channel model for the lowest band, one finds a different scenario. In this case the induced attractive interaction dominates, until it is cancelled by the repulsive background interaction. This means that within the single channel approximation one finds a regime with a normal Fermi-liquid phase in between the BEC/BCS phase and the Mott insulator, which is absent in our phase diagram. This directly indicates that the effect of higher bands is crucial to capture the first-order transition between the superfluid and insulator phase.

Nonzero temperature
Having clarified the ground state phase diagram, we now consider finite temperature. In particular we investigate the critical temperature for the transition to the normal state. Deep in the BEC regime, the critical temperature is constant (T c ≈ 0.21t f ) and completely determined by the properties of the bosons: the bosonic hopping parameter t b , the interbosonic background scattering length a b and the bosonic density. Only very close to resonance the critical temperature suddenly drops (see Fig. 3). This coincides with the magnetic field value for which fermions enter the system. On the BCS side of the resonance, the critical temperature depends on the magnetic field and increases with B (see Fig. 3). This implies that at resonance the critical temperature is minimal. This is in sharp contrast to the situation where no lattice is present, in which case the critical temperature is maximal close to resonance.
This surprising fact can be understood from the behavior of the critical temperature in the single-band attractive Hubbard model [51,52]. In this model, the critical temperature is low both for very large and very small attraction and has a maximum in between. The reason for the low critical temperature at small attraction is the conventional exponential suppression of T c in the BCS regime. For strong attraction the critical temperature decreases again because the fermions start forming bound pairs with a greatly enhanced effective mass. Identifying the resonance position with the case of very large attraction, this explains the low critical temperature at this point. When moving away from resonance, the effective attraction induced by the Feshbach resonance becomes weaker and hence the critical temperature increases again. Far away from resonance one would therefore expect to find a maximum of the critical temperature, after which it decreases again because the BCS regime of weak attraction is entered. However, due to the transition into the Mott insulator phase for the unit filling considered here we cannot see the maximum of the critical temperature. Estimating the induced attractive interaction at the transition point to the Mott insulator by assuming it to be equal to the repulsive background interaction, this indeed gives a value for the induced attractive interaction which is larger (in absolute value) than the position of the maximum in the single-band attractive Hubbard model [51,52].
Our calculation also shows that on both sides of the resonance the ratio c ↑ c ↓ / b † as a function of the temperature for fixed value of the magnetic field is constant. This means that the on-site Bose-Fermi coherence is not affected by the temperature for the low temperatures considered here.

Summary
We have studied ultracold fermionic 40 K atoms in a three-dimensional optical lattice close to a Feshbach resonance. We derived an effective description in terms of a Bose-Fermi Hubbard model, in which the molecular degree of freedom is explicitly present. Our calculations show in agreement with Ref. [35] that the effect of higher bands is crucial for a correct description of the Feshbach physics. We therefore take into account the fermionic occupation of higher bands.
To solve the strongly interacting multi-band problem we decouple the higher bands from the lowest one via a mean-field decoupling and reduce the Hamiltonian to an effective single-band Bose-Fermi Hubbard model which is self-consistently coupled to the higher bands. To solve this resulting model we use GDMFT.
The low temperature physics close to the Feshbach resonance is very rich. Upon changing of the magnetic field, the ratio of fermionic and bosonic densities is changing. Below the resonance the system is mainly occupied by molecular bosons forming a condensate. Close to resonance the number of bosons decreases while the number of fermions increases. The fermions are in the superfluid phase. This resembles the BEC/BCS crossover close to a Feshbach resonance without an optical lattice. In addition, for the unit total filling considered here we found a transition into the fermionic Mott insulating phase when the magnetic field is increased even further. The Mott insulator phase is stabilized by the repulsive fermionic background scattering, which at large magnetic fields overcomes the attractive interaction induced by the Feshbach resonance. The phase transition into the Mott insulator is found to be of first order. We found that higher bands are crucial for a quantitatively correct prediction of the transition point.
We also calculated the critical temperature of the BEC/BCS superfluid phase across the resonance. Below resonance the critical temperature is independent of the magnetic field, until it sharply drops close to resonance. Above the resonance the critical temperature is increasing again, leading to the remarkable result of a minimal critical temperature at resonance.