# Exact quantum Monte Carlo study of one dimensional trapped fermions with attractive contact interactions

## Abstract

Using exact continuous quantum Monte Carlo techniques, we study the zero and finite temperature properties of a system of harmonically trapped one dimensional spin 1/2 fermions with short range interactions. Motivated by experimental searches for modulated Fulde-Ferrel-Larkin-Ovchinikov states, we systematically examine the impact of a spin imbalance on the density profiles. We quantify the accuracy of the Thomas-Fermi approximation, finding that for sufficiently large particle numbers () it quantitatively reproduces most features of the exact density profile. The Thomas-Fermi approximation fails to capture small Friedel-like spin and density oscillations and overestimates the size of the fully paired region in the outer shell of the trap. Based on our results, we suggest a range of experimentally tunable parameters to maximize the visibility of the double shell structure of the system and the Fulde-Ferrel-Larkin-Ovchinikov state in the one dimensional harmonic trap. Furthermore, we analyze the fingerprints of the attractive contact interactions in the features of the momentum and pair momentum distributions.

###### pacs:

03.75.Ss,71.10.Pm,02.70.Ss## I Introduction

Mapping the phase diagram of fermions with attractive interactions and spin imbalance is a challenging problem, both from the experimental and the theoretical side. For an unpolarized system, an effective attractive interaction leads to the formation of Cooper pairs carrying zero momentum.cooper () However in the presence of spin imbalance, where the different spin components have their Fermi energies shifted from one another, the pairing mechanism changes. Depending on parameters, different theories apply. The Sarma or breach-pair theory, which is stable in the presence of long range interactions, suggests a coherent superposition of Bardeen-Cooper-Schrieffer (BCS) pairs and normal fermions,sarma (); liu (); forbes (); gubbels () while in other regimes one finds Fermi surface deformationsmuther (); sedrakian () or a coexistence of normal fluid and superfluid.bedaque () One of the most interesting possibilities is the Fulde-Ferrel-Larkin-Ovchinikov (FFLO) states,fulde (); larkin (); casalbuoni () which reconcile spin imbalance and superconductivity by pairing particles whose momentum does not sum to zero. For example, in the FF state, -spin fermions with momentum are paired with -spin fermions with momentum , resulting in a Cooper pair with momentum . In the LO theory, the Cooper pairs carry no net current, but are in superpositions of momentum states, leading to an order parameter that oscillates in space; the excess spin polarization resides in the nodes of the order parameter. Analogs of the latter state are found in one dimension (1D). We will use an exact quantum Monte-Carlo procedure to study a gas of spin imbalanced harmonically trapped 1D fermions, with a goal of identifying the signatures of FFLO physics which will be seen in cold atom experiments.

Although not completely transparent, the best evidence for the FFLO state in solid-state systems comes from layered organic superconductorssingleton () and heavy-fermion materials.radovan (); miclea () However, recent progress in optical lattice experiments with cold atoms has opened the avenue to study spin imbalanced fermions with attractive interactions in a clean and direct way. zwierlein (); zwierlein2 (); shin (); schunck (); partridge (); partridge2 () Unfortunately, for a bulk three dimensional gas of atoms with short range interactions, the FFLO state is found for a very narrow range of chemical potentials,parish () implying that only a small volume of a trapped atomic cloud would be in this state. Hence attention has been shifting to one dimension, where an analog of the FFLO state has a very large region of stability.liu_hu_drumm_fflo ()

Since there is no true off-diagonal long-range order in a strict 1D system, the 1D FFLO state is characterized by an algebraic quasi-1D order. Although the thermodynamics of a uniform 1D Fermi gas with contact interactions can be found exactly via the Bethe ansatz at zerogaudin (); takahashi () and finite temperature,takahashi_temp (); guan () the characterization of the ground state with spin imbalance is not easy: many correlation functions are unknown due to the complicated structure of the Bethe ansatz. Other properties are well known: for example, Yang firmly established the existence of the FFLO state through a bosonization technique, showing that it emerges via a continuous transition from the BCS state.jang () Further, numerical methods have shown that the FFLO is the ground state of the Hamiltonian on a lattice as soon as the system is polarized batrouni (); rizzi (); luscher ().

