Fermi gas throughout the BCS-BEC crossover: Comparative study of t-matrix approaches with various degrees of self-consistency

TL;DR

This study systematically compares different levels of self-consistency in t-matrix approximations to understand their impact on thermodynamic and dynamical properties of a Fermi gas across the BCS-BEC crossover.

Contribution

It provides a comprehensive, unbiased comparison of all variants of the t-matrix approximation with varying self-consistency degrees in an attractive Fermi gas.

Findings

Different self-consistency levels significantly affect thermodynamic predictions.

The study clarifies which t-matrix variants are more appropriate for describing the normal phase.

Results inform the choice of theoretical approaches for Fermi gases in the BCS-BEC crossover.

Abstract

The diagrammatic t-matrix approximation has often been adopted to describe a dilute Fermi gas. This approximation, originally considered by Galitskii for a repulsive inter-particle interaction (Galitskii-1958), has later been widely utilized for an attractive Fermi gas to describe the BCS-BEC crossover from strongly overlapping Cooper pairs in weak coupling to non-overlapping composite bosons in strong coupling. Several variants of the $t$-matrix approximation have been considered in the literature, which are distinguished by the degree of self-consistency allowed in the building blocks of the diagrammatic structure. Here, we perform a systematic and comparative study of all possible variants on the degree of self-consistency for the $t$-matrix approximation in an attractive Fermi gas, which enables us to confront their outcomes for thermodynamic and dynamical quantities on the same footing in an unbiased way. For definiteness, only the normal phase above the superfluid critical temperature is considered. The dispute that can be raised in this context, about the adequateness of introducing progressive degrees of self-consistency over and above the non-self-consistent $t$-matrix approximation for an attractive Fermi gas, parallels the recent interest in the literature on assessing the importance of various degrees of self-consistency in the context of semiconductors and insulators.

Tables

Table 1. Table 1: Short names adopted for the various t 𝑡 t -matrix approaches, according to the conventions introduced in the expressions ( 4 ) and ( 2 ), with the corresponding key references to these approaches.

$\left(G^{\left(a\right)}{G}^{\left(b\right)}\right)G^{\left(c\right)}$	Reference
$\left(G_0{G}_0\right)G_0$	Ref. PPSC-2002
$\left(G_0{G}_0\right)G$	Ref. Ohashi-2012
$\left(G{G}_0\right)G_0$	Ref. Levin-1997
$\left(GG\right)G_0$	Ref. Micnas-2014
$\left(GG\right)G$	Refs. Haussmann-1993 ; Haussmann-1994

Table 2. Table 2: Values of the coefficient α 𝛼 \alpha of Eq. ( 14 ) obtained within the various t 𝑡 t -matrix approaches, both from analytic calculations (th) and extrapolation of numerical results (extr). One analytic value is not available (na).

$\left(G^{\left(a\right)}{G}^{\left(b\right)}\right)G^{\left(c\right)}$	$\alpha$(th)	$\alpha$(extr)
$\left(G_0{G}_0\right)G_0$	1	1.02
$\left(G_0{G}_0\right)G$	1.07	1.12
$\left(G{G}_0\right)G_0$	-0.5	-0.48
$\left(GG\right)G_0$	na	-2.00
$\left(GG\right)G$	-1	-1.02

Table 3. Table 3: Numerical values of the critical temperature T c subscript 𝑇 𝑐 T_{c} , chemical potential μ c subscript 𝜇 𝑐 \mu_{c} , and contact C c subscript 𝐶 𝑐 C_{c} at T c subscript 𝑇 𝑐 T_{c} are reported for five characteristic couplings in the crossover region about unitarity and for all t 𝑡 t -matrix approaches. A few numerical values are not available (na) for the reasons discussed in the text. For each coupling, reference to different t 𝑡 t -matrix approaches follows the conventions (from top to bottom) of Table  1 .

${\left(k_F{a}_F\right)}^{-1}$	$T_c/T_F$	${\mu }_c/E_F$	$C_c/k_F^4$
	0.08495	0.7997	0.01462
	0.07294	0.7381	0.01139
-1.0	na	na	na
	0.07132	0.8006	0.01930
	0.06776	0.7375	0.01874
	0.1649	0.6846	0.04595
	0.1358	0.6018	0.03076
-0.5	0.1132	0.7230	0.03999
	0.1112	0.7269	0.04099
	0.1077	0.6226	0.03892
	0.2429	0.3655	0.1479
	0.2034	0.3059	0.1036
0.0	0.1709	0.5096	0.1108
	0.1451	0.5734	0.09513
	0.1505	0.4000	0.09170
	0.2588	-0.1890	0.3075
	0.2361	-0.2005	0.2638
0.5	0.1975	0.04935	0.2522
	0.1674	0.2176	0.2148
	0.1870	-0.02339	0.2182
	0.2351	-0.9986	0.4759
	0.2320	-0.9987	0.4666
1.0	0.2097	-0.8042	0.4418
	0.1918	-0.6250	0.4150
	0.2062	-0.8092	0.4342

Equations

$$G(k)$$

$$\Sigma (k)$$

$$\Gamma (Q)$$

$$R_{\mathrm{pp}}(Q)$$

$$G(\mathbf{r},\tau )$$

$$\Gamma (\mathbf{r},\tau )$$

$$\Sigma (r)$$

$$R_{\mathrm{pp}}(r)$$

$$\Delta R_{\mathrm{pp}}(r)=R_{\mathrm{pp}}(r)-R_{\mathrm{pp}}^{(0)}(r)=G^{(\mathrm{a})}(r)G^{(\mathrm{b})}(r)-G_0(r{)}^2,$$

$$\big{[}\Gamma(\mathbf{Q}=0,\Omega_{\nu}=0;T,\mu)\big{]}^{-1}=0.$$

$$n=-2G(\mathbf{r}=0,\tau \to {\beta }^{-};T,\mu )$$

$$T_{c}^{\text{BCS}}=\frac{8e^{\gamma}E_{F}}{\pi e^{2}}\exp\Big{(}\frac{\pi}{2k_{F}a_{F}}\Big{)}$$

$$T_{c}^{\text{BEC}}=\frac{2\pi}{\zeta(3/2)^{2/3}}\frac{\big{(}n_{B}\big{)}^{2/3}}{m_{B}}\simeq 0.218\,E_{F}.$$

$$\frac{T_c-T_c^{\text{BEC}}}{T_c^{\text{BEC}}}=\frac{\alpha }{3\pi }(k_F{a}_F{)}^3,$$

$$\frac{\mu^{\text{Gal}}}{E_{F}}=1+\frac{4}{3\pi}(k_{F}a_{F})+\frac{4}{15\pi^{2}}\big{[}11-2\ln 2\big{]}(k_{F}a_{F})^{2}$$

$$\frac{a_B}{a_F}=\underset{k_F{a}_F\to 0^{+}}{\lim }\frac{3\pi }{4(k_F{a}_F)}\frac{{\mu }_B}{E_F}.$$

$$C_c=\int \frac{d\mathbf{Q}}{(2\pi {)}^3}T\sum _{\nu }\Gamma (\mathbf{Q},{\Omega }_{\nu })e^{i{\Omega }_{\nu }{0}^{+}}=\Gamma (\mathbf{r}=0,\tau \to {\beta }^{-}),$$

$$n_{\sigma }(\mathbf{k})=-G(\mathbf{k},\tau \to {\beta }^{-}).$$

$$\frac{C_{c}}{k_{F}^{4}}=\frac{4(k_{F}a_{F})^{2}}{9\pi^{2}}\bigg{(}1+\frac{12}{35\pi}(11-2\ln 2)(k_{F}a_{F})\bigg{)},$$

$$A(\mathbf{k},\omega )$$

$${\int }_{-\infty }^{+\infty }d\omega A(\mathbf{k},\omega )=1.$$

$$\bar{G}(\mathbf{k},z)=\frac{p_1+p_{2}z+\cdots +p_r{z}^{r-1}}{q_1+q_{2}z+\cdots +q_r{z}^{r-1}+z^r}.$$

$$N(\omega )=\int \frac{d\mathbf{k}}{(2\pi {)}^3}A(\mathbf{k},\omega ).$$

$$\frac{ϵ(k)}{E_F}={\mu }^{\prime }-\sqrt{(k^2-k_L^2{)}^2+{\Delta }^{\prime 2}}.$$

$$G(\mathbf{k},{\omega }_n)\to G(\mathbf{k},\tau )\to G(\mathbf{r},\tau )$$

$$G(\mathbf{k},{\omega }_n)\underset{{\omega }_n\to \infty }{\simeq }{G}_0(\mathbf{k},{\omega }_n)=\frac{1}{i{\omega }_n-{\xi }_{\mathbf{k}}}$$

$$\Delta G(\mathbf{k},{\omega }_n)=G(\mathbf{k},{\omega }_n)-G_0(\mathbf{k},{\omega }_n),$$

$$G_{0}(\mathbf{k},\tau)=e^{-\xi_{\mathbf{k}}\tau}\big{(}f(\xi_{\mathbf{k}})-1\big{)},$$

$$G_0(\mathbf{k},\tau )=G_0^{(\mathrm{n})}(\mathbf{k},\tau )+G_0^{(\mathrm{a})}(\mathbf{k},\tau )$$

$$G_0^{(\mathrm{n})}(\mathbf{k},\tau )$$

$$G_0^{(\mathrm{a})}(\mathbf{k},\tau )$$

$$G_0^{(\mathrm{a})}(\mathbf{r},\tau )=-\frac{e^{\mu \tau }{e}^{-\frac{r^2}{4\tau }}}{8{\pi }^{3/2}{\tau }^{3/2}}$$

$$\tilde{\Delta }G(\mathbf{k},\tau )=G(\mathbf{k},\tau )-G_0^{(\mathrm{a})}(\mathbf{k},\tau ),$$

Full text

I Introduction

The method of functional derivatives provides a general framework to deal with quantum many-particle systems in a non-perturbative fashion. This method, that was originally introduced in condensed matter by Martin and Schwinger MS-1959 and later adopted by Hedin Hedin-1965 (see also Ref. Strinati-RNC-1988 ), starts from the exact equations for the single- and two-particle Green’s functions (namely, the Dyson and Bethe-Salpeter equations, respectively) and introduces approximations only at the level of the kernels of these integral equations.

