Fulde-Ferrell-Larkin-Ovchinnikov superfluidity in one-dimensional optical lattices

Spin-polarized attractive Fermi gases in one-dimensional (1D) optical lattices are expected to be remarkably good candidates for the observation of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase. We model these systems with an attractive Hubbard model with population imbalance. By means of the density-matrix renormalization-group method we compute the pairing correlations as well as the static spin and charge structure factors in the whole range from weak to strong coupling. We demonstrate that pairing correlations exhibit quasi-long range order and oscillations at the wave number expected from FFLO theory. However, we also show by numerically computing the mixed spin-charge static structure factor that charge and spin degrees of freedom appear to be coupled already for small imbalance. We discuss the consequences of this coupling for the observation of the FFLO phase, as well as for the stabilization of the quasi-long range order into long-range order by coupling many identical 1D systems, as in quasi-1D optical lattices.


I. INTRODUCTION
Multi-component attractive fermionic systems with unequal masses, densities or chemical potentials have attracted continued interest for many decades in several fields of physics ranging from high-energy 1,2 to condensed matter 2,3 and, more recently, atomic physics 4,5,6 . The interplay between pairing and density imbalance of the different fermion species leads to a rich scenario, which includes the possibility of various exotic superconducting states 7 . In this context, the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase 8 has recently attracted a great deal of interest from both the experimental and the theoretical community 1,2,3,4,5,6 . In the FFLO phase Cooper pairing occurs between a fermion with momentum k and spin ↑ and a fermion with momentum −k + q (q = 0) and spin ↓. As a result, the superconducting order parameter becomes spatially dependent. Originally, the most favorable systems for the observation of the FFLO phase were predicted to be clean superconducting films in the presence of an in-plane (i.e. Zeeman) magnetic field, above the so-called Clogston-Chandrasekhar limit 9 . Nevertheless, despite the fact that the original prediction dates back to more than thirty years ago, the FFLO phase has been very elusive to detect.
The experimental realization of interacting trapped Fermi gases with population imbalance 4,5 has renewed the hope of observing the FFLO, thus stimulating an intense theoretical activity 6,10 . So far, most of the theoretical analysis has focused on 3D cold atomic systems. However, as in the case of solid-state superconductors, the region of phase diagram where the FFLO phase has been found to be stable is quite small 2,6 . On the other hand, quasi-one dimensional or strongly anisotropic systems (such as coupled chains, heavy-fermion, organic, high-T c , and CeCoIn 5 superconductors) are believed to be good candidates for the realization of the FFLO phase 2,3,11,12 . Since the dimensionality of cold atomic systems can be easily tuned, and indeed cold atoms have already been successfully trapped in 1D geometries 13 , it seems natural to consider these low dimensional systems as the ideal candidates to observe non-homogeneous superconductivity of the FFLO type.
Many important results are available on the properties of spin-polarized 1D Fermi systems with attractive interactions, which have been obtained by different methods and techniques. These include the Bethe-Ansatz solutions of certain exactly solvable models like the Hubbard or Gaudin-Yang models 14,15,16,17 , as well as different types of numerical approaches such as densitymatrix renormalization-group (DMRG) 18,19,20 and quantum Monte Carlo (QMC) 21 , or field theoretical techniques like bosonization 12,18 .
Very recently Orso 14 and Hu et al. 15 have studied the phase diagram of harmonically-trapped 1D polarized Fermi gases by combining the exact solution of the Gaudin-Yang model with a local-density approximation. Mean-field theory has also been applied by Liu et al. 22 , although it is known that it has a number of limitations 23 in 1D, particularly as far as paring correlations are concerned. DMRG has been employed by Feiguin et al. 19 and Tezuka and Ueda 20 , and QMC by Batrouni et al. 21 to investigate the pairing correlations in the spin-polarized ground state of the attractive Hubbard model in the presence of harmonic trapping. Previously, Yang 12 used bosonization to study the pairing correlations and the phase diagram of a single 1D Fermi system as well as an array of weaklycoupled 1D Fermi systems in the presence of a Zeeman field.
Yang's analysis is valid only close to a continuous magnetic-field-driven transition from a uniform BCS phase and an FFLO phase, which he assumed to belong to the commensurate-incommensurate universality class 24 . Another important assumption in Ref. 12 is that charge and spin degrees of freedom are decoupled at low energies for small polarization. However, this scenario does not apply to the Hubbard model away from halffilling 17,25,26 nor to the Gaudin-Yang model 16 , which are the relevant models for current 1D cold atomic systems. As we show in this work, charge and spin degrees of freedom are indeed coupled already for small polarization, which leads to important differences as compared to the scenario described by Yang. We also numerically demonstrate that the pairing correlation function exhibits prominent oscillations with a wave number equal (up to finite size corrections, see below) to the difference of Fermi wave numbers, q FFLO = |k F↑ − k F↓ |, as predicted by FFLO theory in 1D 12 and in agreement with a number of Luttinger-type theorems 27,28 . Thus, the finitewave-number oscillations in the pairing correlation function can be regarded as due to the excess of n ↑ − n ↓ unpaired majority-spin fermions 12 . This is because in 1D the Fermi wave number and the density are proportional to each other: k Fσ = πn σ . On the other hand, this relationship is no longer linear in dimensionality higher than one, and in this case the oscillations in the order parameter are related to the center-of-mass momentum q of the Coopers pairs. These observations seem to indicate that in the 1D case it is not entirely clear whether there is a strict close parallelism with higher dimensional FFLO, and in some respects the system can be also understood as a a coupled Bose-Fermi mixture of spin-singlet pairs (the bosons) and unpaired fermions 15,26 .
The paper is organized as follows. Section II presents the Hamiltonian that we use to describe the system of physical interest, while Section III reports and discusses our main numerical results. Our main conclusions are briefly reported in Section IV.

