Unconventional pairing in few-fermion systems tuned by external confinement

We study the ground-state properties of a two-component one-dimensional system of a few ultra-cold fermions with attractive interactions. We show that, by ramping up an external potential barrier felt by one of the components, it is possible to induce regions of exotic superfluid phases, characterized by a tunable finite net momentum of the Cooper pair, without changing the overall spin populations. We show that these phases, which are the few-body analogs of the celebrated Fulde-Ferrell-Larkin-Ovchinnikov state, can be distinguished by analyzing a specific two-particle correlation encoded in the noise correlation function. Our theoretical results can be addressed in current experiments with cold atoms confined in spin-selective optical traps.


I. INTRODUCTION
Ultra-cold atoms provide an exceptional experimental platform to simulate condensed matter systems in a controlled way [1,2]. One of the most spectacular collective phenomena in solids is superconductivity, where electrons with opposite spin and momenta bind into Cooper pairs, due to an effective attractive interaction mediated by the crystal vibrations. In the presence of a mismatch between the two spin populations, generated for instance by an applied Zeeman field, the conventional Bardeen-Cooper-Schrieffer (BCS) pairing mechanism becomes unstable, as some electrons will inevitably end up without partners. The spin-imbalanced system can nevertheless remain superconducting, at the price of adopting a new pairing mechanism harnessing the excess fermions. A well known example of superconductivity coexisting with a partial spin polarization is the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [3,4], where Cooper pairs condense at a finite momentum. This implies that the associated order parameter becomes spatially modulated, with the excess fermions sitting at the nodes of the wave, where they are less detrimental to superconductivity.
The modulated phase is currently investigated in different physical systems, including one-and twodimensional organic superconductors [5], hybrid structures [6] and quark-gluon plasma [7]. Over the last decade atomic Fermi gases have also emerged as a valid alternative to study this exotic state of matter (for a recent review see [8][9][10]) for a number of reasons: i) the two spin states correspond to two hyperfine levels, whose populations are fixed at the beginning of the experiment via a radio-frequency field, without generating vortices; ii) there is no disorder; iii) fermions can be confined in low-dimensional geometries, where the FFLO state is more robust. In particular, the ground state of one-dimensional (1D) systems with attractive contact interactions is known both theoretically [11][12][13] and nu- * Jacek.Dobrzyniecki@fuw.edu.pl merically [14][15][16][17] to be of the FFLO type, for any finite value of the spin imbalance. The phase diagram of a spin-imbalanced attractive 1D Fermi gas in a harmonic trap has been predicted in [18,19] and experimentally verified in Ref. [20]. In particular the density profiles of the two spin components develop a two-shell structure, with the central part being a FFLO phase, while the wings are either fully paired or fully polarized, depending on the overall spin polarization. These predictions were also confirmed by DMRG studies of the Fermi Hubbard model [21]. To date, obtaining direct experimental evidence of the modulated phase with cold atoms is still challenging. Several detection schemes have been discussed, based on the analysis of collective oscillations [22], the sudden expansion of the gas [23][24][25], interaction quenches [26], noise correlations [27,28], spectroscopy measurements [29][30][31][32] and the coupling to a Bose gas [33].
An interesting problem is to investigate the FFLO pairing in 1D Fermi gases in the presence of spin-dependent external potentials, so that the effective Zeeman field, corresponding to the semi-difference between the local chemical potentials of the two spin components, is no longer uniform throughout the atomic cloud. A natural question then arises: can one tune the external confinement to induce different superconducting phases in the Fermi gas, without changing the overall spin populations? To answer this intriguing question, in this paper we investigate theoretically a 1D spin-1/2 system of a few attractively interacting fermions, confined in a box trap with an additional spin-dependent potential barrier at the trap center. Our main object of interest are the pairing correlations present in this system, which can be analyzed by means of the noise correlation distribution. Depending on the local population imbalance, the system behaves either as a BCS superconductor or a FFLO state, distinguishable by a nonzero center-of-mass-momentum of Cooper pairs. We show that by appropriately tuning the height and width of the potential barrier, it is possible to switch the system between different pairing types, even though the atom numbers of both components remain unchanged.
Although superconducting pairing mechanisms are arXiv:2105.12519v4 [cond-mat.quant-gas] 12 Nov 2021 typically studied for bulk many-body systems, here we explicitly focus on few-body systems far from the thermodynamic limit. In this way our work connects to ongoing experiments on few-body samples [34], such as those being currently undertaken in Jochim's group [35][36][37][38][39]. This work is organized as follows. In Sec. II, we describe the model system under study. In Sec. III we examine the pair correlations that arise in the box trap, without the potential barrier, and show how they can be analyzed through noise correlation distributions. In Sec. IV, we describe the effect of the potential barrier, showing how the dominant net pair momentum changes as the barrier parameters (height and width) are modified. In Sec.V, we describe the particular case where the component split by the potential barrier has an odd number of fermions. Sec. VI contains the conclusions.

