Application of renormalized RPA to polarized Fermi gases

TL;DR

This paper extends the particle-particle RPA method with a renormalization scheme to better predict critical polarization and other properties of polarized Fermi gases near unitarity at zero temperature.

Contribution

It introduces a renormalized RPA approach with self-consistent occupation numbers, improving predictions over standard RPA for polarized Fermi gases.

Findings

Critical polarization is more accurately predicted.

Occupation numbers are determined self-consistently.

The approach highlights some limitations of the method.

Abstract

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.

Tables

Table 1. Table 1: Relative error Δ rel σ superscript subscript Δ rel 𝜎 \Delta_{\text{\tiny{rel}}}^{\sigma} (violation of the Luttinger theorem) for different calculations with a ​ k F ↑ = − 2.5 𝑎 superscript subscript 𝑘 F ↑ 2.5 a\,k_{\text{\tiny{F}}}^{\uparrow}=-2.5 and k F ↓ = k F ↑ / 2 superscript subscript 𝑘 F ↓ superscript subscript 𝑘 F ↑ 2 k_{\text{\tiny{F}}}^{\downarrow}=k_{\text{\tiny{F}}}^{\uparrow}/2

	RPA(1st)	r-RPA(1st)	RPA($\infty$)	r-RPA($\infty$)
${\Delta }_{\text{rel}}^{\uparrow }$	$0$	$0.053\%$	$0.073\%$	$0.013\%$
${\Delta }_{\text{rel}}^{\downarrow }$	$0$	$0.45\%$	$2.9\%$	$2.1\%$

Equations

$$P=\frac{{\rho }_{\uparrow }-{\rho }_{\downarrow }}{{\rho }_{\uparrow }+{\rho }_{\downarrow }},$$

$$\Gamma (k,\omega )=\frac{1}{1/g-J(k,\omega )},$$

$$g=\frac{4\pi a}{m},$$

$$J_{hh}(k,\omega )=-\int \frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}\frac{n^{\uparrow }(|\frac{k}{2}+p|)n^{\downarrow }(|\frac{k}{2}-p|)}{\omega -{\epsilon }_{\frac{k}{2}+p}-{\epsilon }_{\frac{k}{2}-p}-i\eta },$$

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

$$G^{\sigma }=G_0^{\sigma }+G_0^{\sigma }{\Sigma }^{\sigma }{G}_0^{\sigma },$$

$${\Sigma }^{\sigma }(k,\omega )=-i\int \frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}\int \frac{\mathrm{d}{\omega }^{\prime }}{2\pi }{G}_0^{\bar{\sigma }}(p,{\omega }^{\prime })\Gamma (k+p,\omega +{\omega }^{\prime }),$$

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

$$\begin{gathered} n^{\sigma }(k<k_{\text{F}}^{\sigma })=1+\int \frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}{\int }_{{\Omega }_{\text{F}}}^{+\infty }\frac{\mathrm{d}\omega }{\pi }\frac{\theta (k_{\text{F}}^{\bar{\sigma }}-|p-k|)}{(\omega -{\epsilon }_k-{\epsilon }_{p-k}{)}^2} \\ \times \mathrm{Im}\Gamma (p,\omega ), \end{gathered}$$

$$\begin{gathered} n^{\sigma }(k>k_{\text{F}}^{\sigma })=-\int \frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}{\int }_{-\infty }^{{\Omega }_{\text{F}}}\frac{\mathrm{d}\omega }{\pi }\frac{\theta (|p-k|-k_{\text{F}}^{\bar{\sigma }})}{(\omega -{\epsilon }_k-{\epsilon }_{\mathbf{p}-k}{)}^2} \\ \times \mathrm{Im}\Gamma (p,\omega ), \end{gathered}$$

$$n^{\sigma }(k)=Z^{\sigma }\theta (k_{\text{F}}^{\sigma }-k)+n_c^{\sigma }(k).$$

$$C=\underset{k\to \infty }{\lim }{k}^4{n}^{\uparrow }(k)=\underset{k\to \infty }{\lim }{k}^4{n}^{\downarrow }(k)$$

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

$$G=\frac{1}{1/G_0-\Sigma }.$$

