Pseudogap and preformed pairs in the imbalanced Fermi gas in two dimensions

The physics of the pseudogap state is intimately linked with the pairing mechanism that gives rise to superfluidity in quantum gases and to superconductivity in high-Tc cuprates, and therefore, both in quantum gases and superconductors, the pseudogap state and preformed pairs have been under intensive experimental scrutiny. Here, we develop a path integral treatment that provides a divergence-free description of the paired state in two-dimensional Fermi gases. Within this formalism, we derive the pseudogap temperature and the pair fluctuation spectral function, and compare these results with the recent experimental measument of the pairing in the two-dimensional Fermi gas. The removal of the infrared divergence in the number equations is shown both numerically and analytically, through a study of the long-wavelength and low-energy limit of the pair fluctuation density. Besides the pseudogap temperature, also the pair formation temperature and the critical temperature for superfluidity are derived. The latter corresponds to the Berezinski-Kosterlitz-Thouless (BKT) temperature. The pseudogap temperature, which coincides with the pair formation temperature in mean field, is found to be suppressed with respect to the pair formation temperature by fluctuations. This suppression is strongest for large binding energies of the pairs. Finally, we investigate how the pair formation temperature, the pseudogap temperature and the BKT temperature behave as a function of both binding energy and imbalance between the pairing partners in the Fermi gas. This allows to set up phase diagrams for the two-dimensional Fermi gas, in which the superfluid phase, the phase-fluctuating quasicondensate, and the normal state can be identified.


I. INTRODUCTION
Ultracold atomic gases are increasingly used as quantum simulators to probe many-body physics [1]. Recent efforts have focused in particular on understanding superfluidity and Cooper pairing in interacting Fermi systems, including i.a. the effects of varying the interaction strength, introducing population imbalance, and reducing the dimensionality. When these fermionic superfluids are described in the path integral formalism, the thermodynamic potential Ω(T, µ ↑ , µ ↓ ) of the interacting Fermi gas (as a function of temperature T and chemical potentials µ ↑ ,µ ↓ of spin-up and spin-down components) is rewritten as a functional integral over a bosonic field ∆ x,τ that represents the field of the pairs. This field is introduced through the Hubbard-Stratonovich transformation that allows the exact elimination of the fermionic degrees of freedom and results in an action functional for the bosonic field [2]. At first sight, one has only succeeded in rewriting an unsolvable functional integral over fermionic fields by an equally unsolvable functional integral over the bosonic fields. However, the bosonic field lends itself to an obvious simplification when one intuits that a (uniform) Bose-Einstein condensation of pairs is present. In that case, one can surmise that ∆ x,τ ≈ ∆, i.e. all pairs are in the zero-momentum state so the field is a constant in real space. The functional integral can then replaced by its saddle-point value, substituting ∆ x,τ ≈ ∆ and dropping the integrations. The optimal value of ∆ is found by extremizing the action or, equivalently, by minimizing the thermodynamic potential Ω sp (T, µ ↑ , µ ↓ ; ∆) with respect to the saddle point. When ∆ = 0, the (non-interacting) normal Fermi gas is obtained, when ∆ = 0, the saddle-point approximation bears out a Bose-Einstein condensate of pairs.
In two-dimensional superfluid systems, the interpretation of the bosonic field is more subtle. To be precise, the condition for Bose-Einstein condensation (BEC) has been defined by Penrose, Onsager [5] and Yang [6] as the presence of off-diagonal long range order, i.e. lim x−x ′ →∞ ∆ x,τ ∆ x ′ ,τ = 0. At the level of the saddle point, ∆ x,τ ≈ ∆, it is clear that a nonzero ∆ implies BEC. However, as Mermin,Wagner [3] and Hohenberg [4] pointed out, in the two-dimensional system fluctuations play a crucial role: they will prohibit off-diagonal long range order in uniform systems. These fluctuations around the saddle point are commonly taken into account through the Bogoliubov shift ∆ x,τ = ∆ + φ x,τ , we discuss finite-temperature phase diagrams for the imbalanced Fermi gas in 2D. In Sec. V, the theory is applied to the interpretation of the experiment on pairing of cold atoms in 2D. The discussion is followed by conclusions, Sec. VI.
Before going forward in the next section with presenting the functional integral approach in the GPF framework, it is useful to note that the GPF approach that we follow here is not the only way to avoid the divergence problem that occurs in the NSR description for the Fermi gas in two dimensions. The NSR scheme and its modifications are related to the T -matrix perturbation approach, in which the effective interaction between pairs is taken into account diagrammatically. In this context, the divergence problem for a Fermi gas in two dimensions can be remedied by taking into account higher orders of the T -matrix expansion -via an effective interaction between pair fluctuations [29] . This interaction stabilizes the superfluid phase of the 2D fermion system at very low temperatures. However, in 2D the T -matrix method does not predict the universal jump in the superfluid density [30] related to the BKT phase transition. A correct description of the superfluid density becomes possible by explicitly focusing on the phase fluctuations, as in the approach of Refs. [26][27][28]31]. Within that approach, bosonic pair field is gauge transformed ∆ x,τ e iθx,τ and a subsequent gradient expansion of the fluctuation action is performed for phase fluctuations assuming that phase gradients are small. This leads to a quadratic effective action functional of the phase field θ, which, as distinct from the scheme of Ref. [11], has no divergence for ∆ = 0. As far as gradients of the fields are assumed to be small, the resulting effective action is treated as a hydrodynamic action (see, e.g., Ref. [32]). The gradient expansion does not contain the a priori assumption that fluctuations themselves are small. In this connection, the method was categorized in Ref. [26] as non-perturbative. In Ref. [33], the present authors applied the method of Refs. [26][27][28]31] to derive the effective hydrodynamic action for a Fermi gas with a population imbalance. A non-perturbative approach was also the key to develop a description free of infrared and ultraviolet divergences for the 2D Bose gas [34], that successfully describes the crossover between the mean-field regime and the critical fluctuation range corresponding the BKT transition [35].

