Unconventional pairing in few-fermion systems at finite temperature

Attractively interacting two-component mixtures of fermionic particles confined in a one-dimensional harmonic trap are investigated. Properties of balanced and imbalanced systems are systematically explored with the exact diagonalization approach, focusing on the finite-temperature effects. Using single- and two-particle density distributions, specific non-classical pairing correlations are analyzed in terms of the noise correlations—quantity directly accessible in state-of-the-art experiments with ultra-cold atoms. It is shown that along with increasing temperature, any imbalanced system hosting Fulde–Ferrel–Larkin–Ovchinnikov pairs crossovers to a standard Bardeen-Cooper-Schrieffer one characterized by zero net momentum of resulting pairs. By performing calculations for systems with different imbalances, the approximate boundary between the two phases on a phase diagram is determined.

The appearance of non-classical correlations collectively induced by attractive interactions between fermions is one of the most spectacular macroscopic manifestations of quantum mechanics. As justified by Bardeen, Cooper, and Schrieffer in 1 (BCS), it is directly responsible for the remarkable phenomenon of superconductivity-a flow of an electric current without any resistance-discovered by Onnes over a hundred years ago 2 . In the case of homogeneous bulk systems, an existence of pairing induced by attractions is explained on fundamental grounds by showing that for any non-vanishing attractions the system is forced in a nonperturbative manner to rebrand its ground state and minimize energy by forming pairs (called Cooper pairs) composed by oppositespin fermions having exactly opposite momenta 3,4 . Later it was argued that this argumentation is rigorously valid only for perfectly balanced systems, i.e., when both components have exactly the same Fermi spheres. In a more general scenario, when the system is no longer symmetric with respect to the exchange of components, pairing correlations are still favored but resulting pairs may have non-vanishing center-of-mass momentum, as explained independently by Fulde and Ferrel 5 and Larkin and Ovchinnikov 6 (FFLO). This observation triggered many investigations aimed to understand different properties of such systems (see for example [7][8][9][10][11][12][13][14][15][16][17][18] and citations within) and in consequence to unquestionably confirm that such pairs do exist. Unfortunately, indisputable experimental evidence is still lacking. In three-dimensional cases, one of the recognized obstacles is the fact that the final FFLO signal is much less noticeable since any non-zero value of finite net momentum has no unique direction in space. In consequence, the condensation of pairs is spread among different two-particle orbitals and the FFLO phase is stable in a quite small region of the phase diagram 19 . It is argued that systems with reduced dimensionality (particularly one-dimensional systems) may serve as a promising bypassing platforms for FFLO detection [20][21][22][23][24][25][26][27] . For recent reviews see 28,29 . The universal nature of such correlations in different number of spatial dimensions can be seen through the variety of physical systems that are studied. From nuclear matter and neutron stars [30][31][32] , through organic superconductors [33][34][35] and heavy-fermion systems in solid state physics 36-38 , up to ultracold gases 28 .
Significant experimental progress on ultracold atomic systems has opened many interesting and non-obvious paths for exploration. One of them is related to the ability of very precise control of the number of particles in a system. Pioneering experiments of this type [39][40][41] have opened an alternative route to study quantum many-body systems from the side of mesoscopic phenomena 42,43 . Amidst a flurry of different ideas related to such small quantum systems, recently it has been theoretically spotted that small mixtures of attractively interacting fermions may host correlations very similar to pairing known for bulk systems [44][45][46][47][48][49][50] . Depending on inter-component imbalance, they may serve as suitable platforms for observing unconventional pairing mechanisms 51 . Very recently, these theoretical considerations have taken on new importance since innovative methods for detecting and measuring inter-particle correlations in such systems have been just developed [52][53][54] . The majority of previous works have mainly focused on different ground-state properties or some specific dynamical scenarios of such systems with repulsive interactions. In this work, we aim to go beyond these schemes and to extend recent discussion on pairing (the conventional BCS as well as the unconventional FFLO) to the situation when the system is prepared in a mixed state of a given temperature. In this way, we want to establish another theoretical link between few-body problems and finite-temperature results obtained recently for a large number of particles. Particularly, we want to make a correspondence with previously obtained Monte Carlo predictions and solutions based on the Bogoliubov-de Gennes method for harmonically trapped mixtures of fermions [55][56][57] . It is worth noticing that on the repulsive branch of interactions some theoretical analyses of the impact of finite temperatures were presented already [58][59][60][61] . Although finite-temperature effects for attractively interacting few-fermion systems were not deeply studied, they can be important from the experimental perspective, since BCS-like correlation has been already measured 62 .