A systematic method for selecting non-perturbative approximations which satisfy the conservation laws (the so-called “conserving approximations”) has been formulated by Baym and Kadanoff BK-1961 . This method intimately relates the kernel of the Dyson equation (namely, the single-particle self-energy) with the kernel of the Bethe-Salpeter equation. In both cases, these kernels are expressed as functionals of the single-particle propagator that solves the Dyson equation. In this context, the need for a “$\Phi$-derivable” choice of the single-particle self-energy and for the self-consistent solution of both equations has been emphasized Baym-1962 .

Although the use of a fully self-consistent and conserving approximation appears mandatory when dealing with physical problems that involve transport and localization, in other circumstances non-self-consistent or partly self-consistent approximations may provide physically more sensible results with respect to the self-consistent one. An example is provided by the fluctuation exchange (FLEX) approximation introduced for the repulsive Hubbard model (even at half-filling), whereby both the self-consistent and the non-self-consistent versions have been investigated Scalapino-1989 ; Kotliar-2005 . More recently, interest in comparing the results of the self-consistent vs non-self-consistent approaches has arose also in the context of the GW approximation for semiconductors and insulators. In this case, the non-self-consistent calculations turn out to better compare with the experimental values of the band gaps with respect to the self-consistent calculations Kresse-2018 . Similar conclusions have further been drawn in the context of a simpler model Berger-2018 .

In the context of a (dilute) Fermi gas with an attractive inter-particle interaction, the $t$-matrix approximation appears as a natural candidate to describe the system while it evolves throughout the BCS-BEC crossover. The first pioneering approach in this respect goes back to the work by Nozières and Schmitt-Rink (NSR) NSR-1985 , where a simplified version of the non-self-consistent $t$-matrix approximation proved sufficient to highlight the main features of the crossover physics in the normal phase above the superfluid critical temperature. Later on, the NSR approach was extended, either to improve on the treatment of the non-self-consistent $t$-matrix approximation PPSC-2002 , or to include various degrees of self-consistency within this approximation, ranging from partial Levin-1997 ; Ohashi-2012 ; Micnas-2014 to full Haussmann-1993 ; Haussmann-1994 self-consistency. As expected, depending on the degree of self-consistency different numerical results were obtained for various physical quantities, ranging from thermodynamic to dynamic. However, direct comparison among the results obtained with various degrees of self-consistency has been hindered by the (sometimes even drastic) numerical approximations that were introduced in the calculations on top of a specific choice about the degree of self-consistency. For these reasons, it appears that a systematic and direct comparison of the results obtained by adopting various degrees of self-consistency on the $t$-matrix approximation is still lacking, especially when this comparison would be made in an unbiased way by retaining the same level of numerical accuracy for all different approaches.

Primary purpose of this paper is to fill this gap, by undertaking the above systematic study on all the five degrees of self-consistency that have been considered thus far in the literature within the $t$-matrix approximation in the normal phase above the superfluid critical temperature PPSC-2002 ; Levin-1997 ; Ohashi-2012 ; Micnas-2014 ; Haussmann-1993 ; Haussmann-1994 . Although the results of the present study may not lead to definite conclusions, about which one of the above five approximations could account best for the thermodynamic and dynamical properties of a (dilute) Fermi gas undergoing the BCS-BEC crossover, it appears nevertheless interesting and relevant (if not timely) to discover how these properties get modified when passing from one to the other of these five approaches. In addition, the high precision of the numerical calculations that we have implemented has enabled us to accurately check how the two distinct (BCS) and (BEC) limits of the crossover are recovered by the alternative $t$-matrix approaches. Specifically, in the BEC limit the residual interaction among composite bosons extracted from our numerical calculations turns out to be in excellent agreement with the analytic estimates that we also provide (which correct a previous analytic estimate obtained in Refs. Haussmann-1993 ; Haussmann-1994 within the fully self-consistent $t$-matrix approach). In the BCS limit, on the other hand, we have found that a partially self-consistent $t$-matrix approach, developed in Ref. Levin-1997 and often utilized in the literature, breaks down when one avoids using the set of approximations that normally accompany its implementation.

Finally, it should be recalled that, in order to get a refined agreement with Quantum Monte Carlo and experimental data available for a Fermi gas in the unitary region of the crossover intermediate between the BCS and BEC regimes, it may in any case be required to go beyond the $t$-matrix approximation Pisani-2018-I and include on top of it also a class of vertex corrections associated with the Gorkov-Melik-Barkhudarov (GMB) contribution Gorkov-1961 .

The plan of the paper is as follows. Section II sets up the theoretical framework and specifies in details how the alternative $t$-matrix approaches with various degrees of self-consistency need to be handled. Sections III and IV report on the numerical results obtained within the above alternative $t$-matrix approaches, for the thermodynamic and dynamical quantities of interest, respectively. Sections V gives our conclusions. Appendix A gives a detailed account about the numerical procedures we have adopted to achieve (partial or full) self-consistency within the various $t$-matrix approaches, for all relevant sub-units of the many-body structure. Appendix B discusses the optimization procedure that was found necessary to achieve the required convergence toward self-consistency, within some of the above alternative $t$-matrix approaches. Finally, Appendix C compares the results of two partially self-consistent $t$-matrix approaches with that of their approximate treatments usually utilized in the literature.

II Alternative $t$-matrix approaches

In this Section, we set up the theoretical framework for the alternative $t$-matrix approaches that can be used to describe a Fermi gas with an attractive contact interaction throughout the BCS-BEC crossover, in the normal phase at a temperature $T$ above the superfluid critical temperature $T_{c}$. Only a balanced situation, with equal number of spin up and spin down fermions, will be considered in this paper.

A. Basic equations

Within this framework, the basic expressions read (in the following, we set the Planck constant $\hbar$ and Boltzmann constant $k_{B}$ equal to unity):

[TABLE]

Here, $G$ is the single-particle fermionic propagator and $G_{0}$ its non-interacting counterpart, $\Sigma$ is the self-energy, $\Gamma$ is the particle-particle propagator, and $R_{\mathrm{pp}}$ is the (regularized) particle-particle bubble. These quantities are drawn pictorially in Fig. 1. In addition, $k=(\mathbf{k},\omega_{n})$ is a fermionic four-vector with wave vector $\mathbf{k}$ and fermionic Matsubara frequency $\omega_{n}=(2n+1)\pi\,T$ ($n$ integer), and $Q=(\mathbf{Q},\Omega_{\nu})$ a bosonic four-vector with wave vector $\mathbf{Q}$ and bosonic Matsubara frequency $\Omega_{\nu}=2\pi\nu\,T$ ($\nu$ integer). Finally, $m$ is the fermionic mass and $a_{F}$ the scattering length of the associated two-fermion problem.

Note that, in the above equations, the single-particle propagators have been distinguished as $G^{(\mathrm{a})}$, $G^{(\mathrm{b})}$, and $G^{(\mathrm{c})}$, meaning that each of them can be either the dressed $G$ or the bare $G_{0}$, in such a way that different $t$-matrix approaches can be realized by selecting different combinations of these functions. Specifically, the short-hand notation $(G^{(\mathrm{a})}G^{(\mathrm{b})})G^{(\mathrm{c})}$ will be used to identify a given $t$-matrix approach, where the propagators $G^{(\mathrm{a})}$ and $G^{(\mathrm{b})}$ within the parentheses correspond to those entering the particle-particle bubble (4), while the external propagator $G^{(\mathrm{c})}$ enters the self-energy (2). In the following, we shall consider the five combinations of $G^{(\mathrm{a})}$, $G^{(\mathrm{b})}$, and $G^{(\mathrm{c})}$ reported in Table 1, with the corresponding references where the various approaches have been discussed for a Fermi gas with an attractive contact interaction.

B. Routes toward self-consistency

Except for the non-self-consistent $(G_{0}G_{0})G_{0}$ approach, the equations (1)-(4) must be solved in a self-consistent way. To this end, we follow the numerical procedure developed originally in Refs. Haussmann-1993 and Haussmann-1994 . The procedure makes use of the Fourier transforms from the $(\mathbf{k},\omega_{n})$ and $(\mathbf{Q},\Omega_{\nu})$ space to the $(\mathbf{r},\tau)$ space, according to the expressions:

[TABLE]

Analogous transformations hold for $\Sigma(\mathbf{k},\omega_{n})$ and $R_{\mathrm{pp}}(\mathbf{Q},\Omega_{\nu})$. Here, $\tau$ is the imaginary time which varies in the interval $(0,1/T)$. In the $r=(\mathbf{r},\tau)$ space, the equations (2) and (4) acquire the simple form:

[TABLE]

where $\Lambda$ is an appropriate regularization constant that depends on the cutoff in the $\mathbf{k}$-integral of Eq. (4) and becomes infinite together with that cutoff Haussmann-1993 . To avoid dealing directly with $\Lambda$, in $r$-space it is convenient to work in terms of the difference

[TABLE]

where $R^{(0)}_{\mathrm{pp}}(r)$ is the regularized particle-particle bubble built on the non-interacting $G_{0}$. Equations (1), (3), (7), and (9), together with the Fourier transforms (5) and (6), form a complete set of equations that need to be solved self-consistently. The flowchart shown in Fig. 2 summarizes schematically the various routes toward self-consistency to be followed when adopting the alternative $t$-matrix approaches of Table 1.

Although the non-self-consistent $(G_{0}G_{0})G_{0}$ approach does not require any self-consistent cycling, we have reported it in Fig. 2(a) since its calculation has always to be performed at a preliminary level, to the extent that it is also used as input for the self-consistent calculations. In this approach, $R^{(0)}_{\mathrm{pp}}(Q)$ is directly calculated according to Eq. (4) with two non-interacting propagators $G_{0}$ in the place of $G^{(\mathrm{a})}$ and $G^{(\mathrm{b})}$. Figure 2(b) shows the flowchart for the $(G_{0}G_{0})G$ approach (referred to as the “extended $t$-matrix” approach in Ref. Ohashi-2012 ), where self-consistency is present only in the external $G$, while the particle-particle propagator $\Gamma$ coincides with its bare counterpart $\Gamma_{0}$ built on $R^{(0)}(Q)$. Figure 2(c) shows the flowchart for both the $(GG_{0})G_{0}$ and $(GG)G_{0}$ approaches, where self-consistency enters only the particle-particle propagator $\Gamma$, while the external propagator $G^{(\mathrm{c})}$ is a non-interacting $G_{0}$. Finally, Fig. 2(d) shows the flowchart for the fully self-consistent $(GG)G$ approach. Note that, in the flowcharts of Figs. 2(c) and (d), the regularized particle-particle bubble $R_{\mathrm{pp}}(\mathbf{Q},\Omega_{\nu})$ is obtained by Fourier transforming the difference $\Delta R_{\mathrm{pp}}(\mathbf{r},\tau)$ defined in Eq. (9) and then by adding to it the regularized particle-particle bubble $R^{(0)}_{\mathrm{pp}}(\mathbf{Q},\tau)$ of the non-self-consistent approach.

