BCS-BEC Crossover Effects and Pseudogap in Neutron Matter

Due to the large neutron-neutron scattering length, dilute neutron matter resembles the unitary Fermi gas, which lies half-way in the crossover from the BCS phase of weakly coupled Cooper pairs to the Bose-Einstein condensate of dimers. We discuss crossover effects in analogy with the T-matrix theory used in the physics of ultracold atoms, which we generalize to the case of a non-separable finite-range interaction. A problem of the standard Nozieres-Schmitt-Rink approach and different ways to solve it are discussed. It is shown that in the strong-coupling regime, the spectral function exhibits a pseudo-gap at temperatures above the critical temperature T_c. The effect of the correlated density on the density dependence of T_c is found to be rather weak, but a possibly important effect due to the reduced quasiparticle weight is identified.


I. INTRODUCTION
The inner crust of neutron stars is characterized by the presence of a dilute gas of neutrons in between much denser clusters [1]. The superfluidity of this gas is of decisive importance, e.g., to explain pulsar glitches. Nevertheless, in spite of a long effort and the existence of Quantum Monte-Carlo (QMC) calculations of the gap at low densities [2,3], there remain large uncertainties about the density dependence of the pairing gap and the critical temperature of neutron matter.
In reality, the neutron gas in the inner crust is not uniform because of the nuclear clusters. However, if the pairing is strong enough and the Cooper-pair size (coherence length) small enough, i.e., much smaller than the distance between clusters, it should be possible to treat the neutron gas approximately as uniform matter, except near the cluster surface [4]. Furthermore, in the cases where the local-density approximation is not valid, the calculations rely generally on density-dependent effective interactions, fitted such that they reproduce, at the mean-field level, the realistic pairing gaps in uniform matter [4,5]. So, also in this case, the understanding of uniform neutron matter is a necessary prerequisite to describe superfluidity of the inner crust.
The neutron-neutron (nn) interaction has a very large scattering length a −18.7±0.6 fm [6] or −16.59±1.17 fm [7]) while its effective range is only r 0 = 2.83 ± 0.11 fm [7]. Hence, at low densities, when the mean distance between the neutrons is between 3 and 16 fm, the neutron gas resembles a unitary Fermi gas, an idealized system of spin-1/2 Fermions with a contact interaction having infinite scattering length, which was originally introduced as a schematic model of dilute neutron matter [8].
A few years later, the unitary Fermi gas was experimentally realized with ultracold trapped gases of alkali atoms, where the interaction has practically zero range (four orders of magnitude smaller than the interparticle distance) and the scattering length can be tuned with the help of a magnetic field from negative to positive values by passing through infinity at the so-called Feshbach resonance. This allows one to study not only the unitary limit, but the whole crossover from BCS pairing with Cooper pairs in the case of weakly attractive interactions (a < 0) to the Bose-Einstein condensation (BEC) of bound dimers in the case a > 0 [9]. The unitary limit is just a special case in this crossover, corresponding to the situation where the binding energy of the dimer in free space tends to zero.
On the BEC side and in the crossover region, it is crucial for a correct description of the superfluid phase transition that pair correlations persist at temperatures above the critical temperature T c . On the BEC side of the crossover, these correlations can easily be interpreted as dimers which are already formed but not yet condensed. This implies in particular that the BCS relation T c = 0.57 ∆(0) between the critical temperature T c and the zero-temperature pairing gap ∆(0) is not fulfilled in this regime. The simplest theory that correctly interpolates between the BCS and the BEC limits was proposed by Nozières and Schmitt-Rink (NSR) [10].
In nuclear systems, the interaction is of course fixed. Nevertheless, a BCS-BEC crossover is expected to happen in symmetric nuclear matter if one varies the density: at very low density, there will be a BEC of deuterons, which continuously goes over to a BCS superfluid with proton-neutron Cooper pairs at higher density [11]. The critical temperature T c in this crossover was studied in [12][13][14] using an approach similar to the NSR theory. The case of neutron matter is somewhat different, because no bound dineutron state and hence no BEC phase exist. Nevertheless, with varying density, neutron matter passes from the strongly coupled regime close to the unitary limit to the weakly coupled BCS regime. The relevance of BCS-BEC crossover-like physics in dilute neutron matter was pointed out in [15]. Consequently, NSR-like corrections of T c due to pair correlations in the normal phase were studied in [16][17][18][19].
In previous work following the NSR and related approaches, the single-particle self-energy and spectral function were usually not computed. One reason is that in these approaches the density ρ is often obtained as the derivative of the thermodynamic potential with respect to the chemical potential µ [10,19,20]. Computing the thermodynamic potential in the ladder approximation, one automatically obtains an expression which is equivalent to truncating the Dyson equation for the single-particle Green's function at the first order in the self-energy [9]. In this case, when calculating only the density ρ and not the occupation numbers n(k), the somewhat difficult computation of the selfenergy can be avoided. However, some problems of the different approaches show up when one looks at the occupation numbers n(k), and some interesting aspects of the crossover regime, such as the existence of a pseudogap, can only be studied by looking at the spectral function.
In the present work, we revisit this problem, but now we compute the self-energy, which gives us access to the occupation numbers and to the spectral function. In Sec. II, we recall the T matrix formalism and discuss in some detail the way we deal with non-separable interactions. In Sec. III, we will show our results for the spectral function, the pseudogap, the correlated occupation numbers, and the density dependence of the critical temperature. We will in particular discuss the differences between the original NSR approach, its modifications that have been used previously, and our more complete treatment of the self-energy. Some open questions are discussed in Sec. IV, and we conclude in Sec. V.

A. Vertex function
The superfluid phase transition is an instability of the normal phase towards the superfluid one which shows up in the two-particle T matrix (vertex function) Γ. We compute Γ within the ladder approximation, i.e., by resumming ladder diagrams in the medium. In this formalism, the vertex function for total momentum k, in-and outgoing relative momenta q and q , and total energy of the pair ω, is the solution of the following Lippmann-Schwinger like equation corresponding to the Feynman diagrams shown in Fig. 1(a), where v(q, q ) is the matrix element of the s-wave nn interaction. In our calculations, we will use interactions of the V low-k type [21]. At finite temperature, the two-particle Green's function G (2) 0 is written within the Matsubara formalism [22] as where ξ k = 2 k 2 /(2m) − µ is the neutron single-particle energy relative to the chemical potential µ, with the reduced Planck constant and m the neutron mass, β = 1/T is the inverse of the temperature T , Ω N = 2πN T and ω n = (2n + 1)πT are bosonic and fermionic Matsubara frequencies, respectively, and G 0 (k, iω n ) = 1/(iω n − ξ k ) is the free single-particle Green's function. After summation over ω n and analytic continuation to the real ω axis, we obtain the retarded two-particle Green's function where f F (ξ) = 1/(e βξ + 1) denotes the Fermi function. The limit ε → 0 + is implicitly understood. To simplify the writing, we set = m = 1 in the following formulas. In practice, this amounts to measuring energies in units of fm −2 instead of MeV with the conversion factor 1 fm −2 = ( 2 /m) fm −2 = 41.44 MeV. In the present case that ξ k is quadratic in k, the angular integral in Eq. (1) can be easily written with the angle-averaged two-particle Green's functionḠ where θ is the angle between k and p and the angle-averaged Pauli blocking factor is given bȳ In this way, Eq. (1) reduces to a one-dimensional integral equation where denotes the on-shell momentum in the center-of-mass frame. The subscripts l = 0 in Eq. (6) indicate that we use conventions that are common if one works in a partial wave basis, i.e., v = 4πv l=0 and Γ = 4πΓ l=0 .

B. Numerical solution for a non-separable interaction
In order to solve numerically Eq. (6), we will use Weinberg's eigenvalues [16,23]. If the integral equation is schematically written as Γ = V + V GΓ, the idea is to diagonalize the operator V G, i.e., to find the eigenvectors u and eigenvalues η such that V Gu = ηu. Explicitly, one has to solve for each k and ω 2 π pmax 0 dp p 2 v l=0 (q, p)Q (k, p) p 2 0 − p 2 + iε u n (p, k, ω) = η n (k, ω)u n (q, k, ω) .
Here, we have assumed that the matrix elements v(q, q ) fall off fast enough so that in practice the integral can be cut off at some momentum p max . To solve Eq. (8), we discretize the interval from 0 to p max into a grid of n q points q i . Assuming that the eigenfunctions u n can be interpolated between these points with a cubic spline, we can write them as where ϕ i (q) is the piecewise cubic basis function for the B-spline representation [24] of u n , i.e., it is a "hat function" having its maximum at q i and vanishing outside the interval (q i−2 , q i+2 ). The advantage of this method is that, although the kernel of the integral equation may be strongly peaked near the Fermi surface (especially at low temperature), it is not necessary to use a very fine momentum grid q i , because the shape of u(q, k, ω) follows the smooth q dependence of v l=0 (q, p). By injecting the expression (9) into Eq. (8) and taking the values of q on the grid points q i , we get Introducing the matrices we can write Eq. (10) in matrix notation as M c n = η n Ac n or finally which is just an ordinary matrix diagonalization problem. In practice, the imaginary part of the matrix M can be calculated analytically, and the real part is obtained as a principal-value integral.
Let us now express also the vertex function Γ in the basis of the hat functions: Inserting this into Eq. (6) and placing the points q and q on the grid points q i and q j , we get We recognize the matrices M and A defined in Eq. (11). We also define the matrix V ij = v l=0 (q i , q j ), so that Eq. (15) can be written in matrix notation as ACA = V + M CA and therefore with I the identity matrix andṼ = A −1 V (A ) −1 the coefficients one needs to express v l=0 (q, q ) in terms of hat functions analogously to Eq. (14). Now we make use of Weinberg's eigenvalues η n and eigenvectors c n which we determined previously. This allows us to write the matrix A −1 M in the form A −1 M = P DP −1 , with D = diag(η 1 , η 2 , . . . , η nq ) and P = (c 1 , c 2 , . . . , c nq ). For given values of k and ω, the coefficients C ij of the vertex function are given by The vertex function allows us to determine the critical temperature T c . The Thouless criterion [25] tells us that at the critical temperature, a pole appears in the T matrix. This amounts to solving the equation: However, we can also determine T c directly from the Weinberg eigenvalues as explained in [16]. When we look at the coefficients given by Eq. (17), we notice that the pole in the vertex function appears if at least one of the eigenvalues is equal to one. More explicitly, the critical temperature is reached when One can easily see that this is equivalent to writing the BCS gap equation in the limit of a vanishing gap, when the momentum dependent gap becomes proportional to the eigenvector ∆ q ∝ u n (q, 0, 0).

C. Vertex function with a separable interaction
Solving Eq. (6) is much simpler if we consider the case of a separable interaction, i.e., where F is a form factor of the interaction and g the coupling constant. In this case, the vertex function can be written as  with where p 0 is defined in Eq. (7). The V low-k interaction with a cutoff of Λ = 2 fm −1 is reasonably well reproduced with a Gaussian form factor The two parameters to be adjusted are then the coupling constant g and the width of the Gaussian q 0 . One way to obtain these parameters, given in Table I, is to fit the critical temperature curve obtained with the V low-k interaction in the BCS approximation [26], as shown in Fig. 2. In view of the quality of the fit and the gain in simplicity of the calculations, we will perform some of our calculations with this separable interaction instead of the V low-k interaction.
With the parameters given in Table I, one finds a scattering length a = −15.9 fm which is a bit smaller than the real nn scattering length [6,7] but still very large compared to the range of the nn interaction and therefore close to the unitary limit a → ∞. Although the unitary limit is not realized in nature, it is interesting to consider it as a test case for the theory. This can be done by slightly readjusting the coupling constant, see Table I.

D. Self-energy
The Feynman diagram for the self-energy in ladder approximation is represented in Fig. 1(b). The corresponding expression reads [16] Since here in-and outgoing momenta in Γ are equal, we use the notation Γ(q, k, ω) for Γ(q, q, k, ω). The self-energy can be split into two contributions Σ(k, ω) = Σ HF (k) + Σ E (k, ω), where Σ HF is the energy-independent Hartree-Fock (HF) term generated by the first-order (Born) term in Γ, and Σ E is the remaining energy dependent term. After analytic continuation to real energies, the imaginary part of the retarded self-energy can be written as where f B (ξ) = 1/(e βξ − 1) denotes the Bose function. The real part can then be computed as where P denotes the principal value.

E. Occupation numbers
The BCS-BEC crossover (strong-coupling) regime is characterized by the existence of pairing correlations even above the superfluid critical temperature. These modify the occupation numbers which are therefore no longer given by those of an ideal Fermi gas, n 0 (k) = f F (ξ k ). The presence of non-condensed pairs can considerably modify the relationship between the chemical potential µ and the density and thus the density dependence of the critical temperature T c .
The first way to include this effect is the Nozières-Schmitt-Rink (NSR) approach [10,20]. Although originally formulated differently, it amounts to computing the density from a Green's function which is dressed with only one self-energy insertion [9], i.e., Inserting for Σ the spectral representation (right-hand side of Eq. (27) without the principal value and with ω replaced by iω n ), and performing the summations over ω n , one obtains the NSR occupation numbers n NSR (k) = n 0 (k)+n c NSR (k), with By f F (ξ) we denote the derivative df F (ξ)/dξ. However, the occupation numbers obtained within the NSR approach are not satisfactory (see for example Fig. 5 (b) below). The real part of the self-energy produces a shift of the quasiparticle energy which the NSR theory improperly counts as part of the correlation correction. For example, if we consider a constant self-energy (of type Σ HF ), Eq. (29) (of which only the first term remains) is only a perturbative expansion of the occupation numbers n(k) = f F (ξ k ) of an uncorrelated gas of quasiparticles of energy ξ k = ξ k + Σ HF (k).
This problem has been solved for the case of symmetric nuclear matter in Refs. [12][13][14] following an approach presented initially for an electronic system in [27]. In simple words, we must eliminate the shift of the quasiparticle energy from the calculation of the correlations. Within the formalism of self-consistent Green's functions [28], we should compute the T matrix Γ and the self-energy Σ with the dressed Green's function G = 1/(ω − ξ k − Σ). The quasiparticle approximation consists in approximating in Γ and Σ this Green's function by where Σ should be computed self-consistently with G 0 . Thus, the dressed Green's function can be written as which shows that one has to subtract Re Σ(k, ξ k ) from the self-energy. In practice, as our self-energy contains only the s-wave part of the total interaction, we do not have the means to calculate ξ k . We could for example use an energydensity functional of Skyrme or Gogny type to calculate the effective mass and the shift of the chemical potential. In neutron matter, the effective mass is close to the free mass [29]. The rest can be absorbed in a shift in chemical potential which does not impact our results because we do not calculate thermodynamic quantities. For this reason, we directly use ξ k = ξ k . A first approximation presented in [16] consists in subtracting only the HF part from the self-energy. The HF selfenergy is expected to be the dominant term, and this avoids the computation of the complete self-energy. Furthermore, analogous to the approximation (28) used in NSR theory, the dressed Green's function (31) is again truncated to first order. Thus, we have the following expression for the dressed Green's function: Analogously to the NSR case, we get now n s- However, in contrast to Ref. [16], we do compute the self-energy and therefore nothing prevents us from subtracting the total self-energy instead of only the HF part. In this case, we write which leads to the occupation numbers n s-tot (k) = n 0 (k) + n c s-tot (k), with [14] n In the three methods above, the correlated occupation numbers are calculated from the first-order truncated Dyson equation. But we can also consider the complete Dyson equation. In this case, we compute the spectral function as the imaginary part of the Green's function which allows us to calculate the occupation numbers n D with the complete Dyson equation

III. NUMERICAL RESULTS
A. Critical temperature as a function of the chemical potential As pointed out in the end of Sec. II B, the critical temperature T c determined via Eq. (19) as a function of the chemical potential µ is the same as in the BCS approximation. The difference appears only when one computes T c as a function of the density n, since the correlations change the relationship between µ and n.
In Fig. 2, we display T c as a function of µ, computed with different V low-k interactions. Matrix elements of the V low-k interaction are available for arbitrary values of the cutoff Λ. In our calculations, we use matrix elements of Ref. [17]. For each cutoff, the V low-k interaction is constructed such that it reproduces the low-energy nn scattering data, but of course only at momenta below that cutoff. Thus the lower cutoffs make only sense for small values of the chemical potential where the cutoff dependence does not yet set in. For Λ ≥ 2 fm −1 , the T c vs. µ curve remains unchanged over the whole range of chemical potentials and shows the usual behaviour. For the smaller cutoffs, however, we see that the chemical potential must be limited to µ < 5 MeV for Λ = 1 fm −1 , µ < 10 MeV for Λ = 1.5 fm −1 , and so on, which can be summarized as µ[MeV] 4 (Λ[fm −1 ]) 2 . We also display the result obtained with the separable interaction. For µ 35 MeV, it reproduces well the V low-k results obtained with Λ = 2 fm −1 .

