Strong Fulde-Ferrell Larkin-Ovchinnikov pairing fluctuations in polarized Fermi systems

We calculate the pair susceptibility of an attractive spin-polarized Fermi gas in the normal phase, as a function of the pair momentum. Close to unitarity, we find a strong enhancement of Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) pairing fluctuations over an extended region of the temperature-polarization phase diagram, which manifests itself as a pronounced peak in the pair-momentum distribution at a finite pair momentum. This peak should be amenable to experimental observation at achievable temperatures in a box-like trapping potential, as a fingerprint of FFLO pairing. Our calculations rest on a self-consistent t-matrix approach which, for the unitary balanced Fermi gas, has been validated against experimental data for several thermodynamic quantities.

Introduction − One of the main motivations that some time ago [1,2] initiated a vibrant experimental and theoretical research activity on ultracold polarized attractive Fermi gases [3], was the search for exotic superfluid phases, like the FFLO phase predicted many years before by Fulde and Ferrell [4] and Larkin and Ovchinnikov [5] for superconductors in a strong magnetic field. This phase is characterized by pairing at finite center-of-mass momentum of the pairs, with spatially dependent gap parameter ∆(r) = ∆ Q0 e iQ0·r for the Fulde-Ferrell (FF) and ∆(r) = ∆ Q0 cos(Q 0 · r) for the Larkin-Ovchinnikov (LO) solutions. However, despite many efforts in ultracold gases, some (indirect) evidence for the FFLO phase was obtained so far in one dimension only [6]. In the mean time, substantial progress in the quest for the FFLO phase was made in condensed-matter systems, producing solid evidence of an FFLO phase in quasitwo-dimensional organic superconductors [7][8][9], as well as quite robust evidence in the iron-based multi-band superconductor KFe 2 As 2 [10]. In addition, mounting evidence is emerging in cuprate superconductors [11][12][13][14] about the presence of "pair-density-wave" ordering (corresponding to an LO solution in the absence of a magnetic field), which competes with d-wave superconductivity and possibly explains the pseudogap phase observed in these systems as associated with strong pair-density-wave fluctuations [14].
In the present work, we consider the normal (i.e., not superfluid) phase of a polarized Fermi gas, and search for the presence of FFLO pair fluctuations by calculating the pair susceptibility vs pair momentum Q. Close to unitarity [15], we find that the pair susceptibility is considerably enhanced at a finite value Q 0 of Q(= |Q|) over an extended region of the temperature-polarization phase diagram, such that it should be possible to observe this strong tendency towards FFLO ordering at experi-mentally achievable temperatures (say, a few percent of the Fermi temperature) in the normal phase. To this end, we suggest that the projection technique introduced some time ago to detect pair condensation in a strongly interacting Fermi gas [16,17] could as well be used here to measure a "projected" pair-momentum distribution, whereby the occurrence of a pronounced peak at the same Q 0 of the pair susceptibility would provide unambiguous evidence for strong FFLO pair fluctuations. We explicitly calculate this projected pair-momentum distribution, and conclude that its peak at a finite Q 0 should most readily be observed with a box-like trapping potential [18][19][20][21].
Pair susceptibility − We consider a system of spin-1 2 fermions of mass m mutually interacting through a contact interaction, as described by the Hamiltonian (throughout, the reduced Planck constanth and the Boltzmann constant k B are set equal to unity): whereψ σ (r) is a field operator with spin projection σ = (↑, ↓) and v 0 the bare interaction strength (with v 0 → 0 − when the contact interaction is regularized in terms of the two-fermion scattering length a F [22]). Our aim is to calculate the pair susceptibility χ pair (Q), that describes the tendency of the normal Fermi gas towards superfluid ordering with pair center-of-mass momentum Q. This quantity can be obtained as χ pair (Q) = D(Q, Ω ν = 0) [23], where is the Fourier transform of the response function [24,25] at the Matsubara frequency Ω ν = 2πνT (ν integer).
To obtain meaningful results for T c across the whole BCS-BEC crossover, the chemical potentials entering Γ 0 need to take into account the effects of pairing fluctuations in the normal phase [27] (in the balanced case, a single chemical potential µ ↑ = µ ↓ = µ survives). This is achieved by inverting the density equations n σ = dk (2π) 3 T n G σ (k, ω n )e iωn0 + in favor of the chemical potentials, where G σ is the single-particle propagator dressed by the t-matrix self-energy Σ σ , which is in turn constructed by convoluting the ladder series Γ 0 with a bare G 0 σ [28] (σ standing for the opposite of σ). Comparison with experimental data or Quantum Monte Carlo (QMC) results for several thermodynamic quantities in the balanced case [29][30][31][32][33][34], as well as with recent QMC calculations for the polarized case [35], shows, however, that a more reliable diagrammatic approximation is obtained by using a fully self-consistent t-matrix self-energy [36][37][38] where Γ(Q, Ω ν ) is defined by replacing G 0 σ → G σ on the right-hand side of Eq. (4). By this approach, the pair susceptibility of interest is correspondingly given by which retains two-particle diagrams consistent with the fully self-consistent t-matrix self-energy (6) [23]. Our calculations of the pair susceptibility will be based on Eq. (7) and on the self-consistent solution of Eq. (6). Numerical results − In the following, it will be useful to introduce the (effective) Fermi wave vector k F ≡ (3π 2 n) 1/3 defined in terms of the total density n = n ↑ +n ↓ . The dimensionless coupling (k F a F ) −1 then drives the crossover from the BCS and BEC regimes, which correspond to (k F a F ) −1 < ∼ −1 and (k F a F ) −1 > ∼ +1, respectively, while the crossover region in between spans across the unitarity limit (k F a F ) −1 = 0. A dimensionless pair susceptibility is then defined as χ pair (Q)mk F . Figure 1 reports our results for χ pair (Q)mk F for two coupling values in the crossover region, namely, unitar- ity (left panels) and (k F a F ) −1 = −0.5 (right panels). For each coupling, two different polarizations are also considered. From these results one sees that the pair susceptibility gets strongly enhanced about a finite momentum Q 0 as the temperature is progressively lowered, thereby signaling the presence of strong FFLO fluctuations in the normal phase. This peak at Q 0 is rather well pronounced for temperatures , which are well within the current experimental capabilities with ultracold Fermi gases.
To characterize location and strength of FFLO fluctuations in the temperature-polarization phase diagram, a heat map for the peak value of the dimensionless pair susceptibility χ pair (Q)mk F is presented in Fig. 2, for the same couplings (k F a F ) −1 = 0 (a) and -0.5 (b) considered in Fig. 1. Remarkably, for both couplings the region of the phase diagram where the maximum of the pair susceptibility develops at a finite momentum Q 0 is quite sizable, reaching temperatures as high as about 0.15T F and encompassing a wide range of polarization. In heights from Fig. 1, is indicated by the color code and contour lines. Note that, in the region where the peak of χ pair (Q) occurs at Q 0 = 0, the condition Γ(0, 0) −1 = 0 corresponding to a diverging χ pair (0) can be exactly satisfied (green line). This curve yields the transition temperature T c to a uniform polarized superfluid with pairing at Q 0 = 0 until it becomes reentrant. The reentrant part of the curve (green dashed line) is instead expected to be superseded by phase separation or by an FFLO superfluid phase [40].
When Q 0 = 0, on the other hand, the feedback of a diverging χ pair (Q 0 ) on Eq. (6) would yield a diverging self-energy at finite temperatures for all frequencies and momenta (see Refs. [41][42][43][44][45] and [26] for a discussion of a similar phenomenon in related approaches). Accordingly, within the present approach the condition Γ(Q 0 , 0) −1 = 0 (corresponding to a second-order phase transition to the FFLO phase) can be exactly satisfied only at T = 0. Recall, however, that, in the superfluid phase, FFLO fluctuations would turn FFLO ordering from long-range into algebraic [46], in such a way that determining the transition line by diagrammatic methods would be an extremely difficult task (like for the Berezinskii-Kosterlitz-Thouless transition in a 2D superfluid Fermi gas). Nevertheless, this remark does not hinder our investigation, which fo-cuses rather on the effects of the FFLO fluctuations in the normal phase than on the precise determination of the second-order transition line. In this respect, the experimental results for the unitary Fermi gas (symbols in Fig. 2(a)) indicate that at low temperature the transition between normal and superfluid phases becomes actually of first order (with a phase separation between a balanced superfluid and a normal polarized gas). Our calculations show that strong FFLO pairing fluctuations are expected to occur immediately at the right of this phase separation region found in the experiments. Figure 3 complements the results of Fig. 2, by reporting for the same couplings the peak momentum Q 0 as a function of temperature for a set of polarizations. Here, Q 0 is seen to acquire a finite value in a continuous way when entering the FFLO fluctuation region and to further increase as the temperature decreases (except for a slight decrease at the lowest temperatures). Recall in this context that the mean-field solution for the FF phase suggests that at low T (where sharp Fermi surfaces develop) Q 0 should scale with the mismatch k F↑ − k F↓ between the corresponding Fermi wave vectors. This expectation is confirmed by the comparison shown in the insets of Fig. 3, between the polarization dependences of Q 0 (obtained at the low temperature T = 0.02T F ) and the weak-coupling mean-field result 1.2(k F↑ − k F↓ ) for Q 0 at the T = 0 superfluid-normal transition [47].
Proposed experimental test -One expects the pair susceptibility χ pair (Q) to hardly be measured in a direct way. A somewhat related quantity of more direct access to experiments with ultracold gases should be the "projected" pair-momentum distribution n proj pair (Q), which is the momentum distribution of the molecules formed after a rapid sweep of the magnetic field to the BEC side of the crossover. Measurements of the projected pairmomentum distribution have already been successfully applied to detect condensation (or quasi-condensation) of fermionic pairs across the BCS-BEC crossover, both in three [16,17,21,48] and two [49] dimensions, and were also proposed some time ago to detect FFLO superfluidity in trapped Fermi gases [50]. On physical grounds, we expect that, even in the normal phase, strong FFLO pairing fluctuations should result into a peak of n proj pair (Q) at the same finite Q 0 found for χ pair (Q) [51].
To confirm this expectation, we have extended to the normal phase the theoretical approach for the projected pair-momentum distribution introduced in Ref. [57] for the superfluid phase, borrowing also from the formalism recently used in Ref. [58] to address pair correlations in the normal phase of an attractive Fermi gas [23]. The results for n proj pair (Q) are reported in Fig. 4 and show that, in the strong FFLO pairing-fluctuation region, a pronounced peak develops in n proj pair (Q) at a finite momentum which matches the corresponding peak momentum Q 0 found in Fig. 1 for χ pair (Q). The peak of n proj pair (Q) is quite prominent at T = 0.01T F and remains clearly visible up to T ≈ 0.05T F , which is well within the range of temperatures currently attainable with ultracold gases [59]. The results of Fig. 4 are for a uniform system and thus correspond to an idealization of a Fermi gas trapped in a box-like potential. We have performed calculations [23] also for a realistic box-like potential as the one used in Ref. [18], confirming the observability of the predicted peak even in this case. For a gas trapped in a harmonic potential, on the other hand, we have verified that it would be hard to detect a peak of n proj pair (Q) at a finite Q 0 , due to the smearing produced by the trap average [23].
Conclusions -We have uncovered the presence of strong FFLO fluctuations in the normal phase of a polarized Fermi gas, which could experimentally be observable even in a three-dimensional unitary Fermi gas of most interest for the current research in the field. These fluctuations are precursors of an FFLO superfluid phase, which is competing with the phase separation observed so far in experiments [1,39,60]. The outcome of this competition in the T → 0 limit cannot be established from the experimental data, since phase separation could eventually either give the way to an FFLO superfluid phase or continue to suppress it entirely. Nevertheless, irrespective of which one of these two scenarios would actually take place, our investigation has shown that the effects of an underlying FFLO superfluid phase are clearly visible, and should accordingly be sought for, in the normal phase.
Financial support from the Italian MIUR under Project PRIN2017 (20172H2SC4) is acknowledged.

