{'fraction_non_alphanumeric': 0.04230718696623308, 'fraction_numerical': 0.0183506328559541, 'mean_word_length': 4.41218298174923, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase separation of trapped spin-imbalanced Fermi gases in one-dimensional optical lattices

We calculate the density profiles of a trapped spin-imbalanced Fermi gas with attractive interactions in a one-dimensional optical lattice, using both the local density approximation (LDA) and density matrix renormalization group (DMRG) simulations. Based on the exact equation of state obtained by Bethe ansatz, LDA predicts that the gas phase-separates into shells with a partially polarized core and fully paired wings, where the latter occurs below a critical spin polarization. This behavior is also seen in numerically exact DMRG calculations at sufficiently large particle numbers. Unlike the continuum case, we show that the critical polarization is a non monotonic function of the interaction strength and vanishes in the limit of large interactions.

We calculate the density profiles of a trapped spin-imbalanced Fermi gas with attractive interactions in a one-dimensional optical lattice, using both the local density approximation (LDA) and density matrix renormalization group (DMRG) simulations. Based on the exact equation of state obtained by Bethe ansatz, LDA predicts that the gas phase separates into shells with a partially polarized core and fully paired wings, where the latter occurs below a critical spin polarization. This behavior is also seen in numerically exact DMRG calculations at sufficiently large particle numbers. Unlike the continuum case, we show that the critical polarization is a non monotonic function of the interaction strength and vanishes in the limit of large interactions.

I. INTRODUCTION
The physics of population imbalanced Fermi gases in one dimension is currently attracting substantial interest since in these systems, the one-dimensional (1D) counterpart of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [1,2] can be realized. This inhomogeneous superfluid state of fermions was suggested long time ago, yet its unambiguous observation in condensed matter systems has turned out to be quite challenging [3][4][5][6]. Ultracold atomic gases offer a perfect avenue to search for exotic pairing states such as the FFLO state, due to the clean experimental conditions and the tunability of interactions and imbalance [7][8][9]. First experiments for three-dimensional (3D) systems at by Ketterle and co-workers [10,11] and by Hulet and co-workers [12,13] as well as theoretical work for 3D [14] indicate that the chances of observing the FFLO state in 3D are slim since this phase appears to be stable only in a narrow window of the BCS-BEC crossover phase diagram.
A promising direction to detect the FFLO state is to confine atoms in 1D configurations. In this limit the Fermi surface reduces to two points, k F ↑ and k F ↓ , forcing all Cooper pairs to share the same wave-vector Q FFLO = k F,↑ − k F,↓ . For systems of fermions with attractive interactions -usually described by the Hubbard model or the associated continuum model -it has indeed been shown by means of analytical and numerical techniques that any finite imbalance drives the system into the 1D FFLO state, which has algebraically decaying pair-pair correlations in the s-wave channel, modulated with the FFLO wavevector Q FFLO [15][16][17][18][19][20][21][22][23][24].
In the more general case of a two-channel model that also accounts for the formation of composite molecules, the FFLO phase still emerges on the BCS side in the imbalanced case, yet eventually gives room to a Bose-Fermi mixture phase [25,26]. The FFLO phase also exists in mass-imbalanced systems in 1D [27][28][29][30][31], which could be realized in Li-K mixtures (see Ref. [32]). It should further be possible to drive the FFLO state in a globally balanced mixture by using spindependent potentials [33] or spin-dependent optical lattices [34]. Furthermore, FFLO physics has been studied in 3D arrays of coupled 1D systems [23,24,35], the two-leg ladder geometry [36], and the possibility of inhomogeneous pairing states in multi-band systems has been explored in Ref. [37]. Finally, we mention that the spectral function of attractively interacting fermions with spin-imbalance has been calculated and discussed in Ref. [38].
Experimentally, the simplest observables that can be measured in trapped Fermi gases are the density profiles [10][11][12][13]39] of the two spin components [10][11][12][13]39]. On the basis of the exact Bethe ansatz solution for the homogeneous system, it has been independently predicted by Hu et al. [16] and one of the authors [17] that the trapped gas phase-separates into two shells: a partially polarized core, and either fully paired or fully polarized wings. Very recently, this theoretical prediction has been addressed and verified in an experiment with ultracold fermions in 1D wave-guides at Rice [40].
In the case of an optical lattice, the shell structure of spinimbalanced fermions was studied with the density matrix renormalization group (DMRG) method in Refs. [18,19,41] and with quantum Monte Carlo (QMC) in Ref. [20]. In Ref. [18] by some of us, clear evidence for the presence of FFLO correlations of the partially polarized core has been reported on systems with up to 80 fermions. In another DMRG work by Tezuka and Ueda [19], a phase diagram for the shell structure has been proposed, with lower critical polarizations for the emergence of the fully polarized wings. In Ref. [19], the presence of the fully polarized wings was deduced from the presence of certain local maxima in the density difference profile in the outer regions of the trap. However, neither of these studies observed the fully paired wings, which have been seen in QMC calculations, in the continuum case [22].
The purpose of this work is to clarify this open point, namely the presence of fully paired wings in an optical lattice. As a result, numerically exact DMRG simulations indeed show the existence of fully paired wings at low spin polarizations. To that end we perform large scale DMRG [43,44] calculations of the 1D Hubbard model with (effectively) up to 320 fermions in a harmonic trap, which is comparable to the particle numbers typically encountered in experiments [7,45]. We carefully compare the DMRG results with those obtained starting from the exact solution of the homogeneous lattice gas, including the effect of the trap via the local density ap- proximation (LDA), as previously done for the continuum case [16,17].
Our study unveils that the fully paired wings are present below a critical spin polarization, which we calculate as a function of effective density and interaction strength. In contrast with the continuum case, in the strong coupling regime, the fully paired wings shrink in size as the attraction strength increases. As we shall see, this effect is due to the enhanced inertia of the molecules with respect to unpaired fermions.
The outline of this work is as follows. In Sec. II, we shall introduce the model and discuss in detail how we obtain the density profiles of the trapped gas, using either DMRG simulations or LDA applied to the exact equation of state for the homogeneous system. Our main results are contained in Sec. III, where we compare the two methods and discuss density profiles as well as the critical polarization.