Experimentally, an array of quasi one dimensional tubes can be created at the intersection of two orthogonal standing waves generated by laser beams.paredes (); moritz (); moritz2 () Along each tube there would be no lattice, but the particles will be confined by a longitudinal harmonic potential, coming from the spatial profile of the lattice lasers, and possibly from some additional potential. Recently Orsoorso () and Liu at al.liu_hu_drumm_prl () studied the properties of an atomic cloud in such a 1D harmonic trap. By using a Thomas Fermi (TF) approximation based on the Bethe ansatz solution of the homogeneous Hamiltonian, they found that in the trap the system phase separates. The inner shell of the tube is always a partially polarized FFLO state, while the outer shell can be either a fully paired BCS or a fully polarized normal state. Our calculations are designed to go beyond the approximations inherent in those calculations, and accurately quantify the properties of these 1D clouds.

Other researchers have also attempted to investigate the accuracy of the TF approximation. For example, Xianlong and Asgari developed an exchange-correlation functional to perform density functional theory (DFT) calculations with the local density approximation (LDA).xianlong () While revealing, those calculations are only approximate. In a different direction, several researchers have used the density matrix renormalization group (DRMG) to study a related system, where there is a 1D lattice along the tube.feiguin (); tezuka (); capelle () One needs to be cautious about applying those results to the continuous system as the underlying lattice in the model could introduce effects of incommensurability in the density profiles.buchler (); molina () The continuous limit is difficult to access due to the large number of sites required in the DMRG evaluation for an extremely dilute system.

We will also consider how temperature affects the density profiles. Since the superfluid critical temperature is zero in 1D, thermal effects could be quite dramatic, and may not be captured by mean field theories.liu_hu_drumm_finitet () Experiments will inevitably be performed at finite temperature, and it is imperative to know to what extent finite temperature will obscure the clear delineation of phases which are expected in a trapped zero temperature gas.

In this paper, we present exact results for 1D fermions subject to a zero-range attractive interaction and a longitudinal harmonic confinement. We resolve the dependence on the spin imbalance, effective coupling, and temperature, within a continuous quantum Monte Carlo (QMC) framework. We use diffusion Monte Carlo (DMC)reynolds:5593 () for zero temperature ground state calculations and path integral Monte Carlo (PIMC)ceperley_pimc () for finite temperature, which are ideal numerical tools to find the properties of 1D Hamiltonians, as they do not suffer from the “sign problem” which affects QMC simulations in higher dimensions.ceperley91 () In 1D they give unbiased energies, with their accuracy only limited by statistical noise. We correct other potential errors in the DMC algorithm, such as the time step , population bias, and mixed estimate errors, by appropriate extrapolation and forward walking techniques.calandra:1998 () The only systematic bias in the PIMC algorithm is the time step, which we control by working with sufficiently small .

The paper is organized as follows. In Sec. II we outline the model Hamiltonian used in our calculations, and the experimental parameters of the system. Sec. III includes a description of the QMC methods with an emphasis on the features required to simulate 1D fermions interacting with a zero-range potential. In Sec. IV we present our results and compare them with the TF findings. In particular, we focus our attention to the finite size and temperature effects on the phase boundaries. We conclude in Sec. V with some highlights and remarks.

## Ii Model

The Hamiltonian which describes 1D fermions in a trap and interacting with an attractive zero-range potential is

(1) |

where is the total number of particles with mass , is the strength of the contact interaction, and is the harmonic frequency of the longitudinal trap. The effective 1D Hamiltonian in Eq. 1 is a very good representation of a single tube in the optical array provided that and , where is the frequency of the transverse harmonic confinement. The above requirements are met in the experimental setup by properly tuning the height of the lattice depth, which sets the value of . In such a 1D geometry with just one open channel, the interactions can be modeled by as with given byolshanii (); bergeman ()

(2) |

where is the 3D -wave scattering length, is the transverse oscillator length, and , where is the Riemann zeta function. Here and throughout the paper we have used a sign convention where is attractive. An ongoing experiment at Rice Universityrandy () is using of atoms per tube, the actual value depending on the position of the tube in the array. Typical values for the optical lattice are , , while will be tuned by changing . The typical temperatures are around , with the Fermi temperature of fermions with spin . Throughout the paper we will use , i.e. the length will be measured in units of the harmonic oscillator length , and the energy and temperature in units of the level spacing . In these units the energy for is and the binding energy of a pair is .

## Iii Methods

The ground state QMC techniques used in this work are described in Subsec. III.1, while the finite temperature PIMC method is reported in Subsec. III.2. In order to compare our results with the local density approximation, we computed also the TF density profile () and spin density profile () of the system, which are solutions of:

(3) |

where is the chemical potential of the species with spin , and is the local chemical potential found from the Gaudin exact solutiongaudin (); takahashi () of the homogeneous model. The Thomas-Fermi results have been widely discussed elsewhere, and we refer the reader to previous worksorso (); liu_hu_drumm_prl (); liu_hu_drumm_fflo (); xianlong () for complete details. Here, we mention that Eq. 3 can be rescaled into dimensionless variables with coupling strength controlled by . Therefore, the effective coupling constant of the trapped system depends also on the number of particles and scales as .

### iii.1 Ground state QMC

To study the ground state properties of the system, we employed variational and diffusion quantum Monte Carlo methods. Their basic ingredient is the trial wave function , a good approximation to the ground state. We used the form:

(4) |

where is a bosonic function and is the antisymmetrized product of Hermite polynomials of -th order for particles with spin . can be written as determinant of a Van der Monde matrix:

(5) |

The product in the Eq. 5 is convenient, as it can be evaluated in order operations. The Hermite basis set yields an exact variational wavefunction in the limit of vanishing inter-particle interactions. The Slater term determines the nodes of the variational wave function. In 1D, the exact nodes of the ground state are defined by the coalescence conditions between two like-spin particles (i.e. with ). This result holds even for an interacting system, and consequently there is no sign problem in the diffusion Monte Carlo calculations.ceperley91 ()

The bosonic function is real and symmetric under permutation of two fermions, and is the sum of one-body (), two-body (), and three-body () correlations. The one-body part reads:

(6) |

where is the spin of the -th particle, and is written in a Padé form which goes as at large so to assure the right asymptotic behavior given by the harmonic confinement. In the absence of inter-particle interactions, is the exact solution.

The two-body (Jastrow) factor is defined as:

(7) |

with pair functions chosen to fulfill the cusp conditions set by the contact interactions of the Hamiltonian in Eq. 1. In the case of two unlike-spin particles interacting with a delta function potential, the exact solution is given by ; note the cusp at . The exact many-body function will have the same cusp between atoms with different spin. We use the following pair functions:

(8) | |||||

(9) |

where , , and are variational parameters, and is the strength of the contact interaction.

The confining potential leads to inhomogeneities that make it convenient to include an inhomogeneous Jastrow function:

(10) |

where the pair function is modulated by the function which depends on the center of mass of the particles and . The analytic form of is the same as in Eq. 9 also for , while is chosen to be:

(11) |

where controls the inhomogeneity. For one recovers the homogeneous case.

Finally, we employed a three-body term:

(12) |

where also in this case the functional form of is the same as in Eq. 9. We found this term important in the intermediate and strong coupling regime, when the local energy is dominated by the inter-particle potential.

Our variational ansatz includes about 30 parameters , optimized by means of the stochastic reconfiguration algorithm,sorella:024512 (); casula:7110 () which minimizes the total energy of the wave function with an iterative procedure based on the value of the forces, , acting on the parameters at each step. The optimal variational wave function is then projected to the ground state by using the diffusion Monte Carlo method (DMC)reynolds:5593 (), which provides an unbiased estimate of the ground state energy. The physical observables, such as the charge and spin density profiles, are computed with the forward walking propagation of the mixed DMC distribution.liu74 (); calandra:1998 ()

### iii.2 Finite temperature QMC

The finite temperature calculations are performed by using the path integral Monte Carlo (PIMC) methodceperley_pimc () which is based on sampling the thermal density matrix , with , the temperature, by means of a mapping of the quantum system into a classical system of polymers or Feynmanfeynman ()“paths”. and are all-particle configurations. During the Monte Carlo random walk the polymers are distributed according to the probability

(13) |

where the integration from to has been discretized into “time slices” of length , is the action, and . An accurate approximation of for the system under consideration is important to make the PIMC calculations feasible and efficient. In the case of the Hamiltonian in Eq. 1, can be written using a Trotter expansion as a product of four factors: the non-interacting (), the harmonic (), the interacting (), and the fermionic (). The non-interacting part corresponds to the propagator for free particles

(14) |

We use the primitive approximationceperley_pimc () to include the effect of the harmonic external potential:

(15) |

To account for the interactions we make use of the exact action for two particles in free space interacting with an attractive contact potential gaveau (). In terms of the relative distances of the pair (, ), the resulting propagator is:

(16) | |||||

where , , . From this propagator, one obtains the many-body action by means of the pair-product approximationceperley_pimc ()

(17) |

where the product runs over all unlike-spin pairs. This, as the other approximations, becomes exact as .

To include Fermi statistics into the PIMC method ceperley_fn () the paths are constrained by the nodal regions, which act like barriers. Since the nodes are exactly known in 1D, and they do not depend on , the restricted path integral formalism gives the exact result. An approximate fermion action is given by the image approximation:

(18) |

for each spin sector. It includes the nodes between neighboring particles as zero’s of the density matrix, and becomes exact in the limit. We assume that the like spins are ordered for .

The probability distribution in Eq. 13, based on the the density matrix defined above, is sampled via a generalized Metropolis algorithm. The paths are displaced with two types of moves: a multi-level bisection (level two and three have been used here) and a global move, which attempts to displace an entire polymer at once. When the convergence is reached, the thermal expectation values are computed as described in Ref. ceperley_pimc, .

## Iv Results

The QMC results shown in this section characterize the ground state and finite temperature properties of the system. We compute the density , spin density , and local polarization () profiles. The latter quantity is extremely useful, as it locates the paired (), partially polarized (), and fully polarized () shells. For sufficiently high densities, the center of the cloud is always partially polarized, but the edge may be paired or fully polarized. The Thomas-Fermi approximation yields the phase diagram shown in Fig. 1, which shows the expected boundary between these two regimes. We decided to predominantly study the system with , since it is roughly the average number of fermions per tube in typical experiments. As marked on the phase diagram, we present results for a range of interaction strengths and polarizations. In our numerical profiles we are particularly interested in how the visibility of the interface between different regions is influenced by finite particle number and finite temperature. By analyzing the experimental data within a Thomas-Fermi approximations one could use the locations of these boundaries to map out the phase diagram. We find that for and this Thomas-Fermi analysis is reasonable, but for fewer particles or higher temperatures one would need to use a much more sophisticated analysis of the experimental data to learn about the bulk zero temperature phase diagram. We report our results for the ground state and the finite temperature calculations in Secs. IV.1 and IV.2 respectively.

### iv.1 Ground state

In Fig. 2 we report on the density, spin density, and local polarizations of particles with various polarizations and interaction strengths . From the TF results in Fig. 1, one expects that the edge of the cloud will be fully polarized for . Similarly, the edge should be fully paired for . The others: are near-critical, and one expects that the FFLO state should extend nearly to the edge of the cloud. The results in Fig. 2 are consistent with these expectations. For the TF near-critical points, we found that all of them feature fully polarized wings. When the outer shell of the cloud is polarized, the spin density has a peak at the edge, which is replaced by a rapid fall-off when the edge is fully paired. Detecting the fully paired shell will be easiest at strong coupling, where the size of the outer paired region can be large.

Notice that the spin polarization plotted in Fig. 2(c) is much more noisy than for weaker interactions. Indeed, from the computational point of view a problem of ergodicity starts affecting the Monte Carlo (MC) sampling. As the coupling increases, one needs to run longer to equilibrate the sample and have good statistics in the expectation values, because sometimes the evolution gets stuck in a “persistent” configuration. In this regime, it is difficult for unpaired fermions to cross a strongly bound pair due to the Pauli principle. We note that it is possible that the same physics will affect also the dynamics of the atoms in the experiment. The physical origin of the ergodicity problem in the MC evolution, namely the formation of strongly bound pairs, could also cause a slow equilibration in the experimental setup. A coupling between and is a good balance between the visibility of the fully paired region and the ergodicity of the pairs in the tubes.

From the analysis presented above, it seems that the TF mean field phase diagram is a good approximation for the couplings and polarizations studied so far. It is clear that the TF approximation is not able to reproduce the short length-scale charge and spin fluctuations, since the local density approximation assumes a slow varying density profiles, as already pointed out in previous studies.mueller (); xianlong () However, the TF boundaries of the phase separated states are in a quite good agreement with our exact QMC calculations. In order to better understand this agreement, we study the dependence of the local polarization on the number of particles and coupling , and compare our QMC results with the TF predictions. The local polarization for different system sizes (, , ) with and is plotted in Fig. 3. At and , there are qualitative disagreements between the QMC and TF methods, which disappears at . Indeed for small systems our QMC calculations give partially polarized wings, while they are fully paired in the TF approximation. This can be due to a proximity effect not accounted for in the completely local TF approximation. We speculate that the FFLO state partially extends into the wings, making them partially polarized. We cannot confirm this picture since we did not study how the FFLO order parameter depends on position. For the wings become fully polarized, as the effective coupling is reduced, and the polarization profile turns out to be in agreement with the TF results, although the transition from the FFLO to the fully polarized state is much smoother in our QMC calculations. The boundaries between phase separated states are smeared out over a length-scale of order the atomic spacing. For small systems the very concept of phase separation begins to break down.