The system studied
We consider an effectively one-dimensional mixture of two fermionic components σ ∈ {↑, ↓} consisting of N ↑ and N ↓ particles of the same mass m, and trapped in an external harmonic oscillator potential of frequency ω . We assume that on a time scale of considered experiments there are no relevant physical channels causing spin flips, thus the numbers N ↑ and N ↓ are independently conserved quantities. In the case of experiments with ultracold alkaline atoms, this assumption is well-justified. In the second-quantization formalism the Hamiltonian of the system reads: where the fermionic field operator � σ (x) annihilates a particle from a component σ at position x. It obeys standard anti-commutation relations The dimensionless parameter g is related to the one-dimensional interaction strength. In the effective one-dimensional model studied, it can be obtained from the three-dimensional counterpart, by assuming that transverse motion is completely frozen (for example by very strong confinement) and thus it can be integrated out 63 . It means that the value of g depends on the transverse width of the system. This gives one of possible ways to tune its value. Another possibility is granted thank to internal electronic structures of particles, which, although completely irrelevant for spatial dynamics, give a tool for controlling inter-particle scattering via the Feshbach resonances 40,64 . Since the Hamiltonian does commute with number operators N σ = dx� † (x)�(x) we analyze its properties in subspaces of fixed N ↑ and N ↓ . To get better correspondence with standard description in the large system limit, we introduce an intensive quantity characterizing imbalance between components, i.e., the global polarization Note, that the Hamiltonian (1) is fully symmetric under exchange of components. It means that systems with opposite global polarisations ( P and −P ) have exactly the same properties provided that labels ↑ and ↓ are adequately exchanged. Therefore, without losing generality, we can limit ourselves to cases with P ≥ 0.