Pair susceptibility
We consider a system of spin-1 2 fermions of mass m mutually interacting through a contact interaction, as described by the Hamiltonian (with the reduced Planck constanth and the Boltzmann constant k B set equal to unity):Ĥ Here,ψ σ (r) is a field operator with spin projection σ = (↑, ↓) and v 0 is the bare interaction strength (v 0 → 0 − when the contact interaction is regularized in terms of two-fermion scattering length a F [1]).
To quantify the tendency of the normal Fermi gas towards superfluid ordering with pair center-of-mass momentum Q, we consider an associated pair susceptibility χ pair (Q) defined in the following way. We add to the Hamiltonian (8) the symmetry-breaking term whereφ (r) = ∆ (r) +∆ † (r) √ 2 (10) and η(r) is a classical (real) field coupled to the gap op-erator∆(r) = v 0 ψ ↑ (r)ψ ↓ (r) and its hermitian conjugatê ∆ † (r). Within linear-response theory, the pair susceptibility χ pair (Q) is then obtained as the Fourier transform of the local functional derivative Here, Â η = Tr{e −βKÂ } Tr{e −βK } stands for the grand-canonical thermal average of a generic operatorÂ, with the grandcanonical HamiltonianK =Ĥ +Ĥ ext − σ µ σNσ containing the symmetry-breaking term (9), and β = 1/T is the inverse temperature.
When η(r) kept finite (albeit small), the expression (12) implies that a weak probing field of the form η(r) = η cos(Q 0 · r) induces in the normal phase a gap parameter ∆ Q0 (r) proportional to χ pair (Q 0 ) η cos(Q 0 · r), thus signaling that χ pair (Q) quantifies the tendency towards FFLO pairing (with rotational invariance implying that χ pair (Q) is a function of Q = |Q|). Accordingly, when χ pair (Q) is found to diverge at a finite momentum Q 0 , a continuous transition develops from the normal phase to an FFLO phase with pair momentum Q 0 . (Recall in this respect that when the gap parameter ∆ Q0 (r) → 0 from the superfluid phase with η = 0, as expected for a continuous phase transition, the FF and LO solutions become degenerate, such that the transition point is the same for the two phases [2].) More generally, evidence that χ pair (Q) becomes strongly peaked about a finite value Q 0 of Q can be considered as indicating the presence of strong FFLO pairing fluctuations in the normal phase.
The pair susceptibility (12) can be related to an appropriate temperature response function. This is achieved by calculating the functional derivative in Eq. (12) as with η still kept finite, and using the following operator identity e (Â+δÂ)s = eÂ s 1 + s 0 ds ′ e −Âs ′ δÂ e −Âs ′ + · · · (14) to linear order in δÂ. In our case,Â ↔ −K and s ↔ β, such that in Eq. (13) With the definitionφ(r, τ ) = e τKφ (r)e −τK , Eq. (13) becomes eventually: (16) In the normal phase of interest φ(r) η→0 = 0, such that within linear response the local pair susceptibility (12) can be expressed in terms of the temperature response function where T τ is the imaginary-time operator [3]. One obtains: of the main text and Ω ν = 2πν/β (ν integer) is a bosonic Matsubara frequency [3]. To get Eq. (18), homogeneity and isotropy in space and homogeneity in imaginary time have been exploited. In Fourier space, one further obtains that χ pair (Q) = D(Q, Ω ν = 0).

