Quantum Simulation of the Hubbard Model: The Attractive Route

We study the conditions under which, using a canonical transformation, the phases sought after for the repulsive Hubbard model, namely a Mott insulator in the paramagnetic and anti-ferromagnetic phases, and a putative d-wave superfluid can be deduced from observations in an optical lattice loaded with a spin-imbalanced ultra-cold Fermi gas with attractive interactions, thus realizing the attractive Hubbard model. We show that the Mott insulator and antiferromagnetic phase of the repulsive Hubbard model are in fact more easy to observe as a paired, and superfluid phase respectively, in the attractive Hubbard model. The putative d-wave superfluid phase of the repulsive Hubbard model doped away from half-filling is related to a d-wave antiferromagnetic phase for the attractive Hubbard model. We discuss the advantages of this approach to 'quantum simulate' the Hubbard model in an optical lattice over the approach that attempts to directly simulate the doped Hubbard model in the repulsive regime. We also point out a number of technical difficulties of the proposed approach and, in some cases, suggest possible solutions.


I. INTRODUCTION
Understanding the phase diagram of the twodimensional (2D) single-band Hubbard model is considered by many as the 'Holy Grail' of the theory of strongly correlated systems. In the most interesting regime, this model describes a system of spin-1 2 (i.e. two species of) fermions hopping on a 2D square lattice with repulsive (on-site) interactions, and average lattice filling less than one fermion per site. This model has been proposed as the minimal model that explains the observation of dwave superconductivity with fairly high critical temperature in the doped cuprate materials [1] (for a review on doped Mott insulator, see [2]). At half filling, where only one particle per site is allowed, it is by now rather well established that the model is a Mott insulator, which at low temperatures (below a characteristic scale, the Néel temperature, T Néel ) orders anti-ferromagnetically. Away from half filling, the nature of the ground state is a subject of heavy debate. One of the most challenging open issues is whether the Hubbard model on a 2D square lattice would support a d-wave superconducting phase at a (relatively) high temperature. The fact that the Hubbard model can support such an instability has been theoretically proven in double-chain systems coupled by hopping (known as two-leg Hubbard ladders) [3,4,5,6]. However, whether this result extends to the 2D model consisting of an infinite number of coupled chains is still extremely controversial [2,7,8]. At present, neither analytical nor numerical studies are able to settle the issue.
Due to the spectacular advances in the optical manipulation of ultra-cold atomic gases, one very promising route for studying the low temperature phases of the Hubbard model has opened up recently [9]. Indeed, ultra-cold Fermi gases loaded into an optical lattice can be regarded as almost ideal quantum simulators of the Hubbard model, where independent control of the hopping amplitude, t, and the on-site interaction energy, U , are both experimentally available. Exploiting this fact, the Mott insulating phase of the Hubbard model has recently been demonstrated in a 3D cubic optical lattice (where the center is at half-filling) by several experimental groups [10,11]. Many other groups are currently engaged in similar experimental endeavors [12], with the main focus on realizing the repulsive Hubbard model on a 2D square lattice away from half-filling, namely the regime where d-wave superfluidity is speculated to exist. However, one of the main problems that lie ahead in this program has to do with the currently accessible temperatures for the Fermi gases in optical lattices. At present, these temperatures (of the order of a few tenths of the Fermi energy of a non-interacting gas of similar average density) still largely exceed the Néel temperature (T Néel ), thereby washing out any anti-ferromagnetic order in the half-filled system [13].
However, there are other problems with the present approach to simulate the repulsive Hubbard model, which seem not to have received so much attention thus far. One of the most remarkable ones is the difficulty of doping away from half-filling a Fermi gas loaded in an optical lattice. In this case, the situation is very different from doping in solids, mainly because of two reasons: i) The existence of an overall harmonic trapping potential superimposed on the optical lattice potential that makes the system inhomogeneous and tends to favor the maximum site occupancy (i.e. two fermions per site) near the center of the trap. As recently demonstrated experimentally [10,11], for small number of particles this arXiv:0812.4422v2 [cond-mat.str-el] 20 Feb 2009 tendency can be balanced by the onsite repulsion energy, which yields a Mott insulator near the center surrounded by a 'metallic' region where the density of holes is nonuniform. ii) Although the number of available lattice sites can be controlled with accuracy, the total number of atoms is still hard to measure accurately and it is also subject to variations from shot to shot that are inherent to the preparation process [14].
Another source of problems has to do with the need to independently control the on-site interaction, U , and the hopping amplitude, t. As the Néel temperature (T Néel ) (below which the system orders anti-ferromagnetically, and, upon doping, the putative d-wave superfluid may appear), scales as T Néel ∼ t 2 /U , and thus rapidly decreases if the ratio t/U is made very small by increasing the optical lattice depth, it is desirable to have independent control of both t and U . In order to achieve this, the s-wave scattering length, a s (U ∝ a s , roughly speaking) that characterizes the strength and sign of the atom-atom interaction, must be tuned towards a Feshbach resonance where a s → ±∞. Since the current interest is in realizing a Hubbard model with repulsive interactions, the side of the Feshbach resonance where the atom-atom interaction is repulsive (i.e. a s > 0) must be used (see, however, Sect. III for further remarks). On this side of the resonance, there is a weakly bound molecular bound state [15], with which the atoms in the continuum (that is, in the lowest and highest Bloch bands, when loaded in a lattice) have a sizable overlap near the resonance. Thus, at sufficiently low temperatures, Feshbach molecules form resulting from three atom collisions [15]. A collision of one of these molecules with a third atom can cause the molecule to make a transition into a more bound molecular vibrational state. The released energy is taken away by the colliding atom, which therefore causes undesirable heating of the system. Also, the presence of these molecules is not accounted for by the single-band Hubbard model, which is the goal of the quantum simulation. Furthermore, as the scattering length increases, U also increases and become of the order of the separation between Bloch bands, thus leading to the break-down of the single band approximation [16,17,18].
In this article, we propose to explore a different route to simulate the Hubbard model in a regime where the onsite interaction is attractive. Theoretically, the attractive and repulsive regimes are related by a transformation that is well known in the literature of the Hubbard model and it is, for completeness, reviewed in Sect. II. More recently, in the context of cold atomic gases, this transformation has been used by Moreo and Scalapino [19], who pointed out the connection between the Fulde-Ferrel-Larkin-Ovchinikov state in the attractive Hubbard model and a state with stripes in the repulsive model. These authors also briefly considered the relationship between the d-wave superfluid order parameter for the repulsive model and a d-wave anti-ferromagnetic order in the attractive case. For the one-dimensional Hubbard model (note, however, that there is no d-wave superfluid phase in this case), the transformation was also used in Ref. [20] in an analysis of the noise correlations of the attractive Hubbard model with spin imbalance. Moreover, the physics of the attractive Hubbard model has also attracted much interest by itself [21], especially in recent times and in connection with cold atomic gases and the physics of the BEC to BCS crossover [22].
In this article, we explore in depth the possibilities offered by the attractive model to understand the physics of the repulsive Hubbard model. We pay special attention to the effects of the trapping potential as well as the peculiarities of the physical realization of the negative U Hubbard model in optical lattices. We thus argue that the attractive regime presents a number of advantages for the quantum simulation of the Hubbard model in an optical lattice. We also discuss how the (negative-U ) equivalent phases of the paramagnetic Mott insulator, the anti-ferromagnetically ordered Mott insulator, and the putative d-wave superfluid may be observed.
We would like to emphasize that one of the main advantages of the route suggested here is that the equivalence of doping away from half-filling the optical lattice system in the repulsive regime can be achieved by creating a spin-imbalanced gas. Fermi gases with spinimbalanced populations are nowadays routinely created in the laboratory [23], and the magnetization (which is fixed for the duration of the experiment by the preparation method) can be controlled to within a few percent accuracy. Other advantages will be discussed further below. However, our approach does not solve the problem of how to achieve lower temperatures for the fermions on the lattice [13].
The outline of this article is as follows. In Sect. II we discuss the transformation that formally maps the repulsive Hubbard model into an attractive one. We also define the two attractive Hubbard models we shall be concerned with in this paper, together with their corresponding repulsive models. We also describe how various physical operators and order parameters are affected by the mapping. Some important caveats concerning the realization of the attractive Hubbard model are discussed in Sect. III. In Sect. IV we discuss the equivalent state of the paramagnetic phase of the Mott insulator (as well as possible ways of detecting it), whereas in Sect. V we do the same for the equivalent state of the anti-ferromagnetically ordered Mott insulator. The effect of doping with holes, which, as mentioned above, amounts to a spin-imbalanced situation in the attractive case, is analyzed in Sect. VI. Finally, in Sect. VII we offer the conclusions of the present work as well as mentioning some open problems.