We check the agreement between the QMC and TF methods for large number of particles () in a more systematic way by studying the dependence of the polarization on with . These results are plotted in Fig. 4 for , , , and . The TF polarization profile follows quite closely the DMC points. Both theories predict a “transition” from a fully polarized shell to a fully paired one between and . Apart from fluctuations and statistical noise on the points, the only significant difference in the comparison is the size of the fully paired region which is smaller in the QMC results, where the FFLO state seems to extend farther from the center of the trap. This is a common feature of our calculations. The TF approximation always overestimates the size of the fully paired shell. Apparently, this is a drawback of the local density approximation, which makes the effective attractive interaction stronger.xianlong_lda ()

We complete the comparison between the TF and QMC results by also computing the energetics of with different polarizations (, , ) and system sizes (, ), reported in Tab. 1. The agreement in the ground state energies per particle between the two methods is remarkable, even at strong coupling (small ). This means that the total energy is insensitive to the local differences seen in the shell structure of the cloud, interpreted as finite size or edge effects. It also highlights the importance of an accurate treatment of local correlation effects, which can be missed in the TF approximation, and lead to subtle changes in the phase boundaries, whereas the overall trend is correctly captured by the mean field method, even at a quantitative level.

TF energy | DMC energy | ||
---|---|---|---|

20 | 0% | ||

20 | 10% | ||

20 | 20% | ||

200 | 0% | ||

200 | 10% | ||

200 | 20% |

In previous works,batrouni (); tezuka (); feiguin (); rizzi () it has pointed out that the pair momentum distribution function reveals fingerprints of the FFLO state in the center of the trap. The momentum and pair momentum distributions can be extracted in experiments from a statistical analysis of time-of-flight images of the expanding cloud of atoms when the confining potential is removed. Here we study the momentum distribution function defined as

(19) |

where is the one-body density matrix

(20) |

and is the ground state normalized such that

(21) |

In the actual DMC calculations, the density matrix in Eq. 20 is computed as follows

(22) |

where indicates the MC averages over sampled from the DMC mixed distribution. Although the forward walking technique for non-local operators is possible, it has rarely been applied to the evaluation of the momentum distributions.kalos () In our case, we checked that the density matrix computed at the variational MC level is close to the DMC mixed average (MA), due to the accuracy of the trial wave function . Therefore, our best estimate of is a mixed average with a small bias. In Fig. 5 we plot for different couplings and . We found that the long-range tail of the momentum distribution decays as a power law as soon as the attractive contact interaction is switched on, and it is fitted well by . This is the same behavior found by Minguzzi at al.minguzzi:2002 () for the tail of the momentum distribution in a one dimensional gas of hard point-like bosons (Tonks gas) inside a harmonic trap. The power law decay is due to the cusp conditions of the wave function, a consequence of the attractive (or repulsive) contact interactions. The same decay has been shown recently by Santachiara and Calabrese for a one dimensional gas of impenetrable anyons.santachiara () As one can see in Fig. 5, the momentum distribution does not discriminate between different phase separated states, but from the height of its tail one can get information on the effective interaction strength in the trap.

On the other hand, the pair momentum distribution carries fingerprints of the FFLO state. We define the pair two-body density matrix as

(23) | |||||

where the pair is chosen by selecting and finding as the position of the closest unlike-spin particle. () is the center of mass of the () pair. This is an unambiguous and generic way to define the pairs in a continuous system without introducing explicitly a characteristic length. Once is computed, the pair momentum distribution function is evaluated via Eq. 19, with replaced by . Also the normalization is based on the same formula as in Eq. 21, but with replaced by , i.e. the number of pairs in the system. The results for are presented in Fig. 6 for , , and various polarizations. A clearcut signature of the FFLO state is the presence of two peaks at , which signals a pairing with non zero total momentum given by the difference between the Fermi surfaces of up and down spin particles. Since the trap breaks the translational invariance, the Fermi momenta are approximated by , where is the effective spread of the cloud. The location of the peaks in Fig. 6 follows closely the relation .