Connection with many-body diagrammatic theory
There remains to implement the calculation of D(Q, Ω ν = 0) in diagrammatic terms. In this context, the simplest physically meaningful approximation results by summing the series of ladder diagrams in the particleparticle channel, where all rungs contain bare singleparticle propagators. Through a careful handling of the v 0 → 0 limit (where the interaction strength v 0 enters the definition of the gap operator∆(r)), one can show that D(Q, Ω ν ) = Γ 0 (Q, Ω ν ) with the bare pair propagator Γ 0 given by Eq. (4) of the main text (cf. also Ref. [4]).
One knows, however, that an improved description of thermodynamic properties of a Fermi gas with attractive interaction results (at least in the balanced case) within the fully self-consistent t-matrix approximation [5], where the bare Γ 0 is replaced in the single-particle self-energy Σ by the dressed Γ which contains the fully self-consistent single-particle propagators G in the place of the the bare G 0 . It is then natural to replace Γ 0 with Γ also in the expression of χ pair (Q), as it was done in Eq. (7) of the main text.
In this context, the question naturally arises about the need to introduce additional diagrammatic contributions to the temperature response function (17) once the single-particle self-energy Σ is taken within the fully self-consistent t-matrix approximation. Specifically, we are referring to contributions with the topology of the Aslamazov-Larkin (AL) and Maki-Thompson (MT) diagrams, which, at the level of the fully self-consistent tmatrix approximation for the single-particle self-energy here adopted, would influence the (two-particle) response functions through the particle-hole channel. For the temperature response function (17) of interest here, this would be the case if it were calculated in the superfluid phase below the critical temperature T c , where there is no clear distinction between particle-hole and particleparticle channels due to the particle-hole mixing characteristic of the BCS (pairing) theory. However, since we are limiting ourselves to the normal phase above T c , we can readily adapt to the present context the argument described in Appendix A of Ref. [6], and show that above T c the AL and MT contributions cannot affect the temperature response function (17) due to its explicit spin structure and its ultimate particle-particle character.