II. THERMODYNAMIC FUNCTIONS OF THE FERMI GAS IN 2D
A. Gap equation We consider a gas of interacting fermions in 2D, with a contact interaction and with s-wave pairing. In the ultracold regime where only s-wave interactions matter, these interactions only take place between "spin-up" and "spin-down" fermions (in practice, these are usually two different hyperfine states of an atomic species). The thermodynamic functions of the Fermi gas are completely determined by the partition function. Here we will focus on the thermodynamic potential Ω per unit area. The treatment is performed within the path-integral formalism following Ref. [33], building on the original path-integral treatment in Ref. [2] for the case of a balanced three-dimensional Fermi gas. The partition function is represented as the path integral over Grassmann variablesψ σ (x, τ ) , ψ σ (x, τ ), The action functional of interacting fermions is given by the integral where g is the interaction strength and β = 1/(k B T ) is the inverse thermal energy. We choose a system of units where = 1, 2m = 1, and the Fermi wave vector k F ≡ (2πn) 1/2 = 1 with n the total density. Here, we consider also the case when imbalance is present, i.e. the number of spin-up and spin-down atoms are unequal: n ↑ = n ↓ . This in turn implies that the chemical potentials µ ↑ and µ ↓ should be fixed separately. Rather than contemplating the separate components, we will work with the total density n = n ↑ +n ↓ and the density difference δn = n ↑ −n ↓ . Correspondingly, we will use the average chemical potential µ = (µ ↑ + µ ↓ )/2, and the chemical potential difference ζ = (µ ↑ − µ ↓ )/2. Note that the total density is equal to 1/(2π) in our units, so this means that we need to solve the number equation to fix µ (in the non-interacting case, µ = 1 in our units). Only with respect to the imbalance, we have a choice of studying the free energy (and making phase diagrams as a function of δn) or the thermodynamic potential Ω (and making phase diagrams as a function of ζ). The thermodynamic potential is linked to the free energy by the usual Legendre transform and, as mentioned, in our formalism this corresponds to imposing the number equations.
The strength g of the contact interaction is renormalized as in Refs. [36,37] using the binding energy E b for a two-particle bound state, which always exists in 2D [43,44]: with δ a positive infinitesimal number. The BCS regime corresponds to E b /E F ≪ 1, whereas the BEC regime corresponds to the opposite ratio E b /E F ≫ 1. Similarly to Ref. [2], we introduce the pair field ∆ x,τ and perform the Hubbard-Stratonovich transformation, which results in a fermion-boson action quadratic in fermion variables. After integrating out the fermion variables, the following effective bosonic action is obtained as a functional of the Hubbard-Stratonovich pair field ∆ x,τ : where G −1 is the inverse of the Nambu propagator Here, σ j are the Pauli matrices. As far as the effective action S ef f is not a quadratic functional of the Hubbard-Stratonovich pair field, the resulting functional integral over the pair field cannot be calculated analytically exactly. As in the analogous problem in 3D [2,38,39], and as explained in the introduction, we consider approximations provided by an expansion of the effective action S ef f over fluctuations of the pair field ∆ x,τ about its saddle-point value ∆. The phase diagrams of a 2D Fermi gas in the saddle-point approximation have been investigated in Refs. [37,40]. The effective saddle-point action provides the thermodynamic potential per unit area: Here, ξ k = k 2 − µ is the fermion energy, and E k = ξ 2 k + ∆ 2 is the Bogoliubov excitation energy. The gap parameter ∆ is determined from the gap equation generalized to the imbalance case -the minimum condition for the saddle-point thermodynamic potential as a function of the gap parameter ∆ at fixed temperature and chemical potentials: For high temperatures (T > T * c ) or at high levels of imbalance (ζ > ζ c ), thermodynamic potential will have its minimum at ∆ = 0, the unpaired normal state. Following the experimental observation of superfluidity in imbalanced Fermi gases in 3D [41], the phase diagram of the imbalanced Fermi gas has attracted a lot of attention (for a recent review, see Ref. [42]). To find the phase diagrams in 2D, the above gap equation has to be solved in conjunction with the number equations discussed in the remainder of this section.