Most functions appearing in the flowcharts of Fig. 2 suffer from a slowly decaying tail in the variables $(\mathbf{k},\omega_{n})$ or $(\mathbf{Q},\Omega_{\nu})$, which implies a corresponding singular behavior when $(r,\tau)\rightarrow 0$. For this reason, the Fourier transforms should be performed on a logarithmic scale, following the prescriptions given in Refs. Haussmann-1994 and Haussmann-2007 . In addition, one also needs to subtract appropriate semi-analytic expressions from the functions to be numerically Fourier transformed, in order to make their slow decay faster (or, alternatively, their singular behavior weaker). These semi-analytic expressions have to be known also for the transformed representation, so that they can be added back to the functions after having performed the Fourier transform. A detailed account of the semi-analytic expressions used in the numerical calculations is given in Appendix A.

C. Thouless criterion

As already mentioned, all five $t$-matrix approaches considered in this paper will be examined on equal footing throughout the whole BCS-BEC crossover (also with emphasis on the analytic results that can be obtained separately in the BCS and BEC limits). The BCS-BEC crossover is spanned in terms of the (dimensionless) coupling parameter $(k_{F}a_{F})^{-1}$, where $k_{F}=(3\pi^{2}n)^{1/3}$ is the Fermi wave vector associated with the particle density $n$. In practice, the crossover between the BCS and BEC regimes is essentially exhausted within the range $-1\lesssim(k_{F}a_{F})^{-1}\lesssim+1$ across the unitary limit at $(k_{F}a_{F})^{-1}=0$ (for a recent comprehensive account of the BCS-BEC crossover, see Ref. Physics-Reports-2018 ).

In the present paper, we are interested in the normal phase above the the critical temperature $T_{c}$ of the superfluid transition, where the numerical value of $T_{c}$ depends on the specific theoretical approximation one is adopting to describe the Fermi gas. For all five $t$-matrix approaches we are considering, $T_{c}$ is determined by the Thouless criterion Thouless-1960 , in the form:

[TABLE]

This condition has to be supplemented by the density equation to determine the chemical potential $\mu$

[TABLE]

(with $\beta=1/T$ the inverse temperature), where the factor of $2$ accounts for the spin multiplicity. In practice, one fixes the values of the coupling $(k_{F}a_{F})^{-1}$ and of the temperature $T$ to determine $\mu$ from Eq. (11), and then uses this value of $\mu$ to determine $T_{c}$ from Eq. (10). However, for temperatures close to $T_{c}$ this simple iterative procedure may not work properly as far as the $(GG_{0})G_{0}$, $(GG)G_{0}$ and $(GG)G$ approaches are concerned, where difficulties are found in the convergence of the iterative procedure toward self-consistency. To overcome these difficulties, we have found it necessary to introduce two different refinements of the above iterative procedure, the first of which enables us to get close to $T_{c}$ from $T>T_{c}$ while the second one allows us to work exactly at $T=T_{c}$. These refinements are described in detail in Appendix B.

III Numerical results for thermodynamic quantities

Several thermodynamic quantities of interest can be obtained directly in terms of the single-particle fermionic propagator $G$ of Eq. (1) and of the particle-particle propagator $\Gamma$ of Eq. (3). Here, we report on the numerical results obtained within the five alternative $t$-matrix approaches that we are considering, for the critical temperature and chemical potential as well as for the Tan’s contact.

A. Critical temperature

The results for the critical temperature $T_{c}$ obtained from the Thouless criterion (10) are shown in Fig. 3 over a wide coupling range for all the five $t$-matrix approaches reported in Table 1. Several interesting features can be highlighted when comparing the results of the various approaches:

(i) All approaches are seen to interpolate rather well between the BCS and BEC critical temperatures (indicated by black dotted lines in Fig. 3), with the notable exception of the $(GG_{0})G_{0}$ approach which fails to reach the BCS limit since in this case the critical temperature collapses abruptly to zero at coupling $(k_{F}a_{F})^{-1}\simeq-1$. We attribute this failure to the asymmetric treatment of the single-particle propagators $G^{(\mathrm{a})}$ and $G^{(\mathrm{b})}$ in the particle-particle bubble (4) that enters the particle-particle propagator $\Gamma$. The asymmetry generates an artificial imbalance between spin up and down species, which acts to suppress the critical temperature of the superfluid transition (cf. also Appendix C). This feature has apparently passed unnoticed in the literature. It is for this reason that, in what follows, the results of the $(GG_{0})G_{0}$ approach will not be reported for the BCS limit.

(ii) The behavior of the critical temperature in the BCS (weak-coupling) limit $(k_{F}a_{F})^{-1}\ll-1$ is shown as a function of $k_{F}a_{F}$ in the inset of Fig. 3), for all $t$-matrix approaches with the exception of the $(GG_{0})G_{0}$ approach for the reasons discussed in point (i) above. In all cases, the numerical results have been extended toward the extreme BCS limit $k_{F}a_{F}\rightarrow 0^{-}$ through a parabolic extrapolation (dots). Within a numerical error of the order of $1\%$ footnote-1 , in this limit we obtain that the $(GG)G$ and $(GG)G_{0}$ approaches recover the value of the BCS temperature:

[TABLE]

where $\gamma$ is the Euler’s constant. On the other hand, the $(G_{0}G_{0})G$ and $(G_{0}G_{0})G_{0}$ approaches reproduce the BCS critical temperature only within logarithmic accuracy, in the sense that they recover the result (12) with the pre-factor divided by $e^{1/3}$ (see Ref. Pisani-2018-I for a discussion of the origin of this spurious factor in those approaches that do not dress the bare particle-particle propagator $\Gamma_{0}$).

(iii) In the crossover region about unitarity, it turns out that the inclusion of (even a partial degree of) self-consistency in the particle-particle propagator $\Gamma$ of Eq. (3) acts to suppress the maximum of $T_{c}$, which otherwise occurs for the $(G_{0}G_{0})G$ and $(G_{0}G_{0})G_{0}$ approaches.

(iv) In the BEC (strong-coupling) limit $(k_{F}a_{F})^{-1}\gg 1$, all approaches reproduce the value of the critical temperature for a condensate of non-interacting composite bosons made up of fermion pairs (with mass $m_{B}=2m$ and density $n_{B}=n/2$):

[TABLE]

However, alternative $t$-matrix approaches differ in the way the value (13) is reached when $(k_{F}a_{F})^{-1}\gg 1$. Specifically, the sub-leading behavior of $T_{c}$ when approaching $T_{c}^{\text{BEC}}$ can be characterized by the expression

[TABLE]

where the values of the coefficient $\alpha$ for the various $t$-matrix approaches are listed in Table 2.

The theoretical values reported in Table 2 for the $(GG)G$ and $(GG_{0})G_{0}$ approaches are taken from Refs. Haussmann-1993 and Chen-2005 , respectively, while we have calculated independently those for the $(G_{0}G_{0})G_{0}$ and $(G_{0}G_{0})G$ approaches (the corresponding analytic calculations are not reported here owing to their complexity). Note that, for those approaches that include even a partial degree of self-consistency in the particle-particle propagator $\Gamma$, $\alpha<0$ such that the value (13) is approached from below; the opposite occurs for the remaining approaches.

We have also compared the numerical values we have obtained for $T_{c}$ within the various $t$-matrix approaches with other published data for the same quantity. For the $(G_{0}G_{0})G_{0}$, $(G_{0}G_{0})G$, and $(GG)G$ approaches we found good agreement with the data published in Refs. PPSC-2002 , Ohashi-2018 , and Haussmann-1994 , respectively. For the $(GG_{0})G_{0}$ and $(GG)G_{0}$ approaches, on the other hand, direct comparison with previous data throughout the BCS-BEC crossover is not possible, since within these approaches the curves for $T_{c}$ have only been calculated with additional approximations affecting the form of the particle-particle propagator and of the self-energy Levin-2010 Micnas-2014 . In these cases, our data compares well with those of Refs. Levin-2010 and Micnas-2014 only in the strong-coupling regime $(k_{F}a_{F})^{-1}\gtrsim 1$, where the additional approximations introduced in those references remain valid at $T\simeq T_{c}$ (see also Appendix C).

Finally, a comparison with available quantum Monte Carlo (both QMC and DMC) data has been reported in Fig. 3. Note how these data show a steeper coupling dependence with respect to the $t$-matrix calculations. As already mentioned in the Introduction, this steeper dependence can be accounted for by a further inclusion of the GMB vertex corrections Pisani-2018-I .

B. Chemical potential

The corresponding results for the chemical potential $\mu_{c}$ calculated at the critical temperature $T_{c}$ are reported in Fig. 4 as a function of coupling $(k_{F}a_{F})^{-1}$, for each of the five $t$-matrix approaches here considered. The main features that can be identified from this plot are:

(i) In the weak-coupling (BCS) limit $(k_{F}a_{F})^{-1}\ll-1$, the $(GG)G$ and $(G_{0}G_{0})G$ approaches recover the Galitskii’s expression Galitskii-1958

[TABLE]

up to second order in $(k_{F}a_{F})$, while the $(G_{0}G_{0})G_{0}$ and $(G_{0}G_{0})G$ approaches recover this expression up to first order only.