Formalism and calculations for a homogeneous system
The projected pair-momentum distribution can be obtained as follows in terms of the composite-boson propa-gator G B (Q, Ω ν ) [7,8]: Within the self-consistent t-matrix approach, G B (Q, Ω ν ) is, in turn, given by [6] where F j (Q, Ω ν ) (j = 1, 2) are "form factors" defined as Here, φ proj (p) is the molecular wave function onto which the initial correlated pairs are projected during a magnetic sweep, while G σ (p, ω n ) and Γ(Q, Ω ν ) are the self-consistent single-particle Green's functions and the particle-particle (ladder) propagator defined in the main text.
The analysis of Ref. [8] for projection experiments established that projection onto molecules occurs at a coupling (k F a F ) −1 proj on the BEC side of the crossover, in order to optimize the overlap between the initial pair correlations and the molecular wave function. The specific value of the projection coupling (k F a F ) −1 proj depends on the experimental conditions of the magnetic sweep, and it was estimated to be generally in the range 0.5 < ∼ (k F a F ) −1 proj < ∼ 1.5 [8]. Accordingly, for the calculation of n proj pair (Q) shown in Fig. 4 of the main text we have considered a projection coupling (k F a F ) −1 proj = 1 in the middle of the above range. In addition, following again the procedure of Ref. [8], we have taken the molecular wave function φ proj (p) as the normalized two-body bound-state wave function in vacuum at the projection coupling (k F a F ) −1 proj : We have further verified that n proj pair (Q) depends only weakly on the projection coupling (k F a F ) −1 proj . This is evident by comparing the results for n proj pair (Q) shown in the three panels in Figure 5, which correspond to different values of the projection coupling spanning the whole above range.