B. Spectral function and pseudogap
Below the critical temperature T c , the pairing gap corresponds to an energy interval around the Fermi energy (ω = 0) where there are no states, because adding or removing a particle requires the energy to break a pair. Above In figure (a), we can clearly see the position of the peak located at ω = k 2 /2 − µ and a double-peak structure for k kF. Figure  (b) is a zoom on the region with the double peak. T c , such a gap does not exist, but in some cases one observes an energy region of reduced density of states, which is called the pseudogap. The spectral function A(k, ω), Eq. (36), can give give us information on this pseudogap. Figure 3 displays the spectral function computed for µ = 2 MeV at T = 0.53 MeV which is slightly above T c . We can clearly see a double peak for momenta below and around k F , as well as the significant level repulsion between the two peaks near k = k F . Although we are in the normal phase, this behavior reminds the quasiparticle spectrum in BCS theory in the superfluid phase below T c , except that the peaks have a width and the spectral function does not vanish between the peaks.
To quantify the presence of the pseudogap, we define the level density Figure 4 shows the level density calculated for different values of the chemical potential at the corresponding critical temperatures. In all cases, we see a dip in the level density around ω = 0. For µ = 20 MeV, however, this reduction is very weak, because we are no longer in the strong-coupling regime. For µ = 1 MeV the dip looks very small, but one has to remember that also the relevant energy scale given by µ is much smaller. While, to our knowledge, the pseudogap has not yet been discussed for the case of pure neutron matter, we notice that it was theoretically predicted for the case of dilute symmetric nuclear matter, which undergoes a true crossover from a BEC of deuterons to a BCS phase of pn Cooper pairs [30]. In ultracold Fermi gases in the BCS-BEC crossover, some hints for the backbending of the peak energy near k = k F at temperatures above T c were experimentally observed [31]. However, even in the unitary Fermi gas, where pairing correlations are stronger than in neutron matter, it is not clear whether the pseudogap phase exists. Theoretically, it tends to be very visible in non self-consistent T-matrix approaches, such as the present one, but less pronounced in self-consistent Green's functions and quantum Monte-Carlo calculations [32]. A pseudogap can have some effect on the spin susceptibility and the specific heat [33], but in the present case these effects are probably too weak to be observed in neutron stars. Figure 5 shows the shape of the occupation numbers, computed at the critical temperature (in practice, slightly above T c ), for the different approximations defined in Sec. II E. The left panel (a) was computed with µ = 1 MeV, where the neutron matter is in a strong-coupling regime since with T c = 0.27 MeV the ratio T c /µ is not small. We notice that the occupation numbers including the correlations look all very different from the free ones (dotted black line) but they depend strongly on the choice of the approximation which is used (NSR, subtraction of the HF self-energy, subtraction of the total self-energy, or full Dyson equation).

C. Occupation numbers
The right panel (b) was computed with µ = 20 MeV. In this case, the critical temperature of T c = 1.52 MeV is close to its maximum, but the ratio T c /µ is much smaller than for µ = 1 MeV, and we are therefore more in a weak-coupling situation. In this case, one would expect that BCS theory is valid and the occupation numbers above T c are close to the free ones. As we see from the figure, this is indeed the case for all approximations except NSR. The peak in the NSR occupation numbers is clearly unphysical. When we look at the formula (29), we can see that the derivative of the Fermi function will produce a peak near k F . This has nothing to do with correlations but it simply reflects the shift of the quasiparticle energies, treated to first order in a Taylor expansion. This effect is absent in the other approximations where the chemical potential has to be interpreted as an effective one, µ , which already includes the shift, as pointed out below Eq. (31). Therefore, the NSR curve and the other curves are not really comparable as they correspond to different values of the real chemical potentials.