$$\begin{gathered} \mathrm{Im}{\Sigma }_{hh}^{\sigma }(k,\omega )=-{\int }_{p>k_{\text{F}}^{\bar{\sigma }}}\frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}\theta ({\Omega }_{\text{F}}-\omega -{\epsilon }_p) \\ \times \mathrm{Im}\Gamma (k+p,\omega +{\epsilon }_p), \end{gathered}$$

$$\begin{gathered} \mathrm{Im}{\Sigma }_{pp}^{\sigma }(k,\omega )={\int }_{p<k_{\text{F}}^{\bar{\sigma }}}\frac{{\mathrm{d}}^3\mathbf{p}}{(2\pi {)}^3}\theta (\omega +{\epsilon }_p-{\Omega }_{\text{F}}) \\ \times \mathrm{Im}\Gamma (k+p,\omega +{\epsilon }_p) \end{gathered}$$

$$A^{\sigma }(k,\omega )=-\mathrm{Im}\frac{1}{\omega -{ϵ}_k-{\Sigma }^{\sigma }(k,\omega -U^{\sigma })}.$$

$${\int }_{{\epsilon }_{\text{F}}+U^{\sigma }}^{+\infty }\frac{\mathrm{d}\omega }{\pi }{A}^{\sigma }(k,\omega )-{\int }_{-\infty }^{{\epsilon }_{\text{F}}+U^{\sigma }}\frac{\mathrm{d}\omega }{\pi }{A}^{\sigma }(k,\omega )=1.$$

$$n^{\sigma }(k)=\pm {\int }_{-\infty }^{{\epsilon }_{\text{F}}+U^{\sigma }}\frac{\mathrm{d}\omega }{\pi }{A}^{\sigma }(k,\omega ),$$

$$n^{\sigma }(k<k_{\text{F}}^{\sigma })=1-{\int }_{{\epsilon }_{\text{F}}+U^{\sigma }}^{+\infty }\frac{\mathrm{d}\omega }{\pi }{A}^{\sigma }(k,\omega ),$$

$$n^{\sigma }(k>k_{\text{F}}^{\sigma })=-{\int }_{-\infty }^{{\epsilon }_{\text{F}}+U^{\sigma }}\frac{\mathrm{d}\omega }{\pi }{A}^{\sigma }(k,\omega ).$$

$$1/g-\mathrm{Re}J({\mathbf{k}}_{\text{FFLO}},{\Omega }_{\text{F}})=\underset{k}{\max }\left(1/g-\mathrm{Re}J(k,{\Omega }_{\text{F}})\right).$$

$$1/g-\mathrm{Re}J({\mathbf{k}}_{\text{FFLO}},{\Omega }_{\text{F}})=0.$$

$$\begin{gathered} \frac{1}{a_{\text{M}}{k}_{\text{F}}^{\uparrow }}=\frac{4\pi }{m{k}_{\text{F}}^{\uparrow }}J(0,{\epsilon }_{\text{F}};k_{\text{F}}^{\downarrow }=0)=\frac{2}{\pi }-\frac{1}{\sqrt{2}\pi }\ln \frac{\sqrt{2}+1}{\sqrt{2}-1} \\ \simeq 0.24. \end{gathered}$$

Full text

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 Bose-Einstein 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 $\sigma=\uparrow,\downarrow$ 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-Larkin-Ovchinnikov (FFLO) phase Fulde and Ferrell (1964); Larkin and Ovchinnikov (1964) 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\to\infty$, the so-called unitary limit) Partridge et al. (2006); Shin et al. (2006, 2008a) and in the BCS-BEC crossover Shin et al. (2008b); Olsen et al. (2015), 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 Stein et al. (2014) and in the context of QCD with color superconductivity in the case of quark matter which also involves particles of different masses Alford et al. (2008).

Many theoretical studies have already addressed the problem of pairing in polarized Fermi gases, for reviews see Radzihovsky and Sheehy (2010); Chevy and Mora (2010); Gubbels and Stoof (2013). One problem of polarized systems is that the standard Nozières and Schmitt-Rink (NSR) Nozières and Schmitt-Rink (1985) approach, which describes the BCS-BEC crossover at finite temperature, fails in the polarized case Liu and Hu (2006); Parish et al. (2007). Different variants of the NSR approach were proposed to solve this problem Kashimura et al. (2012); Chen et al. (2007); He et al. (2007); Pantel et al. (2014). 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 Urban and Schuck (2014). 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 Pieri et al. (2004), 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) Dukelsky et al. (1998); Schäfer and Schuck (1999); Hirsch et al. (2002), 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) Catara et al. (1996); Delion et al. (2005), 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($\infty$). 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