II. THE MODEL
We consider a two-component mixture with a total of N fermionic atoms loaded in a 1D optical lattice of L sites (the lattice constant is taken to be unity). The fermions are assumed to interact via attractive on-site interactions, whose strength can be tuned e.g. by means of a Feshbach resonance. Sufficiently away from resonance(s), this system is modeled by the attractive Hubbard model: where t is the hopping parameter,ĉ † ℓσ (ĉ ℓσ ) is the creation (destruction) fermion operator in the ℓ-site (ℓ ∈ [1, L]), σ =↑, ↓ the pseudospin-1/2 index (in experiments this labels the two different atomic hyperfine states of the mixture), U > 0 is the strength of the on-site Hubbard attraction,n ℓσ =ĉ † ℓσĉ ℓσ , and N σ = σn ℓσ . In order to simulate the effect of an external trapping potential, open-boundary conditions (OBC) breaking translational symmetry will be used (these are indeed the most suitable conditions for the DMRG treatment of the above model 30 ). Our calculations are performed in the canonical ensemble, and the results apply only to lattices away from half-filling, that is, when N = L. In the calculations, the spin polarization δ = (N ↑ − N ↓ )/(N ↑ + N ↓ ) was varied by decreasing N ↓ while keeping constant the number of "background" up-spin atoms N ↑ , from N ↓ = N ↑ (the unpolarized case, i.e. δ = 0) all the way down to N ↓ = 0 (the fully polarized case, i.e. δ = 1).
In the unpolarized case (δ = 0), all fermions pair into spin singlets due to the attractive on-site interaction. This yields a gap to all spin excitations and therefore spin-spin correlations decay exponentially with distance. Singlet superconducting and charge-density wave correlations exhibit a slower decay (of power-law type in the ground state of a thermodynamically large system), being singlet superconducting correlations the ones that dominate at long distances in systems away from halffiling 29 . The aim of this work is to study the nature of superfluidity for 0 < δ < 1 as a function of the dimensionless ratio U/t (in the fully polarized case, where δ = 1,Ĥ describes a system of N = N ↑ noninteracting fermions). The expectation values ... of all operators below are understood to be taken over the ground state ofĤ.