D. Correlated density
The main objective of the T-matrix theory is to determine the density dependence of the critical temperature. As mentioned in Sec. III A, the relation between T c and µ is the same as in BCS theory, and the effect of the correlations above T c enters only through the µ dependence of the density ρ. The density is calculated directly from the occupation numbers as where the factor of two accounts for the spin degeneracy. As we have seen in Sec. II E, for the approximations that treat the correlations only perturbatively (NSR, HF subtraction or total subtraction), the density is naturally separated into two parts such that ρ = ρ 0 + ρ c , with ρ 0 the free density and ρ c the correlated density. In the case of the full Dyson equation, we directly obtain the total occupation numbers and therefore the total density ρ. In this case, we define the correlated density as ρ c = ρ − ρ 0 . The four panels in Fig. 6 represent the ratios ρ c /ρ 0 , computed at the critical temperature T c (µ), as a function of the chemical potential µ, in the four different approximations: (a) NSR, Eq. (29); (b) HF subtraction, Eq. (33); (c) total subtraction, Eq. (35); and (d) full Dyson equation, Eq. (37). The calculations were done using V low-k interactions with different cutoffs Λ and with the separable interaction. Before we turn to a discussion of the cutoff dependence, let us look at the results obtained with the largest cutoffs Λ = 2 and 3 fm (orange and blue long dashes), which lie practically on top of each other and are very close to the results obtained with the separable interaction (gray solid lines).
We see that at very small chemical potential, the ratio ρ c /ρ 0 rises quickly, because the nn scattering length is so large that one quickly gets into the strong-coupling regime. However, what happens at larger chemical potentials depends on the approximation. Since the nn interaction gets weaker at higher momentum, one would expect that neutron matter returns into the weak-coupling regime at high density, i.e., at high chemical potentials, and in this case, the ratio ρ c /ρ 0 should become small. We see that this is true for all approximations except the NSR one [ Fig. 6(a)]. This shows us once again the problem that exists with the NSR theory, namely that it counts the effect of the shift of the quasiparticle energy as "correlations" (or as "fluctuations" in the terminology of Ref. [19]). Thus, in the weak coupling limit at large values of the chemical potential, the NSR "correlated density" remains very high whereas it should tend towards zero. 1 This problem is solved in the other three methods by the subtraction of the quasiparticle energy shift from the self-energy. However, as mentioned in Sec. III C, some caution should be used when comparing Fig. 6(b-c) with the NSR result in Fig. 6(a), because the chemical potential µ does not have the same meaning in the two cases. We see that among the three approaches employing a subtraction, subtracting only the HF self-energy [ Fig. 6(b)], as it was done in [16], yields a significantly larger correlated density than the subtraction of the total self-energy [ Fig. 6(c)] or the full Dyson equation [Fig. 6(d)]. It is interesting to notice that the last two approximations give almost identical results for ρ c in spite of the different shapes of the corresponding occupation numbers [cf. Fig. 5].
Let us now discuss the results obtained with lower cutoffs Λ = 1 − 1.8 fm −1 (red, green and purple lines in Fig. 6). Keeping in mind the dependence of T c on the cutoff, shown in Fig. 2, we limit these curves to the interval of µ where T c is cutoff independent. Nevertheless, we notice that the correlated densities depend on the cutoff values, but to a greater or lesser extent depending on the approximation that is used. In the NSR case [ Fig. 6(a)], only the Λ = 1 fm −1 result deviates significantly from the results obtained with higher cutoffs. The strongest cutoff dependence, even for Λ = 1.8 fm −1 , can be seen in the calculation using the HF subtraction [ Fig. 6(b)], while the cutoff dependence is hardly visible for calculations using total self-energy subtraction [ Fig. 6(c)] and the full Dyson equation [Fig. 6(d)]. This can be easily understood because the HF self-energy Σ HF , Eq. (25), uses directly the interaction v instead of the T matrix Γ, but the V low-k interaction is made such that Γ in vacuum is independent of the cutoff, which implies a strong cutoff dependence of v and hence of Σ HF . Although the cutoff independence of Γ in vacuum does not necessarily ensure the cutoff independence of Γ in the medium, we see that the results that employ the full Γ, i.e., the subtraction of the total self-energy and the full Dyson equation, satisfy the cutoff independence very well. It is, however, not clear why for Λ = 1 fm −1 , the cutoff dependence of the NSR result sets in already at an unexpectedly small value of µ.