II.1 Recapitulation of pp-RPA at zero temperature

Let’s start by recalling the pp-RPA for a zero-temperature polarized fermion gas. At zero temperature, it is common to fix the densities of up and down populations, denoted $\rho_{\uparrow}$ and $\rho_{\downarrow}$, unlike at finite temperature where one usually fixes the chemical potentials. The polarization $P$ of the gas is then defined as

[TABLE]

assuming $\rho_{\uparrow}>\rho_{\downarrow}$ 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 in-medium 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 $\Gamma$ is written as

[TABLE]

with $J(\operatorname{\textbf{\text{k}}},\,\omega)=J_{hh}(\operatorname{\textbf{\text{k}}},\,\omega)+J_{pp}(\operatorname{\textbf{\text{k}}},\,\omega)$, including both particle-particle (pp) and hole-hole (hh) propagation. The function $J$ needs to be regularized Perali et al. (2002) which gives the following expressions:

[TABLE]

where $\varepsilon_{\operatorname{\textbf{\text{k}}}}=k^{2}/(2m)$, $n^{\sigma}$ are the occupation numbers of spin $\sigma$, and $\bar{n}^{\sigma}=1-n^{\sigma}$. In the case of standard RPA, the expression of $n^{\sigma}$ is $n^{\sigma}(k)=\theta(k_{\text{\tiny{F}}}^{\sigma}-k)$.

To calculate the correlated occupation numbers $n^{\sigma}$ in standard RPA, one uses the Dyson equation truncated at first order [Fig. 1 (c)], i.e.,

[TABLE]

where $G^{\sigma}$ is the dressed Green’s function, $G_{0}^{\sigma}$ the bare Green’s function and $\Sigma^{\sigma}$ the self-energy represented in Fig. 1 (b), which has the form

[TABLE]

where $\bar{\sigma}$ denotes the spin opposite to $\sigma$. This is the same approximation that is used at finite temperature in the NSR theory Nozières and Schmitt-Rink (1985). The general expression for the occupation numbers reads Fetter and Walecka (1975)

[TABLE]

The resulting expressions for the occupation numbers of the holes and particles are respectively Urban and Schuck (2014)

[TABLE]

[TABLE]

where $\Omega_{\text{\tiny{F}}}=(k_{\text{\tiny{F}}}^{\uparrow\,2}+k_{\text{\tiny{F}}}^{\downarrow\,2})/(2m)$ 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^{\sigma}=\lim\limits_{\varepsilon\rightarrow 0}n^{\sigma}(k_{\text{\tiny{F}}}^{\sigma}-\varepsilon)-n^{\sigma}(k_{\text{\tiny{F}}}^{\sigma}+\varepsilon)$, i.e.,

[TABLE]

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) Urban and Schuck (2014).

II.2 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., Dukelsky et al. (1998).

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\to 1$, the r-RPA reduces to the RPA because the up and down occupation numbers tend to $\theta(k_{\text{\tiny{F}}}^{\uparrow}-k)$ and [math], respectively, while the dressed propagator of the down particle (polaron) would remain non-trivial even in this limit.