(ii) In the strong-coupling (BEC) limit $(k_{F}a_{F})^{-1}\gg 1$, all approaches recover the value $\mu_{c}=-\varepsilon_{0}/{2}$ at leading order in $k_{F}a_{F}$, where $\varepsilon_{0}=(ma_{F}^{2})^{-1}$ is the binding energy of the composite bosons. At sub-leading order, on the other hand, the behavior of $\mu_{c}$ depends on the approach. In particular, for the $(G_{0}G_{0})G_{0}$ and $(G_{0}G_{0})G$ approaches the sub-leading term vanishes exponentially when $(k_{F}a_{F})^{-1}\gg 1$, while for the $(GG_{0})G_{0}$, $(GG)G_{0}$, and $(GG)G$ approaches the sub-leading term vanishes linearly in $(k_{F}a_{F})$ (a behavior which has been evidenced in the inset of Fig. 4). This difference is due to the fact that self-consistency in the particle-particle propagator introduces a residual repulsive interaction between composite bosons, as discussed in subsection III-C.

By the present methods, the chemical potential $\mu$ can be calculated not only at $T_{c}$ but also above $T_{c}$. Figure 5 shows $\mu$ as a function of temperature for three characteristic couplings [$(k_{F}a_{F})^{-1}=(-0.5,0,0.5)$] spanning the crossover regime, for all five $t$-matrix approaches we are considering. On the low-$T$ side, the curves terminate at the respective critical temperature $T_{c}$ given in Fig. 3. Note that, for the $(G_{0}G_{0})G_{0}$ and $(G_{0}G_{0})G$ approaches where no degree of self-consistency is introduced in the particle-particle propagator, the curves of $\mu(T)$ show a maximum above $T_{c}$. We have also verified that the curve of Fig. 5(b), which corresponds to the $(GG_{0})G_{0}$ approach at unitarity, compares well with the numerical results reported in Ref. Drummond-2008 . To the best of our knowledge, this is, in fact, the only other reference where a calculation based on the complete $(GG_{0})G_{0}$ approach was performed, without recourse to additional simplifying approximations (albeit in Ref. Drummond-2008 the calculation was limited to the coupling $(k_{F}a_{F})^{-1}=0$ only).

In both Figs. 4 and 5, comparison has also been added with available QMC data, which show overall good agreement with the results of the $t$-matrix approaches.

C. Scattering length of composite bosons

The residual repulsive interaction between composite bosons is characterized by a finite value of the scattering length $a_{B}$. This can, in turn, be determined by comparing the chemical potential for the composite bosons $\mu_{B}=2\mu_{c}+\varepsilon_{0}$ in the strong-coupling limit $(k_{F}a_{F})^{-1}\gg 1$ with the chemical potential $\mu^{0}_{B}=8\pi a_{B}n_{B}/m_{B}$ of a dilute Bose gas at $T_{c}$. This comparison yields:

[TABLE]

The values of $a_{B}$ extrapolated in this way for the $(GG)G$, $(GG)G_{0}$ and $(GG_{0})G_{0}$ approaches are shown in Fig. 6(a). Note, in particular, that the $(GG)G$ fully self-consistent approach yields $a_{B}/a_{F}\simeq 1.16$, which contrasts with the value $a_{B}/a_{F}=2$ obtained analytically in Ref. Haussmann-1993 .

The above limiting numerical values for the bosonic scattering length $a_{B}$, obtained within the $t$-matrix approaches that include the dressed single-particle propagator $G$ in the particle-particle propagator $\Gamma$, can be compared with the corresponding analytic results obtained by expressing the self-consistently dressed $\Gamma$ as an infinite series of diagrams in terms of the bare $G_{0}$ and $\Gamma_{0}$, and then by retaining only the leading-order corrections to the bare particle-particle propagator $\Gamma_{0}$ in the BEC ($\mu/T\to-\infty$) limit.

For the $(GG)G_{0}$ approach, this corresponds to considering the left diagram of Fig. 6(b) (together with the corresponding diagram where is the lower fermionic line to get dressed). Comparison of the analytic evaluation of this diagram in the BEC limit as done in Refs. Haussmann-1993 and Pieri-2000 with the leading self-energy correction $\Sigma_{B}=8\pi a_{B}n_{B}/m_{B}$ for a dilute Bose gas yields the value $a_{B}/a_{F}=2$ (cf. the upper curve of Fig. 6(a)).

For the $(GG_{0})G_{0}$ approach, on the other hand, the absence of the self-energy correction in the lower fermionic line then eliminates the multiplicity of $2$ for the left diagram of Fig. 6(b), yielding $a_{B}/a_{F}=1$ (cf. the lower curve of Fig. 6(a)) Note, however, that this reduction of $a_{B}$ by half in comparison with the more symmetric $(GG)G_{0}$ approach is somewhat artificial, since it corresponds to an incomplete symmetrization of the bosonic interaction vertex.

Finally, for the $(GG)G$ approach also the right diagram of Fig. 6(b) contributes at the leading order, again with a multiplicity of $2$ due to the corresponding dressing of the lower fermionic line. This diagram (which had apparently escaped to the analysis of Ref. Haussmann-1993 ) yields a correction $-0.842a_{F}$ to the bosonic scattering length $a_{B}$ in the BEC limit, thereby resulting altogether in the value $a_{B}/a_{F}=2-0.842=1.158$ in excellent agreement (within $0.2\%$) with that obtained by our numerical extrapolation (cf. the middle curve of Fig. 6(a)). This analytic estimate for the right diagram of Fig. 6(b) can be readily obtained by noting that, in the BEC limit, its leading-order behavior coincides with that of the diagram introduced in Ref. Pisani-2018-I to include the GMB correction for $T_{c}$ throughout the BCS-BEC crossover, whose contribution to $a_{B}$ in the BEC limit was there estimated to be $-0.842a_{F}$. It should also be remarked that this identification between (the numerical values of) the two diagrams holds in the BEC limit only. For this reason, the self-consistent $t$-matrix approximation fails to recover the GMB reduction factor for $T_{c}$ in the BCS limit Pisani-2018-I .

D. Tan’s contact

An important physical quantity that characterizes a Fermi gas with short-range interaction is the Tan’s contact Tan-2008-I ; Tan-2008-II ; Tan-2008-III , which connects two-particle correlations at short distances with thermodynamics. Here, we calculate the contact $C_{c}$ at $T_{c}$ within the various $t$-matrix approaches as the trace of the particle-particle propagator

[TABLE]

according to an expression introduced in Ref. PPS-NP-2009 . The results are shown in Fig. 7 throughout the BCS-BEC crossover. For internal consistency, we have also verified numerically that the values of the contact $C_{c}$ obtained by Eq. (17) coincide with the coefficient of the $\mathbf{k}^{-4}$ tail of the wave-vector distribution (per spin component $\sigma$)

[TABLE]

We have further verified that, at $T_{c}\,(\ll T_{F})$ in weak coupling, the expansion for the contact in powers of $(k_{F}a_{F})$

[TABLE]

that results by taking the derivative with respect to $a_{F}^{-1}$ of the expression for the total energy at $T=0$ obtained by Galitskii in powers of $(k_{F}a_{F})$ Galitskii-1958 , is recovered by the $(G_{0}G_{0})G_{0}$ and $(G_{0}G_{0})G$ approaches up to the leading order and by the $(GG)G$ and $(GG)G_{0}$ approaches up to the next-to-leading order in $(k_{F}a_{F})$.

E. Summary of the main thermodynamic results

The numerical values of the critical temperature $T_{c}$, and of the chemical potential $\mu_{c}$ and the contact $C_{c}$ at $T_{c}$, can be read off directly from Figs. 3, 4, and 7, respectively, for each of the five $t$-matrix approaches that we have considered in this paper. It might be useful, however, to summarize the values of the above physical quantities for a few characteristic couplings in the crossover region of most interest. Accordingly, we report in Table 3 a list of these values for five couplings in the interval $-1\leq(k_{F}a_{F})^{-1}\leq+1$. In this way, differences in the results of the various $t$-matrix approaches can be most readily appreciated.

IV Numerical results for dynamical quantities

The results obtained in Section III, about the single-particle fermionic propagator $G$ of Eq. (1), can now be utilized to obtain a number of spectral features for the attractive Fermi gas in the normal phase throughout the BCS-BEC crossover. To this end, a suitable method is required to perform the analytic continuation from Matsubara to real frequencies, as discussed next.

A. Method for analytic continuation

To assess the dynamical properties of a Fermi gas, the analytic continuation $\bar{G}(\mathbf{k},z)$ of the single-particle fermionic propagator $G(\mathbf{k},\omega_{n})$ is required over the complex $z$-plane. The function $\bar{G}(\mathbf{k},z)$ takes the values $G(\mathbf{k},\omega_{n})$ along the imaginary axis at the Matsubara frequencies $z=i\omega_{n}$, and is analytic everywhere except on the real frequency axis where its imaginary part is discontinuous due to the time-reversal symmetry condition $\bar{G}(\mathbf{k},z^{*})=\bar{G}(\mathbf{k},z)^{*}$. The single-particle spectral function $A(\mathbf{k},\omega)$ is then obtained in terms of the discontinuity of $\bar{G}(\mathbf{k},z)$ across the real frequency axis:

[TABLE]

The positive definite function $A(\mathbf{k},\omega)$ is normalized according to:

[TABLE]

To evaluate the single-particle spectral function (20), a procedure is required to obtain the function $\bar{G}(\mathbf{k},z)$ just above the real frequency axis from the known values $\bar{G}(\mathbf{k},z=i\omega_{n})=G(\mathbf{k},\omega_{n})$ on the imaginary axis. To this end, we make use of the method of Padé approximants Serene-1977 ; Beach-1999 which consists in approximating $\bar{G}(\mathbf{k},z)$ for given $\mathbf{k}$ by a ratio of polynomials, in the form:

[TABLE]

Here, the $2r$ (real) coefficients $\{p_{i},q_{i};i=1,\cdots,r\}$ are determined from the values of $\bar{G}(\mathbf{k},z)$ at $2r$ points on the imaginary axis.

This procedure, however, turns out to be quite sensitive to the presence of numerical uncertainties in the values of the input function $G(\mathbf{k},\omega_{n})$, in such a way that the resulting shape of the function $A(\mathbf{k},\omega)$ may turn out to be rather distorted or even to acquire negative values. To mitigate the occurrence of this sort of problems, we have followed a prescription proposed in Ref. Schott-2016 and averaged over several (typically, 12) runs of analytic continuations with different sets of $2r$ (typically, $50$) points selected on the imaginary axis, consistently discarding those runs that yield $A(\mathbf{k},\omega)<0$ for some $\omega$ intervals. In addition, even though the analytical continuation is performed over the half-plane with $\text{Im}(z)>0$, we have sometimes found it useful to include in the sets of $2r$ points the first few (from one up to three) frequencies on the negative imaginary axis, as also proposed in Ref. Schott-2016 .

