Pairing patterns in one-dimensional spin- and mass-imbalanced Fermi gases

We study spin- and mass-imbalanced mixtures of spin-$\tfrac{1}{2}$ fermions interacting via an attractive contact potential in one spatial dimension. Specifically, we address the influence of unequal particle masses on the pair formation by means of the complex Langevin method. By computing the pair-correlation function and the associated pair-momentum distribution we find that inhomogeneous pairing is present for all studied spin polarizations and mass imbalances. To further characterize the pairing behavior, we analyze the density-density correlations in momentum space, the so-called shot noise, which is experimentally accessible through time-of-flight imaging. At finite spin polarization, the latter is known to show distinct maxima at momentum configurations associated with the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) instability. Besides those maxima, we find that additional features emerge in the noise correlations when mass imbalance is increased, revealing the stability of FFLO-type correlations against mass imbalance and furnishing an experimentally accessible signature to probe this type of pairing.


Introduction
Ultracold quantum gases represent excellent test grounds to challenge our understanding of the physics of strongly interacting Fermi systems. The unprecedented experimental plasticity of these systems allows for the investigation of a rich variety of physical scenarios in a clean and controllable way [1][2][3][4]. An intriguing phenomenon occurs in attractively interacting spin- 1 2 Fermi gases, where singlet pair-formation of spin-up and spin-down particles leads to superfluid behavior at low temperatures. In the case of equally populated spin species, the mechanism is well understood and can be explained via conventional Bardeen-Cooper-Schrieffer (BCS) theory.
Our understanding of fermion pairing can be further tested by considering deviations from the balanced limit of equal Fermi surfaces. A prominent example for such a "deformation" is obtained by introducing a spin imbalance. In this case, superfluidity is generally expected to break down above a critical polarization, commonly referred to as the Chandrasekhar-Clogston limit [5,6]. Below the critical polarization, so-called polarized superfluid states may be stabilized, for which several different mechanisms have been proposed in the literature [7][8][9]. A particularly intriguing scenario was put forward independently by Fulde and Ferrell (FF) and by Larkin and Ovchinnikov (LO), who argued on very general grounds that a spatially inhomogeneous ground state may be energetically most favorable even in infinite systems [10,11]. In this case, the emerging pairs of spin-up and spin-down fermions carry a finite center-of-mass momentum, potentially giving rise to the formation of a spatially varying order parameter, generally referred to as the FFLO state (see, e.g., Refs. [12,13] for an illustration).
Despite the experimental ability to study population-imbalanced Fermi gases [14][15][16][17], a smoking-gun evidence for the spontaneous formation of FFLO-type pairing is still missing. In addition to the increased technical challenge, part of the reason for the absence of a clear signal of inhomogeneous pairing is the "lack of stability" of this inhomogeneous phase: For threedimensional (3D) systems, beyond mean-field calculations suggest that only a thin layer of the parameter space between the homogeneous superfluid and the normal fluid is occupied by an inhomogeneous phase, if at all [13,[18][19][20][21]. For one dimensional (1D) systems, on the other hand, inhomogeneous pairing is expected to exist in a wide range of parameter space [22][23][24][25] and experimental measurements are indeed consistent with FFLO-type pairing [16,17].
On the theory side, continuum 1D systems with zero-range interaction are described by the so-called Gaudin-Yang model which can be solved in closed form via the famous Bethe ansatz (BA) [4,26,27]. While many ground-state properties are accessible within this approach, correlation functions remain notoriously challenging to compute and typically one has to either resort to low-energy effective treatments, such as the Tomonaga-Luttinger-liquid (TLL) description, or use numerical methods, such as the density-matrix renormalization group (DMRG), exact diagonalization (ED) or Monte Carlo (MC) approaches.
In addition to the preparation of systems with two different hyperfine states of the same atom species, recent experiments explored the possibility to realize heteronuclear mixtures such as 6 Li-40 K, 40 K-161 Dy or 6 Li- 53 Cr, which feature a mass imbalance between the spin-up and spin-down fermions [28][29][30][31]. These setups represent a compelling test ground for our picture of the underlying pair formation. However, their theoretical description remains relatively scarce in the literature since the reduced symmetry of the model renders BA studies inapplicable, except for a very few specific values for the mass ratio. Nevertheless, progress on analytic insights has been made for special symmetric configurations in the few-body sector [32][33][34] and numerically the regime was explored by worldline MC [35,36], ED [37] and complex Langevin (CL) studies [38,39]. Beyond the few-body regime, several studies based on effective descriptions have been conducted [40][41][42][43][44] and numerical studies investigated the phase-diagram of the asymmetric 1D Hubbard model (which coincides with the continuum model at low fillings) by means of DMRG [45,46] as well as spin-and mass-imbalanced systems via MC [47], time-evolving block decimation (TEBD) method [48] and DMRG [49] approaches.
The present study aims to complement the above investigations by exploring two-body correlation functions for general spin and mass imbalances for attractive interactions. Our ultimate goal is to arrive at an experimentally accessible quantity to study pairing in spinand mass-imbalanced systems, and to that end we present numerical results for the so-called shot-noise correlation function as obtained from the CL approach [50][51][52]. The latter was recently found to provide a promising way to circumvent the sign problem in nonrelativistic Fermi gases [38,53,54]. It is an additional objective of the present work to further develop the method in this context by presenting the first study of spin-and mass-imbalanced systems in the strict T = 0 limit via a projective formulation, as well as the first CL determination of two-body correlation functions for nonrelativistic systems. The CL method is a versatile tool which, unlike other approaches mentioned above, may be readily extended to higher spatial dimensions. Indeed, this capability was recently demonstrated for 3D fermions at unitarity, where excellent agreement with state-of-the-art numerical and experimental data for the equation of state was reported [55].
The work is organized as follows: In Sec. 2, we introduce the model, define the relevant scales, and briefly discuss some aspects of the CL approach; a detailed discussion of our CL approach to ultracold Fermi gases can be found in Refs. [38,54,56,57]. Numerical results are subsequently presented, where we first validate our numerical study via a comparison of ground-state energies to known exact results obtained with the BA method. We then proceed with the main focus of the present work, namely the computation of two-body correlation functions of the purely spin-polarized case in Sec. 3, and finally of general spin-and massimbalanced systems in Sec. 4. Our conclusions and a brief outlook are presented in Sec. 5.

Model and Scales
We consider systems of two fermionic species interacting via contact interaction in a onedimensional periodic box of extent L. Ignoring interactions between fermions of the same species, the associated Hamilton operator in its second-quantized form readŝ whereψ † σ (ψ σ ) represent creation (annihilation) operators for a fermion of species σ and mass m σ . In the following, units are chosen such that = k B = 1. In the mass-balanced case, we set m ↑ = m ↓ = 1. For our conventions in the case of unequal masses, we refer to Sec. 4. Note that at low densities this system corresponds to the (asymmetric) 1D Hubbard model.
In our numerical studies below, we always work in the canonical ensemble which implies that the particle numbers N σ of the spin-up and spin-down fermions are parameters at our disposal, determining the densities n σ = N σ /L of the two species and their Fermi momenta k σ F = πn σ . Assuming that discretization artifacts in our CL studies on a finite lattice are negligible, physical observables depend on three parameters, namely the mass ratio κ = m ↑ /m ↓ , the spin polarization p, and the dimensionless coupling γ = g/n with n = n ↑ + n ↓ being the total density. The latter can be related to the s-wave scattering length a s via g = 2/a s (see, e.g., Ref. [58]). A meaningful comparison of systems associated with different mass ratios κ requires to identify one scale which is kept fixed to the same value in all systems. In this work, we shall keep the two-body binding energy ε B fixed for g < 0 (attractive case): In practice, this is conveniently achieved by introducing a suitably chosen rescaled coupling constantγ:γ Thus, for fixedγ, the dimensionless binding energy ε B /n 2 = −γ 2 /4 then remains constant when we vary the mass ratio κ. Note that we haveγ = γ in the mass-balanced limit κ = 1.

Outline of the numerical treatment
For our studies of ground-state properties of spin-and mass-imbalanced Fermi gases, we start from the canonical partition function Z which is obtained by projecting an arbitrary trial state |ψ T (typically chosen to be a Slater determinant) onto the ground state in the limit of large imaginary times β: This partition function (not to be confused with its finite temperature counterpart, whose calculation requires summing over a complete set of states) can be rewritten in terms of a Euclidean-time path integral over a Hubbard-Stratonovich (HS) field φ, which is introduced to integrate out the fermions exactly, yielding the usual expression in terms of fermionic determinants This expression is the basis of the well-known class of auxiliary-field methods [59]. The path integral in Z β can now be tackled straightforwardly with conventional stochastic methods, provided that the system is in a mass-and spin-balanced state. In that case, the determinants of the two fermion species are equal and their product yields a nonnegative number. Away from this balanced configuration, a standard MC evaluation of the path integral (6) is complicated by the infamous sign problem [54]. In order to tackle the systems at the heart of this work, we therefore employ the so-called CL approach which allows us to surmount the sign problem by a complexification of the HS field [54,60]. The key ingredient is the use of a complexified version of the Langevin equation to produce a Markov chain of randomly distributed field values φ which then allows us to measure observables in the same fashion as in conventional MC approaches based on, e.g., the Metropolis algorithm. Nevertheless, it should be stressed that the CL approach is still a method under construction and this prescription is not without shortcomings in certain theories, as discussed in Refs. [61][62][63][64]. As remarked above, however, these potential issues seem to be under control for ultracold Fermi gases in a wide parameter range [54], in particular for attractive interactions [36,38]. In practice, the field integration of Eq. (6) is performed by discretizing spacetime (on which the field φ lives) into a spacetime lattice of size N x × N τ , where N x and N τ determine the number of grid points in spatial and imaginary time direction, respectively. Thus, the field integral in Eq. (6) becomes an integral over a finite (albeit large) number of dimensions. In our computations of the pair-correlation functions, and their corresponding pair-momentum distribution functions and ground-state energies, we have fixed the density and polarizations in lattice sizes in the range N x = 120 − 140. Such lattice sizes allow for good resolution of the correlation functions and result in negligible finite-size effects. A discussion of finite-size effects of similar quantities can be found in Ref. [65]. Moreover, we fixed the spatial lattice constant to = 1, such that the first Brillouin zone extends from −π to π. For the imaginarytime axis it was found sufficient to use a discretization step of ∆τ = 0.05. In order to reach the ground state, statistically independent runs were performed for different projection times corresponding to N τ = 50 − 150 (depending on the density). Averaging runs across multiple converged projection times then yield robust estimates of ground-state properties, see, e.g., Ref. [66].
Our sampling strategy consists in integrating the CL equation for the HS field with a suitably chosen step size up to trajectory lengths t L ∼ 1000. Here, t L denotes the fictitious Langevin time that parametrizes the Markov chain and relates to the number of samples N via t L = N ∆t L . The parameter ∆t L denotes the adaptive integration step of the discrete CL equations and has been fixed to an average value of ∆t (0) L = 0.04 which was found to be sufficient for the study of the correlation functions presented below. Moreover, we have used a regulator term of strength ξ = 0.1 in order to stabilize the CL trajectories, as also done in previous studies with this method [38,39,56]. We have checked that a further reduction of this parameter leaves our present results for the correlation functions unaltered, such that an extrapolation to the limit ξ → 0 is not needed. During the evolution of the random process we take "snapshots" of the HS field separated by ∼ 1.0 "Langevin seconds" which is sufficient to yield decorrelated samples. Unless specified otherwise, data points reflect averages over five such CL trajectories associated with randomly chosen initial conditions of the CL equation. This yields roughly 5000 decorrelated samples in total, allowing for a statistical uncertainty of ∼ 1 − 2%. At this point, we would like to emphasize that the error bars presented in this work reflect purely statistical uncertainties and do not contain systematic artifacts associated with the discretizations underlying our numerical solution of the CL equation.
Finally, we note that we have monitored running averages and histograms of computed observables to analyze the reliability of our calculations. For all systems discussed in the present work, we have found well-behaved distributions with well-defined second moments and most importantly without so-called fat tails that could spoil the expectation values, see also Ref. [65] for a related discussion. Further checks, specific to the CL approach, include the investigation of imaginary parts of computed quantities which we have found to be exponentially decaying with increasing t L as expected for physical observables.

Pairing in population-imbalanced Fermi gases
The main focus of this section is the investigation of pairing with the aid of appropriate twobody correlation functions. To this end, we consider systems with a large number of particles on large lattices, i.e., systems sufficiently close to the thermodynamic limit (TL) where finitesize effects should be negligible. Moreover, to be consistent with our previous work [65], we only consider odd numbers of particles for a given species to ensure that we do not explicitly break translational invariance.

Ground-state energy
As a validation of our numerical framework, it is instructive to first look at ground-state energies and compare them with exact results from the literature. In Fig. 1, we show the ground-state energy for systems with fixed polarization p = 0.5 at various densities as a function of the modulus of the coupling strength γ with γ < 0 (attractive case). A comparison with the exact solution [67] reveals significant deviations from the thermodynamic limit (TL) for systems with n = 0.5 and |γ| = 3.0. This is a consequence of the discretization of the spatial volume in our study. In fact, the large density induces a finite interaction range on the lattice which causes a deviation from the pure zero-range limit in the continuum. However, we find convergence to exact results with decreasing density. For dilute systems with n = 0.16, we find excellent agreement across the considered range of interaction strengths, indicating the validity of our numerical approach.

On-site pair-correlation function
A prominent difference between 1D and 3D systems is the effect of long-range fluctuations which prohibits the spontaneous breakdown of a continuous symmetry in the 1D case [68,69]. Instead, quasi long-range order (QLRO) takes place [70][71][72], manifesting itself in a algebraic decay of two-body correlation functions: Here, r = |x − x | and ∆ C denotes the correlation coefficient associated with the underlying ordering mechanism. The dominant instability is the one with the smallest correlation coefficient, i.e., the one whose correlations survive the longest distances. For spin-imbalanced Fermi gases away from half-filling, it has been found that spin-singlet pairing is dominant in the parameter range of interest in this work [23]. For this reason, we compute the on-site pair-correlation function, defined as the overlap of a state with a point-like pair of spin-up and spin-down fermions removed at point x from the ground state and a state with such a pair removed at point x from the ground state:  In addition to the algebraic decay of this correlation function, spin-imbalanced systems are expected to exhibit an additional feature at sufficiently large distances r = |x − x |, namely a spatially oscillating behavior of the form [72]: Here, q = |k ↑ F − k ↓ F | quantifies the population difference and ∆ ↑↓ is the correlation coefficient associated with the pair-correlation function. Since the long-range behavior of ρ ↑↓ is related to the order parameter for the fermion gap, a phenomenological interpretation of the minima of ρ ↑↓ can be given [72]: These minima may be viewed as "pockets" in which the excess fermions (here: excess spin-up fermions for p > 0) tend to reside since the positions of these minima are associated with vanishing (or small) fermion gap.
In Fig. 2, we show our results for the on-site pair-correlation function for a total density n = 0.4 and for different values of the polarization p. In the unpolarized limit, the monotonic decay of the pair-correlation function is apparent, whereas oscillatory behavior is observed for any finite polarization. Moreover, the frequency of the oscillations increases when the polarization increases, as suggested by Eq. (9). We observe the above behavior for all polarizations considered in this work, indicating the absence of a Chandrasekhar-Clogston limit in 1D, in agreement with earlier MC and DMRG studies of the spin-imbalanced 1D Hubbard model [23,73].
To further analyze pairing, it is instructive to study the pair-momentum distribution which is the Fourier transform of Eq. (8): Here, we have introduced the creation and annihilation operators in momentum space: with the the box-normalized single-particle orbitals The pair-momentum distribution quantifies the probability of finding a zero-size (i.e. pointlike) pair composed of one spin-up and one spin-down fermion with a given momentum in the system. Thus, the FFLO state reveals itself through off-center peaks at k = ±q. Of course, from our analysis of the on-site pair-correlation function, the appearance of these off-center peaks does not come unexpected. Indeed, if pairing predominantly occurs at non-zero momenta, the corresponding real-space signal will be modulated by the appropriate wavenumber. In Fig. 3 we show results for the diagonal part of the pair-momentum distribution, i.e., n ↑↓ (k, k), at fixed density n = 0.4 for various spin imbalances. For the spin-balanced case, we observe a pronounced peak in n ↑↓ at k = 0, as it is expected from BCS theory. This is also in agreement with previous MC studies of such systems [65]. The situation changes as soon as we turn on a finite spin imbalance. We then observe peaks at k = 0 as a consequence of pairing across mismatched Fermi points. To illustrate this, we show the pair-momentum distribution for systems with polarizations p = 0.25, 0.5, 0.75 in Fig. 3. The position of the off-center peaks in momentum space follow the prediction for FFLO type pairing. This is highlighted in the inset of the right panel of Fig. 3, where the dependence of the peak positions as a function of the mismatch q = |k ↑ F − k ↓ F | is depicted.

Noise correlations
In addition to the pair correlator, we also discuss the density-density correlation function in momentum space -often referred to as the "atomic shot-noise". It is defined as follows G σσ (k, k ) ≡ ∆n k,σ ∆n k ,σ = n k,σnk ,σ − n k,σ n k ,σ .
Here,n k,σ denotes the number operator in momentum space for the species σ and ∆n k,σ = n k,σ − n k,σ denotes the fluctuation around the thermal average. This quantity has been proposed as a suitable probe for correlations in ultracold quantum gases [74]. Several theoretical studies predict distinct imprints of the underlying ordering mechanism in the shot noise [23,[75][76][77]. Experimentally, it is accessible through the analysis of time-of-flight images of the spatial density profile after ballistic expansion, i.e., after interactions and the trap have been switched off. The so-obtained time-of-flight images are then proportional to the single-particle momentum distribution. A requirement for the experimental measurement is the ability of spin-selective imaging of the sample density which has been applied in several setups so far (see, e.g., Refs. [78][79][80]).
Loosely speaking, the computation of G ↑↓ (k, k ) allows us to study the internal momentum structure of a pair of one spin-up and one spin-down fermion. Indeed, we can combine the momenta k and k such that we obtain k cm = (k + k )/2 and k rel = k − k being the centerof-mass and relative momentum of the two fermions in a pair, respectively. To understand the quantity itself, it is instructive to think of it as the covariance matrix of the singleparticle momentum distributions n σ (k) and n σ (k ): Positive values of G σσ (k, k ) correspond to situations where high (low) values of n σ (k) occur alongside with high (low) values of n σ (k ), whereas negative G σσ (k, k ) corresponds to anti-correlation, i.e. high (low) values of n σ (k)  entail low (high) values of n σ (k ). The two cases can be physically identified with particleparticle (or hole-hole) and particle-hole correlations, respectively. Within this picture, it is also straightforward to see that G σσ (k, k ) vanishes for non-interacting systems, as the two single-particle distributions are then statistically independent. Note that the converse of this statement does not hold. In Fig. 4, we show the density-density correlation function G ↑↓ (k, k ) in the first Brillouin zone at constant total particle number N = 26, for a variety of polarizations at the intermediate coupling γ = −2.0. For all results discussed in the following, we have fixed the spatial lattice size to N x = 64 which we found sufficient to study relevant features (see, e.g., Ref. [23]).
The top-left panel of Fig. 4 shows the unpolarized case, for which the dominant correlations are expected to occur at the momentum space coordinates (±k ↑ F , ∓k ↓ F ), reflecting BCS-type pairing, i.e., pairing at equal and opposite momenta between spin-up and -down particles. Since the Fermi momenta of the two species coincide in this case, the correlation peaks occur on the anti-diagonal k = −k and thus indicate pairing with a center-of-mass momentum of q = k + k = 0.
For non-zero polarizations, a similar -yet distinct -picture emerges: While pairing still appears to happen close the respective Fermi points, the peaks of positive correlation shift "outwards" from the central anti-diagonal line by an offset of |q|. This is a hallmark of FFLOtype correlations which we already discussed in the previous section. Since the positive particleparticle correlations dominate over the negative particle-hole contributions, singlet-pairing clearly dominates other ordering mechanisms, in particular the formation of a charge-density wave.
Another interesting feature is the "checkerboard" pattern, that separates the Brillouin zone into segments of positive and negative correlations at the intersection lines of the Fermi points. Employing the argument above, we expect four particle-hole-like areas with negative correlations alongside with five particle-particle (or hole-hole) regions with positive correlation. This is reflected in our numerical data which is also in excellent agreement with the observations in Ref. [23].
With increasing interaction strength, it becomes energetically more and more favorable to scatter particles from well below the Fermi point which leads to a broadening of the peaks of positive correlation at the opposite Fermi points. In unpolarized systems, this broadening eventually leads to the vanishing of the checkerboard pattern whereas the pattern is extremely stable in the spin-imbalanced case, indicating the "preservation" of the respective Fermi points as expected for systems with gapless excitations [75].

Spin-polarized fermions with unequal masses
Turning now to the spin-and mass-imbalanced case, we note that the reduced symmetry of the Hamiltonian renders the system non-integrable, such that numerical methods have to be employed. Several studies based on DMRG [49], MC [47] and TEBD [48] have explored the phase diagram for either purely mass-imbalanced or spin-imbalanced 1D mixtures before. In this section, we present results for mass ratios which are in the range of potentially realizable values in experiments. This should allow for a check of previously made predictions on the stability of pairing in the presence of a finite mass imbalance in 1D Fermi mixtures.
In the following we set the scales such that the mass m l of the lighter species is fixed to m l = 1. This leaves the mass m h of the heavier fermion species at our disposal, which we characterize by the numerical factor κ as m h = κm l (see also Sec. 2.1 for our conventions). For the balanced case, we thus obtain a reduced mass µ = 1 2 and henceγ = γ, whereas the reduced mass increases for mass ratios κ = 1 from µ = 1 2 up to a maximum of µ = 1 in the limit of κ → ∞. Keepingγ fixed ensures that we take out any trivial dependence on the mass imbalance. Any residual dependence of our results on the mass imbalance should be a true many-body effect.
As mentioned above, we now have to deal with a reduced symmetry of the Hamiltonian. To be specific, because of the additional mass asymmetry, the system is no longer invariant under a flip of the sign of the polarization, p → −p. This requires us to consider two separate scenarios: The case where the majority species is heavier than the minority (heavy-majority case) and the opposite case (heavy-minority case). We keep the convention as in the spin-polarized case, namely N ↑ > N ↓ . This implies that κ > 1 (κ < 1) corresponds to a heavy-majority  (heavy-minority) system.

Pair-momentum distribution
As a first step, we investigate again the diagonal of the pair-momentum distribution n ↑↓ (k, k) as we did in the previous section. In Fig. 5, we show n ↑↓ (k, k) for (N ↑ , N ↓ ) = (17,9) particles on a lattice of N x = 64 sites, corresponding to a density of n ≈ 0.41. Interestingly, we observe almost no variation of n ↑↓ (k, k) as a function of κ for both the heavy-majority case (left panel) and the heavy-minority case (right panel). Most importantly, for all mass imbalances considered in this work, we observe constant positions of the maxima, indicating the stability of the FFLO-type correlations in the presence of a finite mass imbalance. Of course, the robustness of the pair-momentum distribution with respect to a variation of the mass imbalance translates directly into a robustness of the on-site pair-correlation function, which is the Fourier transform of the pair-momentum distribution. For the considered mass ratios, which are in the range of potentially realizable values, we observe FFLO-type correlations, regardless of the spin polarization. At the interaction strength considered here, γ = −2.0, this is in line with previous studies [47][48][49].
For mass imbalances beyond those considered in the present work, it has been observed in the asymmetric 1D Hubbard model that FFLO-type ordering is eventually suppressed. Instead, a charge-density-wave-type phase appears before either phase separation sets in [45] or the system collapses [41,47]. The "critical imbalances" above which these transitions occur depend on the interaction strengths but are expected to be well above the experimentally realizable values for the mass ratio κ. In any case, an exploration of this highly mass-imbalanced regime is beyond the scope of the present study.
For completeness, we comment on the possible occurrence of multi-particle bound states at large mass imbalances. For commensurate spin imbalances (i.e. integer ratios of the spinup and spin-down densities), stable phases -where ordering is associated with the formation of such multi-particle bound states -have been observed [43,44,49,81]. For general spin imbalances away from the commensurate points, however, such phases are unstable and are thus not relevant for the present discussion.

Shot-noise correlations
It may appear counterintuitive that a finite mass imbalance seems to leave the pairing of spin-up and spin-down particles unchanged altogether, as suggested by our results for the pair-momentum distribution. Loosely speaking, by increasing the mass of a fermion species, the spacing between the energy levels decreases (or, in terms of the Hubbard model, the bands flatten with decreased hopping while retaining the same number of states). With such a rescaled energy spectrum, the participation of higher-lying momentum states in scattering processes appears energetically possible and therefore pairing should also take place (far) away from the Fermi points. This is indeed the case. As we shall see below, however, pairs still predominantly form with the FFLO-momentum q = |k ↑ F − k ↓ F | even in the mass-imbalanced case, such that only certain combinations of momenta of the spin-up and spin-down particles away from the Fermi points are also preferred. In essence, the energy in mass-imbalanced systems is still minimized by accommodating the excess particles at the almost gapless points (i.e. at the nodes of the on-site pair-correlation function) which implies pairing at finite k cm . In the following we analyze this aspect with the help of the shot-noise correlator G ↑↓ which shows distinct fingerprints of this mechanism.
Recently, studies based on ED explored shot-noise correlations [82] and a related measure [83] to study the influence of mass imbalance on pair formation in harmonically trapped spin-balanced few-body systems have also been performed. In the present work, however, we address this quantity in the presence of spin imbalance for many-body systems and show that unambiguous features emerge which can potentially be identified in currently available experimental setups.

Heavy-majority systems
First, let us consider the heavy-majority scenario in more detail, for which we show our results in Fig. 6. As discussed above, in the balanced case (top-left panel), pairing predominantly occurs at lines of constant center-of-mass momentum k cm = ±q ≡ ±|k ↑ F − k ↓ F |, which corresponds to k = k ± q and describes the positions of the two maxima in the pair-momentum distribution discussed above. We observe peaks at (±k ↑ F , ∓k ↓ F ) which are associated with the expected FFLO-type pairing in the vicinity of both Fermi points. The dominant peak labeled with "A" in the top-left panel of Fig. 6 is located on the line k cm = −q. 1 With increasing mass ratio κ, additional peaks emerge at (±k ↑ F ∓ 2q, ∓k ↓ F ), indicated by stronger particle-particle correlations (i.e., darker red coloring). For the sake of readability, only the point at (−k ↑ F + 2q, k ↓ F ) is labeled with "B" in the rightmost panel of Fig. 6. The emergence of a peak at the coordinates (k ↑ F −2q, −k ↓ F ) reflects pairing of light fermions "sitting" close to the opposite Fermi point with heavy fermions whose momenta are shifted such that the center-of-mass momentum remains fixed.
The build-up of a second peak is highlighted in detail in the bottom-left panel of with a different internal momentum structure are not formed in the system. The probability to find the latter pairs is just lower.
The bottom-right panel of Fig. 6 illustrates the underlying pairing mechanism: For sufficiently large mass imbalance, it becomes energetically more favorable to pair up with heavy particles from far below the Fermi point while still obeying the constraint |k cm | = |q|. The sketch depicts the perspective of the lighter (spin-down) particle at the Fermi point: In addition to forming a pair with a spin-up particle at k = −k ↑ F (big black dots), there is now a significant probability to pair up with a particle at k ↑ = −k ↑ F + 2q (big blue dots).

Heavy-minority systems
The investigation of the heavy-minority scenario, which is shown in Fig. 7, reveals a similar picture albeit with reversed roles of the majority and minority particles. For the same mass ratios as above, we now observe additional peaks at (±k ↑ F , ±k ↓ F ∓ 2k ↑ F ), already clearly visible as additional "fringes" in the middle and the right panel of Fig. 7. The build-up of a secondary peak becomes again more clearly visible by comparing the shot-noise correlations along the line k cm = −q, see bottom-left panel of Fig. 7. The additional peak "C" emerges at the expense of some weight at the opposite Fermi point of the majority species.
The underlying pairing mechanism is depicted in the bottom-right panel of Fig. 7 where pairing of a light majority particle sitting at its Fermi point is now favored to occur with either a heavy-minority particle sitting at its Fermi point (big black dot) or a fermion far above the Fermi point (big red dot). Similar to the heavy-majority case, the latter option becomes energetically more favorable with increasing mass ratio since higher-lying states can, loosely speaking, be reached more easily. Recall that the difference between the energy levels is proportional to the inverse of the mass of a particle. The corresponding peak in the shot-noise correlation function is labeled with "C" in the top-right panel of Fig. 7.

Discussion & Summary
To summarize, we have investigated the ground-state pairing behavior of non-relativistic spinand mass-imbalanced Fermi gases in 1D, which is expected to be of the FFLO type. To that end, we employed the complex Langevin method (applied to the present system for the first time) and calculated two-body correlation functions, finding remarkable agreement with earlier theoretical predictions and numerical calculations, wherever available. In particular, we have observed hallmark features of the FFLO pairing across all spin polarizations considered in agreement with earlier studies on related systems. Most importantly, we have presented results for the density-density correlation function. The latter shows an exceptionally clean signal for pairing across the mismatched Fermi points and is also accessible in experimental setups. Especially for spin-and mass-imbalanced systems, our study has revealed remarkable features with increasing mass imbalance. In fact, we have found that, even for very large mass imbalances, the dominant correlations are of the FFLO type. However, in addition to the naively expected conventional FFLO-type pairingpeaks appearing at the opposite Fermi points in the density-density correlation function, we have observed that new peaks emerge with increasing mass imbalance. Interestingly, the spectral weight is found to be only shifted along anti-diagonals of constant center-of-mass momentum of the fermion pairs. This is consistent with our finding of a constant peak position in the pair-momentum distribution. In particular, we have found that the aforementioned new peaks build up at the Fermi points of the light fermion species, pushing some of the weight of the heavier component away from its own Fermi point. Phenomenologically, the effect may be understood by recognizing that it is "easier" to move the higher-mass component to a different momentum state since the energy expense for this process decreases with increasing fermion mass.
While the present work focuses on bulk systems, i.e., we have employed periodic boundary conditions, we nevertheless expect it to be of relevance for experiments where the particles reside in non-uniform trapping potentials. In fact, earlier studies of the shot-noise for trapped systems suggest the stability of the distinct pairing patterns in the presence of harmonic (and likely any other) confinement [23], and even in the few-body sector [82,84].
The success of the CL method in computing two-body correlation functions sets the stage for further investigations based on this method. In the future, it will be interesting to study the temperature dependence of the pairing patterns found here. In addition, future efforts will be directed towards systems in higher spatial dimensions, where long-range order is more stable against fluctuations. A particularly exciting avenue concerns the advent of quantum gas microscopes, which give direct access to shot-noise correlation functions for two-dimensional systems in optical lattices and could provide a viable way to study, e.g., the fate of pair