II. THE MODEL
In this work we consider a one-dimensional system of a few fermions of mass m in two different internal states σ ∈ {A, B}, playing the role of effective spins. Motivated by state-of-the-art experiments with ultra-cold atoms in two hyperfine levels, we assume that the particle numbers N A and N B of the two spin components are fixed integers. We assume that atoms are strongly confined along two orthogonal directions, using for instance a tight twodimensional optical lattices, so that their motion along these directions reduces to zero-point oscillations. Under this assumption the system behaves kinematically as 1D and is described by a second-quantized Hamiltonian of the form where the fermionic field operatorΨ σ (x) annihilates a σ-fermion at position x and obeys the conventional fermionic anticommutation relations, . For convenience, we introduced the single-particle density oper- In the following we take the external potential V σ (x) as with V A = V and V B =0, respectively. It means that the particles are confined in an infinite square well of length 2L, and the component A additionally feels a potential barrier of width 2D and height V in the center of the box. From an experimental point of view, a spin-dependent external potential can be achieved e.g. by using a focused laser beam or magnetic field gradient to induce a spatially localized spin-selective energy shift [40][41][42][43]. The inter-particle interactions are modelled as contact interactions between fermions of opposite species with strength g. The interaction strength g is related to the three-dimensional s-wave scattering length [44,45] and can be tuned by magnetic Feshbach resonances [46,47] or by adjusting the confinement in the transverse directions [44].
For convenience, throughout the rest of this paper we employ dimensionless units, i.e., we express all energies, lengths, and momenta in units of 2 /mL 2 , L, and /L, respectively. In these units the interaction strength is expressed in units of 2 /mL. Without losing the generality of the final conclusions, throughout this paper we set the strength of attractive interactions to g = −5. Importantly, we consider 1D systems of few (up to 12) fermions and assume that the external potential has fixed spatial size L. For this reason, our results cannot be straightforwardly extrapolated to the thermodynamic limit, where the size of the system has to be changed together with the number of particles to keep the average density constant.
To numerically obtain the many-body ground state for the given number of particles and external potential configuration, we first solve the corresponding single-particle eigenproblems for each component separately. Then, we use the lowest-energy eigenorbitals to construct the noninteracting many-body Fock basis {|F i }. Specifically, each Fock state |F i is a product of two Slater determinants of N A and N B orbitals, describing the many-body state of A and B component, respectively. The resultant many-body basis, in general, includes all the possible combinations of single-particle orbitals of the N A + N B fermions. Since the many-body basis grows exponentially with the number of particles, we limit the basis to Fock states which have a non-interacting energy below a properly chosen value E max , according to the recipe given in [48]. This procedure is based upon the assumption that very high-energy Fock states will be only negligibly represented in the ground-state wave function of the system. Then the many-body Hamiltonian (1) is expressed as a matrix in the basis {|F i } and diagonalized using the implicitly restarted Arnoldi method [49]. In this way, the ground state |G is found as its decomposition in the basis {|F i } and used for further calculations. In the end, we confirm that the obtained results do not change quantitatively upon increase of the cutoff energy E max . Thus, the method gives practically exact results, in the sense that the ground-state many-body wave function is known almost exactly, i.e., further increase of the Fock basis does not change the results significantly.