B. Single-particle spectral function

The single-particle spectral function $A(k,\omega)$, obtained at unitarity and $T_{c}$ for all $t$-matrix approaches, is shown in Fig. 8 for several values of $k=|\mathbf{k}|$. A notable difference results by comparing the various panels of Fig. 8, between the approaches that dress the fermionic propagator $G^{(\mathrm{c})}$ in Eq. (2) (cf. panels (b) and (e)) and those that do not (cf. panels (a), (c), and (d)). As a matter of fact, in the $(GG)G$ and $(G_{0}G_{0})G$ approaches there is essentially no evidence of a double-peak structure (except at $k\simeq 0$), such that the shape of $A(\mathbf{k},\omega)$ is mostly represented by a single peak that shifts from negative to positive frequencies upon increasing $k$. Yet, this behavior appears not to be consistent with what one would expect for a Fermi liquid, as sometimes claimed instead in the literature Salomon-2011 . This is because the single peak broadens up just at $k\simeq k_{F}$ with a width of the order of $E_{F}$, both features being not consistent with the behavior of a Fermi liquid AGD-1975 ; Nozieres-1964 . On the other hand, the $(GG)G_{0}$, $(GG_{0})G_{0}$ and $(G_{0}G_{0})G_{0}$ approaches (whereby $G^{(c)}$ in Eq. (2) remains $G_{0}$) present a persistent double-peak structure through $k\simeq k_{F}$, with an exchange of weight occurring for increasing $k$ between the peaks at negative and positive frequencies.

C. Single-particle density of states

To avoid reference to a specific wave vector in the the single-particle spectral function, yet maintaining the main features of its frequency dependence, one can integrate $A(\mathbf{k},\omega)$ over all $\mathbf{k}$ and obtain the density of states:

[TABLE]

Figure 9 shows the density of states obtained in this way at $T_{c}$ for the coupling values $(k_{F}a_{F})^{-1}=(-0.5,0.0,0.5)$, within the alternative $t$-matrix approaches. In all cases, a depletion is apparent in the density of states about $\omega=0$. The energy width of this depletion is associated to a pseudo-gap that develops in the normal phase above $T_{c}$ due to pairing fluctuations, as a precursor of the pairing gap that occurs in the superfluid phase below $T_{c}$. Although all $t$-matrix approaches present evidence of a pseudo-gap, its detailed structure depends on the specific approach. Let’s consider, for instance, panel (b) of Fig. 9 for the coupling $(k_{F}a_{F})^{-1}=0$. Here we observe that, similarly to $A(\mathbf{k},\omega)$, also for $N(\omega)$ the approaches can be divided in two classes, namely, those that dress the fermionic propagator $G^{(\mathrm{c})}$ in Eq. (2) and those that do not. The first ones present only a narrow and shallow pseudo-gap feature, while the second ones present a wide and deep pseudo-gap feature (which is especially amplified by the $(G_{0}G_{0})G_{0}$ approach). This is, of course, a direct consequence of the behavior of the single-particle spectral functions shown in Fig. 8, where the single peak in $A(\mathbf{k},\omega)$ that crosses $\omega=0$ tends to partially fill the pseudo-gap region for the $(GG)G$ and $(G_{0}G_{0})G$ approaches, while for the $(GG)G_{0}$, $(GG_{0})G_{0}$ and $(G_{0}G_{0})G_{0}$ approaches a depletion persists about $\omega=0$ owing to the double-peak structure of $A(\mathbf{k},\omega)$. For the weaker coupling $(k_{F}a_{F})^{-1}=-0.5$ (cf. panel (a) of Fig. 9), on the other hand, the density of states has the overall shape $N_{0}(\omega)=m^{3/2}\sqrt{2(\omega+\mu)}/(2\pi^{2})$ of the non-interacting system, with only a rather minor pseudo-gap feature occurring about $\omega=0$. While for this coupling the qualitative behavior of $N(\omega)$ is essentially the same for all approaches, the $(GG)G$ and $(G_{0}G_{0})G$ approaches still have a weaker pseudo-gap behavior than the other approaches. Finally, for the coupling $(k_{F}a_{F})^{-1}=0.5$ (cf. panel (c) of Fig. 9) more marked differences appear among the various approaches. In particular, the $(GG)G_{0}$, $(GG_{0})G_{0}$, and $(G_{0}G_{0})G_{0}$ approaches all present quite a wide and deep pseudo-gap feature, the $(G_{0}G_{0})G$ approach presents a wide but rather shallow pseudo-gap feature, and the $(GG)G$ approach shows almost no evidence of a pseudo-gap (to the extent that the density of states does not even go through a local minimum near $\omega=0$).

D. Pseudo-gap temperature

For given coupling, the evolution of the shape of $N(\omega)$ vs $\omega$ can be followed for increasing temperature starting from $T_{c}$. As an example, Fig. 10 shows this temperature evolution when $(k_{F}a_{F})^{-1}=0$ for each $t$-matrix approach we are considering. Quite generally, the depletion of the pseudo-gap region about $\omega=0$ gradually fades away upon increasing temperature, in such a way that a “crossover” temperature $T^{*}$ can be identified as the highest temperature at which the local minimum of $N(\omega)$ near $\omega=0$ eventually disappears.

The crossover temperature $T^{*}$ obtained in this way throughout the BCS-BEC crossover is reported in Fig. 11 for all $t$-matrix approaches. It turns out that the coupling dependence of $T^{*}$ differs considerably for the various approaches, especially for positive couplings on the BEC side of unitarity. In particular, one notices that $T^{*}$ is considerably suppressed for the $(GG)G$ and $(G_{0}G_{0})G$ approaches with respect to the other approaches. In addition, for the $(GG)G$ and $(G_{0}G_{0})G$ approaches there occurs a coupling interval where $T^{*}$ cannot be defined, because the density of states $N(\omega)$ does not have a local minimum near $\omega=0$ even at $T=T_{c}$. This occurs when $(k_{F}a_{F})^{-1}\simeq 0.3$ for the $(G_{0}G_{0})G$ approach and when $(k_{F}a_{F})^{-1}\simeq 0.6$ for the $(GG)G$ approach. Finally, for all approaches $T^{*}$ begins to increase rapidly with coupling around $(k_{F}a_{F})^{-1}\simeq 0.5-0.7$. This is because this coupling regime is where the actual crossover occurs, from a pseudo-gap phase where the depletion in the density of states is rather shallow and due to a truly many-body effect, to a normal-Bose-gas phase where the depletion in the density of states becomes deep and is just evidence of the two-body binding energy of composite bosons Marsiglio-2015 .

E. Luttinger wave vector

The above crossover, between the pseudo-gap and normal-Bose-gas phases, can be characterized in terms of the Luttinger wave vector $k_{L}$. This wave vector was originally considered in Ref. Perali-2011 within the $(G_{0}G_{0})G_{0}$ approach, as the wave vector $k$ at which the back-bending of the lower branch $\epsilon(k)$ of the single-particle dispersion occurs. In Ref. Perali-2011 , this branch was obtained by following the $k$-dependence of the low-energy peak in the single-particle spectral function $A(\bm{k},\omega)$ (cf. panel (a) of Fig. 8), and then by fitting the dispersion of the lower branch $\epsilon(k)$ obtained in this way through a BCS-like form

[TABLE]

Here, $\mu^{\prime}$ and $\Delta^{\prime}$ are fitting parameters (in units of $E_{F}$), with the energy shift $\mu^{\prime}$ needed to account for the dispersion $\epsilon(k)$ away from the weak-coupling regime, and $k$ and $k_{L}$ are in units of $k_{F}$. On physical grounds, a non-vanishing value of $k_{L}$ signals the presence of an underlying Fermi surface, which endows the system with a persistent fermionic character even in the presence of a strong attractive inter-particle interaction. In these terms, the crossover from the pseudo-gap to the normal-Bose-gas phase is considered complete only when $k_{L}$ reaches zero at some critical coupling, thus signaling the eventual disappearance of the underlying Fermi surface.

The above definition of $k_{L}$, introduced in Ref. Perali-2011 within the $(G_{0}G_{0})G_{0}$ approach, can as well be extended to the $(GG_{0})G_{0}$ and $(GG)G_{0}$ approaches here considered, whereby a lower branch of the dispersion can be clearly identified from the spectra of $A(\bm{k},\omega)$ (cf. panels (c) and (d) of Fig. 8) and a back-bending occurs. However, this definition of $k_{L}$ cannot be transferred to the $(G_{0}G_{0})G$ and $(GG)G$ approaches, where only a single peak appears in the $A(\bm{k},\omega)$ (cf. panels (b) and (e) of Fig. 8) and there is no observable back-bending in the dispersion for most couplings. As a consequence, for the latter two approaches we make use of an alternative operative definition of $k_{L}$ which is closer in spirit to the description one would adopt for a Fermi liquid, and identify $k_{L}$ as the wave vector for which the single peak in $A(\mathbf{k},\omega)$ passes through $\omega=0$. Nevertheless, this definition appears to work properly up to when $k_{L}/k_{F}\gtrsim 0.5$, because at that point a double-peak structure begins to appear in $A(\mathbf{k},\omega)$ even for the $(G_{0}G_{0})G$ and $(GG)G$ approaches. To be able to extend the $k_{L}$-vs-coupling curve up to $k_{L}\rightarrow 0$ when the full collapse of the underlying Fermi surface occurs, for the coupling regime where the double peak occurs we thus found it more appropriate to identify $k_{L}$ as the wave vector at which the peaks in $A(\mathbf{k},\omega)$ at negative and positive frequencies mutually exchange the height of their maxima.

Figure 12 reports the values of the Luttinger wave vector $k_{L}$ as a function of $(k_{F}a_{F})^{-1}$ obtained for all $t$-matrix approaches using the procedures described above. Here, depending on the approach, $k_{L}$ is seen to vanish for couplings values in the rather narrow range between $0.55$ and $0.7$, with the largest critical coupling reached by the $(GG)G$ fully self-consistent approach.