A. Hamiltonian
Our simulations are performed for the 1D Hubbard model with an attractive onsite interaction U < 0: where c † ℓσ creates a fermion with spin σ =↑, ↓ at site l, n ℓσ = c † ℓσ c ℓσ , n ℓ = n ℓ↑ + n ℓ↓ is the local density, and t is the hopping matrix element. We define x = ia, where a is the lattice spacing, set to unity. The local polarization is defined as s i = n i,↑ − n i↓ , the filling factor is n = N/L where N σ = L i=1 n iσ and N = σ N σ . Analogously, we define the densities for spin up and down as n σ = N σ /L.
Our main interest will be in the case of trapped fermions for which we can write the Hamiltonian as where in the second line, V is the trapping potential. To characterize the trapped system, a suitable choice of variable is the total polarization p = (N ↑ − N ↓ )/N and the effective density In the following we fix the energy scale by setting t = 1.

B. Bethe ansatz and Local Density Approximation
The density profiles of the two component Fermi gas can be obtained starting from the exact Bethe ansatz solution of the homogeneous system (1), taking into account the effect of the trap via the local density approximation.

Homogeneous system
The exact ground state energy per site E(n ↑ , n ↓ ) of the gas in the thermodynamic limit is given by ([42, 46, 47] and references therein) where ρ 1 (k) and ρ 2 (k) are spectral functions satisfying the following coupled integral equations: Here K ν (x) = 2νu/(ν 2 u 2 + x 2 ), with u = |U |/4. The integration limits B, Q are related to the total density n = n ↑ +n ↓ and the density difference s = n ↑ − n ↓ by the conditions Equations (4)-(6) are only valid below unit filling, n ≤ 1. For higher concentrations the ground state energy can be obtained through a particle-hole transformation h σ = c † σ and is given by E(n, s) = E(2 − n, s) + U (n − 1). From the energies, we calculate the averaged chemical potential and the effective magnetic field: Interpreting Eq. (7) as a change of variables, one obtains the phase diagram of the Hamiltonian (1), as shown in Fig. 1 for a specific value of the interaction U = −5t (compare [42]). We can distinguish different phases for a low filling n = N/L ≤ 1: (i) the vacuum (V), corresponding to s = n = 0, (ii) the unpolarized phase (ED) characterized by equal densities of majority and minority fermions, (iii) the partially polarized phase (PP), corresponding to 0 < s < max(n, 2 − n); (iv) the fully polarized phase with n < 1 (FP 1 ), and (v) the fully polarized phase with n = 1 (FP 2 ). The remaining phases at higher densities correspond to situations where at least one of the component forms a band insulator.

Trapped system
The density profiles of the gas in a shallow trap can be calculated via the local density approximation, assuming that the system is locally homogeneous: where µ 0 and h 0 are the chemical potential and the magnetic field at the trap center. In the second line of Eq. (8), we have exploited that the external potential is the same for both components, corresponding to vertical trajectories in the phase diagram of Fig. 1. The constants µ 0 , h 0 are fixed by the overall spin populations N σ through the normalization conditions By inverting Eq. (7) we see that the normalization conditions (10) can be written as N = n[µ 0 − V x 2 , h 0 ]dx, with an analogous expression for the density difference. Introducing In this article we mainly consider configurations where the core of the cloud is partially polarized, corresponding to trajectories starting from the PP phase. Then, LDA predicts the existence of two shell structures: a partially polarized (PP) core with either (i) fully polarized (FP 1 ) or (ii) fully paired (ED) wings (see the vertical lines in Fig. 1 for an illustration). The dotted line in Fig. 1 corresponds to the special case h = √ U 2 + 16/2 − 2, where the PP phase extends to the edge of the cloud, in analogy with the continuum case [16,17]. This yields a critical value of the spin polarization p = p c , separating the regime with fully paired (p < p c ) from the one with fully polarized wings (p > p c ). At larger densities, more shell structures are possible, e.g., a band insulator of the majority spins in the center with two surrounding shells [18].

C. DMRG simulations
The density profiles of the gas in the trap can be calculated to a high accuracy using the DMRG algorithm. We use a standard implementation for ground state calculations [44] and up to 600 states. Occasionally, it was necessary to perform an unusually large number of sweeps (sometimes more than 20) to obtain converged results. For a Hamiltonian with sitedependent parameters, it is convenient to grow the system in a linear fashion, say from the left to the right, during the warmup part of the DMRG algorithm. The largest system size used in this work is L = 300 sites for p = 0.5, U = −8t and N = 160 particles. In that case, typical discarded weights in the center of the system are smaller than 10 −10 . In the two following subsection we provide details on the set-up used in the DMRG runs, such as the choice for the trapping potential and a finite-size scaling analysis.

Asymmetric vs symmetric traps
We here consider two types of a harmonic confining potential parameterized by a constant V , namely one that is symmetrically placed in the optical lattice [see Eq. (3)] and one that traps particles in the left part of the chain: We will next show that for the purpose of studying the behavior in the wings, the symmetric Eq. (3) and the asymmetric set-up Eq. (11) yield the same quantitative behavior. Note that in the figures, we shall display the results for the symmetric trap shifted by L/2 lattice sites to the left with respect to Eq. (3). The advantage of using the asymmetric trap is that it allows us to reach much lower densities and larger particle numbers without resorting to large system sizes, and at less computational costs. Note that for the same trapping amplitude V , the particle numbers in the symmetric trap Eq. (3) are twice as large as the ones in the corresponding asymmetric set-up. Figure 2 shows the comparison between the symmetric and the asymmetric trap for the parameters of Ref. . The open boundary that is present when we use the asymmetric set-up clearly induces strong additional oscillations in the density profile, both at p = 0 and p = 0.5. The main point is, however, that away from the left boundary and towards the edge of the cloud, the density and polarization profiles computed from either the symmetric or the asymmetric trap coincide. In the case of p = 0.5, the system enters the fully polarized wings at the same distance from the trap center, independently of the particle number. This justifies the use of the asymmetric trap in our analysis of the shell structure, which we shall use in the remainder of this work.

Scaling analysis in the particle number
As a next step we study the finite-size scaling of the density and polarization profiles. More precisely, we increase the particle number but keep the effective density ρ eff = N √ V fixed [48,49] (the effective density is given for the asymmetric trap unless stated otherwise). Site-dependent quantities should then be plotted versus x/ξ, where ξ = (V /t) −1/2 is the characteristic length scale for a given effective density.
The results of such an analysis are presented in Fig. 3 for the parameters of Fig. 2, with particle numbers N = 40, 80, 160 in the asymmetric trap. As expected, all curves computed for different particle numbers fall on top of each other. We further see that the amplitude of the boundary-induced oscillations becomes much weaker the larger the polarization or the larger N is and therefore, the presence of the oscillations does not pose a problem for our analysis.

A. Fully paired wings
In this section we will focus on N = 160 and N = 80 particles, with a trapping potential of V = 3.125 10 −5 t. The change from fully paired (ED, equal densities) wings to fully polarized (FP 1 ) wings is expected to occur at small polarizations, and we hence concentrate on p < ∼ 0.2. Figure 4(a) contains the results for ρ eff = 20 √ 0.002, i.e., the original density studied in Ref. [18], but with a much larger number of particles. While at p = 0.1, we still clearly observe FP wings, at p = 0.05, the local polarization practically vanishes at the edge of the cloud. Therefore, for these parameters, the critical polarization for the change in the shell structure must be 0.05 < p c < 0.1, and further analysis shows that it is p c ≈ 0.07. Therefore, if there are too few particles in the trap, say N = 40, then the fully paired wings are very difficult to observe since such a p c would correspond to one of the smallest possible polarizations. Moreover, effects due to a finite particle number become more relevant for small N and/or tighter traps, and LDA does not necessarily quantitatively apply in this regime (see also Ref. [22]). In fact, numerically exact DMRG calculations do not yield evidence for fully paired wings at, for instance, N = 40 and V = 0.005t [18].
By going to half as many particles yet while keeping the trapping amplitude constant at V = 3.125 10 −5 t and thus decreasing the effective density, the fully paired wings become better visible, as we demonstrate in Fig. 4(b) for p = 0.05. There, the local polarization vanishes at a distance from the center of the trap close to x/ξ ∼ 0.5, whereas the cloud extends up to x/ξ ∼ 0.6. Clearly, the volume of the fully paired wings increases going to a lower effective density.
Another way of influencing the volume of the fully paired wings is by going to smaller values of the interactions |U |. The respective results for U = −4t and the same density and particle number as in Fig. 4(a)  By contrast, whenever we enter the fully polarized wings, the double occupancy goes to zero. Summarizing, so far, we have illustrated that the fully paired wings indeed show up in 1D optical lattices at small polarizations in shallow traps, similar to the case of the associated continuum model [16,17,22]. Decreasing the effective density stabilizes this shell structure in the sense that, the relative volume of the fully paired wings increases [compare Figs. 4(a) and (b)].
Next, we compare the density profiles calculated by DMRG and by LDA, where the latter makes use of the exact phase diagram for the 1D Hubbard model (Fig. 1, see [42]). This is shown in Fig. 4(b) for U = −8t and a polarization of p = 0.05 at which the wings are fully paired. The agreement between LDA and DMRG is very convincing which is expected in the case of a shallow trap as considered here (V = 3.125 10 −5 t). Most importantly for us here, both methods show the fully paired wings and the cloud's spatial extension for both the density and density difference agree quite well. From the figure, some known shortcomings of LDA are obvious: first, LDA fails to reproduce the correct decrease of the density profile at the edge of the cloud. This implies that within LDA it is possible to simply define the radius of the cloud by its extension, whereas for a finite particle number this quantity fluctuates. Second, the incommensurate oscillations in the density induced by the lattice are not captured by LDA either [50]. While these discrepancies are not crucially relevant to the analysis carried out here, there is another feature that is absent in the LDA, namely the 2Q FFLO oscillation in the density difference. This is obvious since LDA assumes that the gas is locally homogeneous whereas in DMRG calculations and in experiments, one deals with a finite particle number.
The 2Q FFLO modulation is due to the fact that majority fermions reside in the nodes of the quasi-condensate [51], which has previously been demonstrated with DMRG [18,19,26]. In other words, the incommensurate 1D FFLO state is accompanied with a spin-density wave (see Ref. [52]), as is illustrated in the inset of Fig. 4(b). This feature, which is visible in a trap, where the quasi-condensate is pinned, could in principle be observed in experiments with trapped twocomponent Fermi gases [26,51]. This highlights the importance of using exact numerical techniques such as DMRG and QMC that can take the trap into account exactly to complement analytical methods. Nevertheless, LDA combined with the exact solution of the phase diagram [16,17] has been very successful in predicting the shell structures in one dimension, which for the continuum case has been observed experimentally [40].