Many-body spectrum
To get a better understanding of the properties of the system at finite temperature, let us first discuss spectral properties of the many-body Hamiltonian (1) on the attractive branch of interactions. In contrast to the repulsive one 58,65 , this has not been extensively studied in the literature. It is quite natural that in the case studied the spectrum is not bounded from below when attractive g is increasing. In fact, the model (1) is not appropriate for sufficiently large attractions since it completely neglects the possible formation of molecules. However, for finite and reasonable g, the model appropriately captures the system's properties 45,46 and is experimentally relevant 40 .
To present fundamental features of the many-body spectrum of the system studied, we take as an example the system containing N = 8 particles with different imbalances N and we perform numerically exact diagonalization of its many-body Hamiltonian (1) to obtain its eigenstates |i� and corresponding eigenenergies E i (see Appendix for numerical details). Then we plot the lowest energy gaps E i − E 0 as functions of interaction strength g (see Fig. 1a). It is clear that for any attraction the many-body ground state is rather isolated. Indeed, for any number of particles and attractive interactions, the gap to the first excited state is never smaller than the energy of a single excitation in the harmonic confinement. It means that any pairing correlations present in the ground-state of the system probably remain stable against small thermal excitations. It is also clear that along with increasing attractions the degeneracy of the first excited state (existing for vanishing interactions due to the fundamental independence of components) is lifted. Note, however, that only for perfectly balanced systems ( N = 0 ) one of the first excited states rapidly gains energy and very quickly become unimportant for low-energy physics. That suggests that the thermal resistivity of the system could be strongly dependent on the imbalance, thus may be of fundamental importance for experimental discrimination between standard BCS and FFLO pairing hosted respectively by balanced and imbalanced systems.
In the following, we assume that the system studied is prepared in a stationary state of a given absolute temperature T. We do not settle whether thermalization occurred as a result of interactions with an external thermostat or as a consequence of self-thermalization. The latter may be challenging due to the intrinsic properties www.nature.com/scientificreports/ of the system-its lowered dimensionality and mesoscopic character. Nevertheless, we anticipate that the state of the system can be well-described with a density matrix of the Boltzmann form: where P i = Z −1 e −E i /k B T are the Boltzmann distribution probabilities and Z = i e −E i /k B T is the corresponding partition function ( k B is the Boltzmann's constant). Of course, in the limit of vanishing temperature, T → 0 , the state of the system is represented by its many-body ground state |0� . However, even for finite temperatures, due to the mentioned finite gaps to the excited states, the system is quite robust and may remain in the ground state. This property is clearly visible when partitions P i are plotted as functions of temperature (Fig. 1b). Since it will be crucial for the discussion of experimental attainability of pairing phases for realistic systems let us relate this observation to the experiments in Selim Jochim's group. From Fig. 1b it is clear that the lowest excited states start to contribute at a temperature around T ≈ 0.5ℏω/k B . Taking experimentally typical frequency of the trap ω = 2π × 1.488 kHz and lithium 6 Li atom mass m = 1.15 × 10 −26 kg it corresponds to temperature 35.7 nK. Thus, we see that the temperature effects definitely have to be taken into account in the modeling when predicting the system properties.

Two-body correlation and pairing
There are many ambiguities in translating concepts suitable for many-body description into the few-body regime. Mostly, they originate from the fact that in a large number of particles limit one does not care much about their precise counting and therefore performs analysis in the grand canonical ensemble approach. This gives a rise to define mean-field order parameters (pairing fields) as expectation values of appropriately defined annihilation operators 66 . Unfortunately, for small systems we rather focus on properties of the system when the number of particles is well-defined, thus we work in the canonical ensemble framework. Therefore, in the case of mesoscopic systems studied here, instead of considering order parameters, which in fact are not directly measurable, we focus on quantities straightly capturable in experiments that may quantify non-classical pairing correlations in the system. When two-component systems are considered, then the natural quantity describing two-particle correlations collectively induced between particles is the reduced two-particle density matrix. However, from the experimental point of view, measuring its off-diagonal elements (encoding two-body coherence) is not easy if possible at all. Therefore, we focus only on its diagonal part Tr ρ Tn↑ (p 1 )n ↓ (p 2 ) , i.e., the two-particle momentum distribution in the thermal state ρ T . The density operator in the momentum domain n σ (p) is defined The energy gaps between the many-body ground state and excited states for systems with N = 8 particles and different particle imbalances N (polarization P = 0, 0.25, 0.5, 0.75 , respectively). Independently on the case, two of the lowest excited states are always degenerate at g = 0 and the degeneracy is lifted for finite attraction. Note, however, that the lift strongly depends on the polarization. (b) the Boltzmann probabilities P i as functions of rescaled temperature for two extreme imbalances ( P = 0 and P = 0.75 , respectively) and interaction strength g = −3 . Note that for sufficiently low but finite temperatures (due to the finite gap in the spectrum) the system is robust and rather remains in its many-body ground state. All energies and temperatures are expressed in natural units of the problem, i.e., ℏω and ℏω/k B , respectively. www.nature.com/scientificreports/ It is clear that the two-particle momentum distribution defined above encodes not only non-classical correlations forced by interactions but also accidental correlations originating in single-particle distributions 67 . Therefore, to filter out the latter from the description, we go along previous experience [67][68][69][70] and we define the so-called noise correlation encoding pure two-body correlations between particles where we use a short notation �•� T ≡ Tr ρ T • for averaging over the many-body mixed state ρ T of the system of given temperature T. Let us point here, that in the case of ground-state properties, the noise correlation concept has been extremely useful to study two-body correlations in general 71 , but the pairing in particular 27,48,50,72,73 . It turned out that, at least for the ground-state properties, this measurable quantity can be quite easily adopted as a direct indicator of pairing correlations signaling not only standard Cooper pairs but also unconventional pairs with non-vanishing center-of-mass momentum (FFLO) 48 . In the following, we extend this analysis to thermal mixed states of the system.

Balanced systems
Let us start the analysis from the simplest case of balanced system (polarization P = 0 ) containing in total N particles. In Fig. 2 we show the dependence of the noise correlation (4) on the interaction strength and temperature for N = 8 particles. When the system is purely in its many-body ground-state ( T = 0 ), strong anti-correlations between opposite spin momenta are enhanced with increasing interactions. This signals a strong BCS pairing mechanism studied recently 45 . It is also clear that this enhancement is reduced by increasing temperature. One can quantify these competing behaviors by integrating the nose correlation close to its anti-diagonal part. For example, it can be done by calculating the pairing intensity Q T expressed as integration with appropriately localized filtering function Q T = dp 1 dp 2 G T (p 1 , p 2 )F (p 1 + p 2 ), We checked that the final results do not depend on the particular shape of the filter function F , nor on its width κ that should be of reasonable value compared to the resolution of the noise correlation and size of the system. In our calculations we set κ = 0.04 √ ℏmω . In Fig. 3 we display the intensity Q T as a function of temperature for different interaction strengths and as a function of attractions for different temperatures. These results confirm that increasing attractions in the system enhances pairing correlations between opposite spin fermions. At the same time, by increasing temperature, these correlations are suppressed. Note, however, that for some small but finite range of temperatures pairing correlations are almost insensitive to the temperature. This fact is a direct consequence of the finite gap to the excited many-body states mentioned previously.

Polarized systems
The situation is more interesting when imbalanced mixtures are considered. In these cases, although intercomponent attractions force the system to form correlated pairs, their center-of-mass momentum is not vanishing. In the limit of vanishing temperature ( T → 0 ) the situation was studied already with all details in 48 . It was argued that the effect is a direct manifestation of the FFLO mechanism originating in a mismatch of Fermi momenta of interacting components. It was shown that in the case of harmonically trapped particles centerof-mass momentum of the FFLO pair q 0 is related to the difference of quasi-Fermi momenta p F↑ − p F↓ , being defined through the Fermi energies ε Fσ as In the following, we aim to generalize these findings for finite temperatures. The analysis starts by displaying noise correlations for different imbalances and temperatures. In Fig. 4 we show examples for the system containing N = 8 and for quite strong attraction g = −3 . The first row (results for the balanced system) corresponds directly to the middle row in Fig. 2 (note the different range of the scale). The first column, on the other hand, corresponds to the zero-temperature results obtained earlier in 48 . Indeed, we see, that particle imbalance between components leads directly to a specific splitting of the anti-correlation of opposite-spin momenta and force correlated pairs to have finite center-of-mass-momentum. Along with increasing imbalance, we notice that areas of enhanced pairing keep moving away from the anti-diagonal p 1 + p 2 = 0 . When the temperature is increased they are not only substantially diminished but also some tiny reversal shift towards the ani-diagonal is visible. To quantify this quite complicated behavior we generalize the pairing intensity (5) to make it center-of-mass dependent The quantity Q T (q) can be understood as a distribution of the net momentum of correlated pairs. For vanishing temperature, the definition (8) is identical with the corresponding quantity introduced in 48 . Such a construction of Q T (q) gives an insight into the most probable center-of-mass momenta of correlated pairs since they are (8) Q T (q) = dp 1 dp 2 G T (p 1 , p 2 )F (p 1 + p 2 + q). www.nature.com/scientificreports/ encoded in its maxima. Note that due to the general symmetry of the problem studied, Q T (q) ≡ Q T (−q) . Thus in the following, without losing generality, we display results only for non-negative momenta, q ≥ 0. As a working example, in Fig. 5a we show the pairing intensity distribution Q T (q) for the system of N = 8 particles and the highest imbalance N = 6 , obtained for different temperatures T. For zero temperature a clear maximum at q 0 ≈ 2.9 √ mℏω corresponds directly to the difference of quasi-Fermi momenta predicted by (7) and signals strong FFLO pairing with appropriate center-of-mass momentum. When the temperature is increased, the amplitude of the maximum at q 0 is reduced and its position is shifted toward the second maximum at q = 0 representing standard BCS pairs which are also present in the system. Since, the maxima for zero and finite momentum coexist for any temperature, it is hard, if possible at all, to witness correlation of a given type in the system independently. As a consequence, for sufficiently large temperature, the distinction between FFLO and BCS pairs is no longer possible. Especially since both pair formation mechanisms are already very weak at this point and it is difficult to refer to any significant existence of correlated pairs. More importantly, as the temperature increases, there is an initial increase in the amplitude of BCS pairs and a simultaneous decrease for FFLO pairs. Finally, for sufficiently large temperatures both maxima melt to a single, very wide plateau meaning that any pairing correlations in the system are not detectable by the noise correlation means. Between these two extreme cases, i.e. strong FFLO signaled at T = 0 and lack of any pairing at high T limit, we observe that maximum at q 0 drops much faster than the one in the center at q = 0 . In consequence, for some range of temperatures, their amplitudes are comparable and the system undergoes crossover from unconventional FFLO-like pairing to the standard BCS-like mechanism. A very similar effect has been shown in terms of Bogoliubov-de Gennes formalism 57 , and with Monte Carlo calculations for one-dimensional systems of many fermions confined in optical lattice 56 . Our result can be viewed as a few-body precursor of this many-body behavior.
To quantify the mutual relation between the FFLO pairing and the BCS pairing in the system at a given temperature, we introduce the simplest possible measure related to the difference of corresponding maxima We checked that this simple definition gives qualitatively the same predictions as some more sophisticated ones introduced recently 50 . Since a concrete value of ξ(T) has no direct interpretation, in the following we always normalize it to its zero-temperature value, ξ(T)/ξ(0) . With this normalization, value 1 means that FFLO  www.nature.com/scientificreports/ correlations carried by the finite-momentum pairs are the same as at T = 0 . The temperature at which the quantity starts to drop quickly indicates the moment at which the FFLO phase looses its stability.
In Fig. 5b we show how this dependence looks for systems with N = 8 particles and different imbalances N . It is clear that for higher imbalances larger temperature is needed to destroy finite-momentum pair correlations. Dependence of ξ(T)/ξ(0) on temperature supports our previous findings that the transition from FFLO to BCS pairing induced by temperature has clear features of a crossover rather than a sharp transition as predicted for one-dimensional systems of many particles 56 . Alike for Q T in balanced systems, ξ(T)/ξ(0) does not change much for very small temperatures. We define the characteristic temperature T C of the crossover as a point when ξ(T)/ξ(0) drops below 90% (dashed line in Fig. 5b). Such definition captures adequately the moment at which the FFLO-type correlations start decaying and the system is quickly driven towards the BCS regime. We want to emphasize that the chosen threshold 90% is to some extent arbitrary. Therefore the exact value of T C will depend on our criterium. However, we checked that the particular definition does not affect the general connection between polarization and temperature which seems to be universal.
From our numerical results, it follows that the characteristic temperature T C slightly increases with attractions g. What is more, for a fixed interaction strength g we found that increasing the imbalance N makes the tail of ξ(T)/ξ(0) thicker. The reason is simply that the correlations in the system are much more complex than just one of the two types BCS or FFLO, and thus they are beyond noise correlation analysis. Therefore, for sufficiently large imbalances the function ξ(T)/ξ(0) does not drop exactly to zero anymore and after reaching its minimal value starts to grow. However, since the ground state of any imbalanced system is always FFLO-like correlated at T = 0 , our definition of characteristic temperature T C captures appropriately the physics of the crossover even if pure BCS-like correlation at large temperatures cannot be reached for a given interaction g.

Few-body phase diagram
Having all previous results and observations in hand we can repeat all the reasoning for different numbers of particles and different imbalances. In this way, we are able to form a provisional bridge between few-body systems and the many-body counterparts. We assume that in the limit of large systems appropriate description of the system is served by intensive properties which are independent of the size of the system. In the system studied the polarization (2) plays such a role. Therefore, we suspect that for large enough systems we should approach the many-body limit if properties of the system are presented as functions of polarization P and temperature T. It is worth noting that sometimes different many-body properties might saturate unexpectedly quickly (see for example experiments reported in 74 ).
To make the first view on this problem, in Fig. 6 we display characteristic temperatures of the transition T C for systems having different polarizations P. To make it systematic, we performed calculations for all possible distributions of particles among components with up to N = 12 particles in total. Different points correspond to systems with different numbers of particles shown in parenthesis. Note, that there are families of points giving the same polarization P. We find that points having the largest number of particles align on a quite regular and characteristic curve forming a border between FFLO and BCS pairing. We see that the characteristic temperature T C grows with the polarization up to P ≈ 0.5 and then starts to decrease. Due to a numerical limitation of our approach based on the exact diagonalization, we are not able to study systems with P > 0.84 (separated by a horizontal dashed line). At this point, it should be noted also that some points evidently stand out from the border predicted. These points correspond mainly to polaronic systems containing only one particle in a selected component. Of course, these systems, independently on the polarization P and the total number of particles N cannot The phase diagram obtained is surprisingly very similar to the phase diagram obtained earlier with Monte Carlo calculations for many fermions in a one-dimensional optical lattice 56 . Although some quantitative differences are visible, qualitatively the results obtained for these two substantially different systems are compatible.

Conclusions
In this work, using the exact diagonalization approach, we have examined how finite-temperature effects impact two-body correlations in two-component fermionic systems containing few particles with attractive interactions. We focused on correlations encoded in the so-called noise correlation-the quantity which in principle is directly accessible in experiments. We have shown that as the temperature increases, the correlations are initially insensitive (due to the energy gap) and from a certain moment they vanish very quickly.
Importantly, for imbalanced systems hosting the FFLO pairing in their many-body ground state, before pairing correlations are lost, we have observed the transition from the FFLO-like to the BCS-like phase. This transition has features of a crossover rather than a rapid phase transition and can be well characterized by the characteristic temperature T C interpreted as a temperature when the FFLO-like phase loses its stability.
By systematic studies of systems having different numbers of particles, we determine the approximate border in the phase diagram between the two phases. Our predictions qualitatively agree with previous results obtained for larger systems confined in periodic potentials and in harmonic traps. As a side result, we showed that polaronic systems containing a single particle in one of the components have typically substantially different properties and should be studied separately.

Data availability
All data generated or analysed during this study are available from the corresponding author on reasonable request.

Appendix: The numerical method used
Our numerical method is based on a numerically exact diagonalization of the many-body Hamiltonian (1) with applied energetic cut-off of the many-body basis [75][76][77][78] . First we consider a corresponding single-particle problem of a one-dimensional harmonic oscillator whose eigenstates are well-known and represented by wavefunctions where H k (x) is k-th Hermite polynomial and N k = 2 k k! √ πℏ/mω −1/2 is a normalization constant. Corresponding eigenergies are equal E k = ℏω(k + 1/2) . Having these solutions in hand, for a fixed number of particles N ↑ and N ↓ we consider all possible Fock states of the form  www.nature.com/scientificreports/ where α = (i 1 , . . . , i N ↑ ; j 1 , . . . , j N ↑ ) is a set of (sorted in descending order) indices of occupied single-particle orbitals (10) while fermionic operators â i ( b i ) annihilate ↑-particle ( ↓-particle) in corresponding states φ i (x) . A non-interacting energy of a Fock state |F α � is simply a sum of single-particle energies of occupied orbitals By decomposing field operators in the basis of these annihilation operators one can easily represent the many-body Hamiltonian (1) as a matrix in the basis of Fock states {|F α �} . Corresponding matrix elements can be calculated straightforwardly as where interaction integrals are expressed as These integrals can be calculated numerically on a dense spatial grid or analytically by using specific properties of Hermite polynomials 78,79 .
Naturally, in practice, we are not able to consider all possible Fock states since their number grows exponentially with the number of particles and number of single-particle orbitals. Thus some limitation of the Fock basis is required. Since we are interested in the low-lying spectrum of the Hamiltonian (1) and not too far from the non-interacting case ( g = 0 ) we select from a whole basis the states with the lowest non-interacting energy, i.e., only these Fock states whose energies ε α are not larger than some energetic cut-off E C . Since in the case studied single-particle energies are equally distributed and proportional to their occupation indices, this requirement technically means that we pick to the basis only these Fock states |F α � for which A detailed numerical algorithm for quick generation of the Fock basis in this form (also for other confinements) was presented recently 77 . In our calculations, we select the cut-off E C to provide well-converged final results. In practice E C /ℏω − (N ↑ + N ↓ )/2 is never less than 20.
After all these preparations the resulting matrix H αβ is numerically diagonalized using the Arnoldi alghoritm 80 . As a result one obtain the lowest eigenenergies E i and corresponding decomposition coefficients {f (i) α } representing many-body eigenstates |i� in the cropped Fock basis: With this decomposition, one can easily calculate the expectation value of any many-body operator expressed in terms of the field operators � σ (x) by decomposing them according to (13).