B. Gaussian fluctuations
The next-order approximation brings into account fluctuations about the saddle point: We apply the Fourier expansion of the fluctuation coordinates: where L is the linear size of the 2D system, and ω n = 2πn/β (with n = 0, ±1, ±2, . . .) are the bosonic Matsubara frequencies. The quadratic fluctuation contribution to the effective bosonic action is the functional of complex fluctuation coordinates similar to that derived in Ref. [39]: where M j,k (q, iω n ) are the matrix elements of the inverse pair fluctuation propagator. They are determined by the expressions (cf. Ref. [39]): and Here, the following function has been introduced, The integration over fluctuation coordinates gives us the fluctuation contribution Ω f l (T, µ, ζ; ∆) to the total grandcanonical thermodynamic potential Ω per unit area: with

C. Number equations and the GPF approach
The fermion density and the density difference fix the chemical potentials µ and ζ through the derivatives of the total thermodynamic potential per unit area: (remember that in our units n = 1/2π). We can write out these equations by splitting the total thermodynamic potential in saddle point and fluctuation contributions.
We will denote the first and second terms in the right hand side (RHS) of expression (19) for n as n sp and n f l , respectively. Similarly, the terms in the RHS of expression (20) will be denoted by δn sp and δn f l . Note that the thermodynamic potentials obtained from expressions (7) and (15) are expressed as a function not only of T, µ, ζ, but also of ∆. This gap ∆ is not an independent thermodynamic variable, and when considering Ω(T, µ, ζ, ∆) explicitly as a function of also ∆, the implicit dependence of ∆ on the chemical potentials must be taken into account in (19), (20): Note that the gap equation ∂Ω sp /∂∆ = 0 at fixed T, µ, ζ implies that the implicit dependence of ∆ on the chemical potentials will only affect the fluctuation part of the thermodynamic potential in the above equations. Different theories of the BEC-BCS crossover, in any dimension, can be categorized by their choice of number and gap equations. The simplest mean field approach only keeps the terms with Ω sp . The Nozières and Schmitt-Rink approach also includes the second terms in the RHS of expressions (21). Finally, the Gaussian pair fluctuation approach includes also the last term in the RHS of expressions (21). Note that in the literature, there is no common opinion on which approach is best. For example, on the one hand, Randeria et al. [45,46], Hu et al. [13], Keeling et al. [47] state that the derivative over µ must be performed taking into account a variation of the gap determined by the gap equation.
On the other hand, Ohashi et al. [48][49][50], and Strinati et al. [51,52] use the other definition, considering ∆ in the number equations as an independent variable and, therefore, applying the gap equation after taking the derivatives ∂Ω/∂µ. In the papers [48][49][50], it is stated that the last terms in (21) are the higher-order corrections with respect to Gaussian quadratic fluctuations. Keeling et al. [47] correctly argue that both terms in those derivatives are of one and the same order and emphasize that the existence of the second term is crucial in two dimensions. Below, we demonstrate the key significance of taking into account of the last terms (21) for the convergence of fluctuation contributions to the fermion density in 2D.
As stated in the introduction, in order to treat the fluctuations, there is an alternative to the Bogoliubov shift ∆ x,τ = ∆ + φ x,τ , namely the parametrization in amplitude and phase fluctuations ∆ x,τ ≈ ∆ (1 + δ x,τ ) e iθx,τ . When, after this parametrization, the effective action is expanded with respect to δ x,τ and θ x,τ (rather than φ x,τ andφ x,τ ), this leads to another quadratic fluctuation action S ′ f l which differs from expression (11) for S f l only by terms that vanish when applying the gap equation. Correspondingly, the thermodynamic potentials Ω f l and Ω ′ f l provided by those two actions lead to one and the same contribution to the fermion density. Furthermore, keeping in S ′ f l only the leading order long-wavelength and low-energy terms leads to the same effective "hydrodynamic" action as in Ref. [33]. In the particular case of a balanced gas, the effective action of Ref. [33] turns to the result of Refs. [26,31]. This means that the effective action described as the result of the non-perturbative approach in Ref. [26] can be equivalently re-derived within the perturbative NSR-like scheme. Moreover, the present treatment can be considered as an extension of the approach of Refs. [26,31,33] beyond the long-wavelength and low-energy approximation (and to imbalanced 2D gases). The hydrodynamic action is particularly useful in extracting the superfluid density ρ s , by identifying it with the prefactor of the (∇θ x,τ ) 2 /2 term in the expression for S ′ f l . This identification yields straightforwardly [33]: with X(E k ) given by expression (14) and X ′ (E k ) its first derivative, evaluated in E k . Once ∆, µ, ζ are obtained for a given temperature (and interaction strength) by solving the gap and number equations, they can be substituted in this expression to determine whether the system is in the superfluid phase (ρ s = 0) or the normal phase (ρ s = 0). As discussed in the results section, we also use this expression to find the temperature T BKT of the phase transition between those two states. Already we note that ∆ = 0 leads to ρ s = 0, so that T BKT < T * c and the superfluid state requires pair formation, as it should.