V Concluding remarks and perspectives

In this paper, we have performed a systematic theoretical study about several variants of the $t$-matrix approximation for an attractive dilute Fermi gas in its normal phase above the superfluid critical temperature $T_{c}$. This study was extended to the whole BCS-BEC crossover and has regarded both thermodynamic and dynamical quantities that characterize the Fermi gas. Although these variants of the $t$-matrix approximation have already been separately considered in the literature, the novelty here is that all these variants have been treated on equal footing (with the same numerical accuracy also having been pursued for all of them) and in an unbiased way, in order to evidence their individual virtues and shortcoming.

As far as the thermodynamic quantities that we have considered are concerned, from one $t$-matrix approach to the other we have found mostly quantitative differences but similar qualitative trends, apart from the presence vs absence of a maximum for the critical temperature in the intermediate-coupling regime (cf. Fig. 3) and of a residual bosonic interaction affecting the chemical potential in the strong-coupling (BEC) regime (cf. Fig. 4). The most distinctive differences (not only at a quantitative but also at a qualitative level) among the outcomes of the various $t$-matrix approaches have instead been found for the dynamical quantities, specifically, about the occurrence of a one-peak vs two-peak structure in the single-particle spectral function (cf. Fig. 8). This qualitative difference appears relevant, not because it hinges on a dispute about the Fermi liquid vs non-Fermi liquid behavior of an attractive Fermi gas at unitarity Salomon-2011 , but rather because it affects the width (if not the presence itself) of a temperature interval above $T_{c}$ where a pseudo-gap regime would show up, with the simultaneous presence of preformed pairs and of an underlying Fermi surface.

That a fully self-consistent diagrammatic approximation may end up in giving (even considerably) smaller values for the excitation energies with respect to its non-self-consistent version(s) (and possibly “overshoot the mark” when comparing with experimental values) has also been evidenced in other physical contexts. One can specifically refer to the spectra associated with electronic excitations that can be described in terms of the GW approximation, not only for semiconductors and insulators Kresse-2018 but also for complex molecules Rinke-2012 . In these cases, attempts have recently been made to mitigate the effects of self-consistency by introducing vertex corrections on the GW calculations Kutepov-2018 ; Berkelbach-2018 ; footnote-new .

In the present context of the $t$-matrix approximation for an attractive Fermi gas, too, vertex corrections have recently been included on top of a partially self-consistent version of this approximation, ending up with rather good results for the critical temperature $T_{c}$ Pisani-2018-I and the pairing gap at zero temperature Pisani-2018-II throughout the BCS-BEC crossover, when compared with available quantum Monte Carlo calculations and experimental data. It could therefore be interesting to assess whether, including a similar kind of vertex corrections on top of the fully self-consistent $t$-matrix approximation, may result in more favourable conditions for the presence of a pseudo-gap regime about unitarity.

**ACKNOWLEDGMENTS **

M.P. acknowledges L. Fallani for support and discussions. This work was partially supported by the Italian MIUR under Contract PRIN-2015 No. 2015C5SEJJ001.

Appendix A DETAILS OF THE NUMERICAL PROCEDURES FOR ACHIEVING

(PARTIAL OR FULL) SELF-CONSISTENCY

In this Appendix, we present in detail the numerical procedures that are needed to implement the cycles of self-consistency depicted schematically in Fig. 2. For the sake of definiteness, we will specifically consider the fully self-consistent $(GG)G$ approach whose cycle is shown in Fig. 2(d), since all procedures discussed for this approach can as well be applied to the partially self-consistent approaches. We remark that the procedures here presented for the $(GG)G$ approach are in line with those previously suggested in Ref. Haussmann-1994 and partially with those reported in Ref. VanHoucke-2013 .

All expressions reported in this Appendix are given in dimensionless units, such that energies are in units of the Fermi energy $E_{F}=k_{F}^{2}/(2m)$ and wave vectors in units of the Fermi wave vector $k_{F}$. Accordingly, the single-particle fermionic propagator $G(\mathbf{k},\omega_{n})$ is in units of $E_{F}^{-1}$, the fermionic self-energy $\Sigma(\mathbf{k},\omega_{n})$ in units of $E_{F}$, and the particle-particle propagator $\Gamma(\mathbf{Q},\Omega_{\nu})$ in units of $(mk_{F})^{-1}$. In addition, to further shorten the notation, here we use the symbol $v=(k_{F}a_{F})^{-1}$ for the coupling.

1. Transforming from $G(\mathbf{k},\omega_{n})$ to $G(\mathbf{r},\tau)$

The first function to be Fourier transformed in the self-consistent cycle of Fig. 2(d) is the single-particle fermionic propagator $G$. For this function the Fourier transform can be done in two steps, namely,

[TABLE]

with the Fourier transform over the wave vector $\mathbf{k}$ following that over the frequency $\omega_{n}$.

In the first step, to get the Fourier transform over the frequency $\omega_{n}$, we note that the dressed fermionic propagator $G(\mathbf{k},\omega_{n})$ for large frequencies behaves like the free propagator $G_{0}(\mathbf{k},\omega_{n})$:

[TABLE]

where $\xi_{\mathbf{k}}=\mathbf{k}^{2}-\mu$. One can thus calculate numerically the Fourier transform from $\omega_{n}$ to $\tau$ of the difference

[TABLE]

which is easier to obtain than the Fourier transform of $G(\mathbf{k},\omega_{n})$ since the function (27) converges like $\sim\omega_{n}^{-5/2}$ for large frequencies. One then adds to this result the Fourier transform of $G_{0}(\mathbf{k},\omega_{n})$, which is known analytically in the form FW-1971

[TABLE]

where $f(\xi_{\mathbf{k}})=(e^{\beta\xi_{\mathbf{k}}}+1)^{-1}$ is the Fermi distribution function.

In the second step, to get the Fourier transform over the wave vector $\mathbf{k}$, it is further convenient to split the expression (28) in two parts:

[TABLE]

where

[TABLE]

Here, the labels $(\mathrm{n})$ and $(\mathrm{a})$ signify that the term (30) has to be numerically Fourier transformed over $\mathbf{k}$, whereas the term (31) admits an .pdfr transform of the form:

[TABLE]

which in the limit $\tau\rightarrow 0^{+}$ is a representation of the Dirac delta function $\delta^{3}(\mathbf{r})$. This property is related to the anticommutator between fermionic field operators that appears in the fermionic propagator when passing from $\tau=0^{-}$ to $\tau=0^{+}$. The term (32) thus describes the singular behavior when $(\mathbf{r},\tau)\to 0^{+}$, not only for the free propagator $G_{0}$ but also for the dressed propagator $G$. It is then convenient to define a new function

[TABLE]

which can be readily Fourier transformed over $\mathbf{k}$ numerically. The desired Fourier transform $G(\mathbf{r},\tau)$ is eventually obtained by adding $G^{(a)}_{0}(\mathbf{r},\tau)$ of Eq. (32) to the Fourier transform of $\tilde{\Delta}G(\mathbf{k},\tau)$ footnote-2 .

2. Transforming from $\Gamma(\mathbf{Q},\Omega_{\nu})$ to $\Gamma(\mathbf{r},\tau)$

The next function to be Fourier transformed in the self-consistent cycle of Fig. 2(d) is the particle-particle propagator $\Gamma$. Also in this case, the Fourier transform can be done in two steps, namely,

[TABLE]

with the Fourier transform over the wave vector $\mathbf{Q}$ following again that over the frequency $\Omega_{\nu}$.

To perform the first Fourier transform over $\Omega_{\nu}$, we begin by noting that the large-frequency behavior of the particle-particle propagator $\Gamma(\mathbf{Q},\Omega_{\nu})$ coincides with that of its non-self-consistent counterpart taken in the strong-coupling limit $-\beta\mu\gg 1$, which is given by the expression Pieri-2000 :

[TABLE]

However, using $\Gamma_{\mathrm{sc}}(\mathbf{Q},\Omega_{\nu})$ as the reference function to be subtracted in the Fourier transform may lead to problems, because for $\Omega_{\nu}=0$ the function (35) has a pole at $|\mathbf{Q}|=2\sqrt{v^{2}+\mu}$ when $v>0$. To the extent that we are interested only in taking care of the large-frequency behavior of $\Gamma(\mathbf{Q},\Omega_{\nu})$, we are led to introduce the new reference function

[TABLE]

with the zero-frequency term removed from the expression (35). Although the Fourier transform of the function (36) cannot be calculated analytically, it can be computed with limited effort by writing it as an integral over the complex $z$-plane as follows:

[TABLE]

where the contour $\mathcal{C}$ surrounds the poles of the Bose function $b(z)=1/(e^{\beta z}-1)$ on the imaginary axis. Here, the function $\Gamma_{\mathrm{sc}}^{\prime}(\mathbf{Q},z)$ has a branch cut along the negative real axis starting at $z_{\mathrm{c}}=2(\mu-\mathbf{Q}^{2}/4)$ as well as a pole at $z_{\mathrm{p}}=2v^{2}+z_{c}$ when $v>0$. The integral in Eq. (37) then reduces to the calculation of an integral along the branch cut and of the residue of the pole, and can accordingly be split in the following way:

[TABLE]

where the first two terms are contributed by the branch cut and the third term by the pole. Like in Eq. (29), the labels $(\mathrm{n})$ and $(\mathrm{a})$ in the first line of Eq. (A) signify that these contributions are calculated numerically or analytically, respectively.

The first term in Eq. (A) can be cast in the form

[TABLE]

for $z_{\mathrm{c}}<0$ (that is, for $\mu<\mathbf{Q}^{2}/4$), and in the form

[TABLE]

for $z_{\mathrm{c}}\geq 0$ (that is, for $\mu\geq\mathbf{Q}^{2}/4$), where $z_{0}$ is any point along the negative real axis (typically, we have taken $z_{0}=-2T$). The second term in Eq. (A) has instead the semi-analytic form:

[TABLE]

with the coefficient $c(\tau,v)$ given by

[TABLE]

Note that this coefficient is unity for $\tau\to 0^{+}$ or $v=0$.

The term (41) admits also an .pdfr transform from $(\mathbf{Q},\tau)$ to $(\mathbf{r},\tau)$, which is given by:

[TABLE]