II. THE HUBBARD MODEL AND THE PARTICLE-HOLE TRANSFORMATION
The Hamiltonian of the single-band Hubbard model reads: where t is the hopping amplitude and U the one-site interaction. We consider here only the case where the sites i of the lattice constitute an hypercubic lattice (square in two dimensions, and cubic in three dimensions); i, j in the hopping term means that the sum runs over nearest neighbor sites only. In the above equation, n iσ = c † iσ c iσ is the occupancy of spin σ =↑, ↓ fermions at the i-th site. For simplicity and unless otherwise stated, we will refer in the following to the two dimensional case, and thus to a square lattice. All our results are straightforwardly generalizable to the case of any hypercubic lattice.
We have denoted as H ext all the external fields like chemical and trap potentials as well as an external Zeeman (magnetic) field that act upon the system [41]. Their effects will be discussed below. In the grand canonical ensemble, where the total number operator N = i,σ n iσ and the total magnetization, M = i (n i↑ − n i↓ ), and µ along with h are determined by the condition that averages of N and M over the grand canonical ensemble yield the experimentally observed values [42]. However, cold atomic gases are prepared in eigenstates of both N and M , and therefore the relevant ensemble is the canonical instead of the grand canonical. Although it is important to keep this distinction in mind, we expect that for sufficiently large N , the results of both ensembles coincide and thus, we shall use the grand canonical for the calculation of the experimental signatures of the different phases to be described below in Sect. VI.
The Hamiltonian of Eq. (1) has been written in a form such that for a uniform system the ground state will have exactly one particle per site, at any temperature, for U > 0 and H ext = 0. However, note that in real experiments cold-atomic gases are harmonically trapped. Therefore, the most general form of H ext reads: where i is the shift in the local chemical potential caused by the trap. We have added an unimportant constant to the total energy ( = i ( i − µ) ), which will become convenient further below. In current optical lattice experiments, we have i = 1 2 0 i 2 , but more general forms of the trap may become available in the future. In the case relevant to experiments, the Zeeman field h i = h is uniform (and it is used to adjust the total magnetization in the grand canonical ensemble). However, in Eq. (3) we have assumed it to be site-dependent for further convenience.
We next note that, formally, the sign of the interaction term (the term ∝ U in Eq. 1) can be changed by means of the following (particle-hole) transformation on a bipartite lattice such as the 2D square lattice: Note that the transformation leaves the operators of the spin ↑ fermions unchanged. However, it affects the spin ↓ occupation operator: n i↓ ↔ 1 − n i↓ , and thus the sign of the interaction term changes (U ↔ −U ) while hopping term is left invariant (the minus sign in the right-hand side of Eq. (4) takes care of this). However, it will be important for the discussion that follows that the transformation exchanges the roles of h i and ( i −µ) in Eq. (3). Mathematically, If we insist in using the point of view of the canonical ensemble, the transformation (4) implies that N → N = M + N and M → M = N − N , where N is the total number of lattice sites. Thus, if we consider an unmagnetized (i.e. M = 0) system where N = (1 + x)N (x being the doping, where x = 0 corresponds to half-filling in the uniform case), we have that N = N and M = xN . In words, the doped lattice at U > 0 away from half-filling maps onto a U < 0 system at finite magnetization. We emphasize that, as discussed below, the details of the order that the system develops depend not just on usual factors such as the temperature and strength of the trapping and lattice potential, but are also constrained by these globally conserved quantities.
It is also convenient to recall that, in momentum space, the transformation of Eq. (4) becomes: where Q = (π/a, π/a) is the nesting vector with a the lattice spacing. This expression can be used to obtain the way the different order parameters and the corresponding phases transform between the U > 0 and the U < 0 cases. This is shown in table I. The way the transformation affects the different phases expected for the Hubbard model is also illustrated in Fig. 1.
Finally, for the sake of clarity, we will spell out in what follows the two physically different attractive Hubbard models considered in this article, as well as the repulsive Hubbard models onto which, via the particle-hole transformation Eq. 4, they are mapped. The first attractive Hubbard model (from here on called model A1) has the following Hamiltonian (in the grand canonical ensemble): This is the Hamiltonian that describes cold atomic systems trapped in an optical lattice, with the overall harmonic trapping potential i = 1 2 0 i 2 , and a uniform Zeeman field that can be viewed as a knob to tune the spin imbalance (See Sect. III). In this article, we argue that simulating this Hamiltonian has a number of advantages over current attempts at simulating the repulsive Hubbard model in presence of the harmonic trap (from here on called model R2), whose Hamiltonian reads: The other attractive Hubbard model (from here on called model A2) is described by the following Hamiltonian: This is the Hamiltonian that can be obtained from the repulsive model R2 (Eq. 8) via the particle-hole transformation of Eq. (4). Model A2 has an inhomogeneous Zeeman field stemming (via the transformation) from the trapping potential of model R2. In this article, model A2 is used to help us understand e.g. the coexistence of phases in the current experimental regime of the repulsive Hubbard model (R2). Note that the attractive model A1 that we are advocating does not map onto the currently studied repulsive Hubbard model R2, but instead onto the repulsive Hubbard model in a inhomogeneous Zeeman field (called model R1 here): Thus, to summarize, using the mathematical transformation of Eq. (4), we can relate the physics of the attractive Hubbard models to that of the repulsive ones. In particular the observation of one particular phase (e.g.  4) works. In (a) we sketch how it acts on a single site (e.g. red stands for spin up and blue for spin down), and in (b), the way it transforms different types of states on a uniform optical lattice. In addition to the phases depicted above, an anti-ferromagnetic state at U > 0 ordered along the x or y direction corresponds to a superfluid phase of fermion pairs. The bottom diagram illustrates one of the main points of this paper, namely that doping the attractive Hubbard model corresponds to introducing spin inbalance in the repulsive Hubbard model. a d-wave AF phase) at U < 0 would directly imply the existence of the corresponding phase at U > 0 (the putative d-wave superfluid phase). Below we shall see that the realization of some of these phases in the attractive regime requires sometimes less stringent conditions than the corresponding ones in the repulsive regime. Furthermore, as described above, the exchange of roles of the Zeeman field and the chemical potential terms effected by the transformation implies that we can simulate the doping of the Mott insulator by creating a system with a finite magnetization (i.e. a spin-imbalanced system). Also, even if the temperatures that can be achieved in current experiments do not allow for the investigation of the low-temperature ordered states, we may expect that, by the attractive route, some useful insights can be gained into other controversial issues for the high-T c community, such as the nature of the normal state of the 2D repulsive Hubbard model away from half-filling.