D. Pair fluctuation spectral functions
From expression (7) for the saddle-point thermodynamic potential, we derive the following expressions for the saddle-point densities: Similarly, the fluctuation contributions to the fermion densities are determined using (15) and using the GPF approach.
The results can be written as a sum over wavelengths and frequencies of pair fluctuation structure factors: The pair fluctuation structure factors J and K are given by We transform the Matsubara summations in (25) and (26) to the contour integrals in the complex plane as follows: This allows to express the resulting fluctuation contributions to the fermion density through the distribution functions for pair excitations: The fluctuation distribution functions g n (q) and g δn (q) are the integrals over the frequency with the pair fluctuation structure factors: The functions g n (q) and g δn (q) are proportional to the densities of states for the pair fluctuations. The behavior of these functions is crucial for understanding of the pseudogap properties and of different phase transitions in the imbalanced 2D Fermi gas. In the NSR approach, the fluctuation distribution functions have a divergency, and as a consequence no value of the chemical potential µ can be found so that the number equation n sp + n f l = n = 1/(2π) is satisfied. In the GPF approach, the divergency is overcome and the number equation can be satisfied. In order to demonstrate this, we focus in the next section on the long wavelength limit where exact analytic expressions for the distribution functions are obtained.

A. Long-wavelength expansion
In order to investigate the problem of the long-wavelength convergence for the fluctuation contributions to the density, it is necessary to derive analytically the spectrum of low-lying and long-wavelength pair excitations. For this purpose, we expand the matrix elements of the inverse pair fluctuation propagator M j,k (q, z) in powers of (q, z) up to the second-order terms in power of z and q, The derivation of the coefficients is rather tedious. In this connection, here only the final results are represented. The coefficients in the M 1,1 expansion are given by and those in the M 1,2 expansion are given by In these expressions X ′ ,X ′′ and X (3) are the first, second and third derivatives of the function X given by expression (14), with respect to its argument. In the literature, an analogous expansion was performed for 3D at finite temperatures in Refs. [48][49][50] in the strong-coupling limit, and in Ref. [45] at low temperatures. The present expansion is all-coupling and all-temperature, because no restriction is imposed on the thermodynamic parameters.