Fig. 7 displays the dependence on at fixed polarization . At weak coupling () the peaks are short, and in the range the pair momentum distribution is almost flat. On the other hand, at the two peaks at are taller, and the atoms in the cloud have a quite strong FFLO character. In the intermediate case analyzed here, namely , the shows a fluctuating behavior, with peaks at but also at . By computing the Fourier transform (Eq. 19) of constrained in the region with (and plotted in the inset of Fig. 7), we proved that the FFLO peaks come from the inner shell of the cloud, while the central peaks come from the wings. Indeed, around the system undergoes the transition from fully polarized to fully paired (BCS) wings (see Fig. 4). Both the FFLO pairing in the center and the BCS pairing in the outer shell are weak. Going from the center to the edge, the crossover between the two gives rise to this fluctuating pattern in the pair momentum distribution.

### iv.2 Finite temperature

We performed finite temperature PIMC calculations for , , and , with , and . From the DMC analysis and the TF results we know that for the phase separation is between an FFLO center and a fully polarized wing, while for and the outer shell is fully paired. Here we study how the phase separation evolves with the temperature.

In Fig. 8 we plot the local polarization profile for . At a temperature above , the fully polarized shell becomes partially polarized with polarization progressively reduced as the temperature is raised. At the polarization is quite homogeneous throughout the cloud.

The finite temperature results at (), reported in Fig. 9 (10), show a more complicated pattern. Up to a temperature of (), the polarization profile has a dip around (), which evolves into a fully paired region (with zero polarization) below (). When the dip is present, the polarization in the outer part of the cloud increases until a fully polarized state occurs at the very edge of the cloud. This esternal region involves only of the atoms in the cloud, and so it will be hard to detect in the experiment. At higher temperatures, when the dip is gone, the fully polarized edge becomes unstable. For we have a situation analogous to the case with , with a partially polarized state homogeneously spread over the cloud.

A feature common of all couplings is a crossover temperature above which the phase separation disappears and the whole cloud becomes partially polarized. This temperature depends on the polarization, since it scales roughly as , whereas its dependence on the coupling is quite weak, at least until . We found for , and for with . In case of fully polarized wings, one can define another critical temperature based on the stability of the pairs at the edge of the cloud. The fermions dissociate into a partially polarized fluid if the temperature is raised above a threshold proportional to , as it depends on the binding energy of each pair.

As in the DMC calculations, the PIMC simulations suffer from the ergodicity problem mentioned before. This is an issue of the strictly one dimensional Monte Carlo sampling, where the exchange between particles or pairs is suppressed at strong coupling. For this particular system, we found that this issue is less severe in the DMC simulations than in the PIMC ones. Improved moves or working in the grand canonical ensemble could alleviate or solve the problem. Here, however, we could afford simulations with at quite strong coupling, namely up to in the DMC, and in the PIMC framework. Therefore our analysis has not been limited, since we were able to study the intermediate-strong coupling regime, where the fully paired region shows up at not-so-small polarizations.

## V Conclusions

In this paper we have studied the properties of a system of 1D trapped fermions interacting via an attractive contact potential, by performing unbiased DMC and PIMC simulations. We showed that the trap induces the system to phase separate, in accordance to previous calculations done with the TF approximation.orso (); liu_hu_drumm_prl () We compared the TF solution with our exact results, by studying their dependence on the coupling and system size. We found that the DMC phase boundaries are in a good agreement with the TF approximation for large number of particles () even at strong coupling. This is reasonable since the TF approach is applicable when , and the size of the cloud is large compared to the axial oscillator length. However, the local density approximation requires also a slow variation of the mean interparticle spacing, namely , condition which is not fulfilled by the system particularly at the edge of the cloud, where the TF solution slightly overestimates the size of the outer fully paired region. For smaller systems we found that the agreement deteriorates. Indeed, at small a partially polarized state extends until the edge of the cloud even at intermediate coupling and small spin imbalance, in contrast with the TF approximation, which predicts a phase separation with fully paired wings.