III. NUMERICAL RESULTS AND DISCUSSION
Due to the OBC (or, in general, to any external potential that breaks the Bloch translational invariance of the lattice) the spin-resolved site occupation profiles, n ℓσ = n ℓσ , exhibit Friedel oscillations. This is illustrated in Fig. 1. In the unpolarized δ = 0 case the Friedel oscillations in n ↑ are in-phase with those in n ↓ giving rise to large-amplitude atomic-density waves 31,32 in the total site occupation n ℓ↑ + n ℓ↓ . As it is clear from the top panel of Fig. 1, in the general δ = 0 case the total site occupation displays N ↓ maxima associated with the formation of N ↓ spin-singlet pairs that are delocalized over the lattice. In the bottom panel of Fig. 1, we show the local spin polarization n ℓ↑ − n ℓ↓ (which could be measured through phase-sensitive optical imaging 5 ). For small δ (see e.g. the plot for N ↑ = 20 and N ↓ = 16) the local spin polarization displays N ↑ − N ↓ maxima corresponding to the number of fermions that are left unpaired. With increasing δ though the spatial dependence of the local spin polarization becomes more complicated: the amplitude of the oscillations in the bulk becomes indeed smaller, thus making it hard to clearly identify N ↑ − N ↓ maxima. These, however, are not distinctive and unambiguous signals of FFLO pairing. We thus proceed below to present a study of pairing correlations: the model described in Eq. (1), in fact, cannot sustain any true long-range order 29 in 1D, i.e. the ground-state expectation value of the pairing operator ∆ ℓ =ĉ ℓ↓ĉℓ↑ is zero. In the unpolarized case and for an extended system, the correlation function of the pairing operator C ℓℓ ′ = ∆ † ℓ∆ ℓ ′ decays with a power law |ℓ − ℓ ′ | −1/Kρ at large distances, where 1 ≤ K ρ ≤ 2 is an interaction-dependent Luttinger-liquid dimensionless parameter 29 . In the top panel of Fig. 2 we illustrate our DMRG results for the spin-polarization dependence of C ℓℓ ′ =L/2 at U/t = 5, which measures real-space superfluid correlations between the site ℓ ′ = L/2 (the center of the trap) and all the other sites. For δ = 0 the powerlaw decay of the C ℓℓ ′ =L/2 for |ℓ − L/2| ≫ 1 is clearly visible. For finite δ, instead, the superfluid correlator is characterized by a distinctive oscillatory character 12 and a very simple nodal structure with exactly N ↑ − N ↓ zeroes. We have carefully checked that the long-distance decay of C ℓℓ ′ =L/2 is still power-law, signaling quasi-long range superfluid behavior also at finite δ. A careful analysis of the oscillatory character of C ℓℓ ′ =L/2 can be done by means of the Fourier transform of the pairing correlator where ϕ m (ℓ) = [2/(L + 1)] −1/2 sin (k m ℓ) [with k m = πm/(L + 1), m = 1 . . . L] are the eigenstates of the hopping term in Eq. (1). The mode with zero wave number is excluded from the allowed k m values due to the OBC. The lowest energy mode corresponds to k 1 . The diagonal part of the matrix C(k m , k m ′ ) will be simply denoted by C(k m ) ≡ diag{C(k m , k m ′ )} = C(k m , k m ).
In the bottom panel of Fig. 2 we plot the difference ∆C(k m ) = C(k m )− C (0) (k m ) between C(k m ) and its value in the noninteracting gas (i.e. at U/t = 0), C (0) (k m ) 33 . At δ = 0 C(k m ) possesses a very narrow peak at k 1 (see inset in the bottom panel of Fig. 2), signaling quasilong-range superfluid order of the conventional BCS type. For a finite δ, instead, C(k m ) has a local minimum at k 1 and a single well-defined peak appears at a wave number q FFLO = k 1 + |k F↑ − k F↓ |, where k Fσ = πN σ /(L + 1) are the spin-resolved Fermi wave numbers. The peak at q FFLO in the Fourier transform of the pairing correlator, which is a direct consequence of the simple real-space nodal structure illustrated in the top panel of Fig. 2, is a clear-cut signal of FFLO pairing.
The DMRG data shown in Fig. 2 refer only to a single value of U/t = 5. We now turn to illustrate the dependence of the superfluid correlation functions on U/t. In the top panel of Fig. 3 we illustrate the dependence of ∆C(k m ) on the interaction strength U/t for a fixed spin polarization δ = 25%. On decreasing U/t the quasi-longrange FFLO order (i.e. the height of the peak at q FFLO ), which is emphatically strong for large U/t, survives all the way down to the weak coupling regime. This can be quantified better by analyzing the size of the anomaly Γ at k m = q FFLO , which is measured by the difference between left and right (discrete) derivatives of C(k m ) evaluated at q FFLO , In the bottom panel of Fig. 3 we plot Γ as a function of U/t ≤ 5. In this range Γ decreases in a smooth fashion to its noninteracting value (i.e. Γ = 0) as U/t is decreased to zero. In other words, for every finite δ, C(k m ) tends uniformly and smoothly to its noninteracting value C (0) (k m ) as U/t is decreased towards zero. For sufficiently large values of U/t the FFLO phase can also be characterized by the peak visibility defined by This quantity is plotted in an inset to the bottom panel of Fig. 3. Before concluding, we would like to illustrate the behavior of the density-density, spin-spin, and mixed density-spin static structure factors, S nn (k m ), S mm (k m ), and S nm (k m ). These are defined by the sum over all frequencies of the corresponding dynamic structure factors 34 that can be in principle measured through Bragg spectroscopy or Fourier sampling of time-of-flight images 35 . In practice, S nn (k m ), S mm (k m ), and S nm (k m ) are calculated from the following equations: wheren ℓ =n ℓ↑ +n ℓ↓ andm ℓ =n ℓ↑ −n ℓ↓ . In Figs. 4-5 we show the dependence of S nn (k m ), S mm (k m ), and S nm (k m ) on U/t for a slightly asymmetric system with N ↑ = 20 and N ↓ = 18 (δ ∼ 5%). We remind the reader that in the unpolarized δ = 0 case S nn (k m ) has a peak at k m = 2k F↑ = 2k F↓ that signals real-space atomic-density waves 31,32 . In the spin-polarized case, this peak splits into two peaks at 2k F↑ and 2k F↓ . This is clearly visible in Fig. 4 in the static density-density structure factor (top panel), which presents a double-peak structure slightly below k m = 3π/4 (2k F↑ ≈ 2π/3 and 2k F↓ ≈ 3π/5 for the system parameters in this figure). This double-peak structure is not so visible in the magnetic structure factor S mm (k m ), most likely because magnetic correlations near 2k F↑ and 2k F↓ are still quite suppressed by the superfluid correlations, at least in the weakly polarized case (in the unpolarized case they are completely suppressed by the pairing gap). From Fig. 5 we note that S nm (k m ) is non-zero even at small k m , thus indicating that spin and charge degrees of freedom are coupled at long wavelength even for a small imbalance.