B. Structure factor
Let us substitute the expansions (34), (35) to the structure factor J (q, z) in order to obtain its long-wavelength and low-energy form J lw (q, z). The spectrum of pair bosonic excitations is determined by the poles of J lw (q, z). In the long-wavelength and low-energy range, these roots are z = ±ω q with the pair excitation frequency ω q which satisfies the Goldstone theorem, with the parameters (cf. the analogous expansion in the 3D case, Ref. [53]) The parameter v has the dimensionality of velocity and tends to the first sound velocity in the low-temperature limit . As far as the parameter D is proportional to ∆ 2 , the velocity parameter for pair excitations turns to zero at the phase boundary when ∆ = 0. In the BEC limit, when E b ≫ 1, we find that µ → −E b /2 and κ → 1/4. This results in the pair excitation spectrum ω q → q 2 /2 at the BEC side.
In order to calculate the long-wavelength distribution function, we keep the lowest-order terms in powers of z and q 2 in the numerator of J lw (q, z). This gives us the result where the coefficients a q and b q are related to the constants determined above as with the notations Here A µ , B µ , . . . are the derivatives A µ ≡ ∂A/∂µ, etc. The distribution function is calculated setting z = ω + iδ with δ → +0. This gives us the structure factor as a superposition of the delta functions. The distribution function then takes the form For the other distribution function, g δn (q), the derivations are the same, but with a replacement of the derivatives over µ by the corresponding derivatives over ζ.

C. Distribution functions in the paired state
Here we consider the paired state of the quasicondensate in which the gap parameter ∆ = 0. In this case, the gap parameter obeys the gap equation: The strength g of the contact interaction is expressed through the two-particle binding energy E b in 2D by the equation (3). The difference of coefficients A − D is proportional to the left hand side (LHS) of the gap equation. Therefore, as long as the gap equation is satisfied, we obtain D = A. Moreover, because the derivatives of matrix elements are calculated while keeping the gap equation satisfied, we find that A µ = D µ and A ζ = D ζ .This implies in particular that the coefficient α = 0, as is evident from expression (47), and we find: Thus g (lw) n (q) tends to a finite value at q → 0 for ∆ = 0, and behaves as q −2 at q → 0 for ∆ = 0. As a result, the fluctuation contributions n f l and δn f l in 2D are finite at ∆ = 0. They can diverge only at ∆ = 0. This is to be contrasted with the NSR scheme, where the order parameter ∆ is treated as independent variable, and where n f l and δn f l in 2D diverge for all ∆.   1 shows the behavior of the fluctuation distribution functions g n (q) and g δn (q) for different temperatures, at binding energy E b = 0.5 and at the critical value of the chemical potential imbalance ζ = ζ c (E b , T ). The critical value ζ c for a given (E b , T ) is determined as the highest imbalance at which the order parameter ∆ is other than zero. The dashed lines correspond to the NSR scheme and reveal a q −2 long wavelength divergence. The full lines show the results in the GPF scheme, where the long-wavelength divergence is absent. This behavior is seen both for g n (q) and g δn (q). In the limit ∆ → 0, the functions qg n (q) and qg δn (q) become logarithmically divergent. However, the sign of this divergence is opposite to that of the divergence of the functions calculated neglecting the variation of ∆.
For the lower temperature shown in Fig. 1, T /T F = 0.1, g n (q) remains positive, whereas for the higher temperature T /T F = 0.3, there is a sign change in g n (q) as it becomes negative for small q. Regions of negative value for the fluctuation distribution function g δn (q) are expected, as the sign will change depending on which species is the majority species. However g n (q) is expected to remain positive, as it is proportional to the pair fluctuation density of states. The appearance of a long-wavelength instability heralds the breakdown of the paired state. We can track the onset of this instability by studying g n (q → 0) as a function of temperature. At low temperatures, g n (0) is positive, and the fluctuation density function remains positive for all q. At high temperatures g n (0) becomes negative, signalling the long-wavelength instability. We denote the temperature separating the two regions by T p , and find this temperature through solving g n (0) = 0 with respect to temperature. The behavior of g n (0) as a function of temperature is shown in Fig.2, for different values of the binding energy E b and both for balanced and imbalanced systems. The function g n (0) diverges when the temperature achieves the limit T = T * c at which ∆ = 0. This result explicitly follows from the analytic properties of the long-wavelength expansion of the distribution functions as discussed above. The temperature T p at which g n (0) = 0 lies below the temperature T * c where we find ∆ = 0, and above the critical temperature T BKT for superfluidity.
The temperature T p does not correspond to a phase transition, because the gap equation is satisfied with a finite density both below and above T p . Nevertheless, because T p is the temperature at which the fluctuation density of states changes it qualitative behavior, we hypothesize that T p corresponds to a crossover between the normal and pseudogap states. This will be further substantiated by comparing our spectral functions to the experimental ones in sec. V. The joint solution of the gap and number equations within the GPF theory then formally provides a non-superfluid quasicondensate at temperatures below T p . Indeed, for temperatures T BKT < T < T p the phase coherence is destroyed by the phase fluctuations according to the BKT mechanism, resulting in the phase fluctuating quasicondensate discussed by Kagan [7]. Through the interpretation of the spectral function, we will denote this temperature region as the "pseudogap regime". It is worth noting that the total fermion density within the GPF theory is finite at T = T p without the necessity to introduce any cutoff in the integrals over q. In the next section we set up phase diagrams identifying the regions where the superfluid phase and the non-coherent paired phase occur.