III. PAIRING IN THE TRAPPED SYSTEM
For the attractive inter-component interactions (g < 0), fermions of opposite species form strongly correlated pairs. It has been shown that these pairs may display many different features of the Cooper pairs known from the theory of superconductivity [37,[50][51][52][53]. Identifying the type of pairing that arises (i.e. pairing with zero or nonzero net pair momentum) requires a detailed knowledge of the superconducting correlation function. It has been previously shown [27,28,54] that such information can also be obtained from the two-point noise correlation G between the two components, which is directly experimentally accessible from two-body density measurements [55][56][57]. In the momentum domain, the distribution of the noise correlation is given by where the momentum density operatorsπ . The noise correlation is the difference between the two-particle density distribution, and the product of individual single-particle densities. For a non-interacting system, G(p A , p B ) is zero everywhere. It means that G(p A , p B ) expresses the distribution of correlations forced by inter-particle interactions that cannot be captured by a single-particle description, excluding the spurious correlations that arise from single-particle densities.
To demonstrate that different types of correlated pairs are well-captured by the noise correlations, in Fig. 1 we show the distribution G(p A , p B ) in the absence of the potential barrier (V = 0) and different imbalances ∆N = N B − N A . For all balanced systems (∆N = 0), an enhancement of inter-component correlations (G > 0) is visible along the antidiagonal p A = −p B . It means that the probability of finding a pair of fermions with exactly opposite momenta is enhanced, which is a footprint of the standard Cooper-like pairing mechanism with zero pair momentum. For imbalanced systems (∆N > 0), the situation is different. The region of enhanced correlations is split into two ridges, located along the two dashed lines, corresponding to net momenta p A + p B = ±q 0 , with q 0 having a nonzero value, a hallmark of FFLO pairing.
The net FFLO pair momentum q 0 in the box trap is expected to be equal to the difference ∆p F = p F B − p F A between the Fermi momenta p F σ = N σ π/2 of the two spin components [58] Note that q 0 depends only on the population imbalance ∆N and not on the values N A and N B separately, unlike in non-uniform (e.g. harmonically trapped) systems [28]. The relation (4) is confirmed by Fig. 1: for larger particle numbers, the noise correlation enhancement is concentrated into two clear, narrow maxima located at p A ≈ ±N A π/2, p B ≈ ∓N B π/2 (red spots); in contrast, the probability of Cooper pairing between fermions with momenta pointing along the same direction is strongly suppressed (blue spots). Separately, it is worth noting from Fig. 1 that the intensity of the noise correlations diminishes for higher particle numbers, because in 1D systems interactions effects are reduced as the particle density increases [18]. In particular the relevant dimensionless parameter is γ = gm/( 2 n), where n = (N A + N B )/(2L) is the total particle density (notice that we have reintroduced all physical units for better clarity). The parameter γ can range from the weakly interacting mean-field regime (γ 1) to the strongly correlated regime (γ 1). Within our units convention the interaction parameter then reduces to γ = 2g/(N A +N B ). For instance for N A = N B = 4 we find γ = 1.25 for g = −5, implying that the system is in between the weakly and the strongly interacting regimes.
To identify more clearly the most probable net momentum of the pair, q 0 , we use a method previously proposed in [28]. It involves integrating the noise correlation with an appropriate filtering function F(k): For the filtering function, we choose a simple Gaussian function F(k) = (πw) −1/2 exp(−k 2 /2w 2 ). The width parameter w = 0.5 is of the order of the perpendicular width of the enhanced correlation area. We have checked that the form of Q(q) is not significantly affected by small adjustments of w. Note that Q(q) = Q(−q), due to the symmetry of G. Therefore, for simplification, throughout the rest of this paper, we show its values only for positive q. In Fig. 2 we plot the function Q(q) for systems with different particle numbers. The momentum q at which the function Q takes its maximum value can be identified with the most likely net momentum q 0 of the Cooper pairs in the system. Fig. 2a refers to systems with identical population imbalance ∆N = 2, but varying particle numbers N A and N B (they correspond to the third row in Fig. 1). In all these cases, the maximum of Q(q) falls at the same position q π, in agreement with the prediction (4). Conversely, in Fig. 2b we show Q(q) for N A = 2 and different particle imbalances (corresponding to the second column in Fig. 1). The balanced system (∆N = 0) with BCS pairing is characterized by a clear maximum at q 0. As the particle imbalance increases, the maximum occurs at increasingly higher momenta, in each case very close to the predicted value (4), indicated by the vertical dashed lines. This point is further illustrated in the inset in Fig. 2, showing the most likely pair momentum q 0 , defined through the function Q(q), versus the Fermi momentum mismatch ∆p F , for various particle numbers and imbalances (see the numerical data and explanation in Appendix A for details). All data points fall very close to the dashed straight line, corresponding to q 0 = ∆p F as predicted by (4), confirming that the the most likely net momentum of the Cooper pairs basically coincides with the mismatch between the Fermi momenta across a wide variety of system sizes. Additionally, we have verified that this result does not change qualitatively as the interaction strength is varied, although the intensity of the noise correlation distribution reduces by approaching the weakly interacting regime.

IV. ROLE OF THE INTERNAL BARRIER
So far we have assumed that both spin components feel the same external (flat box) potential. Let us now consider the effects of changing the barrier height V felt solely by the component A. As V increases, the A-fermions are progressively pushed towards the lateral wings until the central region is completely emptied. In the high-barrier limit, the A-fermions effectively experience a symmetrical double-well potential with negligible tunneling between the two wells. Since each separate well has width 1 − D, we see from (4) that the most likely net momentum q 0 of Cooper pairs is given by Here ∆N = N B − N A , where N σ is the number of fermions of component σ found within a given well. The expected value of N σ can be determined by integrating the corresponding density profile n σ (x) = n σ (x) over the the well domain.

A. Barrier with varying height
To demonstrate the effect of changing the barrier, let us consider a system with N A = 4 and N B = 6 fermions. In Fig. 3a, we examine its ground state properties for increasing V , assuming that the barrier width is fixed to D = 1/3. We plot the single-particle densities n A (x) and n B (x), along with the corresponding noise correlation distributions G(p A , p B ) and the function Q(q) reflecting the pair momentum distribution. In the homogenous case (V = 0, top row in Fig. 3a), both densities (left column) are roughly evenly distributed throughout the box. Due to the population imbalance, the noise correlations in the system support the FFLO-like pairing. The pair momentum calculated from (4) is q 0 = π, as seen from the locations of the positive maxima in G (middle column). As V increases, the component A is gradually pushed out of the barrier region. Due to the symmetry of the system, the density n A (x) is evenly split between the two side regions. Meanwhile, the density n B (x) is essentially unchanged (except for the slight modifications due to the attractive interaction). The particular choice D = 1/3 means that approximately one-third of the B population (two fermions) is located within each lateral region. As a result, for high V , the population within the lateral regions becomes balanced. This is additionally supported by the fact that densities n A (x) and n B (x) in these regions become almost identical. Thus, in this regime, the pairs created within the lateral regions are standard BCS Cooper-like pairs with zero net momentum. This phenomenological reasoning is supported by the noise correlations -for the large V case (last row) the maxima are found close to the anti-diagonal p A = −p B . The maxima become gradually more indistinct and stretched along the p A direction as V increases, which can be explained by the increasing uncertainty of p A as the A fermions are squeezed into a smaller space.
The transition between FFLO and BCS pairing can be explained in greater detail by inspection of the function Q(q) (right column in Fig. 3a). For V = 0, Q(q) displays a maximum at q 0 = π, exactly as predicted for this imbalanced system from the difference in Fermi momenta. As V increases, this maximum gradually vanishes, while simultaneously another maximum emerges at q 0 = 0. In particular, for V = 30 one can distinguish two separate maxima at the two locations. This indicates that the change between the FFLO and BCS pairings is not a gradual decrease of q 0 , but rather a direct switch between two distinct pairing mechanisms. A separate effect is that for small barrier heights V , the maximum at q 0 = π shifts towards slightly larger momenta, which can be explained by the fact that the momentum of A-fermions slightly increases due to the higher external potential energy.
From an experimental point of view, it is useful to define an additional measurable quantity that indicates whether the ground state of the system displays BCS or FFLO pairing. For this purpose, we define the dimensionless quantity where θ(z) is the Heaviside step function. If the maximum of Q(q) is located at q = 0, ξ is exactly zero, while if the maximum falls at any other position then ξ > 0. The value of ξ can therefore be interpreted as an indicator for FFLO pairing. In Fig. 3b we show in more details the dependence of ξ on the barrier height for the N A = 4, N B = 6 system. As V increases, ξ gradually diminishes and eventually vanishes around V 35, signaling the transition towards the standard BCS pairing.

B. Tuning the barrier width
The FFLO pair momentum q 0 can also be controlled by tuning the width of the barrier. To demonstrate this, in Fig. 4a we show results for a system with N A = 4, N B = 9 particles obtained by progressively increasing the barrier width D, assuming a fixed barrier height V = 150. The latter choice ensures that the middle region is nearly emptied of A-fermions. The single-particle densities n σ (x) (left column) give an approximate view of the changing population difference in the lateral regions.
The first row shows the case of a system without a potential barrier. As seen from the noise correlation distribution G(p A , p B ) and the function Q(q) (middle and right column in Fig. 4a, respectively), in this case FFLO-like pairs are formed with a nonzero momentum, q 0 = 5π/2. The subsequent rows show the effect of changing D to different values. At the values D = 1/9, 3/9, 5/9, the expected value of N B is a clearly defined integer number (four, three, and two, respectively). Meanwhile, N A remains close to two in all cases. For these values of D, the noise correlation distribution and the function Q(q) show that the most probable pair momentum q 0 changes to the value predicted from (6) indicated by vertical dashed lines (q 0 = 9π/4, 3π/2, and 0, respectively). Additionally, we show the cases of intermediate widths (D = 2/9, D = 4/9) which lie in between the above values. In such cases the value N B is non-integer, implying that the ground state is a superposition of different quantum states with different numbers of B-particles in the left and right wells. For instance for D = 2/9 there could be four B-fermions in the left well and three B-fermions in the right one or the other way round.
These results clearly show that adjusting the barrier widths allows tuning the FFLO momentum q 0 , as well as switching from FFLO pairing to BCS pairing. This point is further illustrated in Fig. 4b, showing the behavior of ξ as a function of D when D is smoothly varied. In particular ξ decreases monotonically as D increases, and ultimately vanishes around D 0.4, marking the dominance of BCS pairing.

V. ODD PARTICLE NUMBER NA
So far we have considered systems with an even particle number N A . In those cases, the introduction of the potential barrier leads (on average) to an equal distribution of A-fermions in the two lateral regions. A more complicated situation occurs for systems with odd N A , as the presence of the barrier leads to unequal numbers of A-fermions in the two wells. As a consequence, the pair momentum q 0 can also take distinct values in the two lateral regions. To demonstrate this, we focus on the balanced case with N A = N B = 5 particles. The singleparticle densities, noise correlation and the function Q(q) for increasing values of V and fixed D = 1/5 are shown in is exactly balanced and the system is characterized by BCS-like pairing with net pair momentum q 0 = 0. For very high V the most probable distribution is that of two A-fermions in one well, and three A-fermions in the other. In contrast, the expected number of B-fermions in either of the two wells is two, due to the chosen barrier width. Thus the imbalance ∆N is different in both wells, resulting in a dominance of two different values of net pair momentum, q 0 = 0 and q 0 = 5π/4, corresponding to a BCS and a FFLO state, respectively. This effect is indeed visible in the noise correlation G(p A , p B ) and the function Q(q), where distinct separate maxima are visible at the predicted values of q (last row in Fig. 5).
For completeness, we also investigate the imbalanced case N A = 5, N B = 7, where the system exhibits FFLO pairing already at V = 0. In Fig. 6 we show the singleparticle densities, noise correlations, and the function Q(q) in this imbalanced system for increasing values of V and a fixed barrier width D = 1/7. At V = 0, the system exhibits a FFLO pairing with a preferred net momentum q 0 = π. For very high V the most probable distribution is again that of two A-fermions in one well, and three A-fermions in the other, while the expected number of B fermions in each well is three. In this regime Q(q) exhibits two separate peaks, q 0 = 0 and q 0 = 7π/6. Therefore in the imbalanced odd-N A case, where the system already exhibits FFLO pairing in the absence of a barrier, tuning the barrier height V from zero leads to two distinct effects: i) it can introduce an additional BCS phase in one of the two wells, and ii) it can modify the value of the FFLO momentum q 0 .
At this point it is valuable to clarify a general structure of the many-body ground state in the odd N A case. Due to the spatial left-right symmetry, the many-body ground state of the system can be expressed as a general superposition |G = (|L + |R )/ √ 2, where |L and |R represent many-body states describing configurations with the extra A fermion placed in the left and right well, respectively. This means that it is not possible to distinguish whether it is the left or the right well that contributes to a particular pairing mechanism. In the limit of very large V , the ground state becomes nearly degenerate with the first excited many-body state of the system |G = (|L − |R )/ √ 2, having an essentially different single-particle momentum distribution. As a result, in experimental practice, a system prepared in this regime may end up in either one or any superposition of states |L , |R .
Let us finally investigate the behavior of the system with N A = 5, N B = 7 as the barrier width D is progressively increased from zero, for a fixed barrier height V = 120. The obtained results are displayed in Fig. 7. The distributions Q(q) for no barrier (D = 0) and for D = 1/7 are as shown previously in Fig. 6, with the preferred FFLO momentum being q 0 = π and q 0 = 7π/6, respectively (as marked by vertical dashed lines in the third column). For D = 3/7, the expected value of N B is a clearly defined integer number (N B = 2) and the two momenta are expected to fall at q 0 = 0 and q 0 = 7π/4. However, it can be seen that the position of the maximum in Q(q) visibly deviates from the predicted q 0 in this case, although it is reasonably close to the predicted value. We additionally show the results for two intermediate values of D (D = 1.5/7, D = 2/7). Since in these cases the expected value N B is not a clearly defined integer, the distributions Q(q) in this case are not straightforward to describe.