Calculations for a harmonic trapping potential
We next consider a trapped Fermi gas in a harmonic potential, which can be described by a local-density approximation. Accordingly, wherever they appear, we replace the chemical potentials µ σ of the two σ = (↑, ↓) components by the local chemical potentials µ σ (r) = µ 0σ − V (r), where V (r) = mω 2 0 r 2 /2 is the harmonic trapping potential with frequency ω 0 . By this replacement, the single-particle propagators G σ (k, ω n ; r) and the composite-boson propagator G B (Q, Ω n ; r) become local functions of the position r in the trap through the local chemical potentials µ σ (r). by inverting the number equations while the total projected pair-momentum distribution N proj pair (Q) is obtained in the trap by summing the local projected pair-momentum distribution n proj pair (Q; r) over r N proj pair (Q) = dr n proj pair (Q; r) where n proj pair (Q; r) = − T ν G B (Q, Ω ν ; r)e iΩν 0 + .
To obtain an intensive quantity, one can then divide the total projected pair-momentum distribution of Eq. (24) by the volume N/(k t F ) 3 , where k t F = 2mE t F is the (effective) trap wave vector and E t F = ω 0 (3N ) 1/3 is the (effective) trap Fermi energy, and obtain the projected pair-momentum distribution in the harmonic trap: Figure 6 shows the results of n pair,h proj (Q) obtained for a harmonic trap at unitarity and different temperatures, with global polarization (N ↑ − N ↓ )/N = 0.85 where N = N ↑ +N ↓ . This value corresponds to a local polarization (n ↑ −n ↓ )/n ≃ 0.40 at the trap center, which matches the value of the polarization considered in Fig. 4(a) of the main text. Here, the coupling parameter (k t F a F ) −1 and the pair momentum Q are expressed in terms of the (effective) trap Fermi wave vector k t F , and the temperature T is expressed in terms of the (effective) trap Fermi temperature T t F = E t F . By comparing the results of Fig. 6 for the trapped system with the corresponding results shown Fig. 4(a) in the main text for the uniform system, one is led to conclude that the effect of the harmonic trap is to strongly smear the peak at finite Q which was clearly visible for the uniform system, making it hard to be detected experimentally. This smearing can be justified by extending to the trap the results reported in the inset of Fig. 3(b) of the main text, whereby the local value of the pair momentum Q 0 (r) at which n proj pair (Q; r) is peaked is bound to scale with |k F ↑ (r)−k F ↓ (r)|, which in turn varies along the trap since both local densities and polarization change with r. The net effect is that Eq. (24) effectively averages over locally projected pair-momentum distributions with different peak momenta Q 0 (r), thus smearing out the resulting peak in the projected pair-momentum distribution n pair,h proj (Q) for the harmonic trap with respect to the corresponding peak that would be present for the uniform case.