As pointed out in a number of previous worksbatrouni (); tezuka (); feiguin (); rizzi () the pair momentum distribution reveals the signature of the FFLO inner shell with two symmetric peaks centered at . It could be measured in the experiment by quickly removing the confining lattice and sweeping the magnetic field to the BEC side of the Feshbach resonance. This would bind the pairs into molecules, whose momentum distribution could be measured by means of time-of-flight imaging techniques.greiner () The tails of the one-body momentum distribution strongly depend on the interaction strength and decay as , due to the non analyticity of the exact wave function as a consequence of the delta function interaction. This could be another interesting quantity to measure in the experiment, although it does not carry information on the actual pairing in the system, but only on the effective strength.

The finite temperature PIMC calculations done with , , and three different couplings (, , and ) show a crossover from a phase separated state to a partially polarized state at for , and for . For the strongest coupling analyzed with this method (), the fully paired region disappears already at , leaving a dip in the polarization density profile which becomes shallower and shallower as the temperature increases. This quite low temperature (the experiments will access temperatures down to ) puts some limitation on the visibility of the fully paired phase at weak-intermediate coupling. One has to go to stronger coupling to make the pairs in the outer shell stable up to higher temperatures, being the binding energy proportional to .

We discussed the ergodicity issues in the Monte Carlo sampling at strong coupling, and we identified the problem with the difficulty of unpaired fermions to pass through strongly bound pairs. The same slowing down in the equilibration might occur also in the experiment for ultra thin traps ( and ). On the other hand, at weak coupling there are few particles involved in the fully paired region, which is very narrow, and the signal coming from this phase should be very difficult to detect in the experiment. We suggest that a good spot to detect the transition from the fully paired to the fully polarized state is for effective couplings in the range , where the system features a fully paired region up to of spin imbalance, the size of the two phase separated states is quite large, and the equilibration is not hard to reach.

The information coming from this paper could be combined with the recent work by Parish et al.,huse () who studied the quasi one dimensional effects coming from the array of tubes in the transverse direction by means of a tight-binding model, in order to find out the parameters of the Hamiltonian which maximize the visibility of the FFLO state. In principle, one can extend the present work by explicitly including an array of tubes in the QMC simulations to study quasi one dimensional effects present in the experimental setup.

###### Acknowledgements.

We thank R. Hulet, D. Huse, A. Minguzzi, and D. Trinkle for fruitful discussions. We acknowledge support in the form of the NSF grant DMR-0404853 and from ARO Award W911NF-07-1-0464 with funds from the DARPA OLE Program and computer resources from NCSA .### References

