Ground State Properties of an Asymmetric Hubbard Model for Unbalanced Ultracold Fermionic Quantum Gases

# Ground State Properties of an Asymmetric Hubbard Model for Unbalanced Ultracold Fermionic Quantum Gases

T. Gottwald 1 KOMET 337, Institut für Physik, Johannes Gutenberg-Universität, Staudingerweg 7, 55099 Mainz, Germany
P. G. J. van Dongen 1 KOMET 337, Institut für Physik, Johannes Gutenberg-Universität, Staudingerweg 7, 55099 Mainz, Germany
1 KOMET 337, Institut für Physik, Johannes Gutenberg-Universität, Staudingerweg 7, 55099 Mainz, Germany
Received: date / Revised version: date
###### 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 Zeeman-like magnetic field. In view of the model’s spatial inhomogeneity, we focus in this paper on the solution at Hartree-Fock level. The Hartree-Fock Hamiltonian is diagonalized with particular emphasis on superfluid phases. For the special case of spin-independent hopping we analytically determine the number of solutions of the resulting self-consistency 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 Fermi-mixtures obtained within the local density approximation. In particular, we find a fascinating shell structure, involving normal and superfluid phases. For the general case of spin-dependent 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 Hubbard-type 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 Zeeman-like magnetic field. The symmetric and magnetic field-free version of this Hamiltonian has been intensively discussed in the seminal paper by P. Noizières and S. Schmitt-Rink noizeres:bf94sd12JE . The pseudospin-states 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 spin-dependent hopping (). To deal with the spatial inhomogeneity, caused by the presence of the harmonic trap, we focus in this paper on solutions at Hartree-Fock level; possible extensions and improvements are discussed in section 6. The Hartree-Fock approximation limits predictions to the region of not-too-strong (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

 V(x)=−V0dd∑i=1cos2(πxia), (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 one-band Hubbard model may then be derived from the potential (1). Apart from the on-site interaction term, which is here assumed to be attractive, the Hubbard Hamiltonian contains a nearest-neighbor hopping term. The hopping amplitudes will in general be pseudospin-dependent, , 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 Zeeman-like 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 self-consistency equations for the particle density and the superfluid order parameter. Then, in sections 3 and 4, we discuss the possible solutions of the self-consistency 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 grand-canonical Hamiltonian , describing (possibly spin-dependent) hopping of two different fermionic atomic species in a (possibly unbalanced) mixture, reads:

 K=−∑(ij)σtσic†iσcjσ+U∑ini↑ni↓−μ∑iσniσ+B∑iσσniσ. (2)

The various terms represent the hopping of the fermions, their local Hubbard interaction, the chemical potential and a Zeeman-like magnetic field. Alternatively, the last two terms can of course also be interpreted as a spin-dependent 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 self-consistency equations within the Hartree-Fock 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 spin-exchange transformation and, at half filling on bipartite lattices, also under a general particle-hole transformation , where equals +1 for sites on the A-sublattice and for on the B-sublattice. From these discrete symmetries it can be derived that, at half filling, the chemical potential is given by

 μ=U2. (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 spin-exchange transformation, but it is invariant under a general particle-hole transformation at half filling. Hence, (3) is still valid. If (but possibly ) the Hamiltonian (2) is neither invariant under a general particle-hole transformation nor under a spin-exchange 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 U 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 particle-hole 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 particle-hole 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 Hartree-Fock level

In order to analyze superfluid states at weak coupling, we diagonalize the Hamiltonian (2) within the Hartree-Fock approximation by decoupling the interaction term in the usual manner:

 ni↑ini↓i→n↑ini↓i+n↓ini↑i−n↑in↓i+Δc†i↑ic†i↓i+Δ∗ci↓ici↑i−|Δ|2, (4)

where and are assumed to be translationally invariant thermal averages. The resulting mean-field Hamiltonian can as usual be diagonalized by a Bogoliubov-transformation. One finds

 KHF = ∑k{E(εki)−sign[E(εki)]√[E(εki)]2+U2|Δ|2} (5)

where is the number operator of the Bogoliubov quasiparticles. Moreover, we defined and is the number of lattice sites. The spin-averaged non-superfluid particle energies and the quasiparticle energies are defined as:

 E(ε)= −(t↑+t↓)ε−μ+U(n↑+n↓)2 (6a) Eσi(ε)=sign(E(ε))√E2(ε)+U2|Δ|2++σ((t↓i−t↑i)ε+B+U(n↓−n↑)2), (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 self-consistency equations have to be satisfied in the thermodynamic limit:

 n↑= 12⎡⎢ ⎢⎣1+∫∞−∞dενd(ε)∣∣E(ε)∣∣sign(−E↑(ε))√E2(ε)+U2|Δ|2⎤⎥ ⎥⎦ (7a) n↓= 12⎡⎢ ⎢⎣1+∫∞−∞dενd(ε)∣∣E(ε)∣∣sign(−E↓(ε))√E2(ε)+U2|Δ|2⎤⎥ ⎥⎦ (7b) Δ=UΔ2∫∞−∞dενd(ε)sign(E(ε))√E2(ε)+U2|Δ|2××[Θ(−E↑(ε))+Θ(−E↓(ε))−1], (7c)

where is the -dimensional density of states of the interaction-free Hubbard model. Since particle-hole 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 self-consistency equations for symmetric hopping (\boldmatht↑=t↓≡t)

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 :

 ||I(n+δn)−In||∞≤q||δn||∞ (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 Taylor-expansion 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 Taylor-expansion yields as a criterion for the validity of (8) with :

 |U|t<2νd(εmax) (9)

and in the second case it yields:

 |U|t<12νd(εmax). (10)

Here is defined as the energy, for which the non-interacting 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 non-superfluid phase. For and symmetric hopping amplitudes it is useful to distinguish between , where and Equation (7c) reduces to

 0=|U|2∫∞−∞dενd(ε)1√E2(ε)+U2|Δ|2−1, (11)

and , where (7c) may be rewritten as:

 0=|U|2∫∞−∞dενd(ε)Θ(E↑(ε)⋅E↓(ε))√E2(ε)+U2|Δ|2−1. (12)

Since, from section 3.1, we know that there is only one solution for and , we are able to do a one-dimensional 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 non-trivial solutions on either side of . Note that the modulus of is rigorously bounded from above by .

### 3.3 Possible phases occurring at Hartree-Fock level

Since, for , the magnetization is rigorously zero () and there are at most three non-equivalent solutions of the coupled Equations (7a) - (7c), we find that (for ) there can in principle be up to three competing phases:

• a non-superfluid magnetized phase ()

• a superfluid magnetized phase ()

• a superfluid non-magnetized phase ()  .

For it is well known that the superfluid non-magnetized solution has the lowest grand potential. Since this is not necessarily true for , phase transitions will generally be of first order, also at Hartree-Fock 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 self-consistency 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 self-consistently 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 magneto-optical trap. In Hubbard model language this potential may be described as:

 HTrap=∑iσ(i⋅Ω⋅i)niσi, (13)

where is a real, symmetric, positive definite matrix. We assume the trapping potential to be spin-independent. The trapping potential is treated within the local density approximation (LDA), meaning that local quantities at site are determined from the Hartree-Fock 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 Hartree-Fock-LDA 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 non-magnetized 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 Zeeman-like 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 non-superfluid 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 shell-structure in the trap as an experimental signature for superfluidity in unbalanced Fermi-mixtures. The observable can be resolved spatially with the use of in-situ imaging patridge:pairing ; zwierlein:fermionic .

#### Comparison with optical lattice-free 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 lattice-free 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 one-band situation discussed in this paper, the quasimomentum is bounded from above by , whereas in optical lattice-free systems there is no upper bound for the components of the momentum vector. We have shown that superfluid magnetized phases occurring in optical lattice-free 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 lattice-free 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 lattice-free systems are discussed. At present it seems not possible to combine this formalism, used for spatially homogeneous (or periodic) potentials, with the LDA-approach, used in this paper for spatially inhomogeneous systems, since the FFLO-ground 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 potential-free systems.

In recent literature iskin:bgu93dbE , trapping potential-free 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 non-magnetized phase and a normal phase, we can formulate predictions for systems with fixed particle numbers and instead of fixed conjugate variables and : For non-magnetized systems the system is completely superfluid. For a population imbalance, , the following case differentiation has to be done:

1. 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.

2. 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 Hartree-Fock 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 U 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 half-filling, since in this case charge density wave (CDW) states are known to be degenerate with s-wave superfluidity if . This is immediately clear from a special particle-hole transformation from the negative- to the positive- Hubbard model, both at half-filling and , since under this transformation s-wave 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 half-filling, however, s-wave superfluidity is stabilized with respect to the CDW phase.

For moderately asymmetric hopping () the situation is similar: Exactly at half-filling the CDW phase is energetically even somewhat lower than s-wave superfluidity, so that one expects the CDW state to dominate in the phase diagram in a narrow but finite strip around half-filling. We conclude that the negative- Hubbard model, considered here, predicts phases near half-filling, 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 half-filling, so that the superfluid phases, considered here, can safely be assumed to be stable koppe:ground .

### 5.2 Hartree-Fock density of states

The absence of symmetry in the hopping amplitudes () clearly also influences the (Hartree-Fock) density of states, which is defined as:

 νHFd,σ(E)=∫∞−∞dενd(ε)δ(E−Eσ(ε)). (14)

The asymmetric hopping now causes additional singularities (jumps) in the spin-dependent densities of states. If we assume , the following changes occur:

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

• The DOS retains its inverse square-root 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 three-dimensional 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 self-consistency equations (7a) - (7c). The spin-dependent 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 spin-independent 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 spin-dependent) hopping amplitudes and a Zeeman-like 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 Hartree-Fock 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 spin-dependent hopping and we calculated the (spin-dependent) 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 spin-independent hopping.

Finally, we discuss possible extensions of our results. In this paper, we restricted consideration to weak-coupling solutions at Hartree-Fock-LDA level. For our purposes, this restriction is meaningful, since the Hartree-Fock-LDA 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 magneto-optical 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 LDA-approximation and calculate the spatially inhomogeneous Hartree-Fock solution for an optical lattice in a parabolical trap. However, even for a two-dimensional lattice without spin-asymmetries ( 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 Hartree-Fock 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 self-consistently for a three-dimensional 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 Hartree-Fock 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 Hartree-Fock results.

Clearly it would, for the description of non-perturbative and strong-coupling 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 three-dimensional 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 symmetry-broken low-temperature phase, the DMFT has been demonstrated to be well applicable, e.g., to the description of spatially inhomogeneous Mott metal-insulator transitions helmes:bf5sc8H .

For the time being, we conclude that we discovered several fascinating new phenomena in Hubbard-type 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 weak-coupling 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. Schmitt-Rink, J. Low Temp. Phys. 59, 195 (1985)
• (14) C. Pethick, H. Smith, Bose-Einstein 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:cond-mat/07091669v1 (2007)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters