Enhanced visibility of the Fulde-Ferrel-Larkin-Ovchinnikov state in one dimensional Bose-Fermi mixtures near the immiscibility point

Based on the matrix product states method, we investigate numerically the ground state properties of one-dimensional mixtures of repulsive bosons and spin-imbalanced attractive fermions, the latter being in the Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) state, where Cooper pairs condense at a finite momentum $k=k_{FFLO}$. We find that the visibility of such a state is dramatically enhanced as the strength of the boson-fermion repulsion is increased and the system becomes close to the phase-separation point. In particular, self-induced oscillations with wave-vector $2k_{FFLO}$ appear in both the fermion total density and the boson density profiles, leaving sharp fingerprints in the corresponding static structure factors. We show that these features remain well visible in cold atoms samples trapped longitudinally by a flat-bottom potential.

Based on the matrix product states method, we investigate numerically the ground state properties of one-dimensional mixtures of repulsive bosons and spin-imbalanced attractive fermions, the latter being in the Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) state, where Cooper pairs condense at a finite momentum k = kF F LO . We find that the visibility of such a state is dramatically enhanced as the strength of the boson-fermion repulsion is increased and the system becomes close to the phase-separation point. In particular, self-induced oscillations with wave-vector 2kF F LO appear in both the fermion total density and the boson density profiles, leaving sharp fingerprints in the corresponding static structure factors. We show that these features remain well visible in cold atoms samples trapped longitudinally by a flat-bottom potential.
According to the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity, electrons with opposite spin bind into bosonic pairs, which then condense in the state of zero center-of-mass momentum, leading to macroscopic phase coherence and vanishing electrical resistance. An intriguing question is the possible coexistence of superconductivity with a spin imbalance, which destabilizes the BCS mechanism. Several exotic superfluid states have been proposed theoretically, including the breached pair or Sarma state [1][2][3], states with deformed Fermi surfaces [4,5] and Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) state [6,7], to name a few.
The FFLO state is characterized by the condensation of Cooper pairs at finite momentum k F F LO , corresponding in real space to a spatially modulated order parameter. As a consequence, excess fermions, which are detrimental to superconductivity, are stored preferentially at the nodes of the pairing field, leading to an oscillation in the spin density with wave-vector 2k F F LO (corresponding to two nodes per wavelength). The FFLO state is currently being investigated in a variety of physical systems, including layered organic [8], heavy fermion [9] and iron-based [10] superconductors, hybrid superconducting-ferromagnetic structures [11] and quark-gluon plasma [12]. To date, its experimental evidence relies mostly on thermodynamic measurements, although recent NMR spectra of organic superconductors are consistent with a periodic modulation of the spin density [8].
In this Letter, we unveil a general mechanism to significantly enhance the visibility of the FFLO state by placing the spin-imbalanced Fermi system in contact with a superfluid Bose gas and driving the mixture close to the phase-separation point. Specifically, we find that robust self-induced density modulations with wave-vector 2k F F LO appear both in the boson and in the fermion total density profiles, leading to sharp kinks in the corresponding static structure factors (much sharper than the original kink in the magnetic response).
Let us start by discussing the ground state properties of the homogeneous system (effects of the trap will be presented in the second part of the Letter). We describe the Bose-Fermi mixture by the following lattice Hamiltonian: where the first two terms represent the Fermi Hubbard model, where c † iσ (n iσ ) is the local creation (density) operators for fermions with spin component σ =↑, ↓, t f is the tunneling rate of fermions between neighboring sites arXiv:1911.03448v1 [cond-mat.quant-gas] 8 Nov 2019 and U f (U f < 0) is the strength of the attractive interaction between fermions with opposite spin. The third and fourth terms yield the Bose-Hubbard model, where b † i (n bi ) is the bosonic creation (density) operator at site i, while t b and U b (U b > 0) are the corresponding tunneling rate and the onsite repulsion strength. Bosons and fermions are coupled by repulsive contact interactions of strength U bf (U bf > 0), as described by the last term of the model (1). For definiteness, in the following we assume that bosons and fermions have equal tunneling rates and fix the energy scale by setting t f = t b = 1.
Our numerical results are based on the Density Matrix Renormalization Group (DMRG) method, expressed in terms of Matrix Product States (MPS) (for a review see Ref. [68]). Specifically, we use the mps-optim code of the ALPS library [69]. We consider a chain of L = 120 sites with open boundary conditions, containing N ↑ = 40 spinup fermions, N ↓ = 32 spin-down fermions and N b = 60 bosonic atoms (the choice of half filling for bosons is not crucial for our results). We set U b = 4 and U f = −2 while varying the Bose-Fermi coupling U bf . To ensure proper convergence, we allow a maximum occupancy of four bosons per-site along with bond dimension up to 4000 and 80 sweeps. For the above sets of parameters we find that phase separation effects start to appear around U bf = 3.6 [70].
In order to characterize the 1D FFLO phase, we study the pair momentum distribution (PMD) where is the superconducting correlation function in the singlet channel. Figure 1 shows the PMD for U bf = 0 (green circles) and for U bf = 3.4 (red diamonds), where the system is at the verge of phase separation. The 1D FFLO state is signaled by a sharp peak in the PMD at a finite momentum where k F ↑ = πN ↑ /L and k F ↓ = πN ↓ /L are the Fermi momenta of the majority and minority spin components. We see from Fig.1 that close to phase separation the FFLO peak remains well visible and is even slightly enhanced. The main effect of the repulsive boson-fermion interaction is the appearance of Cooper pairs with large momentum, close to the edge of the Brillouin zone. Since k n pair k = i n i↑ n i↓ , this results in an increase of the number of doubly occupied sites from 17.4 to 23.3. We have verified numerically that this behavior is not intrinsic to the FFLO state, as it also occurs for equal spin populations, N ↑ = N ↓ . The PMD could be measured by projecting the Cooper pairs into deep molecular states and performing time-of-flight experiments [71], as previously done for three-dimensional Fermi superfluids, although interactions effects during the expansion can complicate the picture.
In coordinates space, the FFLO state appears as a selfgenerated spatial modulation of the superconducting correlation function, Eq.(3), superimposed to the algebraic decay typical of 1D systems. Our numerical results for P pair ij are displayed in the inset of Fig.1 for the aforementioned values of U bf . Interestingly, as U bf increases, we see that this quantity smoothens out and becomes more symmetric with respect to the center of the chain, which might favors the observation of the exotic superfluid with interferometric techniques.
Let us now investigate the effect of the Bose-Fermi interactions on the density profiles. With the spectacular recent advances in quantum gas microscopy [72][73][74], it is now possible to measure these local observables with high resolution. In Fig.2 we plot the distributions of the spin density n i↑ − n i↓ (red circles), the total fermionic density n i↑ + n i↓ (blue squares) and the boson density n ib (green diamond) in the absence of coupling (panel a) and close to phase separation (panel b). Although in 1D systems true long range order is absent due to the strong quantum fluctuations, the crystalline structure in the local spin density can appear in chains of finite size.
We see in Fig.2 (a) that for U bf = 0 the FFLO oscillation is barely visible, as the attraction between fermions is relatively weak, U f = −2. Indeed, previous DMRG [20]  The three data curves correspond to the spin density (red circles), total fermionic density (blue squares) and boson density (green diamonds), respectively. The particle numbers are the same as in Fig.1. The panel c shows the amplitudes of the FFLO oscillations in the three density profiles as a function of 1/L for U bf = 3.4 and for three different values of the system size, L = 60, 90 and 120 at fixed particle densities. The intraspecies interaction strengths are U f = −2 and U b = 4.
and quantum Monte Carlo [37] studies have observed a clear periodic modulation only for relatively strong attractions between fermions, say |U f | 5. A Fourier analysis of the data [70] shows that the fermion total density displays oscillations with wave-vectors 2k F ↓ and 2k F ↑ , while the 2k F F LO modulation is nearly absent. In contrast, the boson density displays oscillations with dominant wave-vector k = π.
As the strength of the Bose-Fermi repulsion increases, we see in Fig.2 (b) that a clear periodic structure with the expected wavelength, 2π/(2k F F LO ) = 15, progressively emerges in the spin-density profile (notice indeed that this quantity has a minimum when the absolute value of the superconducting correlation in Fig.1 reaches a local  maximum). Remarkably, when the system becomes sufficiently close to the immiscibility point, the FFLO order is imprinted in both the total fermionic density and the boson density profiles, as shown in Fig.2 (panel b).
We see in Fig.2 (b) that the oscillation in the boson density is out of phase by a factor of π with respect to the other two density profiles, due to the repulsive Bose-Fermi interaction. In contrast to the uncoupled case (panel a), we also notice that at the edge of the system there are more fermions than bosons. Hence, in the middle of the chain the total fermionic (bosonic) density profile oscillates around a lower (higher) value.
Since the modulations displayed in Fig.2 (b)signal the FFLO state, we expect that their amplitudes A vanish in the limit of infinite chains. In order to verify this point, we have performed analog calculations for system size L = 60 and L = 90, keeping the densities N σ /L and N b /L unchanged. We compute A by fitting the numerical data in the central region of the chain (45 < i < 75) with a function n(x) = A cos(2k F F LO x + φ). The result is shown in Fig.2 (c) as a function of 1/L. All the data curves are well fitted by straight lines with approximately zero intercept, thus confirming our claim.
Next, we discuss the fingerprints of the FFLO state in the static structure factors for the aforementioned observables. These quantities can be accessed experimentally in cold atoms samples via Bragg scattering [75] or quantum polarization spectroscopy [37]. They are defined as where for fermions we distinguish between the spin response, corresponding to the operator O = m = n ↑ − n ↓ , and the total density response, where O = n = n ↑ + n ↓ . Fig.3 displays the momentum dependence of the three static structure factors for three different values of the Bose-Fermi coupling, U bf = 0 (green circles), 3.0 (blue squares) and 3.4 (red diamonds). For U bf = 0 we see that both S f m and S f n exhibit similar shapes with kinks at k = 2k F ↓ and k = 2k F ↑ (dashed lines). This can be understood by noticing that in a noninteracting Fermi gas the two static structure factors coincide, S m = S n , and are given by [37] S m (k) = |k|/π for 0 < |k| < The inclusion of a moderate attraction between fermions slightly smears the two kinks and decreases (increases) the overall scale of the magnetic (total density) structure factor, as indicated by the green curves in the panels (a) and (b). In contrast, the density response of the Bose gas is smooth and increases monotonously as the momentum increases, as shown by the panel(c) of the same figure. As U bf increases, we see that a kink progressively stands out in the spin response at k = 2k F F LO , signaling the FFLO state. At the same time both the total density and the boson density responses become approximately flat for k > 2k F ↓ . By further approaching the immiscibility point, these two quantities develop a sharp kink at k = 2k F F LO , as shown in Fig.3 (b) and (c). In particular, we see that the two kinks are significantly more pronounced than the corresponding one in the magnetic response, providing a further tool to detect the FFLO state. We emphasize that the results shown in Fig.3 have negligible finite size effects [70]. This is due to the fact that the static structure factors measure density correlations at different sites, which remain finite in the thermodynamic limit.
The effects of a trapping potential acting on bosons and fermions can be taken into account through the generalized Hamiltonian where V and p are positive numbers. Since the FFLO wave-vector is fixed by the value of the local spin density, the latter should stay approximately constant over a wide region of the trap for the corresponding density modulations to be observable. Hence a confinement sharper than harmonic, p > 2, is generally required. Flat-bottom potential can be realized optically using a digital micromirror device (DMD) (for instance p = 16 in Ref. [76]). In Fig. 4 we display the calculated density profiles for a mixture with N ↑ = 20, N ↓ = 14 and N b = 42 in a trap with p = 12 and L = 94, using the same values for the interactions strengths as in Fig.2 (b). We see that the characteristic FFLO modulations in the density profiles of bosons and fermions are well visible in the central region of the chain. We have also checked [70] that the kinks in the static structure factors S f n and S b remain very sharp in the trapped case, although the wings of the mixture in Fig. 4 are not in the FFLO phase.
In summary, we have unveiled a general mechanism to strongly enhance the visibility of the 1D FFLO state by coupling the spin-imbalanced Fermi gas to a Bose superfluid. As the mixture is brought close to the phaseseparation limit, strong FFLO modulations appear in the density profiles of the two species, resulting in sharp peaks in the corresponding static structure factors. We have shown that these features remain well visible when atoms are confined longitudinally by flat-bottom traps, already for moderate values of the trap power p. Our results open new prospects to directly observe the exotic superfluid phase in cold atoms systems via high resolution in-situ imaging and Bragg spectroscopy.
While the density modulations in Fig.2 (b) vanish for 1D mixtures of infinite size, a weak interchain hopping is expected to establish true long range (FFLO) order in higher dimensions. This, in turn, could drive the formation of bosonic supersolid phases in Bose-Fermi mixtures near phase separation.
We acknowledge C. Salomon and F. Chevy for fruitful discussions. This work was supported by ANR (grant SpifBox). M.S. acknowledges funding from MULTIPLY fellowships under the Marie Sk lodowska-Curie COFUND Action (grant agreement No. 713694).
In this Supplementary Material we provide additional information about our results. In Sec.I we identify a generic set of values of the model parameters, where the Bose-Fermi mixture is close to the phase-separation limit (but still homogeneous). In Sec.II we provide a detailed Fourier analysis of the oscillations in the density profiles in order to identify the dominant wave vectors. In Sec.III we present additional results for the static structure factors, including the finite-size scaling analysis for homogeneous systems and the effect of a flat-bottom trapping potential in the axial direction.

