## Application of renormalized RPA to polarized Fermi gases

David Durel ∗ and Michael Urban † Paris-Sud 11,

Institut de Physique Nucl´aire, e CNRS/IN2P3, Univ. Universit´ e Paris-Saclay, F-91406 Orsay Cedex, France

We consider a spin imbalanced Fermi gas at zero temperature in the normal phase on the BCS side of the BCS-BEC crossover and around unitarity. We compute the critical polarization for pairing, the correlated occupation numbers and the contact in an extension of particle-particle RPA (also called non self-consistent T -matrix approach or ladder approximation). The so-called renormalized RPA consists in computing the T matrix with self-consistently determined occupation numbers. The occupation numbers are determined either by keeping the self-energy only to first order or by resumming the Dyson equation. In this way, the result for the critical polarization, strongly overestimated in standard RPA, is clearly improved. We also discuss some problems of this approach.

## I. INTRODUCTION

Cold atoms have allowed for the first realization of the crossover from superfluidity of Cooper pairs as described by the Bardeen-Cooper-Schrieffer (BCS) theory to BoseEinstein condensation (BEC) of dimers. This is possible because the scattering length a , i.e., the interaction strength, can be changed thanks to the Feshbach resonance. At zero temperature, a Fermi gas with two spin states σ = ↑ ↓ , of equal masses and populations is superfluid at any interaction strength. This is not the case for a polarized gas, in which pairing disappears beyond some critical polarization P c whose value depends on the interaction strength. In the polarized case, other forms of pairing may exist such as the Fulde-Ferrell-LarkinOvchinnikov (FFLO) phase [1, 2] characterized by an oscillating order parameter corresponding to pairs with a non-zero total momentum. So far, in experiments on polarized Fermi gases performed on resonance ( a → ∞ , the so-called unitary limit) [3-5] and in the BCS-BEC crossover [6, 7], the FFLO phase was not seen, but a phase separation into a paired and an unpaired state was found. In fact, the FFLO phase may be difficult to see in a harmonic trap and the use of flat traps may clarify this question in the future. The problem of pairing in asymmetric systems has also been studied, e.g., in the case of proton-neutron pairing in asymmetric nuclear matter [8] and in the context of QCD with color superconductivity in the case of quark matter which also involves particles of different masses [9].