E. Density dependence of the critical temperature
Having computed the critical temperature T c and the total density ρ as functions of the chemical potential µ, we can now compute T c as a function ρ. To emphasize the low-density behavior, where strong-coupling effects are most pronounced, we show in Fig. 7 the critical temperature T c as a function of the Fermi momentum k F , which is related to ρ by k F = (3π 2 ρ) 1/3 .
Compared to the BCS result (black curve), which is computed with the uncorrelated density ρ 0 , the curves that account for the correlations in the normal phase are more or less strongly shifted to the right, depending on the correlated density ρ c . The fact that the NSR correlated density does not tend towards zero in the high-density limit causes the NSR curve (red) to be shifted even at the highest densities where T c is very low, as it was also found in Ref. [19]. In the other treatments of the correlations, the curves tend towards the BCS one at high density, which is the expected result since this limit corresponds to the weak-coupling region where the BCS theory should be valid. The correlated densities calculated with the subtraction of the total self-energy (purple) and with the full Dyson equation (orange) being almost identical, it is difficult to distinguish these two curves in Fig. 7.
In order to better see where the weak-and strong coupling regimes are located, we plot in Fig. 8(a) the ratio T c /E F as a function of k F , where E F = 2 k 2 F /(2m) denotes the Fermi energy. The BCS theory should hold in the weak-coupling limit T c /E F 1. Therefore, largest deviation from the BCS solid curve (blue) should be seen around k F ≈ 0.3 fm −1 , where the ratio T c /E F has its maximum. Again, as discussed above, the NSR dotted curve (red) is the only one that does not tend towards the BCS curve at high density. In the case of HF subtraction (green short dashed curve), used in Ref. [16], a significant reduction of T c compared to the BCS one is seen in the interval k F ≈ 0.1 − 0.8 fm −1 (ρ ≈ 5 · 10 −5 − 2 · 10 −2 fm −3 ), and the reduction can be up to ∼ 30 %. With the subtraction of the total self-energy (purple short dashed curve) or the full Dyson equation (orange long dashed curve), a significant deviation from the BCS result is observed only in the smaller interval k F ≈ 0.1 − 0.6 fm −1 (i.e., ρ ∼ 5 · 10 −5 − 7 · 10 −3 fm −3 ) and the reduction is at most ∼ 12 %.
The small correlated densities that we get with the subtraction of the total self-energy or with the Dyson equation are quite astonishing. In ultracold atoms (i.e., with a contact interaction) in the unitary limit, the experimental result is T c /E F 0.16 [34,35], while the BCS theory predicts T c /E F 0.5. A large part of the observed suppression is believed to come from the correlated pairs above T c , for instance with NSR one finds T c /E F 0.23 [20], i.e., a suppression by more than 50 % compared to BCS. The remaining reduction could be, e.g., due to screening effects [36], but within the self-consistent Green's function formalism one finds T c /E F 0.16 even though screening is not included [28].
Let us therefore test our methods by verifying that we reproduce the known results. To that end, we repeat the calculations with the modified value of the coupling constant g of the separable interaction corresponding to the unitary limit given in Table I. Figure 8(b) shows us the ratio T c /E F as a function of k F for this case. We notice that T c /E F remains finite even in the limit k F → 0. This is a peculiarity of the unitary limit, because the dimensionless parameter k F a remains infinite independently of the value of k F . However, the finite range r 0 of the interaction becomes negligible in this limit because the relevant dimensionless combination k F r 0 tends to zero. Hence, in the limit k F → 0, our results reproduce those obtained in ultracold atoms with a contact interaction: In BCS, we find T c /E F 0.5 and in the correlated cases (NSR, HF subtraction, total self-energy subtraction, and full Dyson equation) we find T c /E F 0.23 which is compatible with the results obtained in [20] for NSR and in [37] for the subtraction of the total self-energy (denoted there Zimmermann-Stolz (ZS) scheme after Ref. [27]).
For finite values of k F , the interaction gets weaker because of its momentum dependence, and we expect that, as a function of k F , we should find a similar behavior as if one increases the parameter −1/(k F a) in the BEC-BCS crossover of cold atoms. All curves except the NSR one tend towards the BCS result. We also find that the behavior of the curves for NSR and total self-energy subtraction is similar to that obtained in cold atoms, namely that the NSR curve is strictly increasing with increasing interaction on the BCS side (since the maximum lies on the BEC side [20]), while the curve for total self-energy subtraction has its maximum on the BCS side [37].
From this comparison we conclude that the smallness of the correlation correction that we find in Fig. 8(a) is a consequence of the finite scattering length a = −15.9 fm of the realistic separable interaction. Although this value is quite large compared to other nuclear scales, the parameter k F a gets only large when the finite range of the interaction already starts to weaken it, so that the BCS-BEC crossover effects remain rather weak.

IV. OPEN QUESTIONS
A. Problem of the subtraction In Sec. II E, we argued that in a non-self-consistent treatment of the propagators, we have to subtract from the self-energy the mean-field like quantity Σ(k, ξ k ). In Fig. 9, we represent the momentum dependence of this shift in the case of the separable interaction. By defining k µ = √ 2mµ/ which roughly indicates the position of the Fermi surface, we notice a big jump in this region. It is clear that a momentum dependent mean field as computed with, e.g., a Skyrme or Gogny interaction, would never have such a shape. We suspect that by doing this subtraction, while computing the internal propagators in Γ and Σ with the free dispersion relation, we remove parts of the physical effects caused by the self-energy. It seems probable that this subtraction, as it was implicitly also used in [12-14, 27, 37], reduces the pseudogap and the correlated densities. We notice that in the literature there are other possible ways to do this subtraction. One of them, used in [38,39], consists in neglecting the momentum dependence and taking it at k = k µ , i.e., subtracting a constant Σ 0 = Re Σ(k µ , 0). Looking at the graph, we see that this is problematic, too, because the subtraction is precisely fixed in the zone of the jump where the function varies rapidly, so that the value of the subtraction depends very sensitively on the approximation. Qualitatively, we think that it would be more appropriate to subtract a smooth curve as schematically represented in Fig. 9 by the red dotted line.