III. REALIZATION OF THE ATTRACTIVE HUBBARD MODEL
In principle, since the onsite interaction energy U is naïvely proportional to the atomic scattering length in free space [18], a s , the U < 0 regime can be accessed by sweeping the magnetic field to the side of an inter-species Feshbach resonance [15] where a s < 0. In the literature of the BCS to BEC crossover [15,23,24], this side is known as the "BCS side" of the resonance.
However, the above point of view entirely neglects the subtleties of the scattering problem on a lattice potential, as it turns out that U is not a linear function of a s , in the general case [25,26,27]. The details of the dependence of the zero momentum scattering amplitude, f 0 (a s ) (and U ∝ f 0 ), on the atomic scattering length, a s , are determined by the dimensionality and other parameters of the lattice [25,26,27]. However, all these results share one common feature, namely the existence of a particular length scale, l * > 0, such that for a s = −l * the scattering amplitude f 0 exhibits a (geometric) resonance [25,26,27]: f 0 ∝ a s /(1 + a s /l * ). Indeed, the resonance can be approached either by changing the lattice parameters, or by changing a s through a Feshbach resonance. We shall focus on the latter case here. To realize the attractive Hubbard model, U < 0, we require tuning the scattering length to the attractive side (a s < 0) but such that |a s | < l * . As we approach the Feshbach resonance and f 0 diverges to −∞, crossing the geometric resonance beyond l * leads to the interaction becoming effectively repulsive, and as a consequence, close to the geometric resonance, a weakly bound state appears. The existence of this bound state has been discussed in the literature [25,26], and if the temperature is sufficiently low (compared to the binding energy of the bound state), many bound states will be created even if the effective interaction between the atoms in the lattice is repulsive. This regime is clearly not described by the repulsive Hubbard model, because it does not take into account the bound states. The situation is similar (although probably less harmful for the system [25]) to the one encountered as a s → +∞, on the so-called BEC side of the Feshbach resonance. In this regime, the scattering amplitude corresponds to that of a repulsive effective interaction, which may lead us to think that the system is described by the repulsive Hubbard model, except crucially, for the existence of the lattice molecular bound states. But indeed, as described in the introduction, the existence of Feshbach molecules [15] leads to inelastic losses.
To summarize, the attractive regime can be reached by making the scattering length a s negative, but not beyond a certain limit where for a s = −l * , (l * depending on the lattice dimensionality and other parameters [25,26,27]) the scattering amplitude has a geometric resonance.

IV. ANALOG OF THE MOTT INSULATOR IN THE ATTRACTIVE REGIME
Let us first start by looking at the U < 0 system with a balanced population of spin up and down fermions (Model A1, Eq. 7 with h = 0) and at temperatures T where |U | can be made such that T < |U | (however, in this section T is assumed to be large compared with t 2 /|U |). Using the transformation, Eq. (4), this corresponds to a half-filled system for U > 0 (model R1, Eq. 10 with h = 0). This is the situation where it is known that a Mott insulator appears in model R1. For U > 0 the existence of the Mott insulator means that the states with more than one particle per site are strongly disfavored due to the large on-site repulsion, U . The corresponding situation for U < 0 is that the only allowed states for every lattice site are either zero or doubly occupied states, as shown in Fig.1. The existence of the Mott insulator in the repulsive regime thus corresponds in the attractive regime to having all fermions form pairs [43]. For sufficiently large |U | these pairs are tightly bound, which means that their existence can be probed by sending photons to photo-associate them into dimers, which are no longer trapped, and therefore can be detected as a loss of atoms from the lattice [28].
It is worth noticing that this system of pairs exhibits a pairing gap (∼ |U |, for large |U |) to spin excitations. Therefore, a measurement of the single particle properties, like the single-particle spectral function which is accessible by photoemission-like spectroscopy proposed in [29] and recently applied to ultra-cold Fermi gases by the JILA group [30], should be able to detect the pairing gap. However, since in model A1 the harmonic trap only couples to the atomic density (cf. Eq. 7), it does not lead to the breaking of the pairs and, therefore, it does not take the system out of the subspace where n i = 0 or n i = 2 (which, by virtue of Eq. 4 corresponds to halffilling, i.e. n i = 1). Indeed, the role of the trap is to lift the large energy degeneracy (for t = 0) of the states in this subspace. In other words, in absence of the trapping potential in model A1, and since the chemical potential µ couples to the total particle number, all states with the same total number of fermions and only doubly occupied or empty sites are degenerate. The trap breaks this degeneracy and selects as the ground state of model A1 the state where all pairs uniformly occupy the lattice sites at the trap bottom. In the strong coupling limit, this state can be regarded as a 'band insulator' of the pairs (see Fig. 3).
A complementary way of arriving at the same conclu-sion relies on the transformation of Eq. (4). After the transformation, model A1 at h = 0 is mapped to model R1, that is a repulsive Hubbard model but in a inhomogeneous Zeeman field (note that h = 0 in this model too, but in this case it couples to the total density). Although the Zeeman field affects the magnetic ordering by ferromagnetically polarizing the fermions at the center of the trap, it does not lead to doubly occupied sites, and therefore it does not affect the characteristic incompressibility of the Mott insulator, which depends on the existence of a gap (∼ U , for large U ) to all density excitations. It is worth to contrast the situation described in previous paragraphs with the one found in current experiments, which are performed in the repulsive regime of the Hubbard model (model R2, Eq. 8). In such a system, one needs to adjust the chemical potential µ (that is, the number of fermions, N ) in order to have a halffilled lattice with one fermion per site at the trap bottom. Otherwise, at too large a µ/U (i.e. large N ), the system energetically prefers to pay the energy cost of having doubly occupied sites near the bottom of the trap, rather than accumulating them far from the center, where the trapping energy is very large. Thus, the lattice at the center of the trap ceases to be in the Mott insulating phase, becoming a band insulator. On the other hand, for µ too small (i.e. small N ) the optical lattice is not uniformly occupied by one fermion site. Testing for the existence of the Mott insulator at U > 0 thus requires checking for the absence of doubly occupied sites, which has been already achieved experimentally by the Zürich and Mainz groups [10,11]. However, testing the absence of holes (which may appear due to thermal or quantum fluctuations, especially as N or U decrase) is a more difficult task. In this regard, in the attractive regime, only the absence of singly occupied states needs to be tested. Such measurement, which implies recording the spatial distribution of lattice sites with different occupations may be accessible through spectroscopic techniques similar to those employed to observe the 'wedding-cake' structure of the Bosonic Mott insulator [31,32].
In Figs. 3,2 we show illustrative plots (in the large |U |/t limit and for T < |U |) of the experimentally measurable/measured density profiles n(r) = n ↑ (r)+n ↓ (r), both in the attractive regime with no spin imbalance (Fig. 3 for model A1 at h = 0) and in the repulsive regime, the half-filled lattice (Fig. 2 for model R2).
[44] We have used data similar to those in current experiments, e.g. the Zurich group [10] (see figure captions).
It is noticeable that for U > 0 (model R2), for the given number of particles vs. trap energy 0 = mω 2 a 2 (m is the mass of a single fermion and ω is the trapping frequency), there is single occupancy in the central region, which is the Mott insulator region. At temperatures that are low compared to U , but higher than any second order (exchange, etc.) processes at the energy scale ∼ t 2 /U (i.e. at currently achievable temperatures), the ↑ and ↓ fermions are equally likely to be found on any site of the Mott-insulating region, as there is no magnetic ordering. The latter can only emerge as the temperature is lowered below the Néel scale, T N eel ∼ t 2 /U . As we move towards the edge of the sample, the site occupancy deviates noticeably from one fermion per site, and a "metallic" shell appears. Its width depends on the temperature, as well as t, U and the trapping energy. On the other hand, for the U < 0 case of model A1 at h = 0, the center of the lattice is filled with fermion pairs resulting from the attractive on-site interaction. Thus, in model A1 (Fig. 3) the effect of temperature (in the limit where one can disregard the second order effects due to the hopping), is to create a finite density of empty sites near the edge of the 'band insulator' of pairs. This explains the deviation of the site occupancy from n(r) = 2 observed in Fig. 3, as the distance to the center, r, increases. We also notice that, because the trap is not in competition with the pairing gap as described above, the stability of the band insulator state is not threatened by the trap (contrary to the Mott insulator at U > 0), and thus the fermion numbers in this regime (U < 0) can be larger than those at U > 0, which is also illustrated by the sizes used to generate the figures. However, for model R2, the trap, as described above, tends to favor the band insulator in the middle for large enough N . When the total number of fermions N becomes smaller than a critical value (see next section), consideration of the effects of hopping will be required.

V. ANTIFERROMAGNETIC ORDER FOR U > 0 AND ITS ANALOG IN ATTRACTIVE CASE
As discussed in the previous section, for large on-site attraction, the only two possible states for a single site are empty sites and doubly occupied. Singly occupied sites are separated from these states by the (large) energy gap ∼ |U |. Thus all the fermions are paired and can be regarded as hard core bosonic entities hopping from site to site with amplitude ∼ t 2 /|U |. To see this, recall that, in the repulsive case [33] and ignoring for the moment the overall harmonic trapping, the low energy effective model for the U > 0 half-filled Hubbard model (with equal proportion of the two species) is a spin-1/2 nearestneighbor Heisenberg model H eff (U > 0) = J i,j S i · S j , J = 4t 2 /|U |. This model transforms for U < 0, upon applying (4), into a half-filled lattice described by the following effective Hamiltonian in terms of hard-core bosons (b i = c i↓ c i↑ ): where i, j means that the sum runs over nearest neighbor sites only. For the U > 0 case, the ground state of the Heisenberg model is a (s-wave) antiferromagnet (sAF). In the absence of any terms in the Hamiltonian that distinguish spin up and down fermions, that is, for a spin-isotropic Hamiltonian, the staggered magnetization can point in any of the spin directions, x, y, or z. For the U < 0 case, the resulting system, Eq. (11), is thus a system of hard-core bosons (arising from the pairing of two different spin fermions) hopping on the lattice with the kinetic energy ∼ J, and experiencing nearest neighbor repulsion interaction also ∼ J. Several phases can be realized in such a system. The resulting nearest neighbor repulsion and kinetic energy in Eq. (11) favors an alternating pattern of empty and doubly occupied states, i.e., the checkerboard state in 2D also known as the "charge density wave" (CDW) state. This alternating pattern of empty and doubly occupied states for U < 0 is the analogue of the antiferromagnetic (along zdirection) state for U > 0 as shown in Fig.1 and Table I. Another alternative ground state is simply an s-wave superfluid of these bosonic pairs (sSF). In fact, the degeneracy for ordering along any direction x, y, z of the sAF for the U > 0 case maps to a degeneracy between the CDW state and the SF state for the U < 0 case, which one can see formally from the order parameter mapping, using the transformation Eq. (6): the checkerboard order parameter ∆ CDW maps to the antiferromagnetic order parameter in the z spin direction M z stag , while the SF order parameter ∆ (s) , ∆ (s) † maps to the AF order along −, + direction, Table I. However, the presence of the harmonic trap on the attractive side leads to interesting effects. As can be seen from Eq. (3), in model A1 Eq. (7), the trap acts as a chemical potential for the pairs in the U < 0 case, and thus favors a completely filled center of pairs (which is just a band insulator in the center, cf. Fig. 3), rather than either a CDW or superfluid state. In fact, in the limit where tunneling between sites is suppressed, the ground state in the trap potential is this band insulator, and only the kinetic energy and the nearest neighbor repulsion of order J prevent this state to occur. Mapping back to the U > 0 case, one sees that the trap transforms into a Zeeman field along the z direction to become model R1 Eq. (10). This Zeeman field will lift the degeneracy between the various magnetic states and then polarize ferromagnetically the spins rather than favor an antiferromagnetic order. This competition is summarized in Fig. 4.
This effective Zeeman field in model R1 goes from large and positive in the center of the trap to large and negative in the periphery. In the center of the trap one would thus have all the spins polarized up. This phase corresponds via the transformation, to a pair on each site and thus to the band insulator of pairs in model A1. Whether the effective field in the center of the trap is enough to polarize fully the spin depends on the number of particles in the trap and will be discussed below.
Further from the center, the Zeeman field decreases and becomes negative. One has thus in a certain radius a shell where the spins are not fully polarized. In this region the spins preserve an antiferromagnetic order in the direction perpendicular to the Zeeman field albeit with a reduced amplitude (i.e. the AF order parameters M (s)± Q are non-zero). In other words, in this shell, the system of model R1 thus possesses antiferromagnetic order in the x, y direction in spin space. After the transformation this x, y-antiferromagnetic phase for U > 0 maps for U < 0 to an s-wave superfluid phase, as shown in Fig. 4. We thus see that looking for antiferromagnetism in the repulsive Hubbard model amounts to probing for superfluidity in the U < 0 one. Then beyond a certain radius the trap prevents the pair to exist, and the system is empty (see Fig. 3). In the U > 0 system, this corresponds to a Zeeman field that is large enough to fully polarize to down the system of spins [46].
In summary, on the U < 0 side the trap does not really spoil the search for the analogue of the antiferromagnetic phase. The corresponding phase is now a simple s-wave superfluid that can be probed by similar techniques that have already been used in the continuum. The extra phases induced by the trap only potentially reduce the spatial extent of the superfluid shell, but since they both corresponds to insulating regions (either band insulator or empty region), they should not spoil the observation of the superfluidity. To minimize the central band insulator region for the attractive model A1 (and hence maximise the superfluid signal), the chemical potential at the trap center must not be too strong, so that when tunneling is turned on, the resulting pair hopping and pair repulsion energy scale J (see Eq. 11) can overcome the trapping potential to delocalize the band insulator. First, in the strong coupling limit, consider the case where the whole trap consists of mostly the band insulator. We can estimate the instability of the band insulator to this scale J: taking a pair at the edge of the band insulator system (R = R 2 ) to a site just beyond the edge (by a lattice spacing a) costs an energy ∼ 0 (R 2 /a), and this is to be balanced against J = 4t 2 /|U |. For a large enough system, π(R 2 /a) 2 = N/2 where N/2 is the total number of pairs. Thus the band insulator state will become unstable when roughly, t 2 /|U | > 0 N/2π. However, for a full SF state in the central region of the trap, one needs to move outward more than just the pairs at the boundary of the band insulator. A similar estimation indicates that the √ N in the above criterion is replaced by linear in N . Thus for a given trap energy 0 , there is an upper limit to total number of fermions that can be loaded into the trap. For example, using data similar to those in the 40 K experiments of the Zurich group, at an optical lattice depth of 5 recoil energy for a laser wavelength of 825 nm, a mean trapping frequency of 80 Hz, and an (independently tuned) ratio |U |/t = 8, the total number of fermions has to be smaller than ∼ 500 to have a pure superfluid core. Remembering that this is only a rough estimate, this is nevertheless a rather small number to achieve under current conditions [34]. This critical number can be larger if we allow for a central region of band insulator in addition to a superfluid shell. Depending on the sentivity to detect the superfluid shell, and thus the number of atoms that one needs to have in the superfluid, the number of fermions can be larger than the above estimate. The situation is better for the lighter atom 6 Li. Using the same numbers as above, but with a laser wavelength of 1064 nm, the upper critical number becomes ∼ 7000, which should be feasible in current experiments.
As for experimental signature, the observation of the coherence peak in the momentum distribution should signify the onset of the BEC of pairs, and superfluidity can be proven when vortices are observed when rotating the trap and optical lattice (although to date, vortices have been seen in a rotating bosonic BEC superfluid in an optical lattice [35]). These are simpler probes than say, using noise correlation [36,37,38] to deduce the broken translation symmetry of the AF state in the U > 0 case [39].

VI. EFFECT OF DOPING
So far, we have shown that the spin-balanced population for U < 0 already presents several advantages to tackle the Mott and AF physics. But one of its main advantages is the possibility to effectively "dope" the U > 0 system by looking at spin imbalanced U < 0 systems. This then allows to settle experimentally the still controversial issue of the presence or absence of the d-wave superfluid (dSF) in the repulsive Hubbard model doped away from half-filling. As mentioned in the Introduction, this doping may be difficult to do directly in the U > 0 system due to the presence of the overall harmonic trap. Via the transformation, the U > 0 model away from halffilling maps to the U < 0 model with an effective Zeeman field. In the context of cold atom experiments, this corresponds to a (fixed) imbalance of spin-up versus spin-down fermions, which can readily be achieved to an accuracy of a few percent currently [40] (see [23] for a review).
We now examine some of the observables in that case, and in particular what would be the consequences of the existence of a d-wave superfluid phase in the repulsive Hubbard model for the U < 0 phase. We perform this analysis for the homogeneous system and will discuss the possible effects of the trap at the end of this section.

A. Transformation of the operators
Under the transformation, Eq. (6), the superfluid order parameter ∆ (α) † (the label α = s, d indicates the swave or d-wave symmetry of the order parameter) maps to the commensurate antiferromagnetic order parameter where previously, φ (s) k = 1 is the s-wave form factor and φ Thus, if there is a regime of dSC in the U > 0 Hubbard model, then correspondingly, there is a regime of dAF in the U < 0 model. This is the analogue of the well-known mapping at precisely half-filling of the Hubbard model between the ground states of sSF or sCDW at U < 0 and the sAF at U > 0 (see Fig.1 and Table I).

B. Momentum distribution and noise correlation: a mean field calculation
In this section, we briefly outline a simple calculation to illustrate experimental signatures of dAF states that may exist in the U < 0 system. We need to emphasize from the outset that since there are no microscopic analytical calculation of the dSF nor the dAF states in the 2D Hubbard model, we will instead use a toy mean field model that does give rise to such states, in order to calculate some expected responses such as the momentum distribution and noise correlation [36,37,38] . While such a mean field model misses out correlations and quantum fluctuations, the objective here is to demonstrate that the symmetry of the CDW or SF order parameters has very definite signatures in noise correlation experiments [39]. Thus, we simply assume that the U < 0 Hubbard model acquires a mean field (MF) form, which in momentum space becomes: where we have introduced (by hand) the mean field order parameter Eq. (12), and as before, α labels the s or d x 2 −y 2 order parameter. As usual, via a global gauge transformation, the order parameter can be chosen to be real: M Our main interest is in the phases of the U > 0 Hubbard mode without Zeeman field away from half-filling. This then corresponds to the U < 0 mean field model above (Eq. 12) at µ = 0 and finite h. In principle, the nesting wavevector Q should be a variational parameter to be determined from the particular band structure etc. However, the following mean field theory only makes sense if Q = (π/a, π/a) is a commensurate wavevector: the interaction above couples k with k + Q, and this in turn couples to k + 2Q which is the same as k only for commensurate Q.
This MF Hamiltonian is diagonalized by a Bogoliubov rotation: Then, provided the MF Hamiltonian becomes with the "magnetic band" energies: The ground state is then made up by filling these bands up to the respective Fermi surface k F ασ with the occupation number n α k = AF Q |α † k α k |AF Q = Θ(k F α − k) (and similarly for the β band): Note that this wavefunction does not have definite number N ↑ and N ↓ for each species: this is a consequence of the quasi-particles α and β carrying indefinite spin, an analogue of the textbook number non-conserving Bogoliubov quasi-particles or the BCS wavefunction or the BCS mean field Hamiltonian. In turn, here for the AF, we have indefinite spin (but definite charge) quasi-particles because we have assumed the mean field decoupling to be in the S + -axis in spin space. It is straightforward to show (at least when h = 0) that the same results can be gotten for a spin-conserving mean field wavefunction (eg. when the mean field decoupling is in the z spin axis). On the square lattice, since kσ + k+Qσ = −2µ σ , taking the same chemical potential and bare dispersion for the two spin species, E α,β k = −µ ± Ω k . Hence for the s-wave case there is always a band-gap of size ≥ 2g Q M Q , and the minimum gap occurs at k↑ − k+Q↓ = 2h. Thus the half-filled lattice (1 fermion per site) with N ↑ = N ↓ has the β-band completely filled up and an empty αband: this is a magnetic insulator since adding one more ↑ fermion has c † k↑ → cos(θ k ) α † k creating an α-particle, but this costs at least the band-gap energy (and similarly for adding a ↓ spin particle). For the d x 2 −y 2 -wave case, since the form factor φ k = cos k x − cos k y has nodes at k x = ±k y , the band gap also vanishes at these positions, giving rise to a pseudo-gap only. At less (more) than half-filling, the "magnetic Fermi surface" is within the β-(α-) band and there is no energy gap to adding an extra fermion: the system is "metallic". We should be careful about the meaning of E α,β k : this is not the physical excitation energy (unlike in the BCS MF theory): in particular, the spin-wave (Goldstone mode) spectrum is missing.
Using the ground state Eq. (16), the T = 0 momentum distribution for each spin component is: In experiments, assuming we can image perpendicular to the plane of the lattice, this result has to be convolved with the Wannier function w(k) (in momentum space) for the optical lattice: n σ k=mR/τ ∝ |w(k)| 2 c † kσ c kσ . (When there are more than 1 plane of the lattice, we also need to integrate over the planes.) The T = 0 noise correlation is proportional to the connected correlation function: (Again, for the experimentally measured noise correlation, this result has to be multiplied by |w(k)| 2 |w(k )| 2 .) We plot in Fig. 5 the noise correlation for the 2D d x 2 −y 2wave AF, at µ = 0. Noise correlation should clearly distinguish between s-wave and d x 2 −y 2 AF, thanks to the appearance of nodes in the ±(π/a, π/a) directions.

