Ground State Properties of an Asymmetric Hubbard Model for Unbalanced Ultracold Fermionic Quantum Gases
Abstract
In order to describe unbalanced ultracold fermionic quantum gases on optical lattices in a harmonic trap, we investigate an attractive () asymmetric () Hubbard model with a Zeemanlike magnetic field. In view of the model’s spatial inhomogeneity, we focus in this paper on the solution at HartreeFock level. The HartreeFock Hamiltonian is diagonalized with particular emphasis on superfluid phases. For the special case of spinindependent hopping we analytically determine the number of solutions of the resulting selfconsistency equations and the nature of the possible ground states at weak coupling. We present the phase diagram of the homogeneous system and numerical results for unbalanced Fermimixtures obtained within the local density approximation. In particular, we find a fascinating shell structure, involving normal and superfluid phases. For the general case of spindependent hopping we calculate the density of states and the possible superfluid phases in the ground state. In particular, we find a new magnetized superfluid phase.
pacs:
03.75.HhStatic properties of condensates; thermodynamical, statistical, and structural properties and 39.25.+kAtom manipulation (scanning probe microscopy, laser cooling, etc.) (see also 82.37.Gk STM and AFM manipulations of a single molecule in physical chemistry and chemical physics; for atom manipulation in nanofabrication and processing, see 81.16.Ta) and 71.10.FdLattice fermion models (Hubbard model, etc.)1 Introduction
Experimental background
Experiments on unbalanced ultracold fermionic quantum gases have become of significant interest during the last year patridge:pairing ; zwierlein:fermionic . While these experiments were performed on quantum gases in the absence of an optical lattice, it will only be a matter of time before a lattice structure is superimposed on the trap potential. Since the experiments without a lattice have already been theoretically analyzed stoof:deformation ; stoof:sarma , in this paper we do the next step and investigate systems with an additional optical lattice. Such lattice systems can be described with the help of Hubbardtype models zoller:cold ; greiner:mott ; bloch:interfrence ; oosten:quantum ; bloch:nature ; bloch:formation ; hofstetter:hightemp ; maelo:superfluidity . Hence, since we are interested in superfluid phases with a possible imbalance in the pseudospins’ occupation numbers we analyze an appropriate Hubbard model with an attractive interaction () and a Zeemanlike magnetic field. The symmetric and magnetic fieldfree version of this Hamiltonian has been intensively discussed in the seminal paper by P. Noizières and S. SchmittRink noizeres:bf94sd12JE . The pseudospinstates in this Hamiltonian represent different hyperfine states of one and the same atomic species or, alternatively, fermionic atoms with different masses. Accordingly, we also allow for a possible asymmetry in the model, i.e., for a possibly spindependent hopping (). To deal with the spatial inhomogeneity, caused by the presence of the harmonic trap, we focus in this paper on solutions at HartreeFock level; possible extensions and improvements are discussed in section 6. The HartreeFock approximation limits predictions to the region of nottoostrong (attractive) coupling between the fermions. In this paper, therefore, we restrict consideration to the moderate weak coupling regime. Clearly this regime can be realized experimentally, since interactions are tunable in ultracold quantum gases with the help of Feshbach resonances (see greiner:mott ; oosten:quantum ; pethick:bose ).
Model Hamiltonian
A class of models, which is well suited for describing ultracold quantum gases on optical lattices, are models with a local Hubbard interaction. Here we describe only the basic ingredients of the Hamiltionian. Details of the appropriate model and its properties are discussed in section 2.
In ultracold quantum gases on optical lattices, the Hubbard interaction is generated as follows. The periodic potential induced by the coherent laser beam has the form
(1) 
where is proportional to the laser beam’s intensity, is the lattice constant and is the laser beam’s wavelength. In the tight binding limit, , only the local overlap matrix element of the interaction contributes to the Hamiltonian. At sufficiently low temperatures (, where is the band gap between the lowest two Bloch bands) an effective oneband Hubbard model may then be derived from the potential (1). Apart from the onsite interaction term, which is here assumed to be attractive, the Hubbard Hamiltonian contains a nearestneighbor hopping term. The hopping amplitudes will in general be pseudospindependent, , since the pseudospin index labels atoms in different hyperfine states or with different masses (see giamarchi:twocomp ). Furthermore, the particle density in the grand canonical ensemble is tuned with the help of a chemical potential, and the imbalance in the occupation numbers of the pseudospin species is regulated by a Zeemanlike magnetic field.
Structure of the paper
The plan of the paper is as follows. First, in section 2, we describe the Hamiltionian, its symmetry properties and the relevant selfconsistency equations for the particle density and the superfluid order parameter. Then, in sections 3 and 4, we discuss the possible solutions of the selfconsistency relations for symmetric hopping, , and the corresponding numerical results. Our most important result here is the finding of an interesting shell structure for the spatial distribution of the quantum gas in the trap. Results for the general case of asymmetric hopping, , are discussed in section 5. Here we find a new magnetized superfluid phase. We end (in section 6) with a brief summary and an outlook.
2 Properties of the Hamiltonian
The full grandcanonical Hamiltonian , describing (possibly spindependent) hopping of two different fermionic atomic species in a (possibly unbalanced) mixture, reads:
(2) 
The various terms represent the hopping of the fermions, their local Hubbard interaction, the chemical potential and a Zeemanlike magnetic field. Alternatively, the last two terms can of course also be interpreted as a spindependent chemical potential . The summation label in the kinetic energy indicates that all nearest neighbor sites and are summed over, so that every bond is counted twice. Throughout, we assume that the lattice has (hyper)cubical structure, although many symmetry properties are more generally valid for bipartite lattices. In the following, we first discuss some symmetry properties of the Hamiltonian (2), and then we derive the selfconsistency equations within the HartreeFock approximation.
2.1 Symmetry properties
Chemical potential at half filling.
The discrete symmetries of the standard Hubbard model ( and ) are well known from the literature. In this case, the grand canonical Hamiltonian is invariant, e.g., under a spinexchange transformation and, at half filling on bipartite lattices, also under a general particlehole transformation , where equals +1 for sites on the Asublattice and for on the Bsublattice. From these discrete symmetries it can be derived that, at half filling, the chemical potential is given by
(3) 
In the general case (i.e., both and ) these symmetries do not exist, so that at half filling the chemical potential is a more complicated function of the system parameters. However, if either or , equation (3) is still valid: For (but possibly ), the Hamiltonian is not invariant under a spinexchange transformation, but it is invariant under a general particlehole transformation at half filling. Hence, (3) is still valid. If (but possibly ) the Hamiltonian (2) is neither invariant under a general particlehole transformation nor under a spinexchange transformation; nevertheless it is mapped onto itself at half filling by consecutively performing both transformations, so that (3) is also satisfied in this case.
Relation to the repulsive model.
In order to understand the properties of the model (2) at and arbitrary filling, it is often helpful to map the Hamiltonian to a canonically equivalent model at in a magnetic field. With the use of a special particlehole transformation, and , the Hamiltonian (3) with and parameters is mapped onto a formally identical Hamiltonian with and the (in general different) parameters . This is especially useful for , since the new chemical potential is then given by . The main advantage of the mapping is that (for symmetric hopping amplitudes) many properties of the repulsive model are well known from the literature koppe:ground . Furthermore, with the help of a special particlehole transformation, superfluid states in the attractive model are mapped onto antiferromagnetic states with a staggered magnetization in the plane in the repulsive model, while charge density wave states (CDW) are mapped onto antiferromagnetic states with a staggered magnetization in direction.
2.2 Diagonalization of the Hamiltonian at HartreeFock level
In order to analyze superfluid states at weak coupling, we diagonalize the Hamiltonian (2) within the HartreeFock approximation by decoupling the interaction term in the usual manner:
(4) 
where and are assumed to be translationally invariant thermal averages. The resulting meanfield Hamiltonian can as usual be diagonalized by a Bogoliubovtransformation. One finds
(5)  
where is the number operator of the Bogoliubov quasiparticles. Moreover, we defined and is the number of lattice sites. The spinaveraged nonsuperfluid particle energies and the quasiparticle energies are defined as:
(6a)  
(6b) 
where is alternatively interpreted as or in (6b). The three quantities , and can also be determined from their definitions as thermal averages. Suppressing all details, we just mention as a result that, at zero temperature, the following three selfconsistency equations have to be satisfied in the thermodynamic limit:
(7a)  
(7b)  
(7c) 
where is the dimensional density of states of the interactionfree Hubbard model. Since particlehole symmetry is lost if both and (see section 2.1), we consider in the following only special cases with either or .
3 Properties of the selfconsistency equations for symmetric hopping ()
Since the coupled equations (7a)  (7c) may have more than one solution, we show analytically that the number of solutions at fixed is between one and three. If more than one solution is found, the one with the lowest grand potential is thermodynamically stable.
3.1 Particle numbers at fixed
Integral equations as a mapping.
In order to determine the number of solutions of the coupled equations (7a)  (7c), we first consider Equations (7a) and (7b) separately at fixed order parameter . We show that, at fixed , there is a unique solution of (7a) and (7b) by using Banach’s fixed point theorem. As vectors we use , where denotes different points in the parameter space, i.e., different values of . The integrations of (7a) and (7b) can be interpreted as a mapping from the parameter space onto itself. Furthermore, in order to apply Banach’s theorem, we consider the maximum metric . To prove uniqueness it is sufficient to show that a number exists, so that for all :
(8) 
holds, where is a small variation in parameter space. Hence, it suffices to show that the mapping is a contraction. The solution of (7a) and (7b) is then given by the unique fixed point of the equation .
Convergence at small interaction strength.
For small variations it can be shown that the mapping is a contraction for all points in parameter space. In order to prove that Equation (8) is indeed satisfied at every point , we use a Taylorexpansion of in the parameter and some estimates, whose validity depends upon the relative interaction being weak. In doing so, we must distinguish two cases, namely that the functions and change sign at the same value or at different values of . In the first case the Taylorexpansion yields as a criterion for the validity of (8) with :
(9) 
and in the second case it yields:
(10) 
Here is defined as the energy, for which the noninteracting density of states is maximal. Since the condition (10) is more restrictive than (9), Equation (10) guarantees that the mapping is a contraction in the whole parameter space at sufficiently weak coupling. Note that this proof does not work if the density of states diverges at a van Hove singularity. For a simple cubic lattice, Equation (10) yields for to be a contraction.
3.2 Variation of the order parameter
In order to solve the coupled Equations (7a)  (7c) we also have to discuss the properties of (7c) explicitly. Equation (7c) is trivially solved by , which corresponds to a nonsuperfluid phase. For and symmetric hopping amplitudes it is useful to distinguish between , where and Equation (7c) reduces to
(11) 
and , where (7c) may be rewritten as:
(12) 
Since, from section 3.1, we know that there is only one solution for and , we are able to do a onedimensional parameter scan through all possible values in order to find the solutions of (7c), using Equations (12) for and (11) for . Since the integral on the right in (11) decreases with increasing if , there is at most one solution in this interval. Similarly, there is also at most one solution for , since the integral on the right in (12) in general increases with in this case, so that the residual is monotonic on both sides of . As a consequence, there exist at most three solutions for , namely the trivial solution and possibly two nontrivial solutions on either side of . Note that the modulus of is rigorously bounded from above by .
3.3 Possible phases occurring at HartreeFock level
Since, for , the magnetization is rigorously zero () and there are at most three nonequivalent solutions of the coupled Equations (7a)  (7c), we find that (for ) there can in principle be up to three competing phases:

a nonsuperfluid magnetized phase ()

a superfluid magnetized phase ()

a superfluid nonmagnetized phase () .
For it is well known that the superfluid nonmagnetized solution has the lowest grand potential. Since this is not necessarily true for , phase transitions will generally be of first order, also at HartreeFock level. We should also mention that there is no thermodynamically stable superfluid magnetized phase at and , since that type of solution is always a local maximum of the grand potential which can be shown with elementary mathematical methods.
4 Numerical results for symmetric hopping
Numerical method.
With the knowledge of the properties of the selfconsistency equations (7a)  (7c), established in section 3, it is possible to do accurate numerical calculations for the phases occurring in the ground state. Specifically, we first determine all possible solutions of Equations (7a)  (7c), and then identify the thermodynamically stable solution by comparing the respective numerical values of the grand potentials. The various solutions of the selfconsistency equations are determined by performing a scan over values between and in order to find the roots of (7c). The parameters and are varied selfconsistently at each step by successive approximation; this method is guaranteed to work on account of Banach’s fixed point theorem (see section 3.1). The phase diagram of the homogeneous system is presented in Figure 1. The parameters , and are presented as a function of and at fixed and . At the system is always superfluid and the parameters , and do not depend on at fixed until a critical magnetic field , where superfluidity breks down and magnetization occurs. Note that this first order phase transition generally takes place at . The superfluid order parameter is of course maximal at half filling ().
Numerical results within the LDA.
In order to formulate predictions for experiments on ultracold quantum gases we also have to treat the parabolic potential arising from the magnetooptical trap. In Hubbard model language this potential may be described as:
(13) 
where is a real, symmetric, positive definite matrix. We assume the trapping potential to be spinindependent. The trapping potential is treated within the local density approximation (LDA), meaning that local quantities at site are determined from the HartreeFock selfconsistency equations corresponding to the effective chemical potential ; the global chemical potential is then tuned so as to reproduce the desired total number of particles in the trap.
The numerical results, obtained from the HartreeFockLDA approximation, are presented in Figures 2  4, which show the spatial distribution of the occupation number , the spatial distribution of polarization and the spatial distribution of the superfluid order parameter’s absolute value , respectively. In accordance with the ellipsoidal geometry of the trap, we find a fascinating and intuitively extremely plausible shell structure for the spatial distribution of the various densities and the local order parameter. In Figure 2 there is a normal phase in the center of the trap surrounded by a superfluid nonmagnetized shell and a normal phase in the outer region of the trap. Apparently, the normal phase in the center of the trap is stabilized by the Zeemanlike magnetic field and reacts very sensitively to changes in this field. This can be seen from Figure 3, where the (absolute value of the) magnetic field is slightly lower than in Figure 2; also the global particle number (chemical potential ) is slightly different. These changes suffice to cause the disappearance of the nonsuperfluid core. By enhancing the stiffness of the trap potential in the  and directions or slightly decreasing the interaction strength , it is possible to suppress superfluidity alltogether: There is no superfluid phase in Figure 4.
Experimental consequences.
In experiments, the superfluid order parameter cannot be observed directly. On the basis of the results, presented above, it seems that little information of use concerning superfluidity can be deduced from the total particle number . However, the other observable considered above, i.e., the local polarization , might in fact be quite helpful in identifying the shellstructure in the trap as an experimental signature for superfluidity in unbalanced Fermimixtures. The observable can be resolved spatially with the use of insitu imaging patridge:pairing ; zwierlein:fermionic .
Comparison with optical latticefree systems.
Systems without optical lattices have been intensely discussed in the literature, both from a theoretical stoof:deformation ; stoof:sarma as well as from an experimental point of view patridge:pairing ; zwierlein:fermionic . Important similarities between systems with and without optical lattices are, first, the occurrence of superfluidity below a critical temperature and secondly, the observation of a shell structure containing normal and superfluid rings. Moreover, there are also several differences between these classes of systems, which will now be highlightned.
In optical latticefree systems it is found, in contrast to our results, that, if superfluidity occurs, the atoms in the center of the trap are always superfluid. This difference arises from the fact that, in the oneband situation discussed in this paper, the quasimomentum is bounded from above by , whereas in optical latticefree systems there is no upper bound for the components of the momentum vector. We have shown that superfluid magnetized phases occurring in optical latticefree systems as demonstrated, e.g., by Parish et al. parish:finite:6Ge ; parish:polarized:v8E do not occur in systems, considered in this work (see section 3.3).
Deformations of the shell structure in the optical latticefree system, discussed by H. Stoof in stoof:deformation , seem to arise from a surface tension between a superfluid core and a normal state, which is not discussed in our work, where we consider systems described by discrete lattice vectors.
Moreover, in recent publications sheehy:b84F3oa FFLO phases in latticefree systems are discussed. At present it seems not possible to combine this formalism, used for spatially homogeneous (or periodic) potentials, with the LDAapproach, used in this paper for spatially inhomogeneous systems, since the FFLOground state is spatially inhomogeneous by itself. For this reason, and since FFLO phases were not yet observed experimentally, we disregard FFLO phases in this paper.
Trapping potentialfree systems.
In recent literature iskin:bgu93dbE , trapping potentialfree systems have been discussed in the weak as well as in the strong coupling regime. In the weak coupling regime we obtain the phase diagram in the grand canonical ensemble, presented in Figure 1. Since we find a first order phase transition between a superfluid nonmagnetized phase and a normal phase, we can formulate predictions for systems with fixed particle numbers and instead of fixed conjugate variables and : For nonmagnetized systems the system is completely superfluid. For a population imbalance, , the following case differentiation has to be done:

If in the normal phase (i.e. assuming ) the particle numbers and can be realized (within the grand canonical ensemble) by choosing definite values and for the chemical potential and the magnetic field, respectively, then this normal phase is also the ground state.

If in the normal phase there exists no such combination of ()values, then phase separation occurs in the ground state. The superfluid and the normal fraction have to be determined by a Maxwell construction.
As mentioned in section 4 a superfluid magnetized phase does not exist, at least not in the weak coupling limit.
5 Results for asymmetric hopping
In this section we discuss the properties of the HartreeFock Hamiltonian for superfluid states with asymmetric hopping (). For clarity and convenience, we focus on translationally invariant solutions () for .
5.1 Charge density wave states and the repulsive model
As emphasized before, we are interested in (and, hence, focussing on) superfluid phases in the negative Hubbard model, disregarding possible competing phases, which seem to be irrelevant for experiments on ultracold quantum gases (see for example hofstetter:hightemp ; maelo:superfluidity ; orso:sound ; yamada:novel ; shly:superfluid ; reddy:asymmetric ). The issue of possible competing phases becomes particularly relevant for and halffilling, since in this case charge density wave (CDW) states are known to be degenerate with swave superfluidity if . This is immediately clear from a special particlehole transformation from the negative to the positive Hubbard model, both at halffilling and , since under this transformation swave and CDW states are mapped onto the various directions of the staggered magnetization, which are energetically degenerate due to the Hubbard model’s symmetry in the spin sector. Away from halffilling, however, swave superfluidity is stabilized with respect to the CDW phase.
For moderately asymmetric hopping () the situation is similar: Exactly at halffilling the CDW phase is energetically even somewhat lower than swave superfluidity, so that one expects the CDW state to dominate in the phase diagram in a narrow but finite strip around halffilling. We conclude that the negative Hubbard model, considered here, predicts phases near halffilling, which at present seem to be of little experimental interest. For this reason we concentrate in the following on band fillings, which deviate significantly from halffilling, so that the superfluid phases, considered here, can safely be assumed to be stable koppe:ground .
5.2 HartreeFock density of states
The absence of symmetry in the hopping amplitudes () clearly also influences the (HartreeFock) density of states, which is defined as:
(14) 
The asymmetric hopping now causes additional singularities (jumps) in the spindependent densities of states. If we assume , the following changes occur:

The inverse squareroot divergence of the density of states at the superfluid gap in is replaced by a discontinuity (jump).

The DOS retains its inverse squareroot divergence but now displays two additional discontinuities within the band (one below and one above the gap).

The superfluid gap in for is smaller than the corresponding value for , while the gap remains unchanged in .
The normalization of is obviously not affected.
The densities of states for both spin species, as determined from Equations (14) and (6b), are shown in Figures 5 and 6 for the parameter values and a threedimensional simple cubic lattice. These parameter values are only used to illustrate the mentioned effects on the DOS; they do not necessarily represent a solution of the selfconsistency equations (7a)  (7c). The spindependent size of the gap, as seen in Figures 5 and 6, can of course also be observed experimentally.
5.3 Phases in the ground state
In this section we present the results of the numerical solutions of (7a)  (7c) for translationally invariant simple cubic lattices at fixed particle numbers and . The proof presented in section 3.1 is also valid for asymmetric hopping (), albeit for smaller values of the dimensionless interaction parameter than for systems with symmetric hopping. This is due to the fact that the quasiparticle energy function in (6b) is not a monotonic function in for the case of asymmetric hopping. Numerical calculations show that, for and , Equation (7c) has at most one solution with . If two solutions (one with and one with ) exist, we find that the superfluid phase always minimizes the free energy.
In Figures 7 and 8 we illustrate the behavior of the chemical potential , the magnetization and the absolute value of the superfluid order parameter as a function of the interaction strength at different fixed occupation numbers . We find that superfluid solutions exist only for attractive interaction strengths above a certain minimal value, . In contrast to systems with spinindependent hopping, we now also find thermodynamically stable magnetized superfluid solutions. Note that the trasitions occurring here are in general of second order.
6 Summary and Discussion
To summarize, we analyzed a generalized attractive Hubbard model with (possibly spindependent) hopping amplitudes and a Zeemanlike magnetic field, which is relevant for ultracold quantum gases. In view of the spatial inhomogeneity due to the presence of a trap, we studied the model at HartreeFock level. For the special case of symmetric hopping we analytically demonstrated the uniqueness of the solution of Equations (7a) and (7b) at fixed order parameter . Hence we are certain to have detected all solutions of (7a)  (7c). We found that phase transitions at are generally of first order. Within the local density approximation we showed that phase separation with an interesting shell structure occur if a parabolic potential is added to the Hamiltonian (2). For general spindependent hopping and we calculated the (spindependent) density of states. Interestingly, we found that, sufficiently far away from half filling, a superfluid magnetized phase occurs, which was absent in the phase diagram for spinindependent hopping.
Finally, we discuss possible extensions of our results. In this paper, we restricted consideration to weakcoupling solutions at HartreeFockLDA level. For our purposes, this restriction is meaningful, since the HartreeFockLDA approximation is expected to yield qualitatively correct results in the weak coupling regime, provided spatial gradients are not too large; moreover, any other method would be virtually impracticable in a spatially inhomogeneous system with a magnetooptical trap. Nevertheless, with an eye on the future, it seems worthwhile to review possible alternatives.
First of all, it would probably be possible to drop the LDAapproximation and calculate the spatially inhomogeneous HartreeFock solution for an optical lattice in a parabolical trap. However, even for a twodimensional lattice without spinasymmetries ( and ), this is known to be computationally fairly expensive for realistic lattice sizes (with typically particles) freimuth:super .
Secondly, it would of course be important to go beyond the HartreeFock approximation and include correlation effects. This could be done in principle by taking into account the dominant quantum fluctuations, which are described by second order Feynman diagrams in a perturbation expansion. In principle, such an expansion could be formulated also for spatially inhomogeneous systems with symmetry breaking, but the computational expense of solving these equations selfconsistently for a threedimensional system of reasonable size would be tremendous. Here we just recall that it is known for translationally invariant systems dongen:futr93Ea ; dongen:bgE58saB ; dongen:gt8Yaqv4 in three or more dimensions, that the leading effect of quantum fluctuations is to renormalize HartreeFock results and reduce all energy scales (critical temperatures, gaps) by factors of typically three to four. For the calculations of the present paper this suggests that inclusion of second order diagrams would quantitatively but not qualitatively change the HartreeFock results.
Clearly it would, for the description of nonperturbative and strongcoupling effects, also be desirable to apply simulation methods (in particular Quantum Monte Carlo techniques) to the Hubbard models studied in this paper. However, a simulation of a threedimensional system of physically interesting size at sufficiently low temperatures seems at present unfeasible. It also seems unlikely that approximations like Dynamical Mean Field Theory (DMFT) or its cluster extensions, which are known to be extremely helpful in translationally invariant systems in and , are applicable to spatially inhomogeneous systems in a trap at ultralow temperatures, where superfluidity occurs. However, away from the symmetrybroken lowtemperature phase, the DMFT has been demonstrated to be well applicable, e.g., to the description of spatially inhomogeneous Mott metalinsulator transitions helmes:bf5sc8H .
For the time being, we conclude that we discovered several fascinating new phenomena in Hubbardtype models for asymmetric fermionic ultracold quantum gases. Specifically we found phase separation with an interesting shell structure for systems in a trap and also (in a translationally invariant model) a new superfluid magnetized phase. It would be important if the calculations presented here could be extended beyond the weakcoupling regime, using different methods, along the lines pointed out above. We hope and expect that our results can soon be compared to experiment.
References
 (1) G. B. Partridge et al., Science 311, 503 (2006)
 (2) M. Zwierlein et al., Science 311, 492 (2006)
 (3) G. B. Partridge et al., Phys. Rev. Lett. 97, 190407 (2006)
 (4) K. B. Gubbels et al., Phys. Rev. Lett. 97, 210402 (2006)
 (5) D. Jaksch et al., Phys. Rev. Lett. 81, 3108 (1998)
 (6) M. Greiner, Ph.D. thesis, LMU München (2003)
 (7) F. Gerbier et al., Phys. Rev. A 72, 053606 (2005)
 (8) D. van Oosten, Ph.D. thesis, Universiteit Utrecht, Utrecht (2004)
 (9) I. Bloch, Nature Physics 1, 23 (2005)
 (10) S. Fölling et al., Phys. Rev. Lett. 97, 060403 (2006)
 (11) W. Hofstetter et al., Phys. Rev. Lett. 89, 220407 (2002)
 (12) M. Iskin and C. A. R. Sá de Melo, Phys. Rev. B 72, 224513 (2005)
 (13) P. Noizières and S. SchmittRink, J. Low Temp. Phys. 59, 195 (1985)
 (14) C. Pethick, H. Smith, BoseEinstein Condensation in Dilute Gases (Cambridge University Press, Cambridge, 2002)
 (15) M. A. Cazalilla et al., Phys. Rev. Lett. 95, 226402 (2005)
 (16) K. Dichtel et al., Z. Physik 246, 248 (1971)
 (17) Meera M. Parish et al., Nature Physics 3, 124 (2007)
 (18) Meera M. Parish et al., Phys. Rev. Lett. 98, 160402 (2007)
 (19) D. E. Sheehy and L. Radzihovsky, Phys. Rev. Lett. 96, 060401 (2006)
 (20) M. Iskin and C. A. R. Sá de Melo, Phys. Rev. Lett. 99, 080403 (2007)
 (21) L. P. Pitaevskii et al., Phys. Rev. A 71, 053602 (2005)
 (22) M. Machida et al., Phys. Rev. Lett. 93, 200402 (2004)
 (23) G. Orso and G. V. Shlyapnikov, Phys. Rev. Lett. 95, 260402 (2005)
 (24) J. Carlson and S. Reddy, Phys. Rev. Lett. 95, 060401 (2005)
 (25) F. Freimuth (2005)
 (26) P. G. J. van Dongen, Phys. Rev. Lett. 67, 757 (1991)
 (27) P. G. J. van Dongen, Phys. Rev. Lett. 74, 182 (1995)
 (28) T. Schauerte and P. G. J. van Dongen, Phys. Rev. B 65, 081105 (2002)
 (29) R. W. Helmes, T. A. Costi and A. Rosch, arXiv:condmat/07091669v1 (2007)