Nuclear Pairing from Two-body Microscopic Forces: Analysis of the Cooper Pair Wavefunctions

In a recent paper we studied the behavior of the pairing gaps $\Delta_F$ as a function of the Fermi momentum $k_F$ for neutron and nuclear matter in all relevant angular momentum channels where superfluidity is believed to naturally emerge. The calculations employed realistic chiral nucleon-nucleon potentials with the inclusion of three-body forces and self-energy effects. In this contribution, after a detailed description of the numerical method we employed in the solution of the BCS equations, we will show a preliminary analysis of the Cooper pair wavefunctions.


I. Introduction
The goal of this article is to solve the BCS equations starting from modern nucleon-nucleon potentials (at N3LO in the chiral expansion [2]) and to perform preliminary calculations of the Cooper pair wavefunctions as a first step towards an analysis of the occurrence of BCS-BEC crossover in nuclear systems [4].

II. Khodel's method
In this section we explain the method suggested in Ref. [5] to solve the BCS equations by partial-wave decomposition that has been proven to be stable even for small values of the gap and to require only the initial assumption of a scale factor δ (results, of course, will be δ-independent). The BCS equation reads in terms of the NN potential V(k, k ) = k |V| k as follows with E(k) 2 = ξ(k) 2 + |∆(k)| 2 and where ξ(k) = ε(k) − µ, ε(k) denotes the singleparticle energy and µ is the chemical potential. We can decompose both the interaction and the gap function where Y lm (k) denotes the spherical harmonics, l and m are the quantum numbers associated with the orbital angular momentum and its projection along the z axis and P l (k ·k ) refers to the Legendre polynomials. After performing an angle-average approximation we have the following equation for any value of l where Λ = 1 + (l − l )/2, j refers to the total angular momentum (J = l + S) quantum number including spin S and now E(k) 2 = ξ(k) 2 + ∑ jl ∆ j l (k) 2 . Gaps with different l and j are coupled due to the energy denominator but we assume that different components of the interaction mainly act on non-overlapping intervals in density.
We define an auxiliary potential W according to vanishes on the Fermi surface. The coupled gap equations can be rewritten as where dτ = k 2 dk/π and the coefficients D ll satisfy The gap is defined as follows where The property that W ll (k, k ) vanishes on the Fermi surface ensures a very weak dependence of χ l 1 l 2 l (k) on the exact value of the gap so that, in first approximation, it is possible to rewrite the previous equation (9) as We use this equation to evaluate χ l 1 l 2 l (k) initially by matrix inversion, then we use this function to self-consistently evaluate D ll . Finally, we solve the system given by Eqs. (7)-(9) in a self-consistent procedure as shown in Fig. 1 (left panel). We always assume µ = ε F and adopt the relativistic version of the single-particle where M N is the nucleon mass.

Numerical analysis
As a benchmark, we will consider the 1 S 0 pairing gap in neutron matter, but the same conclusions can be drawn for all interaction channels. Concerning the singlet channel, at the two-body level, we find a good agreement with the gap computed from well known realistic potentials like the CD-Bonn or Nijmegen interactions [6], except for larger densities where the N3LO gap exhibits a higher value (phase shifts from the chiral N3LO potential exhibit more attraction than the CD-Bonn potential for high momenta [7]). We tested Khodel's method [5] against the variation of the following three parameters: n gauss (number of Gauss integration points), Λ k (cutoff for integrals in momentum space, see Eq. (4)) and δ (the scale factor). In Fig. 1 (right side) we summarise our results. In the upper panel (a) we calculated ∆ F for different values of the momentum cutoff (using n gauss = 200 and δ = 1 × 10 −10 MeV) where in the second panel (b) we varied n gauss (keeping Λ k = 4.5 fm −1 and δ = 1 × 10 −10 MeV) and in the lower panel (c) we changed δ (with n gauss = 200 and Λ k = 4.5 fm −1 ) by orders of magnitude.
Our conclusion is that the method proposed by Khodel [5] is a very stable procedure to study nuclear superfluidity if a reasonable number of Gaussian points (≥ 100) and a realistic momentum cutoff (≥ 4 fm −1 ) are employed.