One can show that this is the leading term of the singular behavior of the full $\Gamma(\mathbf{r},\tau)$ for $(\mathbf{r},\tau)\rightarrow 0^{+}$, not only in the strong-coupling regime but also throughout the BCS-BEC crossover. Finally, the third term in Eq. (A) is the residue of the pole at $z=z_{\mathrm{p}}$, given by:

[TABLE]

One can also show that the divergence of this term for $z_{\mathrm{p}}=0$ (that is, for $|\mathbf{Q}|=2\sqrt{v^{2}+\mu}$) is exactly compensated by the fourth term of (A), in such a way that $\Gamma^{\prime}_{\mathrm{sc}}(\mathbf{Q},\tau)$ is always a smooth function of $\mathbf{Q}$.

At this point, we calculate numerically the Fourier transform from $(\mathbf{Q},\Omega_{\nu})$ to $(\mathbf{Q},\tau)$ of the difference

[TABLE]

to obtain $\Delta\Gamma(\mathbf{Q},\tau)$, and then add to it $\Gamma^{\prime}_{\mathrm{sc}}(\mathbf{Q},\tau)$ given by Eq. (A) to obtain $\Gamma(\mathbf{Q},\tau)$.

Finally, for the remaining Fourier transform from $(\mathbf{Q},\tau)$ to $(\mathbf{r},\tau)$, we can make use of $\Gamma_{\mathrm{sc}}^{(\mathrm{a})}(\mathbf{Q},\tau)$ defined by Eqs. (41) and (42) to subtract the leading term of the singular behavior when $(\mathbf{r},\tau)\to 0^{+}$, similarly to what we did in Eq. (33). Accordingly, we define the difference function:

[TABLE]

and Fourier transform it from $\mathbf{Q}$ to $\mathbf{r}$ to obtain $\tilde{\Delta}\Gamma(\mathbf{r},\tau)$. The desired function $\Gamma(\mathbf{r},\tau)$ eventually results by adding to $\tilde{\Delta}\Gamma(\mathbf{r},\tau)$ the semi-analytic expression (43) for $\Gamma_{\mathrm{sc}}^{(\mathrm{a})}(\mathbf{r},\tau)$ footnote-3 .

3. Transforming from $\Sigma(\mathbf{r},\tau)$ to $\Sigma(\mathbf{k},\omega_{n})$

The last function to be Fourier transformed in the self-consistent cycle of Fig. 2(d) is the self-energy $\Sigma$. Also in this case, the Fourier transform is done in two steps, namely,

[TABLE]

where now one first transforms from $\mathbf{r}$ to $\mathbf{k}$ and then from $\tau$ to $\omega_{n}$, in the reversed direction to what was done for $G$ (cf. Eq. (25)) and for $\Gamma$ (cf. Eq. (34)). This is because from Eq. (7), once $G(\mathbf{r},\tau)$ and $\Gamma(\mathbf{r},\tau)$ are know, $\Sigma(\mathbf{r},\tau)$ is also known.

From Eq. (7) one can also determine the singular behavior of $\Sigma(\mathbf{r},\tau)$ for $(\mathbf{r},\tau)\rightarrow 0^{+}$ in terms of those of $G(\mathbf{r},\tau)$ [cf. the discussion after Eq. (32)] and of $\Gamma(\mathbf{r},\tau)$ [cf. the discussion after Eq. (43)]. To this end, we rewrite Eq. (7) in the alternative forms:

[TABLE]

(where the factor of two originates by having expressed the particle-particle propagator $\Gamma$ in terms of dimensionless quantities). In the second line we have used the spatial isotropy and the anti-periodicity in $\tau$ of the fermionic propagator, to write $G(-\mathbf{r},-\tau)=-G(\mathbf{r},\beta-\tau)$. As discussed previously, both $G(\mathbf{r},\tau)$ and $\Gamma(\mathbf{r},\tau)$ are strongly peaked for $(\mathbf{r},\tau)\rightarrow 0^{+}$. We then expect the singular behavior of $\Sigma(\mathbf{r},\tau)$ to be captured by the following alternative expressions:

[TABLE]

in the limit $(\mathbf{r},\tau)\rightarrow(0,0^{+})$ and

[TABLE]

in the limit $(\mathbf{r},\tau)\rightarrow(0,\beta^{-})$. Out of the above terms, (50) is the dominant one because $\Gamma(\mathbf{r},\tau)$ [cf. Eq. (43)] is more strongly peaked than $G(\mathbf{r},\tau)$ [cf. Eq. (32)] in the limit $\tau\rightarrow 0^{+}$.

Accordingly, in order to perform the Fourier transform of $\Sigma(\mathbf{r},\tau)$ from $\mathbf{r}$ to $\mathbf{k}$, we consider the difference

[TABLE]

with the term (50) only. In addition, to the extent that in Eq. (52) we are interested in the leading behavior of $\Sigma^{(+)}$ for $(\mathbf{r},\tau)\rightarrow(0,0^{+})$, we can approximate the particle-particle propagator $\Gamma(\mathbf{r},\tau)$ in Eq. (50) by the analytic result (43) and write

[TABLE]

where we have consistently set $c(\tau,v)\to 1$ and $e^{\mu\tau}\rightarrow 1$ and used the result $n=-2\,G(\mathbf{r}=0,\beta^{-})$ in terms of the fermionic density $n$ (in dimensionless units). [For the $t$-matrix approaches that use an external $G_{0}$ in the place of $G$, the expression (53) needs to be modified by replacing $n$ by the free-fermion density $n_{0}=-2\,G_{0}(\mathbf{r}=0,\beta^{-})$.]

The expression (53) has the further advantage that its Fourier transform from $\tau$ to $\omega_{n}$ can be obtained analytically in terms of the error function AS-1972 , in the form:

[TABLE]

This expression correctly reproduces the leading $\omega_{n}^{-1/2}$ behavior of $\Sigma(\mathbf{k},\omega_{n})$ for large frequencies. [Similarly, $\Sigma^{(-)}$ of Eq. (51) can be shown to behave like $\sim\omega_{n}^{-1}$ for large frequencies VanHoucke-2013 , thereby confirming that this is a sub-leading contribution also in the frequency domain.]

Once the difference (52) has been Fourier transformed numerically (first from $\mathbf{r}$ to $\mathbf{k}$ and then from $\tau$ to $\omega_{n}$) according to the above prescriptions to obtain $\Delta\Sigma(\mathbf{k},\omega_{n})$, one can eventually add to it the analytic expression (54) and obtain the desired function $\Sigma(\mathbf{k},\omega_{n})$.

Appendix B OPTIMIZING THE CONVERGENCE TOWARDS SELF-CONSISTENCY

In this Appendix, we discuss the procedures that we have adopted to achieve optimal convergence toward self-consistency within the $t$-matrix approaches considered in this paper. A judicious optimization in achieving convergence is, in fact, especially required for those approaches that implement partial or total self-consistency in the particle-particle propagator $\Gamma$, to the extent that these approaches become intrinsically unstable when the temperature approaches $T_{c}$ when one uses the straightforward iterative procedure sketched in Section II (cf. Fig. 2 therein). The origin of this problem can be highlighted through the following analytic considerations.

1. General considerations on the iterative procedure

For definiteness, we shall consider in detail the fully self-consistent $(GG)G$ approach. The equations (1)-(4), that need to be self-consistently solved, can be written in a compact way as a functional equation for the self-energy $\Sigma$, in the form:

[TABLE]

This is because, once $\Sigma$ is know, $G$ and $\Gamma$ in Eqs. (1)-(4) can also be readily obtained. Let us then see what happens when trying to solve Eq. (55) iteratively.

Suppose that \Sigma^{\mathrm{SC}}(k)=\mathcal{F}\big{[}\Sigma^{\mathrm{SC}}(p)\big{]}\!(k) is the self-consistent ($\mathrm{SC}$) solution to Eq. (55) which the iterative method is expected to reach. At a generic iteration step $\mathrm{(i)}$ toward self-consistency, the self-energy $\Sigma^{\mathrm{(i)}}(k)$ deviates from $\Sigma^{\mathrm{SC}}(k)$ by some amount $\delta\Sigma^{\mathrm{(i)}}(k)$:

[TABLE]

Provided one is close enough to the self-consistent solution, Eq. (55) can be linearized about $\Sigma^{\mathrm{SC}}$, to write:

[TABLE]

Here, the subscript $\mathrm{SC}$ indicates that the functional derivative is taken at self-consistency, and the integral over $p$ contains both an integral over the wave vector $\mathbf{p}$ and a sum over the Matsubara frequency $\omega_{n}$. This provides a relation between the distance $\delta\Sigma$ from the self-consistent solution between the steps $\mathrm{(i-1)}$ and $\mathrm{(i)}$:

[TABLE]

When $\Sigma(k)$ is calculated on a $k$-grid of points, like it is done in our numerical calculations, $\delta\Sigma$ can be regarded as a vector which is acted upon by the functional derivative matrix $[\delta\mathcal{F}/\delta\Sigma]_{\mathrm{SC}}$. The convergence of the iterative procedure, from $\mathrm{(i-1)}$ to $\mathrm{(i)}$ and so on, is then governed by the behavior of this functional derivative matrix.

To find an explicit expression for $[\delta\mathcal{F}/\delta\Sigma]_{\mathrm{SC}}$, we rewrite it in the form:

[TABLE]

Here, the factor on the right is given by

[TABLE]

with the use of Eq. (1), while the factor on the left can be calculated by recalling that (cf. Eqs. (2) and (55))

[TABLE]

yielding

[TABLE]

In this expression, the functional derivative of the particle-particle propagator $\Gamma(Q)$ can be obtained from Eq. (3)

[TABLE]

where the functional derivative of the regularized particle-particle bubble $R_{\mathrm{pp}}(Q)$ is obtained from Eq. (4)

[TABLE]

Grouping all the above results together, Eq. (59) becomes eventually:

[TABLE]

When this result is used in Eq. (58), $\delta\Sigma^{(i)}$ therein can be conveniently split in two terms

[TABLE]

where we have defined

[TABLE]

Suppose now that $T=T_{c}$. The Thouless criterion (10) then implies that the particle-particle propagator $\Gamma(Q)$ has a pole for $Q=0$, such that one expects $\Gamma(Q)$ to remain strongly peaked in the vicinity of $Q=0$. The expressions (67) and (68) can be simplified accordingly, by setting to zero the arguments of the particle-particle propagators in the smooth functions that multiply them. For the term (67) we thus have:

[TABLE]

where $C$ is the Tan’s contact according to Eq. (17). This term poses no problem to the convergence, since the quantity that multiplies $\delta\Sigma^{\mathrm{(i-1)}}(-k)$ is finite. For the term (68) we instead obtain:

[TABLE]

where the factor

[TABLE]

is infrared divergent at $T=T_{c}$ even in three dimensions. This is because, for the term with $\Omega_{\nu}=0$ therein, $\Gamma(\mathbf{Q},\Omega_{\nu}=0)\sim\mathbf{Q}^{-2}$ when $\mathbf{Q}\rightarrow 0$ at $T=T_{c}$ PPSC-2002 .

This divergence represents a problem for the convergence of the iterative algorithm, because it implies that, no matter how close the step $\mathrm{(i-1)}$ might be to the self-consistent solution, the step $\mathrm{(i)}$ is bound to run infinitely away from it. This problem can affect the convergence of the iterative algorithm also for temperatures $T\gtrsim T_{c}$, depending on how much $\Gamma(Q)$ is peaked about $Q=0$.

The only (partially) self-consistent $t$-matrix approach not affected by this problem is the $(G_{0}G_{0})G$ approach. This is because in this approach, by construction, the particle-particle propagator $\Gamma$ coincides with the bare $\Gamma_{0}$, whereby $\delta\Gamma_{0}/\delta G=0$. As a consequence, the second term on the right-hand side of Eq. (62) identically vanishes and with it the divergent factor in Eq. (70).

We pass now to show how this problem can be overcome in practice, both for $T>T_{c}$ and $T\rightarrow T_{c}$.

2. Improved method for $T>T_{c}$

Although the iterative approach (58) to solve Eq. (55) cannot converge at exactly $T=T_{c}$, there exists a simple method to make it converge for $T\gtrsim T_{c}$, with the factor (71) keeping a finite (albeit large) value. This method, which has already been used for the self-consistent calculations of electronic structures in atoms and molecules Ferreira-1980 , consists in redefining the iterative steps in terms of the weighted sum:

[TABLE]

where the weight factor $\alpha$ ranges between [math] and $1$. This method is found to reduce the effects of the divergence due to the term (70), thereby making the iterative process to converge even sufficiently close to $T_{c}$. Nevertheless, the method fails upon approaching $T_{c}$, because smaller and smaller values of $\alpha$ are needed for attaining convergence. A smaller value of $\alpha$, in turn, implies that more iterative steps are required for convergence, since the contribution of a single step, too, becomes smaller and smaller. In practice, we have found that this method can conveniently been used down to temperatures for which $(T-T_{c})/T_{c}$ is of order $1\%$, before it becomes numerically too demanding.

3. Improved method for $T=T_{c}$

Alternatively, exactly at $T=T_{c}$ one can rely on a different method that avoids the convergence problem discussed above. This method can be summarized as follows:

(i) One begins by fixing the value $n$ of the density, in terms of which one obtains the Fermi wave vector $k_{F}=(3\pi^{2}n)^{1/3}$ and the Fermi energy $E_{F}=k_{F}^{2}/(2m)$. One also fixes a guess value $T_{\mathrm{g}}/T_{F}$ for the temperature in units of the Fermi energy, as well of the ratio $\mu/T_{\mathrm{g}}$ between the chemical potential and the guess temperature $T_{\mathrm{g}}$ (recall that we have set the Boltzmann constant $k_{B}$ equal to unity throughout).

(ii) Next, one replaces the particle-particle propagator of Eq. (3) with the following expression

[TABLE]

in such a way that $\tilde{\Gamma}(Q=0)^{-1}=0$ by construction. This implies that the Thouless criterion (10) is always satisfied by the modified particle-particle propagator (73), no matter what was the initial guess for $T_{\mathrm{g}}/E_{F}$.

(iii) At this point one can proceed and perform the iterative procedure toward self-consistency on the set of equations (1)-(4), with $\Gamma(Q)$ therein replaced by $\tilde{\Gamma}(Q)$. It will be shown below that this replacement avoids the occurrence of the infrared divergence that plagues instead the expression (70) obtained in terms of the original $\Gamma(Q)$.

(iv) Once self-consistency has been achieved with this modified set of equations, one can obtain the modified density as $\tilde{n}=-2\,\tilde{G}(\mathbf{r}=0,\beta^{-})$, and the modified scattering length $\tilde{a}_{F}$ from the expression

[TABLE]

where $\tilde{R}_{\mathrm{pp}}$ is obtained from Eq. (4) with $\tilde{G}$ replacing $G$. The result (74) follows directly from the Thouless criterion corresponding to the modified density $\tilde{n}$.

(v) Finally, in terms of $\tilde{n}$ one obtains the modified Fermi wave vector $\tilde{k}_{F}=(3\pi^{2}\tilde{n})^{1/3}$ and the modified Fermi energy $\tilde{E}_{F}=\tilde{k}_{F}^{2}/(2m)$. The desired value of the critical temperature is then obtained by

[TABLE]

in terms of the initial guess $T_{\mathrm{g}}/E_{F}$. From Eq. (74) the corresponding coupling value is given by:

[TABLE]

There remains to explain the reason why this method is not plagued by the convergence problem discussed above for the iterative procedure. The point is that, using the modified particle-particle propagator (73) in the place of the original one, one also modifies the structure of the functional derivative in Eq. (58). Specifically, for the functional derivative of $\tilde{\Gamma}(Q)$ with respect to $\tilde{G}(p)$ one obtains (cf. Eq. (63)):

[TABLE]

The corresponding variation of the self-energy related to this functional derivative then becomes (cf. Eq. (68)):

[TABLE]

Comparing this result with Eq. (68), one notes that the singular behavior of $\tilde{\Gamma}(Q)^{2}$ for $Q\rightarrow 0$ is now suppressed by the presence of the factor $[\tilde{G}(Q-p)-\tilde{G}(-p)]$. This feature makes it possible to reach convergence exactly at $T=T_{c}$ with a limited number of iterations, without the need for the weighted sum of Eq. (72).

Appendix C COMPARISON WITH THE PSEUDO-GAP APPROXIMATION AT $T_{c}$

The $(GG_{0})G_{0}$ and $(GG)G_{0}$ $t$-matrix approaches have almost invariably been implemented in the literature using a set of approximations (sometimes referred to as the “pseudo-gap approximation”), which considerably simplify the numerical calculations. Specifically, close to $T_{c}$ where the particle-particle propagator $\Gamma(Q)$ is strongly peaked about $Q=0$, the fermionic self-energy (2) has been approximated as follows Maly-1999 ; Chen-2005 :

[TABLE]

Due to this approximation for the self-energy, the dressed single-particle propagator $G$ (and thus also the equation for the particle number) coincides with that of BCS theory, with the pseudo-gap energy $\Delta_{\rm pg}$ now playing the role in the normal phase of the BCS gap $\Delta$ in the superfluid phase. In addition, the Thouless criterion (10) for the dressed $\Gamma$ coincides with the BCS gap equation (again with the replacement $\Delta\rightarrow\Delta_{\rm pg}$) for the $(GG_{0})G_{0}$ approach Maly-1999 ; Chen-2005 , or with the BCS gap equation plus additional corrections (which become anyway negligible both in the BCS and BEC limits) for the $(GG)G_{0}$ approach Micnas-2014 . A further approximation, which is usually adopted within the pseudo-gap approximation when calculating $\Delta_{\rm pg}$ by means of Eq. (79), is the use of an expansion of $\Gamma({\bf Q},\Omega_{\nu})$ for small values of ${\bf Q}$ and $\Omega_{\nu}$.

It should, however, be remarked that the pseudo-gap approximation (79) appears justified only in the strong-coupling (BEC) regime $(k_{F}a_{F})^{-1}\gtrsim+1$, where the large fermionic energy scale $|\mu|$ (with $\mu<0$) dominates over the bosonic energy scales and the approximation (79) becomes fully correct. This can also be verified numerically as shown in Fig. 13, where a comparison is presented for the calculation of $T_{c}$ between the complete $(GG_{0})G_{0}$ and $(GG)G_{0}$ approaches (full lines) and their corresponding pseudo-gap approximations (dashed lines). This comparison shows that a good agreement between the complete and approximate calculations occurs only for $(k_{F}a_{F})^{-1}\gtrsim 1$, while significant deviations result both at intermediate and weak couplings.

From Fig. 13 one also notices that, in the weak-coupling (BCS) regime $(k_{F}a_{F})^{-1}\lesssim-1$, the curves for $T_{c}$ obtained within the pseudo-gap approximation converge rapidly to the corresponding BCS curve for $T_{c}$. This is because $\Delta_{\rm pg}(T)$ is bounded by the value $\Delta_{0}$ of the BCS gap at $T=0$ Chen-2005 , which in turn vanishes exponentially in the weak-coupling limit. This implies that the approximate self-energy (79), too, vanishes exponentially, in such a way that the BCS result for $T_{c}$ is recovered. Without the use of the approximation (79), in weak coupling the self-energy would instead approach the value $\Sigma\simeq 2\pi a_{F}n/m$ associated with a mean-field shift.

For the complete $(GG_{0})G_{0}$ approach, whereby this shift appears in just one of the two single-particle propagators that enter the particle-particle bubble in $\Gamma$, the equation for $T_{c}$ corresponds to that of a Fermi system in the presence of a chemical potential imbalance $\delta\mu$. This equation is known to have no solution when this imbalance about exceeds the BCS gap $\Delta_{0}$ of the balanced system at $T=0$ Sarma-1963 . Given the exponential dependence of $\Delta_{0}$ on coupling, to be contrasted with the linear dependence $\delta\mu\simeq 2\pi a_{F}n/m$ associated with the mean-field shift, the condition $\delta\mu>\Delta_{0}$ is readily met in the weak-coupling regime. This explains why, in the weak-coupling regime, the $(GG_{0})G_{0}$ approach does not admit solution for $T_{c}$, as it was already noted in the discussion of Fig. 3. For the complete $(GG)G_{0}$ approach, on the other hand, the differences with respect to its pseudo-gap approximation are overall less pronounced, albeit still significant.