B. Critical polarization
We now present our theoretical predictions for the critical polarization p c associated to the change in the shell structure of the density profiles, from fully paired (p < p c ) to fully polarized (p > p c ) wings. For a shallow trap, this quantity depends only on the values of the interaction U and the effective density ρ eff = N √ V . In Fig. 5 we plot the critical polarization as a function of |U | comparing DMRG (squares) with LDA (circles), for V = 3.125 10 −5 t and two different particle numbers N = 160 and N = 80. From the DMRG simulations, we estimate p c as the point at which, by visual inspection, we can clearly resolve regions in the wings with s i = n i . Due to the finite particle number, these results carry uncertainties of δp c ≈ ±0.0125, 0.03 for N = 80, 160, respectively.
The agreement between the DMRG results and LDA is obviously very good. For |U | < ∼ 5t, the critical polarization is an monotonically increasing function of the interaction, in analogy with the continuum model [16,17]. However, for stronger interactions the critical polarization in a lattice reaches a maximum value and then decreases at large |U |, in contrast with the continuum case [16,17].
The physical explanation is that in this regime the kinetic energy of the pairs is suppressed, since their effective mass becomes increasingly large (∼ |U |/t 2 ) compared to that of single fermions (1/t). As a consequence, the energy needed to break a pair is given by |U |, independently of the particle density. This means that the boundary between the fully paired (ED) and the partially polarized phase in Fig. 1 (the line connecting the points A,B, and C counterclockwise) reduces to the dotted vertical line for U → −∞, implying the absence of fully paired wings in a trap. In Fig. 6 we plot LDA results for the critical polarization versus the effective density for different values of the inter-action. We see that the polarization increases by decreasing the effective density and saturates at a maximum value of p max c (U ). Then, if one sends U to zero as well, this maximum will approach p max c (U ) → 0.2 as is illustrated in the inset of Fig. 6, consistent with the findings reported for the continuum limit [16,17].

IV. SUMMARY AND DISCUSSION
In summary, we investigated the density profiles of trapped spin-imbalanced Fermi gases in a 1D optical lattice. We demonstrated that, when performed on sufficiently large systems, numerically exact DMRG simulations are fully consistent with LDA results. In particular, we reported numerical evidence for the fully paired wings that were predicted to exist at small polarizations [16,17]. However, differently from the continuum case, we found that the upper critical polarization p c for the appearance of the fully paired wings is a non monotonic function of the interaction strength U at a fixed filling. This behavior is a direct consequence of the increased inertia of the pairs in the strong coupling regime (|U | ≫ 4t).