Experimental signatures of the FFLO phase
The most direct way to detect FFLO pairing would be to measure the pairing correlation function C ℓℓ ′ . We would like to remark here that this correlation function is, at least in principle, measurable via interferometric schemes 36 in which two atomic wave packets are coherently extracted from the gas at different positions and then are mixed by a matter-wave beam splitter. The atom counting statistics in the beam splitter output channels has been shown 36 to reflect the spatial dependence of C ℓℓ ′ . The oscillations of the pairing correlations will also leave a detectable signature in the noise correlations 37,38 , G ↑↓ (k, k ′ ) = n k,↑nk ′ ,↓ − n k,↑ n k ′ ,↓ , wheren k,σ measures the number of fermions with momentum k and spin σ in a time-of-flight experiment. With increasing spinpolarization, in fact, the peak at k = −k ′ = (k, 0, 0) 39,40 [here (1, 0, 0) is the direction along the axis of the 1D system] will shift to a finite relative momentum (see e.g. the work by Yang in Ref. 10 and the very recent DMRG calculation by Lüscher et al. 41 ).
However, it is worth pointing out that the strength of noise signal in a strictly-1D system will be strongly affected by finite-size and temperature effects. This is because in 1D the order is not long-range but quasi-long range, and therefore the slowest decay exhibited by correlations (like the pairing correlations) is a power law. Thus, in order to enhance the strength of the experimental signal for FFLO, it would be desirable to couple many 1D systems, as in a tight 2D optical lattice (arrays of "atomic quantum wires") 42 , so that the quasilong range FFLO order can become true long-range order. The phase diagram of many coupled 1D systems has been worked out in Ref. 12, where the author found that, at small polarization δ, true long-range 3D FFLO order will occur when the Luttinger liquid parameter for the charge excitations, K ρ , is larger than 3/2. In such a case the low temperature properties of the system are dominated by hopping of pairs (i.e. Josephson coupling) rather than by single-particle hopping (the latter would turn the system into an anisotropic Fermi liquid, which could in turn become unstable to the FFLO state under appropriate conditions 43 ).
However, the analysis of Ref. 12 assumed the transition from the unpolarized to the polarized case to belong the the commensurate-incommensurate universality class. This assumption implicitly neglects the coupling between charge and spin degrees of freedom at low energies, which is known to modify the behavior of physical observables at the transition 25,26 . In this work this coupling has been demonstrated to exist also at long wavelengths by an explicit numerical evaluation of the mixed static structure factor S nm (k m ) in a weakly polarized system (see Fig. 5). Thus, the phase diagram depicted in Fig. 1 of Ref. 12 seems not appropriate for coupled 1D Hubbard (or Gaudin-Yang 16 ) models, and the transition to long-range order will not take place in general for K ρ = 3/2 and may in general depend on the system parameters (i.e. U/t and the lattice filling for the Hubbard model).
To the best of our knowledge, a quantitative phase diagram of coupled 1D systems lacking spin-charge separation has not yet been calculated. Furthermore, it is worth noticing that in cold atomic systems with short range interactions, another important factor must be taken into account, namely the relative strength of the pair hopping when compared to the single-particle hopping. The strength of the latter is given by t ⊥ (t ⊥ ≪ ε F , where ε F is the Fermi energy, for the analysis based on coupled Tomonaga-Luttinger liquids to hold), where t ⊥ is the hopping amplitude between two neighboring 1D systems. However, in absence of long-range interactions pair hopping can be only generated by (virtual) single-particle tunneling events, which at lowest order, yield a (Josephson) coupling strength of order t 2 ⊥ /∆ σ , where ∆ σ is the spin gap. The phase diagram predicted in Ref. 12 is the result of a calculation which only compares the scaling dimensions of the pair hopping and single-particle hopping operators, and thus does not take into account the microscopic details of the coupling between 1D systems.
Thus, the stabilization of FFLO long-range order by a weak coupling between 1D systems in the FFLO phase (the extreme anisotropic limit that could not be accessed by the authors of Ref. 43) does not seem easily achievable. In turn, the most likely scenario for arbitrary polarization is that the single-particle hopping will control the physics at low temperatures, and the system will behave as a spin-polarized normal Fermi liquid, which, in turn, could become unstable towards 3D FFLO ordering 43 under appropriate conditions.

IV. CONCLUSIONS
In summary, we have shown how ultracold spinpolarized two-component Fermi gases confined in 1D optical lattices are FFLO superfluids whose pairing correlation functions are characterized by a power-law decay and a simple nodal structure. However, we have also shown that charge and spin degrees of freedom appear to be coupled already for a small value of the spin polarization. Finally, we have commented on the impact of this coupling on the detectability of true long-range FFLO order arising from Josephson coupling between 1D systems.