III. Cooper pair wavefunctions
Recently it has been shown that, at low density, nuclear matter can be subjected to a phase transition, belonging to the BCS-BEC crossover phenomenon [8], caused by the collapse of the Cooper pairs in deuteron-like particles. In this phase, the nucleons forming the pairs are strongly correlated; this correlation gives an extra binding energy to the nuclear equation of state and could lead to identify new properties at low density. Because BCS equations are still valid in the BEC regime [9], it is useful to study the evolution of pairs of correlated nucleons as a function of density and spatial coordinate, starting from the solution of Eqs. (7-9). The Cooper pair wavefunction is defined as follows where |Φ 0 is the BCS ground state, ψ † is the particle creation operator and C, C are normalization factors that can be fixed by imposing normalization conditions.
In the following discussion we will indicate as ρ(r) the probability density to find the nucleons forming the pair at a distance r, P(r) = r 0 dr ρ (r ) where ρ(r) = |Ψ(r)| 2 r 2 in the singlet channel. Observing that the pairing gap in the S channel is larger then the one in the D channel we assumed, in the SD channel, |Ψ(r)| 2 ≈ |Ψ l=0 (r)| 2 and we have approximated ∆(k) with ∆ l=0 (k). Analogously, in the PF channel, ∆(k) ∆ l=1 (k).
As suggested by Matsuo [10] the spatial behaviour of the Cooper pair varies strongly with the Fermi momentum k F (or density ρ 0 = g/(2π) 3 k 3 F , where g = 1 for neutron matter and g = 2 in symmetric nuclear matter), in particular the ρ-profile shows strong variations. The weak-coupling limit is known to lead to an exponential falloff convoluted with an oscillation suggesting a large correlation length. On the other hand, a pronounced peak with small oscillations could be interpreted as a sign of a transition to a different regime, i.e. the so called BCS-BEC crossover. Using phenomenological pairing interactions, Matsuo [10] suggested that in the singlet channel, over a wide range of densities (ρ/ρ 0 10 −4 − 0.5), the spatial dineutron correlation is strong and a possible crossover region could be found in the density region ρ/ρ 0 10 −4 − 10 −1 .
In the following paragraphs we summarise our results obtained with microscopic two-body NN forces.
Singlet case ( 1 S 0 in neutron matter) In Fig. 2 we show in the upper panel (a) the pairing gap ∆ F as a function of the Fermi momentum k F [1] for the singlet channel in neutron matter, where in the lower panel (b) the probability density ρ is plotted as a function of the coordinate r and Fermi momentum k F . For the sake of clarity we also show ρ(r) for few selected values of k F in the graph labelled by (c). Results of Ref. [10] are substantially confirmed: at densities below k F 0.8 fm −1 the system belongs to a crossover region, while at higher densities is in the pure BCS region. Triplet cases ( 3 SD 1 in nuclear matter and 3 PF 2 in neutron matter) In Figs. 3 and 4 we show the same informations as before, respectively for the 3 SD 1 and 3 PF 2 gaps. In the first case, our results seems to support the previous indication that the Cooper pair wave function merges smoothly into a "deuteronlike" wavefunction as density decreases [11] and for higher densities the system belongs to a BCS-BEC crossover region. In the latter case, it is very clear that the possibility of a BCS-BEC crossover is ruled out because ρ(r) behaves exactly as expected in a BCS weak coupling regime.
Of course, the previous results represent only a very preliminary analysis that need to be refined and improved, i.e. including three-body forces and self-energy effects [4].

IV. Conclusions
We have presented preliminary calculations of the Cooper pair wavefunctions in infinite matter employing realistic two-body nuclear forces derived within the framework of chiral effective field theory. The BCS gap equation is solved employing Khodel's method, which is found to be numerically stable with respect to variations of Gaussian integration points, momentum cutoffs and the scale factor.     Fig. 2 for the 3 PF 2 channel. In the upper panel (a) we include as benchmarks some well known results (we refer the reader to Ref. [6] for more details).