B. Effect of the quasiparticle weight
In Ref. [40], a new effect was discussed that strongly reduces the critical temperature, namely the reduced quasiparticle weight. It is well known that the Cooper instability comes from the jump of the Pauli factor 1 − f F (ξ k/2+p ) − f F (ξ k/2−p ) in the numerator of Eq. (3) at p = k µ for k = 0. However, this Pauli factor is obtained with uncorrelated Green's functions having a quasiparticle weight of unity, Z = 1. In Ref. [40], it was argued that the amplitude of the jump should be multiplied by the quasiparticle weight Z < 1, which in that work was determined from a Brückner calculation.
A more self-consistent method to include this effect would be the so-called renormalized RPA (r-RPA) [41], where, in the particle-particle channel, the Pauli factor is replaced by with n(k) the correlated occupation numbers, which should be calculated self-consistently. At zero-temperature, the correlated occupation numbers have a jump of height Z at k = k F and therefore the height of the jump of the Pauli factor is also reduced by a factor of Z. We recently applied this method to the case of strongly polarized Fermi gases at zero temperature [42] where it reduces the critical polarization, which plays a similar role as the critical temperature in the present case. However, here it turns out that it is not possible to realize a self-consistent calculation of r-RPA type. When we insert the correlated occupation numbers in the Pauli factorQ(k, p), it does no longer vanish at p = k µ , leading to a non-vanishing imaginary part of the vertex function Γ at ω = 0. This makes it impossible to calculate the self-energy, as one can see from Eq. (26): Im Γ must vanish where the Bose function f B diverges.
Despite this problem, let us estimate the order of magnitude of the effect. To that end, we replace in the angle averaged Pauli factorQ the free occupation numbers by the correlated ones according to Eq. (40), calculated with the full Dyson equation with the separable interaction. For this choice, in certain regions (in the not too strongly coupled regime), the imaginary part of Im J(k = 0, ω) [cf. Eq. (22)] changes sign almost at ω = 0 so that we can try to estimate the critical temperature from Re J(k = 0, ω = 0 ; T = T c ) = 1/g [cf. Eq. (21)]. This is illustrated in Fig. 10. Notice that the occupation numbers n(k) used in Eq. (40) are computed with the ordinary J and not self-consistently, because, as mentioned above, we cannot compute the self-energy with the modified J. Hence, our modified J, denoted J (1) , can be regarded as a first step of a self-consistent iteration. Furthermore, as it is impossible to compute the correlated occupation numbers below T BCS c , it is necessary to extrapolate J (1) to the new critical temperature T (1) c . In the example µ = 20 MeV presented in Fig. 10, we find that at temperatures not too close to T BCS c this effect reduces Re J(0, 0) by about 10 %, reducing T c by about 30 %. Although this is only a rough estimate, it indicates that this effect might be important.

V. CONCLUSIONS
The in-medium T-matrix formalism makes it possible to describe pair correlations in neutron matter. To model the nn interaction, we use the effective low-momentum interaction V low-k . We presented a numerical method to solve the integral equation for the vertex function Γ. It is also possible to replace V low-k by a separable interaction which gives very similar results and largely facilitates the numerical computations.
Because of the large nn scattering length, it is believed that dilute neutron matter is close to the unitary Fermi gas and that BCS-BEC crossover physics plays an important role. We pointed out that the standard NSR theory, which is widely used for the BCS-BEC crossover, has the shortcoming not to tend towards the weak-coupling BCS theory at high density. It leads to unphysical occupation numbers and to an incorrect shift of the T c vs. ρ curve, because it attributes the mean-field like shift of the quasiparticle energy to the correlations.
We compared different prescriptions to avoid these deficiencies by subtracting this shift from the dynamical (i.e., energy dependent) part of the self-energy. Subtracting only the Hartree-Fock (HF) part, as in Ref. [16], is already enough to make the correlations tend towards zero in the weak-coupling limit. To go beyond this approximation, we subtract the total on-shell self-energy, as implicitly done within the Zimmermann-Stolz scheme [27] which was previously used for the description of the BCS-BEC crossover in dilute symmetric nuclear matter [12][13][14]. This prescription avoids the unphysical dependence of the HF subtraction on the cutoff of the V low-k interaction, but has the effect of further reducing the correlations and therefore the crossover effects. With the explicit calculation of the self-energy, we are also able to calculate the correlations through the full Dyson equation and not just its truncation to the first order as it is customary to do. This method gives us access to the spectral function and allows us to highlight the existence of a pseudo-gap in neutron matter above T c . As discussed in Sec. IV A, it is possible that we underestimate the pseudogap as a consequence of the subtraction method.
At a quantitative level, it seems that the suppression of the critical temperature for a given density, due to the density of preformed pairs, is too weak to have much effect on neutron-star observables. However, a more important effect is expected to result from the quasiparticle weight Z < 1, as also pointed out in Ref. [40]. In the weak-coupling regime we estimate that this effect may reduce T c by ∼ 30 %, but our present theory does not allow us to compute it self-consistently or in the strong-coupling regime.
Finally, a large uncertainty comes from screening of the bare nn interaction in the medium [17,18,40] which is not included in the present study. These effects depend sensitively on the approximations that are used, on the Landau parameters and on the momentum dependence of the effective interaction used in the residual particle-hole interaction. From the results of Ref. [18], it seems most likely that in the low-density region corresponding to the strong-coupling regime, screening suppresses the critical temperature by about 30 − 40 %, thereby reducing further the crossover effects.