VI. CONCLUSION
Few-body cold atom systems represent an intriguing platform to study pairing phenomena. Here we have investigated a one-dimensional system of few attractively interacting spin-1/2 fermions confined in a flat box trap, through the numerical study of the ground state density profiles and noise correlations. We show that by ramping up a central potential barrier felt by one of the two components, and thus restricting pair formation to regions outside the barrier, the system can undergo different pairing mechanisms without changing the overall spin populations. Specifically, solely by adjusting the barrier height and width, the particles in the two wells can form either BCS-like pairs with zero center-of-mass momentum, or FFLO-like pairs with a tunable finite momentum. Moreover, we found that for odd particle numbers both BCS and FFLO type correlations can coexist in different spatial regions of the system, even in the absence of an overall spin imbalance, provided the barrier parameters are tailored appropriately. Our theoretical results are relevant for current experiments using quasi-1D atomic samples with few particles per tube. In these systems the noise correlations have been measured with great accuracy, while the spindependent external potential can be tailored with stateof-the-art optical techniques. Our work therefore provides a promising route to investigate the FFLO pairing mechanism starting from experiments with small-size cold atoms systems acting as superconducting grains. The present study can also be generalized to higher dimensions, assuming that the transverse confinement is reduced so that some of the excited single-particle states in the transverse directions become populated. For instance, one can study how the intensity of the noise correlation function depends on the numbers of particles in the system, for a fixed system size; we expect that, contrary to the 1D case discussed here, the intensity increases as the atom density increases, because interaction effects become stronger. It would also be interesting to understand whether the signatures of FFLO pairing in the noise correlation function remain visible in higher dimensions, where the exotic superfluid is known to be less robust.
Another promising direction, which directly connects with π-phases [59] and hybrid Josephson junctions [6], is to consider that both left and right side of the barrier represent bulk 1D systems with uniform densities (far from the barrier region). This corresponds to taking the thermodynamic limit N σ → +∞ and L → +∞, with N σ /L and the barrier width D being finite. It is worth stressing that such a limit is out of reach for the method used in this work, because the computational effort required to obtain convergent results grows too fast with system size and particle numbers. Other numerical approaches are more suitable to extract the ground state properties for bulk systems, including Quantum Monte Carlo or the Density Matrix Renormalization Group (DMRG). We leave this program for future studies.
All numerical data presented in this paper is freely available online [60]. NA NB ∆N ∆pF /π q0/π TABLE I. The most probable net pair momentum q0 for systems with different particle numbers NA, NB and imbalances ∆N = NB − NA. The system parameters are g = −5 and V = 0. The pair momentum q0 is found as the location of the maximum of function Q(q), defined in (5). It is compared to the theoretically predicted Fermi momentum difference ∆pF = ∆N π/2. It is seen that q0 ≈ ∆pF in all cases, as predicted. Not shown are the entries for NA = NB, for which in all tested cases q0 = ∆pF = 0.