C. Experimental observation and complications
We thus see that the negative U side offers the great advantage to access directly the doped regime, without having to suffer directly from the presence of the trap. We have presented in the previous section some possible experimental signatures that a d-wave superfluid phase would give when suitably transformed to the negative side. Of course, even if the situation is potentially improved by the transformation to the negative U side, probing the doped phase is a considerable challenge. Some of the limitations are obvious. The most immediate one is the effect of the temperature. Indeed, already for the antiferromagnetic phase, the temperature is in competition with an energy of order J 4t 2 /U the kinetic energy of the pairs. When looking at the doping effects one has to face even smaller energies. Lowering the temperature is thus a must.
The second limitation is again the trap. Although the mapping to the negative side avoids the direct effect of the trap on the doped holes, the trap still has a potentially indirect effect. Indeed as we saw in the section on antiferromagnetism, the trap will act, for the U > 0 side (model R1), as an effective Zeeman field, and lead to two shell regions, with fully polarized spins up or down (Fig. 4). When holes are introduced into the system of model R1, they will thus have the possibility to go in one of these two polarized regions or in the antiferromagnetic region, which corresponds to the shell where the effective Zeeman field is not strong enough to fully polarize the system. Where the holes go will depend on how much kinetic and interaction energy they can gain in the three different phases since they are not sensitive directly to the presence of the Zeeman field. A naive calculation neglecting the presence of the energy scale J, leads directly to the holes going in the interface between the fully polarized up and fully polarized down regions. In the U < 0 language, i.e. for model A1, this corresponds to the excess of one spin species going to the edge of the band insulator region. This seems to suggest, that when the energy scale J is put back in the problem, the holes will indeed go into the antiferromagnetic region, and thus can lead potentially to the d-wave superconducting phase there. This is however a delicate question since the kinetic energy of a hole in a ferromagnetic environment is in principle higher due to the absence of frustration of the antiferromagnetic order upon hole motion. This important question thus fully deserves more analytical and numerical study. It remains however academic until serious progress on the temperature issue have been made.

VII. SUMMARY, CONCLUSION
In this article, we have explored in depth the possibilities offered by quantum simulating the attractive model (A1) in cold atoms in optical lattices, to understand the physics of the repulsive Hubbard model, via a well known canonical particle-hole transformation. We have argued that there are certain advantages in doing experiments in the attractive regime.
For the undoped case the attractive side replaces the Mott phase and the antiferromagnetic phase by a phase composed only of pairs (Mott insulator for U > 0) and that undergoes a superfluid transition (antiferromagnet for U > 0). The trap which exists in any realistic experiment does not really affect the observation of these two phases since it can only add a core of band insulator at the center, thereby not spoiling the observation of the superfluid. The attractive side also offers the advantage of only having to test for the pairing for the observation of the band insulator that must be simpler than testing for the absence of doubly occupied and empty sites on the repulsive side.
Another key advantage of the attractive side is the relative ease in controlling the spin population imbalance. Via the canonical transformation, this corresponds to doping away from half filling for the repulsive side, which is hard to achieve because of the presence of the harmonic trapping potential that moves the holes away from the central region of the trap. This doping in the repulsive regime is needed to quantum simulate and explore the possibility of a d-wave superfluid, the fundamental question of great relevance to the cuprate high temperature superconductors. This question can be answered instead in the attractive regime by exploring the presence or absence of a d-wave antiferromagnet. We indicate in this paper several ways to probe for the existence of such a phase. We have also pointed out a number of technical difficulties of the proposed approach and, in some cases, suggest possible solutions. that the effect of the physical magnetic field is purely to change the Zeeman energy of the hyperfine states of the two species of fermions considered here. In particular, since cold atoms are neutral, there are no orbital effects on the superfluid, such as the Meissner effect for a superconductor. Henceforth, we shall call the physical magnetic field the Zeeman field, to emphasize this.
[42] Note that here because the interacting term contains a term proportional to N , the definition of µ is somewhat non-standard.
[43] By pairs it is meant in this article a many-body bound state of the same nature as a Cooper pair. These pairs should be distinguished from molecular bound states arising from poles in the scattering amplitude of the twobody problem.
The total number of fermions N = P i n i↑ + n i↓ , which fixes the chemical potential µ(T, U ). We have also found that, at low temperatures, the density profiles for U < 0 can be obtained from a mapping to a free fermion model, details of this model and its consequences will be reported elsewhere [27].
[45] Another way to see the same correspondence is using wavefunctions. A cartoon ground state for AF ordering in the x direction is of the form: |AFx ∼ Πj[c † j↑ + (−1) j c † j↓ ]|vacc , where |vacc is the vaccum state for the c-fermions. Defining d † j↑ = c † j↑ and d j↓ = (−1) j c † j↓ , with |vacc = Πjd † j↓ |vac d , then, |AFx −→ Πj[d † j↑ d † j↓ + 1]|vac d , which is a coherent state of pairs b † j = d † j↑ d † j↓ .
[46] Note that for small U the interplay of the inhomogeneous magnetic field and the kinetic energy leads to interesting effects for the intermediate phase, due to the difference between the radial and angular directions for the hopping of a particle. These effects will be discussed elsewhere [27]