Many theoretical studies have already addressed the problem of pairing in polarized Fermi gases, for reviews see [10-12]. One problem of polarized systems is that the standard Nozi`res and Schmitt-Rink (NSR) [13] ape proach, which describes the BCS-BEC crossover at finite temperature, fails in the polarized case [14, 15]. Different variants of the NSR approach were proposed to solve

∗ Electronic address: dureldavid@ipno.in2p3.fr

† Electronic address: urban@ipno.in2p3.fr

this problem [16-19]. All these approaches have in common to be based on the ladder approximation for the in-medium T matrix and to include the self-energy in a somewhat more self-consistent way than it is done in the NSR theory.

In this paper, we consider only the case of zero temperature and the normal phase. This implies that the polarization of the gas must exceed the critical polarization P c below which superfluidity sets in. The ladder approximation becomes then equivalent to what is known in nuclear physics as the particle-particle Random-Phase Approximation (pp-RPA) which was applied to polarized Fermi gases in [20]. The pp-RPA gives satisfactory results on the BCS side for not too strong interaction, but at unitarity it overestimates strongly the critical polarization. Within this formalism, the onset of superfluidity appears as an instability, but it is not possible to describe the superfluid phase. Maybe one could generalize the formalism to the superfluid phase using Gor'kov Green's functions, similarly to what was done in the unpolarized case at finite temperature, e.g., in [21], but this is beyond the scope of the present study.

From a general perspective, the RPA describes correlations in the medium, e.g., the correlation energy or correlated occupation numbers, starting from an uncorrelated ground state. In this sense, the RPA is not fully consistent. Different extensions of RPA were developed to take into account these correlations in a more consistent way. One of them is the self-consistent RPA (SCRPA) [22-24], which is based on the correlated RPA ground state. In practice, however, the SCRPA is very difficult to implement except in simple toy models. An approximation to the SCRPA is the renormalized RPA (r-RPA) [25, 26], which instead of using the correlated ground state uses only the correlated occupation numbers in the calculation. The inclusion of the correlated occupation numbers can be expected to solve, at least partially, the problem of the critical polarization in the polarized Fermi gas at zero temperature.

Our article is organized as follows. In Sec. II, we recall briefly the pp-RPA formalism and describe the basic idea of the self-consistent processing of the occupation numbers. We use two different methods to obtain the

occupation numbers, one based on the Dyson equation truncated at first order that we call RPA(1st), and the other using the full Dyson equation that we call RPA( ∞ ). We will discuss the results obtained for the occupation numbers and the contact. In Sec. III, we present the calculation of the critical polarization, first in the case of RPA and then in the case of the r-RPA. Finally, we conclude in Sec. IV.

## II. SELF-CONSISTENT OCCUPATION NUMBERS

## A. Recapitulation of pp-RPA at zero temperature

Let's start by recalling the pp-RPA for a zerotemperature polarized fermion gas. At zero temperature, it is common to fix the densities of up and down populations, denoted ρ ↑ and ρ ↓ , unlike at finite temperature where one usually fixes the chemical potentials. The polarization P of the gas is then defined as

$$P = \frac { \rho _ { \uparrow } - \rho _ { \downarrow } } { \rho _ { \uparrow } + \rho _ { \downarrow } } \,, \quad \quad \quad ( 1 ) \quad \begin{matrix} \text{ro} \\ \text{in} \end{matrix}$$

assuming ρ ↑ &gt; ρ ↓ as convention. Since the formalism breaks down at the superfluid phase transition (see Sec. III), we limit ourselves to sufficiently strong polarizations where the system remains normal fluid even at zero temperature.

The pp-RPA is based on the formalism of the inmedium T matrix which is written as the sum of ladder diagrams that we can see in Fig. 1 (a).

In the case of cold atoms, it is appropriate to take a contact interaction with coupling constant g . The vertex function Γ is written as

$$\Gamma ( { \mathbf k }, \omega ) = \frac { 1 } { 1 / g - J ( { \mathbf k }, \omega ) } \,, \quad \quad ( 2 ) \quad \text{apf} \\. \quad \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

with J ( k , ω ) = J hh ( k , ω ) + J pp ( k , ω ), including both particle-particle (pp) and hole-hole (hh) propagation. The function J needs to be regularized [27] which gives the following expressions:

$$g = \frac { 4 \pi a } { m } \,, \quad \ \ \ \ ( 3 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$J _ { h h } ( { \mathbf k }, \omega ) & = - \int \frac { d ^ { 3 } { \mathbf p } } { ( 2 \pi ) ^ { 3 } } \frac { n ^ { \uparrow } ( | \frac { k } { 2 } + { \mathbf p } | ) \, n ^ { \downarrow } ( | \frac { k } { 2 } - { \mathbf p } | ) } { \omega - \varepsilon _ { \frac { k } { 2 } + { \mathbf p } } - \varepsilon _ { \frac { k } { 2 } - { \mathbf p } } - \mathbf i \, \eta } \,, \quad ( 4 ) \quad \ n ^ { \sigma } ( \quad \\ & \quad \lambda \quad / = \uparrow / \mathbf k \quad - \vdots = \downarrow / \mathbf k \quad - \vdots \quad \ \, \cdots$$

$$J _ { p p } ( k, \omega ) = \int \frac { \mathrm d ^ { 3 } \mathbf p } { ( 2 \pi ) ^ { 3 } } \left ( \frac { \bar { n } ^ { \uparrow } ( | \frac { k } { 2 } + \mathbf p | ) \, \bar { n } ^ { \downarrow } ( | \frac { k } { 2 } - \mathbf p | ) } { \omega - \varepsilon _ { \frac { k } { 2 } + \mathbf p } - \varepsilon _ { \frac { k } { 2 } - \mathbf p } + \mathbf i \, \eta } } + \frac { m } { p ^ { 2 } } \right ) \,,$$

where ε k = k / 2 (2 m ), n σ are the occupation numbers of spin σ , and ¯ n σ = 1 -n σ . In the case of standard RPA, the expression of n σ is n σ ( k ) = θ k ( σ F -k ).

To calculate the correlated occupation numbers n σ in standard RPA, one uses the Dyson equation truncated

FIG. 1: Representation in terms of Feynman diagrams of (a) the vertex function Γ, (b) the self-energy Σ, and the occupation numbers n σ calculated with (c) the truncated Dyson equation [RPA(1st)] and (d) the full Dyson equation [RPA( ∞ )]. Inserting self-consistent occupation numbers into the two-particle propagator J within the r-RPA corresponds roughly to an approximate way of dressing the thick blue lines in diagram (a), but there is no one-to-one correspondence in terms of diagrams.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-f7131c6e-a4e7-4aa5-8620-b09f28a6be85.jpg)

at first order [Fig. 1 (c)], i.e.,

$$G ^ { \sigma } = G _ { 0 } ^ { \sigma } + G _ { 0 } ^ { \sigma } \Sigma ^ { \sigma } G _ { 0 } ^ { \sigma } \,,$$

where G σ is the dressed Green's function, G σ 0 the bare Green's function and Σ σ the self-energy represented in Fig. 1 (b), which has the form

$$\text{take a} \quad & \Sigma ^ { \sigma } ( { \mathbf k }, \omega ) = - \, \mathfrak { i } \int \frac { d ^ { 3 } { \mathbf p } } { ( 2 \pi ) ^ { 3 } } \int \frac { d \omega ^ { \prime } } { 2 \pi } G _ { 0 } ^ { \bar { \sigma } } ( { \mathbf p }, \omega ^ { \prime } ) \, \Gamma ( { \mathbf k } + { \mathbf p }, \omega + \omega ^ { \prime } ) \,, \\ & \text{--} \text{bounds} = \text{d-most} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \max \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfra { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { w } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + \, \mathfrak { h } \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min \, \text{and} \, + w \, \min.$$

where ¯ denotes the spin opposite to σ σ . This is the same approximation that is used at finite temperature in the NSR theory [13]. The general expression for the occupation numbers reads [28]

$$n ^ { \sigma } ( k ) = - \, \mathfrak { i } \int \frac { \mathrm d \omega } { 2 \pi } \, \mathrm e ^ { \mathrm i \omega \eta } \, G ^ { \sigma } ( \mathbf k, \omega ) \,. \quad \quad ( 8 )$$

The resulting expressions for the occupation numbers of the holes and particles are respectively [20]

$$& \L ^ { - \nu } _ { \nu } \\ & ( 4 ) & \quad n ^ { \sigma } ( k < k ^ { \sigma } _ { \nu } ) = 1 + \int \frac { d ^ { 3 } { \mathbf p } } { ( 2 \pi ) ^ { 3 } } \int _ { \Omega _ { \nu } } ^ { + \infty } \frac { d \omega } { \pi } \, \frac { \theta ( k ^ { \bar { \sigma } } _ { \nu } - | { \mathbf p } - { \mathbf k } | ) } { ( \omega - \varepsilon _ { \mathbf k } - \varepsilon _ { \mathbf p } - { \mathbf k } ) ^ { 2 } } \\ & \quad \times \mathcal { I } m \Gamma ( { \mathbf p }, \, \omega ) \,, \quad ( 9 ) \\ & \Big. \,,$$

$$& \int _ { ( 5 ) } ^ { \,, } \quad n ^ { \sigma } ( k > k _ { _ { F } } ^ { \sigma } ) = - \int \frac { \mathrm d ^ { 3 } \mathbf p } { ( 2 \pi ) ^ { 3 } } \, \int _ { - \infty } ^ { \Omega _ { F } } \frac { \mathrm d \omega } { \pi } \, \frac { \theta ( | \mathbf p - \mathbf k | - k _ { _ { F } } ^ { \bar { \sigma } } ) } { ( \omega - \varepsilon _ { \mathbf k } - \varepsilon _ { \mathbf p - \mathbf k } ) ^ { 2 } } \\ & \quad \times \mathcal { I } \mathrm m \Gamma ( \mathbf p, \omega ) \,, \quad ( 1 0 ) \\ & \quad \text{$\mathbb{m}$} \, \circ \, - \, / \mathcal { I } \uparrow 2 \, + \, \mathcal { I } \downarrow 2 \, / / \mathcal { I } \downarrow \, \colon \, \colon \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \colon \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$} \, \text{$\mathbb{m}$,}$$

where Ω F = ( k ↑ 2 F + k ↓ 2 F ) / (2 m ) is the energy separating the two-particle and two-hole continua. The occupation numbers that we obtain can be decomposed

into a continuous part n c and a step of height Z σ = lim ε → 0 n σ ( k σ F -ε ) -n σ ( k σ F + ), i.e., ε

$$n ^ { \sigma } ( k ) = Z ^ { \sigma } \, \theta ( k _ { F } ^ { \sigma } - k ) + n _ { c } ^ { \sigma } ( k ) \,. \quad \ \ ( 1 1 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

On the one hand, the procedure of keeping in Eq. (6) only the first order in the self-energy makes it possible to satisfy the Luttinger theorem, i.e., that the correlations do not modify the densities of the gas. On the other hand, it leads to a pathology when the interaction becomes too strong, the height of the step Z becoming negative (and occupation numbers can become negative, too) [20].

## B. Renormalized pp-RPA

The idea of renormalized RPA (r-RPA) is to reuse the pp-RPA formalism presented in the preceding subsection but in which the occupation numbers in Eqs. (4) and (5) are no longer Heaviside functions but the correlated occupation numbers writted in Eq. (11). This approximation can be justified using the equation-of-motion method for the two-particle Green's function, see, e.g., [22].

In some sense, this prescription can be viewed as an approximation to the two-particle Green's function J one would obtain by dressing the single-particle propagators appearing in the ladder diagrams in Fig. 1 (a). However, there are important differences. For instance, in the limit P → 1, the r-RPA reduces to the RPA because the up and down occupation numbers tend to θ k ( ↑ F -k ) and 0, respectively, while the dressed propagator of the down particle (polaron) would remain non-trivial even in this limit.

By including the correlations ( Z σ &lt; 1) in the calculation of the two-particle propagator J , the logarithmic singularity of R e J , responsible for the instability of the normal phase [28], will be reduced. This can be seen in Fig. 2, where we show J ( k , ω ) for small non-vanishing k for better visibility. This softening of the singularity will allow the normal phase to remain stable at lower polarization.

With this formalism, the results of the pp-RPA can be considered as the first iteration of the self-consistent calculation. To carry out the iteration, the correlated occupation numbers are calculated according to Eqs. (9) and (10) and are reinjected into the functions J , Eqs. (4) and (5). This procedure is repeated until convergence is reached. One important thing to notice is that the converged result is independent of the initial occupation numbers. The comparison between the RPA and r-RPA occupation numbers is shown in Fig. 3 for different values of the interaction strength. For weak interactions, i.e., | ak ↑ F | &lt; 1, the r-RPA provides virtually no correction to the RPA. In contrast, for stronger interactions, the r-RPA reduces the value of Z more strongly than the RPA calculation. Therefore the self-consistent treatment does not cure the pathology of RPA that the value of Z becomes negative for interactions that are too strong;

FIG. 2: Imaginary and real part of J as function of ω for k = 0 01 . k ↑ F and k ↓ F = k ↑ F / 2. The blue solid lines represent J with uncorrelated occupation numbers and the red dashdotted lines are obtained with self-consistent occupation numbers for ak ↑ F = -2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-5f09d3e7-2c4a-4e3d-a66a-d138f4840a6a.jpg)

on the contrary, the negative step appears already for weaker interactions (for instance, in Fig. 3(b), the step of n ↓ ( k ) is negative in r-RPA as can be seen in the zoom).

Let us note that it follows from the spectral representation of the two-particle Green function J (see chapter 15.2 of [29]) that in principle the two-hole continuum of I m J hh should be restricted to energies below Ω F and the two-particle continuum of I m J pp to energies above Ω , F as it is the case in standard pp-RPA. In the selfconsistent treatment, i.e., by including the correlations in the calculation of J , we see that the two-particle continuum I m J pp extends into the two-hole continuum and vice versa as shown in Fig. 4. This is a general problem of the r-RPA approach.

## C. Contact

A very interesting property of the occupation numbers is the asymptotic behavior of the momentum distribution tails. In the case of a contact interaction, the asymptotic behavior ( k glyph[greatermuch] k σ F ) follows the power law n σ ( k ) ∼ C/k 4 . The coefficient

$$C = \lim _ { k \to \infty } k ^ { 4 } n ^ { \uparrow } ( k ) = \lim _ { k \to \infty } k ^ { 4 } n ^ { \downarrow } ( k ) \quad \ \ ( 1 2 )$$

is independent of the spin and, in the notation of Tan [30, 31], it is called contact . This relationship has been restated in a field theory context in Ref. [32]. 1 The

1 In Ref. [32] a different notation is used where C denotes Tan's contact integrated over the volume. Here we adopt Tan's nota-

FIG. 3: Up and down occupation numbers for the polarization k ↓ F = k ↑ F / 2 ( P = 0 778) for the RPA (blue solid line) and the . r-RPA (at convergence, red dashed line). (a) ak ↑ F = -2, (b) ak ↑ F = -3. The inset shows a zoom on the unphysical negative step ( Z ↓ &lt; 0) of the down occupation numbers.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-b531b73a-d895-4a49-b758-a7b64cd0b442.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-4f42bd61-782c-4a0f-bf48-8128e27dca46.jpg)

value of the contact C is related to different thermodynamic properties of the Fermi gas [31, 33]. For instance, it determines the dependence of the energy density E/V on the interaction strength as [31]

$$\frac { d ( E / V ) } { d ( - 1 / a ) } = \frac { C } { 4 \pi m } \,. \quad \quad \quad ( 1 3 )$$

Figure 5 shows the dependence of the contact C on the interaction and polarization parameters for the RPA and the r-RPA. Note that the value of the contact is almost identical within the RPA and the r-RPA, we do not know if the small difference between these curves is only due to the numerical precision or not. In the limit of weak interaction (large -1 / a k ( ↑ F )), the contact approaches the perturbative result C = 16 π 2 a 2 ρ ρ ↑ ↓ [34], as can be seen in Fig. 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-232a782a-7bfe-4932-af2e-04526611c9e0.jpg)

FIG. 4: Imaginary part of J as function of ω for k = 0 and k ↓ F = k ↑ F / 2. The blue line is the imaginary part of J for holes and the red dashed line for particles. As can be seen, each curve extends into the energy zone of the other.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-9424f1bf-7bcd-400c-a8aa-842d11676a9e.jpg)

-

FIG. 5: Contact C (a) as a function of the interaction for an asymmetry of k ↓ F = k ↑ F / 2; (b) as a function of the polarization for interactions of ak ↑ F = -1 and ak ↑ F = -2. The blue solid lines are the RPA results and the red dashed lines are for the r-RPA. The black dashed-dotted line represents the perturbative expression for the contact in the weakly interacting limit.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-1e387aec-7704-4526-b7b5-fd7b0d98e76f.jpg)

## D. Occupation numbers calculated with the full Dyson equation

The problem of the negative step that appears in the occupation numbers calculated with the standard RPA when the interaction becomes too strong is not improved by the self-consistent treatment. One way to cure this pathology is to use the complete Dyson equation instead of the truncated version (6) presented in Sec. II A, i.e., to dress the Green function as

$$G = \frac { 1 } { 1 / G _ { 0 } - \Sigma } \,. \quad & \quad \text{th} \\ \text{areme the accumulation numbers calculated}$$

In terms of diagrams, the occupation numbers calculated with the complete Dyson equation are represented in Fig. 1 (d).

When we consider the complete equation (14) to calculate the occupation numbers, the self-energy Σ, Eq. (7), must be explicitly calculated, which is not the case when using the first-order truncated Dyson equation (6) as done before. The expressions for the imaginary part of the self-energy Σ = Σ hh +Σ pp are given by [20]

$$\mathcal { I } m \, \Sigma _ { h h } ^ { \sigma } ( \mathbf k, \, \omega ) = - \int _ { p > k _ { F } ^ { \sigma } } \frac { \mathrm d ^ { 3 } \mathbf p } { ( 2 \pi ) ^ { 3 } } \, \theta ( \Omega _ { F } - \omega - \varepsilon _ { \mathbf p } ) \\ & \times \mathcal { I } m \, \Gamma ( \mathbf k + \mathbf p, \, \omega + \varepsilon _ { \mathbf p } ) \,, \quad ( 1 5 )$$

$$\mathcal { I } m \Sigma _ { p p } ^ { \sigma } ( \mathbf k, \, \omega ) = \int _ { p < k _ { F } ^ { \sigma } } \frac { d ^ { 3 } \mathbf p } { ( 2 \pi ) ^ { 3 } } \, \theta ( \omega + \varepsilon _ { \mathbf p } - \Omega _ { F } ) \, \begin{array} { c } \alpha \, \cos \, \binom { \alpha } { \text{calcul} } \\ \text{tion} \, ( \, \text{num} \, ( \, \text{mat} \, \text{in} \, \text{s} \, \text{dores} \, \text{s} \end{array} )$$

and the corresponding real parts are calculated with a dispersion relation. Then we calculate the spectral function

$$A ^ { \sigma } ( { \mathbf k }, \, \omega ) = - \mathcal { I } m \, \frac { 1 } { \omega - \epsilon _ { \mathbf k } - \Sigma ^ { \sigma } ( { \mathbf k }, \, \omega - U ^ { \sigma } ) } \,. \quad ( 1 7 ) \quad \text{as in} \\ \intertext { w e h o n i n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n dt h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d t h o n d \cdot }$$

We have introduced the quantity U σ = R e Σ σ ( k , ε F F ) to take into account the shift of the Fermi energy caused by the real part of Σ. A useful property of the spectral function is

$$\int _ { \varepsilon _ { F } + U ^ { \sigma } } ^ { + \infty } \frac { d \omega } { \pi } \, A ^ { \sigma } ( { \mathbf k }, \omega ) - \int _ { - \infty } ^ { \varepsilon _ { F } + U ^ { \sigma } } \frac { d \omega } { \pi } \, A ^ { \sigma } ( { \mathbf k }, \omega ) = 1 \,. \, ( 1 8 ) & & \frac { \cdot } { \infty }$$

The occupation numbers are obtained from

$$n ^ { \sigma } ( k ) = \pm \int _ { - \infty } ^ { \varepsilon _ { F } + U ^ { \sigma } } \frac { d \omega } { \pi } \, A ^ { \sigma } ( { \mathbf k }, \omega ) \,, \quad \ \ ( 1 9 )$$

with + for holes and -for particles. To avoid having to integrate the peak present in the spectral function when calculating the occupation numbers of the holes, we use Eq. (18) that normalizes the spectral function and integrate over the complementary interval where there is no peak. Finally this leads to the formulas

$$n ^ { \sigma } ( k < k _ { _ { F } } ^ { \sigma } ) = 1 - \int _ { \varepsilon _ { F } + U ^ { \sigma } } ^ { + \infty } \frac { \mathrm d \omega } { \pi } \, A ^ { \sigma } ( \mathbf k, \omega ) \,, \quad ( 2 0 ) \quad \text{ with }$$

$$n ^ { \sigma } ( k > k _ { _ { F } } ^ { \sigma } ) = - \int _ { - \infty } ^ { \varepsilon _ { F } + U ^ { \sigma } } \frac { \mathrm d \omega } { \pi } \, A ^ { \sigma } ( \mathbf k, \, \omega ) \,. \quad ( 2 1 )$$

By comparing the results of this method, denoted RPA( ∞ ), and those of the standard RPA (RPA(1st)), we see in Fig. 6 that the occupation numbers within the RPA( ∞ ) are less modified than those within RPA(1st). In particular, the negative step disappears even at the strongest interactions. Roughly speaking, since Z 1st glyph[similarequal] 1 + dΣ / d ω and Z ∞ glyph[similarequal] 1 / (1 -dΣ d / ω ) the step heights of the two methods are related by Z 1st glyph[similarequal] 2 -1 /Z ∞ . For weak interactions, both methods give similar results.

FIG. 6: The solid blue lines represent the occupation numbers calculated with RPA(1st), i.e., the truncated Dyson equation (6) and the dashed red lines represent the occupation numbers calculated with RPA( ∞ ), i.e., the complete Dyson equation (14). The polarization is fixed at k ↓ F = k ↑ F / 2, interaction strength is ak ↑ F = -4. The Z ↓ calculated with RPA( ∞ ) decreases much less than with the RPA(1st).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-d935e2af-6919-4be1-ab27-5e69c8b5258f.jpg)

To set up self-consistency, we adopt the same approach as in Sec. II B. Figure 7 shows the difference of the occupation numbers calculated with the r-RPA(1st) and rRPA( ∞ ). Again, in the occupation numbers calculated

FIG. 7: Same as Fig. 6 but using self-consistent occupation numbers (r-RPA). Z ↓ calculated with r-RPA( ∞ ) is positive whereas the one calculated with r-RPA(1st) is negative.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-304e173b-d8ac-4976-92cc-d37cf67db416.jpg)

with the full Dyson equation, the problem of the negative

TABLE I: Relative error ∆ σ rel (violation of the Luttinger theorem) for different calculations with ak ↑ F = -2 5 . and k ↓ F = k ↑ F / 2

|         |   RPA(1st) | r-RPA(1st)   | RPA( ∞ )   | r-RPA( ∞ )   |
|---------|------------|--------------|------------|--------------|
| ∆ ↑ rel |          0 | 0 . 053 %    | 0 . 073 %  | 0 . 013 %    |
| ∆ ↓ rel |          0 | 0 . 45 %     | 2 . 9 %    | 2 . 1 %      |

step (and even negative occupation numbers) present in the r-RPA(1st) has disappeared.

However there is a price to pay. If we define ρ L σ = k σ 3 F / (6 π 2 ), the Luttinger theorem [35] states that the density ρ σ calculated as the integral of the correlated occupation numbers satisfies the relation ρ σ = ρ L σ . As shown in [20] this is exactly fulfilled within RPA(1st), but it is no longer true if the occupation numbers are calculated with the full Dyson equation (RPA( ∞ )). More quantitatively, we define the relative error ∆ σ rel = | ρ L σ -ρ σ | /ρ σ . Table I shows the violation of the Luttinger theorem for one specific example ( ak ↑ F = -2 5, . k ↓ F = k / ↑ F 2). We see that only the RPA(1st) satisfies the Luttinger theorem exactly. However, the r-RPA(1st) violates the Luttinger theorem only very slightly and the error observed can be due to the accumulation of numerical errors. The RPA( ∞ ), on the contrary, clearly violates the Luttinger theorem and therefore the r-RPA( ∞ ), too. At stronger interactions, the violation within RPA( ∞ ) and r-RPA( ∞ ) can be much worse, e.g., near the critical polarization at the unitary limit (see Sec. III C).

## III. CRITICAL POLARIZATION

## A. FFLO transition

The Thouless criterion [36] states that the superfluid transition occurs when a pole appears in the T matrix at ω = Ω , i.e., 1 F /g -J ( k , ω = Ω ) = 0. F As long as the condition 1 /g -J ( k , Ω ) F &lt; 0 is fulfilled for all k , we are in the normal phase as shown in Fig. 8 as the blue dashed curve.

glyph[negationslash]

In the case of a non-polarized gas at finite temperature, approaching T c from above, the instability of the normal phase sets in first at k = 0 . However, in the case of a polarized gas, the difference between the Fermi levels favors the creation of pairs with non-zero total momentum, resulting in the emergence of a new type of superfluidity called FFLO phase [1, 2]. Now we look for the appearance of the pole when approaching P c from above. Considering 1 /g -J ( k , Ω ) as function of F k for different values of P , as shown in Fig. 8, we notice that the pole at k = 0 appears at a higher polarization than the one at k = 0 . The value of k at which the pole appears first,

FIG. 8: Typical behavior of 1 /g -J ( k , Ω ) (computed within F RPA for ak ↑ F = -1 913) as function of . k for five polarizations P . At the three highest polarizations associated with lower three (blue, black and purple) curves, the system is in the normal phase. The green solid line is at the critical polarization P c where the function vanishes for a value of k which corresponds to k FFLO . The red dash-dotted line is obtained for the value of the polarization which gives the pole at k = 0 (BCS superfluidity), but it lies already in the FFLO superfluid region. For the definitions of P L and a M , see text.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-ed236374-3407-459e-93a7-bfe50eeb8e14.jpg)

coming from strong polarization, is denoted k FFLO :

$$\text{At} \quad \ 1 / g - \mathcal { R e } \, J ( \mathbf k _ { F F L O }, \Omega _ { F } ) = \max _ { k } \left ( 1 / g - \mathcal { R e } \, J ( \mathbf k, \Omega _ { F } ) \right ) \,. \\ \text{Do} \quad \ \ \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \ \. \. \.$$

To determine the critical polarization, we must therefore determine the value of P such that

$$1 / g - \mathcal { R e } \, J ( \mathbf k _ { F F L O }, \Omega _ { F } ) = 0 \,. \quad \quad ( 2 3 )$$

At the critical polarization, there is an instability of the system towards a formation of pairs with momentum k FFLO . Therefore k FFLO corresponds more or less to the wave vector of the order-parameter oscillations in the FFLO phase. However, our theory does not tell us whether the paired phase that will be formed corresponds to a Fulde-Ferrel (FF) state with just one wave vector [2], a Larkin-Ovchinnikov (LO) state with spatial modulations of the order parameter but without a varying phase [1], or even more complicated states with a crystal-like structure.

The qualitative picture described above and illustrated in Fig. 8 remains the same in a large range of interaction strengths 1 / ak ( F ). Increasing 1 / ak ( F ), i.e., 1 /g , simply shifts the curves upwards and thereby increases the critical polarization, which in turn leads to a higher value of k FFLO . The dependence of k FFLO computed within RPA as a function of the ratio k /k ↓ F ↑ F (which is directly related to P , small values of k /k ↓ F ↑ F corresponding to large P and vice versa) is displayed in Fig. 9 as the solid blue line.

At some positive value of 1 /g , corresponding to a strong critical polarization, the picture changes. As it can be seen from the double-dashed black curve in Fig. 8,

1

FIG. 9: Value of k FFLO as a function of k ↓ F /k ↑ F within RPA (solid blue line) and within r-RPA( ∞ ) (dashed red line). In RPA, it is independent of the strength of the interaction. For the r-RPA( ∞ ), we have used ak ↑ F for each k ↓ F /k ↑ F such that P = P c . The FFLO phase disappears for k ↓ F /k ↑ F &lt; 0 224 in . the RPA, and for k ↓ F /k ↑ F &lt; . 0 435 in the r-RPA( ∞ ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-3cccde5b-ad2b-4f45-a765-0680563f2d7a.jpg)

glyph[negationslash]

at some strong polarization denoted P L , the local maximum at k = k FFLO has become so flat that its value coincides with another local maximum that has built up at k = 0. Hence, at P = P L , the global maximum of 1 /g -J changes from the one at k = 0 to the one at k = 0, and as a consequence, k FFLO jumps from a finite value to zero, as can be seen in Fig. 9. This corresponds to the Lifshitz point L in the schematic phase diagram of [37].

## B. Implementation of self-consistency

As explained in the preceding subsection, for fixed values of k ↑ F and a , the method to find the critical polarization is to determine the zero of 1 /g -J ( k FFLO , Ω ) as a F function of P , i.e., in practice as a function of k ↓ F . This curve is shown in Fig. 10. The value of k ↓ F corresponding to P c will be denoted k ↓ F,c . For values of k ↓ F lower than k ↓ F,c , we are in the normal phase, otherwise we are in the superfluid phase which cannot be described with our theory.

However, it is impossible to compute this curve for the r-RPA up to the corresponding k ↓ F,c (for given a and k ↑ F ) if one initializes the self-consistent iteration with the uncorrelated occupation numbers (step functions). The reason becomes clear from Fig. 10: with increasing correlations, P c becomes smaller, i.e., k ↓ F,c becomes larger. Therefore, k ↓ F,c for r-RPA lies in the superfluid area of RPA and already the first iteration step cannot be performed. To avoid going through the superfluid zone, we need to find initial conditions that maximize the correlations and at the same time k ↓ F,c so that it decreases with each iteration. Fortunately the result is independent of the initial

FIG. 10: Schematic behavior of the maximum value of the function 1 /g -R e J ( k , Ω ) depending on F k ↓ F for fixed k ↑ F . The blue solid line is obtained from the RPA and the red dashed line from the r-RPA. As can be seen, the critical polarization obtained from r-RPA is in the superfluid phase of the RPA.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-95c02049-3d5a-4330-8294-d1ee69751edb.jpg)

occupation numbers. The occupation numbers for the initialization are constructed from those calculated with the RPA by artificially increasing as much as possible the correlated part δn σ ( k ) = n σ ( k ) -θ k ( σ F -k ) such that Z ↓ = 0 (cf. Fig. 11).

FIG. 11: Up and down occupation numbers for the polarization k ↓ F = k ↑ F / 2 ( P = 0 778) and . ak ↑ F = -2. Solid blue lines represent the occupation numbers used for the initialization of the self-consistent iteration. Dashed red lines represent the converged occupation numbers within the r-RPA(1st) calculation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-5e6689c8-d1a6-4e83-a118-f6f63353f8c6.jpg)

We then use the procedure described in the preceding section to find the k FFLO . In the r-RPA( ∞ ), we observe a more and more important deviation from the RPA as the polarization increases as can be seen by comparing the solid blue and the dashed red lines in Fig. 9. Also the critical value of k ↓ F changes (where k FFLO disappears). In RPA, we find this point at ( k /k ↓ F ↑ F ) L = 0 224 and slightly . on the BEC side close to unitarity at 1 / a k ( L ↑ F ) glyph[similarequal] 0 210, . in agreement with the disappearence of the FFLO phase

in mean-field theory, see Fig. 6(a) of [34]. In r-RPA( ∞ ), we find ( k /k ↓ F ↑ F ) L = 0 435 at 1 . / a k ( L ↑ F ) glyph[similarequal] 0 131 which is . closer to the unitarity compared to the RPA result. It is due to the fact that the maximum that gives the k FFLO becomes flatter with the self-consistency. Therefore, it drops more quickly below the global maximum at k = 0 and the FFLO phase is lost.

## C. Phase diagram

Following the procedure discussed in the preceding subsections, we are only able to detect the instability of the normal phase corresponding to a second-order phase transition to the superfluid phase. However, if there is a first-order phase transition from the normal to the superfluid phase with a coexistence region extending to a higher polarization than our critical polarization P c , we are not able to see it, because this would require to compute the energy of the superfluid phase. Keeping this word of caution in mind, we show in Fig. 12 the phase diagram giving the critical polarization as a function of the interaction strength.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00530v2-figure-3d6144b1-daed-489d-a660-22cd55a7cce0.jpg)

-

FIG. 12: Critical polarization from the Thouless criterion as function of -1 ( / a k ↑ F ). Thick lines correspond to P c calculated with the occupation numbers and thin lines to P L c that one would obtain if the Luttinger theorem was satisfied. The blue solid lines represent the RPA calculation. The red lines are for r-RPA(1st), in the range where they are dotted it gives a negative step ( Z &lt; 0). The green dashed lines are for r-RPA( ∞ ). The black cross is the theoretical prediction extracted from Fig. 1 of Ref. [38]. We also show experimental results for the limit of the coexistence region corresponding to a first-order phase transition towards an unpolarized superfluid phase; the data are from Refs. [5, 6] (Shin et al.) and [7] (Olsen et al.).

Notice that on the BEC side of the crossover, beyond some interaction -1 / a ( M F k ↑ ) &lt; 0, the polarized normal phase (except at P = 1 where the system is non-interacting) does not exist any more. This point, denoted M, corresponds to the polaron-to-molecule tran- sition discussed in Refs. [39, 40] where few down particles added to a fully polarized system of up particles do not form a normal fluid Fermi sea of polarons any more, but a BEC of molecules with total momentum k = 0. In the present framework, this transition happens when the lowest curve shown in Fig. 8, corresponding to P = 1, is shifted upwards to the horizontal axis, i.e., at

$$\begin{array} { c } \dots \quad \text{$\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\ \quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad$\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $ \quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\quad $\ \end{array} \quad \begin{array} { c } 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ \text{ing} & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \quad & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ \text{ing} & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \quad & & & & \quad & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ \text{of} & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & && & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & - & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & = & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &$$

This value has to be compared with the exact one, 1 / a ( M F k ↑ ) glyph[similarequal] 0 9 . [39]. The origin of this discrepancy is that the ladder approximation includes only 2-particle (2p) but no 3-particle-1-hole (3p1h) and more complicated states [40]. Since at P c = 1 no down particles are present any more, the self-consistency does not have any effect and all our curves, whether they are calculated with RPA or r-RPA, end in the same point M.

Next to the point M, we have the transition towards the k = 0 superfluid phase. Only to the right of the Lifshitz point L mentioned in Sec. III A, visible as kink in the curves pointed by the black arrows in Fig. 12, the transition is of the FFLO type.

We see that, as expected, the critical polarization in r-RPA is always lower than in RPA. Notice that the rRPA(1st) (red lines) breaks down at interactions stronger than 1 / a k ( ↑ F ) glyph[similarequal] -0 675 where a negative step appears in . the occupation numbers. Nevertheless, it is still possible to reach convergence in some range and these results are shown as red dotted lines in Fig. 12. In contrast to the r-RPA(1st), the r-RPA( ∞ ), shown as the green dashed lines, allows us to describe the whole range of interactions up to the point M, since a negative step does not appear. At weak interaction, the computation with r-RPA( ∞ ) gives the same result as r-RPA(1st) which is obvious since the occupation numbers are practically the same in this case.

The critical polarization within RPA is actually the same as within mean field including the possibility of the FFLO phase but no phase separation [41] (see also Fig. 6(a) of [34]). 2 In particular, in the unitary limit, the RPA predicts P c = 0 834. This is clearly higher than . the critical polarization P c = 0 562 obtained . from an energy-density functional (ASLDA) fitted to QuantumMonte-Carlo (QMC) results [38], marked by the cross in Fig. 12. Thus the r-RPA( ∞ ) result of P c = 0 543 . (thick dashed line) would be a significant improvement compared to the RPA. However this low P c is to some extent a consequence of the violation of the Luttinger theorem. In fact, if one computes P c under the assumption that the Luttinger theorem is satisfied, i.e.,

2 In [41] and in Fig. 6(a) of [34], P c is shown as a function of 1 ( / a k F ) = (1 + P ) 1 / 3 / a k ( ↑ F ).

P L c = ( k ↑ 3 F -k ↓ 3 F ) / k ( ↑ 3 F + k ↓ 3 F ), one obtains only a weaker reduction from 0 834 to 0 709 (thin dashed line). . .

Since in RPA(1st) the Luttinger theorem is exactly fulfilled, we have P L c (RPA(1st)) = P c (RPA(1st)) (thin blue line). Since the RPA( ∞ ) uses the same uncorrelated occupation numbers as the RPA(1st), it is clear that also P L c (RPA( ∞ )) = P c (RPA(1st)). The reduction of P c (RPA( ∞ )) (thick blue line) near unitarity is therefore only due to the violation of the Luttinger theorem. The violation of the Luttinger theorem in RPA( ∞ ) and r-RPA( ∞ ) gets stronger and stronger as we approach the unitary limit. It becomes visible in the region of interaction and polarization parameters where the RPA(1st) and r-RPA(1st) give a negative step. On the contrary, for the r-RPA(1st), we have P c (r-RPA(1st)) glyph[similarequal] P L c (r-RPA(1st)) for all values of interaction, therefore we tend to suspect that the Luttinger theorem is fulfilled for the r-RPA(1st) and that the small discrepancy comes from the accumulation of numerical errors. Finally, treating the RPA with the Dyson equation avoids occupation numbers with negative steps but there are still unphysical features.

Finally, let us say a few words about the phase separation observed in [5-7]. These experiments indicate that at the polarizations shown by the circles in Fig. 12, there is a first-order phase transition towards an unpolarized superfluid, in good agreement with QMC results which did not include the possibility of FFLO-type phases [42]. The experimental critical polarization for phase separation is much lower than the one obtained in mean field [41] and leaves some room for a possible continuous (second order) transition towards the FFLO phase before phase separation happens, which corresponds qualitatively to the scenario suggested in [38] for the case of unitarity.

## IV. CONCLUSIONS

We have implemented a self-consistent calculation of the occupation numbers in the formalism of the r-RPA. This is an approximation to the so-called self-consistent RPA (by RPA we mean the RPA in the particle-particle channel, pp-RPA). On the one hand, in the context of polarized Fermi gases, the r-RPA is interesting because it reduces the critical polarization at T = 0 that is strongly overestimated with RPA. The r-RPA is therefore an improvement in the description of this system. On the other hand, the formalism does not cure certain pathologies specific to the RPA. One of the problems of standard RPA [RPA(1st)] is that in the strongly coupled regime the quasiparticle Z factor at the Fermi surface (i.e., the height of the step in the occupation numbers) becomes

[1] P. Fulde and R. A. Ferrell, Physical Review 135 , A550 (1964).

negative. A way to cure this pathology is to consider the complete Dyson equation [RPA( ∞ )] instead of the truncated one which is used in RPA(1st) (and also in the NSR approach at finite temperature). This method does not conserve the number of particles (Luttinger theorem) but it gives at least physical occupation numbers for arbitrary strength of interaction.

Other self-consistent calculations exist to treat inmedium correlations in cold atoms. These numerically demanding methods are based on the self-consistent Green's functions also called Luttinger-Ward formalism at finite temperature. These methods were applied to study the finite temperature BCS-BEC crossover, in particular the phase diagram of non-polarized [43, 44] and polarized [45] gases. In [46], it was shown that within the Luttinger-Ward formalism, the Luttinger theorem is exactly fulfilled with the full Dyson equation.

The r-RPA method allows us to see the evolution of occupation numbers by including correlations in the ground state. An interesting aspect is that the tails of the occupation numbers, whose information is contained in the contact term, remains almost unchanged with the inclusion of self-consistency.

We discussed the transition to the superfluid phase as an instability in the T matrix (Thouless criterion). Except at very strong polarizations (beyond the Lifshitz point), the transition is predicted to be of the FFLO type and our r-RPA result for the corresponding critical polarization in the unitary limit is close to the one of Ref. [38].

Although the FFLO phase has attracted a lot of attention (also in the context of imbalanced gases with two components of different masses [47]), its existence is still an open question. In particular, there is a competition between phase separation (as seen in experiments) and the formation of the FFLO phase [37, 38, 48, 49]. With the present theory, which is limited to the normal phase, we cannot address the question of phase separation because it requires the calculation of the energy in the superfluid phase. The FFLO phase and phase separation can perhaps be reconciled by interpreting the FFLO phase (to be precise, the LO phase) as a periodic 'microphase separation' [10]. Furthermore, as it was pointed out, e.g., in [10], the harmonic potential of the trap makes it difficult to see the FFLO phase in experiments. Hopefully, future experiments with flat traps as they are being built will clarify this question.

## Acknowledgments

We thank Peter Schuck for useful discussions and for comments on the manuscript.

- 47 , 1136 (1964).
- [3] G. B. Partridge, W. Li, R. I. Kamar, Y.-a. Liao, and R. G. Hulet, Nature 311 , 503 (2006).
- [4] Y.-i. Shin, M. W. Zwierlein, C. H. Schunck, A. Schirotzek, and W. Ketterle, Phys. Rev. Lett. 97 , 049901 (2006).
- [5] Y.-i. Shin, C. H. Schunck, A. Schirotzek, and W. Ketterle, Nature 451 , 689 (2008).
- [6] Y.-i. Shin, A. Schirotzek, C. H. Schunck, and W. Ketterle, Phys. Rev. Lett. 101 , 070404 (2008).
- [7] B. A. Olsen, M. C. Revelle, J. A. Fry, D. E. Sheehy, and R. G. Hulet, Phys. Rev. A 92 , 063616 (2015).
- [8] M. Stein, A. Sedrakian, X.-G. Huang, and J. W. Clark, Phys. Rev. C 90 , 065804 (2014).
- [9] M. G. Alford, K. Rajagopal, T. Schaefer, and A. Schmitt, Reviews of Modern Physics 80 , 1455 (2008).
- [10] L. Radzihovsky and D. E. Sheehy, Rep. Prog. Phys. 73 , 076501 (2010).
- [11] F. Chevy and C. Mora, Reports on Progress in Physics 73 , 112401 (2010).
- [12] K. Gubbels and H. Stoof, Physics Reports 525 , 255 (2013), ISSN 0370-1573, imbalanced Fermi Gases at Unitarity.
- [13] P. Nozi`res and S. Schmitt-Rink, Journal of Low Teme perature Physics 59 , 195 (1985).
- [14] X.-J. Liu and H. Hu, Europhys. Lett. 75 , 364 (2006).
- [15] M. M. Parish, A. Marchetti, F. M. andLamacraft, and B. D. Simons, Nat. Phys. 3 , 124 (2007).
- [16] T. Kashimura, R. Watanabe, and Y. Ohashi, Phys. Rev. A 86 , 043622 (2012).
- [17] Q. Chen, Y. He, C.-C. Chien, and K. Levin, Phys. Rev. B 75 , 014521 (2007).
- [18] Y. He, C.-C. Chien, Q. Chen, and K. Levin, Phys. Rev. A 75 , 021602 (2007).
- [19] P.-A. Pantel, D. Davesne, and M. Urban, Phys. Rev. A 90 , 053629 (2014).
- [20] M. Urban and P. Schuck, Phys. Rev. A 90 , 023632 (2014).
- [21] P. Pieri, L. Pisani, and G. C. Strinati, Phys. Rev. B 70 , 094508 (2004).
- [22] J. Dukelsky, G. R¨ opke, and P. Schuck, Nucl. Phys. A 628 , 17 (1998).
- [23] S. Sch¨fer and P. Schuck, Phys. Rev. B a 59 , 1712 (1999).
- [24] J. G. Hirsch, A. Mariano, J. Dukelsky, and P. Schuck, Annals of Physics 296 , 187 (2002).
- [25] F. Catara, G. Piccitto, M. Sambataro, and N. Van Giai, Phys. Rev. B 54 , 17536 (1996).
- [26] D. S. Delion, P. Schuck, and J. Dukelsky, Phys. Rev. C 72 , 064305 (2005).
- [27] A. Perali, P. Pieri, G. C. Strinati, and C. Castellani, Phys. Rev. B 66 , 024510 (2002).
- [28] A. Fetter and J. Walecka, Quantum Theory of Manyparticle Systems (McGraw-Hill, New York, 1975).
- [29] J.-P. Blaizot and G. Ripka, Quantum Theory of Finite Systems (MIT Press, 1986).
- [30] S. Tan, Annals of Physics 323 , 2952 (2008).
- [31] S. Tan, Annals of Physics 323 , 2971 (2008).
- [32] E. Braaten and L. Platter, Phys. Rev. Lett. 100 , 205301 (2008).
- [33] S. Tan, Annals of Physics 323 , 2987 (2008).
- [34]
- G. Calvanese Strinati, P. Pieri, G. R¨pke, P. Schuck, and o M. Urban, Physics Reports 738 , 1 (2018).
- [35] J. M. Luttinger, Physical Review 119 , 1153 (1960).
- [36] D. J. Thouless, Annals of Physics 10 , 553 (1960).
- [37] D. T. Son and M. A. Stephanov, Phys. Rev. A 74 , 013614 (2006).
- [38] A. Bulgac and M. M. Forbes, Phys. Rev. Lett. 101 , 215301 (2008).
- [39] N. V. Prokof'ev and B. V. Svistunov, Phys. Rev. B 77 , 125101 (2008).
- [40] M. Punk, P. T. Dumitrescu, and W. Zwerger, Phys. Rev. A 80 , 053605 (2009).
- [41] C.-H. Pao and S.-K. Yip, Journal of Physics: Conference Series 150 , 032078 (2009).
- [42] S. Pilati and S. Giorgini, Phys. Rev. Lett. 100 , 030401 (2008).
- [43] R. Haussmann, Phys. Rev. B 94 , 12975 (1994).
- [44] R. Haussmann, W. Rantner, S. Cerrito, and W. Zwerger, Phys. Rev. A 75 , 023610 (2007).
- [45] B. Frank, J. Lang, and W. Zwerger, J. Exp. Theor. Phys. 125 , 812 (2018).
- [46] P. Pieri and G. C. Strinati, Eur. Phys. J. B 90 , 68 (2017).
- [47] M. M. Forbes, E. Gubankova, W. V. Liu, and F. Wilczek, Phys. Rev. Lett. 94 , 017001 (2005).
- [48] P. F. Bedaque, H. Caldas, and G. Rupak, Phys. Rev. Lett. 91 , 247002 (2003).
- [49] W. V. Liu and F. Wilczek, Phys. Rev. Lett. 90 , 047002 (2003).
David Durel Michael Urban
Application of renormalized RPA to polarized Fermi gases
We consider a spin imbalanced Fermi gas at zero temperature in the normal
phase on the BCS side of the BCS-BEC crossover and around unitarity. We compute
the critical polarization for pairing, the correlated occupation numbers and
the contact in an extension of particle-particle RPA (also called non
self-consistent \textit{T}-matrix approach or ladder approximation). The
so-called renormalized RPA consists in computing the \textit{T} matrix with
self-consistently determined occupation numbers. The occupation numbers are
determined either by keeping the self-energy only to first order or by
resumming the Dyson equation. In this way, the result for the critical
polarization, strongly overestimated in standard RPA, is clearly improved. We
also discuss some problems of this approach.