By including the correlations ($Z^{\sigma}<1$) in the calculation of the two-particle propagator $J$, the logarithmic singularity of $\operatorname{\mathcal{R}e}J$, responsible for the instability of the normal phase Fetter and Walecka (1975), will be reduced. This can be seen in Fig. 2, where we show $J(\operatorname{\textbf{\text{k}}},\omega)$ for small non-vanishing $\operatorname{\textbf{\text{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., $|a\,k_{\text{\tiny{F}}}^{\uparrow}|<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; on the contrary, the negative step appears already for weaker interactions (for instance, in Fig. 3(b), the step of $n^{\downarrow}(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 Blaizot and Ripka (1986)) that in principle the two-hole continuum of $\operatorname{\mathcal{I}m}J_{hh}$ should be restricted to energies below $\Omega_{\text{\tiny{F}}}$ and the two-particle continuum of $\operatorname{\mathcal{I}m}J_{pp}$ to energies above $\Omega_{\text{\tiny{F}}}$, as it is the case in standard pp-RPA. In the self-consistent treatment, i.e., by including the correlations in the calculation of $J$, we see that the two-particle continuum $\operatorname{\mathcal{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.

II.3 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\gg k_{\text{\tiny{F}}}^{\sigma}$) follows the power law $n^{\sigma}(k)\sim C/k^{4}$. The coefficient

[TABLE]

is independent of the spin and, in the notation of Tan Tan (2008a, b), it is called contact. This relationship has been restated in a field theory context in Ref. Braaten and Platter (2008).111In Ref. Braaten and Platter (2008) a different notation is used where $C$ denotes Tan’s contact integrated over the volume. Here we adopt Tan’s notation, which is more convenient for uniform systems. The value of the contact $C$ is related to different thermodynamic properties of the Fermi gas Tan (2008b, c). For instance, it determines the dependence of the energy density $E/V$ on the interaction strength as Tan (2008b)

[TABLE]

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_{\text{\tiny{F}}}^{\uparrow}))$, the contact approaches the perturbative result $C=16\pi^{2}\,a^{2}\,\rho^{\uparrow}\rho^{\downarrow}$ Calvanese Strinati et al. (2018), as can be seen in Fig. 5.

II.4 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.1, i.e., to dress the Green function as

[TABLE]

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 $\Sigma$, 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 $\Sigma=\Sigma_{hh}+\Sigma_{pp}$ are given by Urban and Schuck (2014)

[TABLE]

[TABLE]

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

[TABLE]

We have introduced the quantity $U^{\sigma}=\operatorname{\mathcal{R}e}\Sigma^{\sigma}(k_{\text{\tiny{F}}},\,\varepsilon_{\text{\tiny{F}}})$ to take into account the shift of the Fermi energy caused by the real part of $\Sigma$. A useful property of the spectral function is

[TABLE]

The occupation numbers are obtained from

[TABLE]

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

[TABLE]

[TABLE]

By comparing the results of this method, denoted RPA($\infty$), and those of the standard RPA (RPA(1st)), we see in Fig. 6 that the occupation numbers within the RPA($\infty$) are less modified than those within RPA(1st). In particular, the negative step disappears even at the strongest interactions. Roughly speaking, since $Z_{1\text{st}}\simeq 1+\mathrm{d}\Sigma/\mathrm{d}\omega$ and $Z_{\infty}\simeq 1/(1-\mathrm{d}\Sigma/\mathrm{d}\omega)$ the step heights of the two methods are related by $Z_{1\text{st}}\simeq 2-1/Z_{\infty}$. For weak interactions, both methods give similar results.

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

Again, in the occupation numbers calculated with the full Dyson equation, the problem of the negative step (and even negative occupation numbers) present in the r-RPA(1st) has disappeared.

However there is a price to pay. If we define $\rho_{\sigma}^{\text{\tiny{L}}}=k_{\text{\tiny{F}}}^{\sigma\,3}/(6\pi^{2})$, the Luttinger theorem Luttinger (1960) states that the density $\rho_{\sigma}$ calculated as the integral of the correlated occupation numbers satisfies the relation $\rho_{\sigma}=\rho_{\sigma}^{\text{\tiny{L}}}$. As shown in Urban and Schuck (2014) 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($\infty$)). More quantitatively, we define the relative error $\Delta_{\text{\tiny{rel}}}^{\sigma}=|\rho_{\sigma}^{\text{\tiny{L}}}-\rho_{\sigma}|/\rho_{\sigma}$. Table 1 shows the violation of the Luttinger theorem for one specific example ($a\,k_{\text{\tiny{F}}}^{\uparrow}=-2.5$, $k_{\text{\tiny{F}}}^{\downarrow}=k_{\text{\tiny{F}}}^{\uparrow}/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($\infty$), on the contrary, clearly violates the Luttinger theorem and therefore the r-RPA($\infty$), too. At stronger interactions, the violation within RPA($\infty$) and r-RPA($\infty$) can be much worse, e.g., near the critical polarization at the unitary limit (see Sec. III.3).

III Critical polarization

III.1 FFLO transition

The Thouless criterion Thouless (1960) states that the superfluid transition occurs when a pole appears in the T matrix at $\omega=\Omega_{\text{\tiny{F}}}$, i.e., $1/g-J(\operatorname{\textbf{\text{k}}},\omega=\Omega_{\text{\tiny{F}}})=0$. As long as the condition $1/g-J(\operatorname{\textbf{\text{k}}},\Omega_{\text{\tiny{F}}})<0$ is fulfilled for all $\operatorname{\textbf{\text{k}}}$, we are in the normal phase as shown in Fig. 8 as the blue dashed curve.

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 $\operatorname{\textbf{\text{k}}}=\boldsymbol{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 Fulde and Ferrell (1964); Larkin and Ovchinnikov (1964). Now we look for the appearance of the pole when approaching $P_{c}$ from above. Considering $1/g-J(\operatorname{\textbf{\text{k}}},\Omega_{\text{\tiny{F}}})$ as function of $k$ for different values of $P$, as shown in Fig. 8, we notice that the pole at $\operatorname{\textbf{\text{k}}}\neq\boldsymbol{0}$ appears at a higher polarization than the one at $\operatorname{\textbf{\text{k}}}=\boldsymbol{0}$. The value of $k$ at which the pole appears first, coming from strong polarization, is denoted $k_{\text{\tiny{FFLO}}}$:

[TABLE]

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

[TABLE]

At the critical polarization, there is an instability of the system towards a formation of pairs with momentum $k_{\text{\tiny{FFLO}}}$. Therefore $k_{\text{\tiny{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 Larkin and Ovchinnikov (1964), a Larkin-Ovchinnikov (LO) state with spatial modulations of the order parameter but without a varying phase Fulde and Ferrell (1964), 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_{\text{\tiny{F}}})$. Increasing $1/(ak_{\text{\tiny{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_{\text{\tiny{FFLO}}}$. The dependence of $k_{\text{\tiny{FFLO}}}$ computed within RPA as a function of the ratio $k_{\text{\tiny{F}}}^{\downarrow}/k_{\text{\tiny{F}}}^{\uparrow}$ (which is directly related to $P$, small values of $k_{\text{\tiny{F}}}^{\downarrow}/k_{\text{\tiny{F}}}^{\uparrow}$ 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, at some strong polarization denoted $P_{\text{\tiny{L}}}$, the local maximum at $k=k_{\text{\tiny{FFLO}}}$ has become so flat that its value coincides with another local maximum that has built up at $k=0$. Hence, at $P=P_{\text{\tiny{L}}}$, the global maximum of $1/g-J$ changes from the one at $k\neq 0$ to the one at $k=0$, and as a consequence, $k_{\text{\tiny{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 Son and Stephanov (2006).

III.2 Implementation of self-consistency

As explained in the preceding subsection, for fixed values of $k_{\text{\tiny{F}}}^{\uparrow}$ and $a$, the method to find the critical polarization is to determine the zero of $1/g-J(\operatorname{\textbf{\text{k}}}_{\text{\tiny{FFLO}}},\Omega_{\text{\tiny{F}}})$ as a function of $P$, i.e., in practice as a function of $k_{\text{\tiny{F}}}^{\downarrow}$. This curve is shown in Fig. 10.

The value of $k_{\text{\tiny{F}}}^{\downarrow}$ corresponding to $P_{c}$ will be denoted $k_{\text{\tiny{F,c}}}^{\downarrow}$. For values of $k_{\text{\tiny{F}}}^{\downarrow}$ lower than $k_{\text{\tiny{F,c}}}^{\downarrow}$, 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_{\text{\tiny{F,c}}}^{\downarrow}$ (for given $a$ and $k_{\text{\tiny{F}}}^{\uparrow}$) 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_{\text{\tiny{F,c}}}^{\downarrow}$ becomes larger. Therefore, $k_{\text{\tiny{F,c}}}^{\downarrow}$ 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_{\text{\tiny F,c}}^{\downarrow}$ so that it decreases with each iteration. Fortunately the result is independent of the initial 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 $\delta n^{\sigma}(k)=n^{\sigma}(k)-\theta(k_{\text{\tiny{F}}}^{\sigma}-k)$ such that $Z^{\downarrow}=0$ (cf. Fig. 11).

We then use the procedure described in the preceding section to find the $k_{\text{\tiny{FFLO}}}$. In the r-RPA($\infty$), 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_{\text{\tiny{F}}}^{\downarrow}$ changes (where $k_{\text{\tiny{FFLO}}}$ disappears). In RPA, we find this point at $(k_{\text{\tiny{F}}}^{\downarrow}/k_{\text{\tiny{F}}}^{\uparrow})_{\text{\tiny{L}}}=0.224$ and slightly on the BEC side close to unitarity at $1/(a_{\text{\tiny{L}}}k_{\text{\tiny{F}}}^{\uparrow})\simeq 0.210$, in agreement with the disappearence of the FFLO phase in mean-field theory, see Fig. 6(a) of Calvanese Strinati et al. (2018). In r-RPA($\infty$), we find $(k_{\text{\tiny{F}}}^{\downarrow}/k_{\text{\tiny{F}}}^{\uparrow})_{\text{\tiny{L}}}=0.435$ at $1/(a_{\text{\tiny{L}}}k_{\text{\tiny{F}}}^{\uparrow})\simeq 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_{\text{\tiny{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.

III.3 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.

Notice that on the BEC side of the crossover, beyond some interaction $-1/(a_{\text{\tiny{M}}}k_{\text{\tiny{F}}}^{\uparrow})<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 transition discussed in Refs. Prokof’ev and Svistunov (2008); Punk et al. (2009) 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

[TABLE]

This value has to be compared with the exact one, $1/(a_{\text{\tiny{M}}}k_{\text{\tiny{F}}}^{\uparrow})\simeq 0.9$ Prokof’ev and Svistunov (2008). 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 Punk et al. (2009). 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.1, 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 r-RPA(1st) (red lines) breaks down at interactions stronger than $1/(a\,\smash{k_{\text{\tiny{F}}}^{\uparrow}})\simeq-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($\infty$), 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($\infty$) 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 Pao and Yip (2009) (see also Fig. 6(a) of Calvanese Strinati et al. (2018)).222In Pao and Yip (2009) and in Fig. 6(a) of Calvanese Strinati et al. (2018), $P_{c}$ is shown as a function of $1/(a\,k_{\text{\tiny{F}}})=(1+P)^{1/3}/(a\,k_{\text{\tiny{F}}}^{\uparrow})$. 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 Quantum-Monte-Carlo (QMC) results Bulgac and Forbes (2008), marked by the cross in Fig. 12. Thus the r-RPA($\infty$) 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., $P_{c}^{\text{\tiny{L}}}=(k_{\text{\tiny{F}}}^{\uparrow 3}-k_{\text{\tiny{F}}}^{\downarrow 3})/(k_{\text{\tiny{F}}}^{\uparrow 3}+k_{\text{\tiny{F}}}^{\downarrow 3})$, 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_{c}^{\text{\tiny{L}}}(\text{RPA(1st)})=P_{c}(\text{RPA(1st)})$ (thin blue line). Since the RPA($\infty$) uses the same uncorrelated occupation numbers as the RPA(1st), it is clear that also P_{c}^{\text{\tiny{L}}}(\text{RPA(\infty)})=P_{c}(\text{RPA(1st)}). The reduction of P_{c}(\text{RPA(\infty)}) (thick blue line) near unitarity is therefore only due to the violation of the Luttinger theorem. The violation of the Luttinger theorem in RPA($\infty$) and r-RPA($\infty$) 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}(\text{r-RPA(1st)})\simeq P_{c}^{\text{\tiny{L}}}(\text{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 Shin et al. (2008a, b); Olsen et al. (2015). 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 Pilati and Giorgini (2008). The experimental critical polarization for phase separation is much lower than the one obtained in mean field Pao and Yip (2009) 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 Bulgac and Forbes (2008) 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 negative. A way to cure this pathology is to consider the complete Dyson equation [RPA($\infty$)] 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 in-medium 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 Haussmann (1994); Haussmann et al. (2007) and polarized Frank et al. (2018) gases. In Pieri and Strinati (2017), 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 LABEL:Bulgac2008.

Although the FFLO phase has attracted a lot of attention (also in the context of imbalanced gases with two components of different masses Forbes et al. (2005)), 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 Bulgac and Forbes (2008); Bedaque et al. (2003); Liu and Wilczek (2003); Son and Stephanov (2006). 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 “micro-phase separation” Radzihovsky and Sheehy (2010). Furthermore, as it was pointed out, e.g., in Radzihovsky and Sheehy (2010), 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.