IV. PHASE DIAGRAMS
In order to get the complete set of equations for phase diagrams, the number equations (21), (21) and the generalized gap equation (8) are solved jointly with the equation for the BKT transition temperature T BKT determined by [30] where ρ s is the superfluid pair density given by Eq. (22). To investigate the phase transitions for the Fermi gas in 2D for different binding energies, we have calculated the critical temperatures of the BKT phase transition T BKT and the critical temperature T p below which the phase fluctuating quasicondensate is formed, as a function of the binding energy E b . Because the fluctuation contribution to the density is finite at T p and at T BKT , these temperatures can be self-consistently determined from the joint solution of the gap and number equations with the complete thermodynamic potential Ω = Ω sp + Ω f luct . The phase diagrams in Fig. 3 show the critical temperatures for cold fermions in 2D as a function of the binding energy E b for the balanced case (panel a) and for the chemical potential imbalance ζ = 0.5 (panel b). The formation of the superfluid state is indicated by the critical temperature T BKT of the BKT phase transition. The pseudogap temperature T p is the upper bound for the existence of the phase fluctuating quasicondensate described in the previous section. We also show the mean-field temperature for pair formation, T * c , obtained by solving gap and number equations with Ω = Ω sp . According to Ref. [54] (for 3D), the unitary gas can exist in the normal state with pairing correlations called preformed pairs which survive at temperatures up to this T * c . The critical temperatures for the balanced case, Fig. 3 (a), were calculated in Ref. [55]. Here, they are reproduced in order to compare them with those for a nonzero imbalance. The population imbalance brings new features to the phase diagram: a phase separation region (between T * c1 and T * c2 ) and a minimum binding energy E b,cr required for superfluidity. For the balanced Fermi gas the superfluid state exists for any value of the binding energy: the BKT critical temperature as well as other critical temperatures gradually decrease with decreasing E b , remaining always finite. However, when ζ = 0 a minimum value of the binding energy E b,cr is required for superfluidity to exist even at T = 0. As shown in Fig. 3, the pseudogap temperature T p does not grow unboundedly when increasing the binding energy E b . For ζ = 0.5 it achieves its maximum at around E b ≈ 4 and then slowly decreases tending to a finite value. Consequently, in the strong-coupling regime the pseudogap temperature is suppressed with respect to the mean-field prediction, where it is often identified with our pair formation temperature T * c , as in [54]. This behavior is qualitatively similar to that for the critical temperature T c as a function of 1/a s for the cold fermions in 3D obtained first in Ref. [2] accounting for the Gaussian fluctuations.
The critical temperatures T * c1 and T * c2 coincide with each other in the balanced case, and they can be different in the imbalanced case: the area between T * c1 and T * c2 is the "phase-separated state". In the phase-separated state, uniform phases are not possible. The temperatures T * c1 and T * c2 were already calculated in Ref. [33]. The temperature T p is determined for the state with ∆ = 0. Therefore a non-zero imbalance does not lead to a splitting of this critical temperature. However, a tricritical point appears at T p = T BKT in the phase diagram joining three regions: the superfluid state, the pseudogap regime and the normal state. This tricritical point is rather conventional as far as the pseudogap temperature indicates a crossover rather than a sharp transition. At zero imbalance, T p > T BKT , and the phase coherence in the range T BKT < T < T p is destroyed by phase fluctuations that lead to a phase fluctuating quasicondensate. However, at nonzero imbalance, there is a region where pseudogap temperature crosses the BKT temperature for superfluidity. This result is interesting in connection with recent experiments on high-T c superconductors [56], that show a crossing of the zero-field superconducting transition temperature and the temperature indicating the opening of the pseudogap in overdoped La 2−x Sr x CuO 4 . The crossing of pseudogap temperature and BKT temperature is also seen in Fig. 4, showing the phase diagram in the variables (T, ζ), for the binding energy E b /E F = 0.04. Here, the same critical temperatures and phase regions are identified as in Fig. 3(b). Increasing imbalance is not only detrimental to the superfluid phase, it also suppresses the pseudogap regime.