Calculations for a realistic box-like trapping potential
It is relevant to calculate the projected pair-momentum distribution within a local-density approximation for a realistic box-like trapping potential. To this end, we consider the same kind of cylinder-shaped potential utilized in Ref. [9]. The trapping potential in cylindrical coordinates is given by V (ρ, z)/E F (0) = (ρ/R) p + α z (z/L) 2 (in units of the (effective) Fermi energy E F (0) = [3π 2 n(0)] 2/3 /(2m), where n(0) is the total density at the trap center). Here, V (ρ, z) is the sum of a radial powerlaw potential (∼ ρ p ), which describes the confinement of the ring beam, and of a weak axial harmonic potential (with α z ≪ 1), which is needed to obtain the momentum distributions of atoms and pairs in a quarter period timeof-flight expansion [9]. A hard-wall potential at z = ±L is also included to describe the light sheets acting as the end caps for the axial trapping. In Ref. [9] it was estimated that p ≃ 16 for the power-law exponent and α z ≤ 0.05 for the axial harmonic confinement parameter, while L ≃ R.
To perform a direct comparison with the uniform case, the (intensive) projected pair-momentum distribution for the box-like trap is obtained as n proj,b pair (Q) = 1 V 0 dr n proj pair (Q; r) , where the local n proj pair (Q; r) is defined like in Eq. (25) and V 0 = N/n(0) is the volume that would be occupied by the N = N ↑ + N ↓ particles for a uniform system with density equal to the density n(0) at the trap center. This is because, for a completely uniform trap such that n proj pair (Q; r) = n proj pair (Q) independent of r, the integration  Figure 7 shows the results for the projected pairmomentum distribution in a box-like trap with p = 16 and L = R, corresponding to different values of the axial harmonic confinement parameter α z ≤ 0.05 like in Ref. [9]. For comparison, the corresponding result for the uniform system is also shown (cf. Fig. 4(a) of the main text). The calculations are performed by taking the local temperature and polarization at the trap center to coincide with those of the uniform system. From this comparison one concludes that the peak originating at finite momentum in the projected pair -momentum distribution for a uniform system still appears to be quite prominent even when the effects of a realistic box-like trapping potential are taken into account (and this is especially true when the axial harmonic trapping potential is kept sufficiently weak like in the experiment of Ref. [9]). The favorable outcome of this specific test considerably reinforces our expectation as discussed in the main text, that the presence of strong FFLO fluctuations in the normal phase of a polarized Fermi gas can be uncovered by measurements of the projected pair-momentum distribution under realistic experimental conditions.
[1] P. Pieri and G. C. Strinati, Strong-coupling limit in the evolution from BCS superconductivity to Bose-Einstein