I. BOSE-FERMI MIXTURE NEAR PHASE SEPARATION
As explained in the main text, we fix the values of the intraspecies interactions strengths and we increase the value of U bf until phase separation effects start appearing in the calculated density profiles. For this search we use shorter chains of L = 60 sites, with re-scaled particle numbers N ↑ = 20, N ↓ = 16 and N b = 30. The obtained results are shown in Fig.1.
At U bf = 3.4 (panel a), we see that the density profiles of bosons and fermions display FFLO oscillations around a well defined value, implying that the mixture is homogeneous, except for a small region near the edge of the chain, due to the choice of open boundary conditions. Similar results are found also for U bf = 3.5 (panel b), although the inhomogeneities of the local densities at the periphery are somewhat more pronounced.
For U bf = 3.6, however, there is a drastic change in the density profiles of the two species, as shown in Fig.1(c). In particular we see that inhomogeneous regions appear in the bulk of the mixture, signaling the proximity to the phase separation regime.
Since we are interested in mixtures close to the immiscibility point, but still homogeneous, in the main text we present results for U bf = 3.4.

II. FOURIER ANALYSIS OF THE DENSITY PROFILES
In this section we present a Fourier analysis of the density profiles of bosons and fermions plotted in Fig.2(a) and 2(b) of the main text. In order to avoid boundary effects, we consider only the central region of the chain, corresponding to m = 44 sites (approximately one third of the total length). For each type of profiles (spin density, fermion total density and boson density) we proceed as follows. Let x i denotes the i-th element of the truncated data set. To avoid uninteresting peaks in the Fourier spectrum we first normalize our data by setting where µ and σ denote, respectively, the mean and the standard deviation of the data. We then compute the In each column, the three panels from top to bottom, display results for the spin density, total fermion density, and boson density. The three dashed lines from left to right, refer to 2kF F LO , 2k F ↓ and 2k F ↑ momenta, respectively.
discrete Fourier transform of the data, characterized by the set of m coefficients X r defined as where r is an integer. Since our data are real, it is enough to study the first half of the Fourier coefficients. Each integer value r is associated to a bin in Fourier space centered at wave-vector k r = 2πr/m and of size 2π/m. Finally we plot the absolute value of the normalized coefficients X r /m, as a function of k r . Figs.2(a)-(c) show the obtained results for the density profiles in the absence of the Bose-Fermi coupling, U bf = 0, corresponding to Fig.2(a) of the main text. While in the aforementioned figure it is difficult to identify the underlying wave-vectors of the oscillations, the picture considerably simplifies in the Fourier space. The spin-density (panel a) shows three distinct peaks, corresponding to 2k F F LO = 0.133π, 2k F ↓ = 0.533π and 2k F ↑ = 0.667π, as given in the main text.
The latter two peaks appear in the spectrum of the total fermionic density also (panel b), while the FFLO peak is barely visible. The oscillations in the boson density profile display instead dominant wave-vector at 2πN b /L = π, as shown in Fig.2 (panel c). All the above results are consistent with the Luttinger liquid theory.
Let us now present the Fourier analysis of the local density distributions in Fig.2(b) of the main text, corresponding to U bf = 3.4. This figure shows that all density profiles exhibit a periodic pattern with the same wave-vector 2k F F LO , signaling the FFLO state. This translates into a strong peak in the Fourier spectrum of the corresponding data, as shown in Figs.2(d)-(f).
Importantly, the 2k F ↑ response in the fermion density distributions (panel d and e) is absent close to the phase separation limit. Similar considerations hold for the boson density distribution (panel f), where the peak at k = π is washed out by the boson-fermion repulsion.
We emphasize that the outcomes of the above Fourier analysis of the density profiles are fully consistent with the results for the corresponding static structure factors (SSF) presented in Fig.3 of the main text. In particular the peaks of the Fourier transform observed in the local densities appear as kinks, namely discontinuities of the derivative as a function of momentum, in the density correlation functions.

III. STATIC STRUCTURE FACTORS
A. Finite-size scaling analysis In Fig. 2(c) in the main text we have shown that for U bf = 3.4 the amplitudes of the FFLO modulations of the local densities of bosons and fermions scale as the inverse of the system size L. This reflects the quasilong range order of the 1D FFLO state. To complement these results, in this section we present the corresponding scaling analysis of the SSFs of the density correlations, for the same set of parameters.
In Fig.3 we display the momentum dependence of the SSFs for the fermion spin density (S f m ), fermion total density (S f n ) and the boson density (S b ) for L = 60 (red circles), 90 (blue squares) and 120 (green diamonds) sites. We see that finite-size effects in the SSFs are very small, as anticipated in the main text. that atoms are confined longitudinally by a smooth boxlike potential with moderate power p = 12. The figure shows that the central part of the system is an FFLO state, as witnessed by the characteristic modulations of the density profiles. Notice in this respect that the value of the FFLO momentum is fixed by the average value of the spin density in the middle of the chain, yielding 2k F F LO = 0.218. Fig.4 displays the corresponding results for the SSFs, obtained using the same set of parameters values as in Fig.4 of the main text. We see that the SSFs of the bosonic and the total fermionic densities still exhibit sharp kinks, as stated in the manuscript.