V. COMPARISON WITH EXPERIMENT
In the experiment [17] on pairing of cold fermions in two dimensions, the single-particle spectral function A (q, ω) is measured for different values of the wave number q. The spectral function exhibits peaks whose positions indicate the energies of the pair excitations. In the strong-coupling regime, these energies are close to the pair binding energy E b . However, as stated in the paper, some discrepancies remain between the peak positions observed in the experiment and those predicted by the mean-field theory. The deviation "could stem from beyond mean-field effects provoked by our two-dimensional geometry and interaction energy shifts" [17].
In the GPF approach, the pair fluctuation contribution of the fermion density is expressed through the integral (32), where the structure factor J (q, ω) describes the spectrum of the pair excitations of the fermion system. Thus there should be a correspondence of the peaks of the structure factor J (q, ω) with the peaks of the spectral function A (q, ω). In this connection, we compare the positions of the peaks of the spectral function measured in Ref. [17] with those of the structure factor calculated within the GPF approach. The results are shown in Fig. 5 for q = 0 and T /T F = 0.27, where k B T F ≡ E F . The 2D scattering length a 2D is related to the binding energy E b as a 2D = / √ mE b . The value of the Fermi wave vector taken from Ref. [17] is k F = 8.1 µm −1 . When using the mass of the fermion atom m ≈ 39.964 u, we found that the frequency ν F ≡ E F / (2π ) corresponding to the Fermi energy is ν F = 8.2967 kHz.
For the visualization of the peaks of the structure factor, we have used J (q, ω + iγ) with a finite damping parameter γ (as in Refs. [38,39], where this parameter was introduced to facilitate the numeric calculations). Here, the value The parameters of the state (the chemical potential µ and the gap parameter ∆) are determined for each plot from the joint solution of the gap and number equations. In the number equation, the Gaussian fluctuations are included within the GPF formalism. The GPF method provides a finite (convergent) pair fluctuation contribution for any finite ∆ without any cutoff for the pair momentum. This is to be contrasted with the standard NSR scheme which leads to a divergence of the fluctuation contribution at any ∆. Therefore the standard NSR scheme cannot be used for the description of the pseudogap state, whereas the GPF approach can describe this regime.
FIG. 5: Full dots: measured energy distribution curve A (q = 0, ω) for ln(kF a2D) = 0.8 from Ref. [17]. The red solid line is the fit by elementary functions to the experimental data performed in Ref. [17]. The black solid line: the structure factor J (q = 0, ω) calculated in the present work within the GPF approach.
In Fig. 5, the high peak at ω = 0 in our results has no relation to the energies of the pair excitations: it is an intrinsic feature of the structure factor. The other peak of our structure factor at ω < 0 is positioned remarkably close to the measured peak of the spectral function attributed to the pair excitation energy in Ref. [17], especially for the relative high coupling strength at ln (k F a 2D ) = 0. A possible reason for the remaining difference between the peak positions of the calculated structure factor J (q, ω) and the measured peak spectra can be the experimental uncertainty in the determination of the Fermi wave vector, which can be slightly different from the reported value k F = 8.1 µm −1 . Another possible source of the remaining difference is the similar experimental uncertainty on ln(k F a 2D ). In particular, this uncertainty can be provided by the facts that the Fermi wave vector determined in Ref. [17] is a trap-averaged rather than local quantity. It should be noted that the structure factor J (q, ω) calculated with the mean-field values for µ and ∆ leads to a large discrepancy between the peaks of J (q, ω) and those of the measured spectral function. This confirms the importance of including fluctuations through the GPF approach for the description of the pseudogap state of cold fermions in 2D. In Ref. [17], the pairing crossover temperature T * and the pseudogap pairing temperature T * B < T * have been introduced. The temperature T * coincides with the mean-field transition temperature T * c . The temperature T * B , as stated in Ref. [17], indicates the formation of pairs, and has the same physical meaning as the temperature T p obtained in our study. As far as the transition between the normal and paired states is a crossover rather than a true phase transition, the pairing temperatures T * B and T p only approximately indicate the formation of a paired state. In Fig. 6, the pseudogap pairing temperature T p is compared with the experimental data for T * B . The dotted curve shows the scaled mean-field transition temperature from Ref. [17]; the scaling indicates that the experimental result for T * B is a factor 0.36 smaller than the mean-field prediction. We see that, in contrast to the mean-field result, the value of T p obtained in the present treatment lies in the same range as the experimentally determined temperature T * B . This coincidence is worth remarking. However, the conclusions from the latter comparison of two temperatures need care, because the temperature in [17] is measured in the weakly interacting regime and hence it may differ from the actual temperature in the strongly interacting regime.
It is stated in Ref. [17] that the discrepancy between the mean-field and experimental pairing temperatures could suggest that the appearance of a back-bending feature in the spectral function [20], which has been interpreted as a signature for many-body pairing, is only a qualitative evidence. However, the present results show that the fluctuations can drastically reduce the pairing temperature T p with respect to T * c . Thus there is no discrepancy between experiment and theory when taking account the fluctuations.