- L. N. Cooper, Phys. Rev. 104, 1189 (1956).
- G. Sarma, J. Phys. Chem. Solids 24, 1029 (1963).
- W. V. Liu and F. Wilczek, Phys. Rev. Lett. 90, 047002 (2003).
- M. M. Forbes, E. Gubankova, W. V. Liu, and F. Wilczek, Phys. Rev. Lett. 94, 017001 (2005).
- K. B. Gubbels, M. W. J. Romans, and H. T. C. Stoof, Phys. Rev. Lett. 97, 210402 (2006).
- H. Müther and A. Sedrakian, Phys. Rev. Lett. 88, 252503 (2002).
- A. Sedrakian, J. Mur-Petit, A. Polls, and H. Müther, Phys. Rev. A 72, 013613 (2005).
- P. F. Bedaque, H. Caldas, and G. Rupak, Phys. Rev. Lett. 91, 247002 (2003).
- P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964).
- A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP 20, 762 (1965).
- R. Casalbuoni and G. Nardulli, Rev. Mod. Phys. 76, 263 (2004).
- J. Singleton et al., J. Phys. Condens. Matter 12, L641 (2000).
- H. A. Radovan et al., Nature 425, 51 (2003).
- C. F. Miclea et al., Phys. Rev. Lett. 96, 117001 (2006).
- M. W. Zwierlein, A. Schirotzek, C. H. Schunck, and W. Ketterle, Science 311, 5760 (2006).
- M. W. Zwierlein, C. H. Schunck, and A. S. W. Ketterle, Nature 442, 54 (2006).
- Y. Shin, M. W. Zwierlein, C. H. Schunck, A. Schirotzek, and W. Ketterle, Phys. Rev. Lett. 97, 030401 (2006).
- C. H. Schunck, Y. Shin, A. Schirotzek, M. W. Zwierlein, and W. Ketterle, Science 316, 867 (2007), http://www.sciencemag.org/cgi/reprint/316/5826/867.pdf.
- G. B. Partridge, W. Li, R. I. Kamar, Y.-a. Liao, and R. G. Hulet, Science 311, 503 (2006), http://www.sciencemag.org/cgi/reprint/311/5760/503.pdf.
- G. B. Partridge et al., Phys. Rev. Lett. 97, 190407 (2006).
- M. M. Parish, F. M. Marchetti, A. Lamacraft, and B. D. Simons, Nature Physics 3, 124 (2007).
- X.-J. Liu, H. Hu, and P. D. Drummond, Phys. Rev. A 76, 043605 (2007).
- M. Gaudin, Phys. Lett. 24A, 55 (1967).
- M. Takahashi, Prog. Theor. Phys. 44, 348 (1970).
- M. Takahashi, Prog. Theor. Phys. 46, 1388 (1971).
- X. W. Guan, M. T. Batchelor, C. Lee, and M. Bortz, Phys. Rev. B 76, 085120 (2007).
- K. Yang, Phys. Rev. B 63, 140511(R) (2001).
- G. G. Batrouni, M. H. Huntley, V. G. Rousseau, and R. T. Scalettar, Phys. Rev. Lett. 100, 116405 (2008).
- M. Rizzi et al., arXiv:cond-mat/0712.3364v1.
- A. Lüscher, R. M. Noack, and A. M. Läuchli, arXiv:cond-mat/0712.1808v2.
- B. Paredes et al., Nature 429, 277 (2004).
- H. Moritz, T. Stöferle, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 91, 250402 (2003).
- H. Moritz, T. Stöferle, K. Günter, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 94, 210401 (2005).
- G. Orso, Phys. Rev. Lett. 98, 070402 (2007).
- H. Hu, X.-J. Liu, and P. D. Drummond, Phys. Rev. Lett. 98, 070403 (2007).
- G. Xianlong and R. Asgari, Phys. Rev. A 77, 033604 (2008).
- A. E. Feiguin and F. Heidrich-Meisner, Phys. Rev. B 76, 220508 (2007).
- M. Tezuka and M. Ueda, Phys. Rev. Lett. 100, 110403 (2008).
- G. Xianlong et al., Phys. Rev. Lett. 98, 030404 (2007).
- H. P. Büchler, G. Blatter, and W. Zwerger, Phys. Rev. Lett. 90, 130401 (2003).
- R. A. Molina, J. Dukelsky, and P. Schmitteckert, Phys. Rev. Lett. 99, 080404 (2007).
- X.-J. Liu, H. Hu, and P. D. Drummond, arXiv:cond-mat/0711.4390v1.
- P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, J. Chem. Phys. 77, 5593 (1982).
- D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
- D. M. Ceperley, J. Stat. Phys. 63, 1237 (1991).
- M. Calandra Buonaura and S. Sorella, Phys. Rev. B 57, 11446 (1998).
- M. Olshanii, Phys. Rev. Lett. 81, 938 (1998).
- T. Bergeman, M. G. Moore, and M. Olshanii, Phys. Rev. Lett. 91, 163201 (2003).
- R. G. Hulet, private communication.
- S. Sorella, Phys. Rev. B 64, 024512 (2001).
- M. Casula, C. Attaccalite, and S. Sorella, J. Chem. Phys. 121, 7110 (2004).
- K. S. Liu, M. H. Kalos, and G. V. Chester, Phys. Rev. A 10, 303 (1974).
- R. P. Feynman, Statistical Mechanics (Benjamin, New York, 1972).
- B. Gaveau and L. S. Schulman, J. Phys. A 19, 1833 (1986).
- D. M. Ceperley, Phys. Rev. Lett. 69, 331 (1992).
- E. J. Mueller, Phys. Rev. B 72, 075322 (2005).
- G. Xianlong, M. Polini, R. Asgari, and M. P. Tosi, Phys. Rev. A 73, 033609 (2006).
- M. H. Kalos, J. Comput. Phys. 1, 257 (1966).
- M. Minguzzi, P. Vignolo, and M. P. Tosi, Phys. Lett. A 294, 222 (2002).
- R. Santachiara and P. Calabrese, arXiv:cond-mat/0802.1913.
- M. Greiner, C. A. Regal, J. T. Stewart, and D. S. Jin, Phys. Rev. Lett. 94, 110401 (2005).
- M. M. Parish, S. K. Baur, E. J. Mueller, and D. A. Huse, Phys. Rev. Lett. 99, 250403 (2007).