VI. CONCLUSIONS
The T -matrix approach straightforwardly applied to cold fermions in two dimensions leads to a divergent fermion density for any finite temperature. We have shown in the present work that taking into account the variation of the order parameter in the number equations, as suggested in the GPF approach [13][14][15], provides a divergence-free description of the paired state in two dimensions. This was shown both through numerical calculations and through an analytic expansion at long wavelengths and low energies, where the divergency occurs in the standard Nozières & Schmitt-Rink approach. The formalism allows to study the effects of the fluctuations both at zero and at finite temperatures, and we find that fluctuations affect the critical binding energy to obtain pairing and superfluidity in the presence of imbalance. Moreover, the formalism also gives access to the density of states of the pair fluctuations, from which we have defined a pseudogap temperature T p as the temperature where an instability appears in the pair fluctuation density. The pseudogap temperature defined in this way agrees with the measured values of the pseudogap temperature in 2D fermi gases. Also the location of the peaks in the spectral functions for pair fluctuations is shown to agree with the experimental observations. The pseudogap temperature T p , along with the critical temperature T BKT for superfluidity and the pair formation temperature T * c , have been calculated as a function of binding energy, temperature and imbalance, from which we obtain the phase diagram as shown in Figs. 3 and 4. Whereas in meanfield the pseudogap temperature is usually identified with the pair formation temperature, we find that the inclusion of fluctuations beyond mean field strongly suppresses the pseudogap temperature with respect to the mean-field pair formation temperature. Moreover, in the presence of imbalance, the pseudogap temperature may cross the BKT temperature for superfluidity. The results obtained here in the context of superfluid quantum gases shed new light on the study of the pseudogap phase in layered high-temperature superconductors, where the question of the crossing of the pseudogap temperature with the superconducting temperature, and the presence of preformed pairs, remains an open question.