$^1S_0$ pairing in neutron matter
Abstract
We report calculations of the superfluid pairing gap in neutron matter for the components of the Reid softcore and the Argonne twonucleon interactions. Groundstate calculations have been carried out using the central part of the operatorbasis representation of these interactions to determine optimal JastrowFeenberg correlations and corresponding effective pairing interactions within the correlatedbasis formalism (CBF), the required matrix elements in the correlated basis being evaluated by Fermi hypernettedchain techniques. Different implementations of the FermiHypernetted Chain EulerLagrange method (FHNCEL) agree at the percent level up to nuclear matter saturation density. For the assumed interactions, which are realistic within the low density range involved in neutron pairing, we did not find a dimerization instability arising from divergence of the inmedium scattering length, as was reported recently for simple squarewell and LennardJones potential models (Phys. Rev. A 92, 023640 (2015)).
Keywords:
Superfluidity, Quantum Fluids, Neutron Matterupdate
1 Introduction
1.1 Adaptation of BCS theory to nuclear systems
The nature and role of fermionic pairing and superfluidity in nuclei and nuclear matter became a subject of great interest shortly after publication of the landmark paper by Bardeen, Cooper, and Schrieffer (BCS) establishing the physical basis of superconductivity in metals ^{1, 2}. Bohr, Mottelson, and Pines ^{3} were quick to recognize implications of this development for a deeper understanding of nuclear phenomena, relating it to evidence for a characteristic energy gap between the ground state and the first intrinsic excitation in a certain class of nuclei.
Concurrently, there was growing interest among nuclear theorists in what could be learned from the quantum manybody problem of infinite nuclear matter composed of nucleons interacting through the best nucleonnucleon (NN) potentials available at the time. Cooper, Mills, and Sessler^{4} (CMS) were the first to apply BCS theory to such a system. They encountered two obstacles when attempting to solve the BCS equation for the superfluid energy gap as a function of momentum .
To understand what they faced, it is necessary to consider the BCS gap equation, written in the generic form
(1) 
where defines the pairing matrix elements of the bare twobody potential , while
(2) 
represents the (gapped) quasiparticle energy in the superfluid state, with an “appropriate” singleparticle energy related to the normal state. Given the original BCS trial ground state
(3) 
(but written slightly differently in terms of Bogoliubov amplitudes , satisfying the normalizing condition ), the expression (1) of the gap equation can be derived from the EulerLagrange variational principle following exactly the same path as in the 1957 BCS paper ^{1} and in Schrieffer’s book ^{5}. As the BCS state does not have a definite particle number, the chemical potential (determined from the number density) is introduced as a Lagrange parameter to accommodate the constraint that the particle number is conserved on average.
Of the two problems CMS faced in implementing BCS theory for nuclear matter, they managed to solve what appeared to be the more difficult one, and finessed the other. During this same period in the midtolate 1950s, it had become apparent that an acceptable model of the NN interaction, fitted to the available NN scattering data and the deuteron, must possess a strong inner repulsion, most commonly taken to be a hard core. This precluded solving the BCS gap equation as formulated in momentum space, because the necessary pairing matrix elements of the NN potential would be undefined. However, CMS recognized that the BCS gap equation could be transformed to coordinate space to yield a nonlinear but Schrödingerlike equation for an underlying twobody problem. The analog of the wave function for the separation vector is the pairing function , which may also be regarded as the superfluid order parameter. Basically, is the Fourier transform of the product of Bogoliubov amplitudes, or equivalently of . Therefore the problem created by the hard core of the NN potential could be solved, for the same reason that the Schrödinger equation for a hardsphere scattering potential has a solution.
The second problem confronting CMS was what to take for the singleparticle energy in the expression for . There is first a subtlety relating to that should be exposed, for the record. The above derivation leads to the actual expression
(4) 
This contains the Fermisurface smearing factor represented by , and hence requires a solution of the pairing problem before can be evaluated. In practice, this factor is almost always replaced by the Fermi step, converting to a standard HartreeFock singleparticle energy. It is argued, in most cases safely, that this can be done because the gap is much smaller than the Fermi energy, thus decoupling from the rest of the gap problem.
The primary issue raised by the expression (4) is not at all subtle. If the bare NN interaction contains a hard core, the HartreeFock matrix elements it contains are infinite; nor would the results for be sensible if the interaction remains finite, but features an internal repulsion strong enough to achieve empirical saturation of nuclear matter. CMS were forced to finesse this second problem; they imposed an effectivemass spectrum . With this step, the problem was welldefined and in principal soluble for ; however, for a time only the existence of a superfluid solution was established ^{6}, due to the limited computational resources of that period.
In summary, the nature of the BCS theory of superfluidity is such that its application to nuclear systems is practical, in particular for the hypothetical system of infinite nuclear matter and certain nucleonic subsystems existing in neutron stars. However, due to the presence of a strong shortrange repulsion in the bare NN interaction, one must make a reasonable, but ad hoc, assumption for the normalstate singleparticle energy. The theory has the capacity to generate twobody correlations that can accommodate even the effects of a hard core, although the problem must then be solved in coordinate space. Solution of the problem in momentum space, i.e. the original gap equation (1), does in fact become possible if the NN interaction, even though strongly repulsive at short distance, has a Fourier transform. (For some interactions including the Reid softcore potential ^{7, 8}, numerical solution can present some technical difficulty; this can be avoided by applying the separation approach developed in Ref. 9).
Yet the status of nuclear BCS as described remains unsatisfactory for potentials with repulsive cores. This issue was naturally addressed by the introduction of JastrowFeenberg correlation factors ^{10, 11, 12}. Clusterexpansion techniques were applied to evaluate the required expectation values^{13}. The corresponding gap equations were studied and procedures for their solution explored, with applications not only to isospinsymmetric nuclear matter (equal numbers of neutrons and protons) ^{10}, but also pure neutron matter and stable nucleonic matter relevant to neutronstar interiors ^{11, 14}. In the mid1970s, major advances in microscopic quantum manybody theory involving correlated basis functions (CBF) were made with the replacement of cluster expansions by Fermi HypernettedChain (FHNC) diagramresummation techniques ^{13, 15}, facilitating accurate evaluation of expectation values and matrix elements of observables in a correlated basis and culminating in a framework for EulerLagrange optimization of JastrowFeenberg correlations. When implemented in a BCS extension, these advances have made possible the development of a rigorous correlated BCS (CBCS) theory ^{16} (see also Ref. 17) that respects the U(1) symmetrybreaking aspect of the superfluid state – i.e., nonconservation of particle number. Some earlier applications of CBCS theory to nuclear systems, and especially neutronstar matter, may be found in Refs. 8, 16. A recent indepth study of correlations in the lowdensity Fermi gas^{18}, with emphasis on the presence of Cooper pairing and dimerization, documents the power of the EulerLagrange FHNC approach adopted in the present work, especially when coupled with CBF perturbation theory.
1.2 Extension of BCS asymptotics to structured interactions and inmedium effects
Having derived the equivalent of the gap equation (1), BCS went on to simplify the pairing interaction in a way suitable for electron liquids in solids, arriving at the iconic asymptotic result
(5) 
for the value of the energy gap in terms of a cutoff and the coupling constant of the attractive pairing interaction , with denoting the density of states around the Fermi surface. It is important to recognize that this result, being restricted to the weakcoupling regime , is not at all appropriate for nuclear problems. In nuclear systems, the bare twobody interaction is strong, and strongly nonmonotonic in coordinate space. Two parameters are not sufficient to characterize the asymptotic behavior of the gap at relevant densities. See Refs. 9, 19 for extensive analysis and computational exploration of this important distinction. The latter reference includes an asymptotic study in which the pairing interaction is characterized by an additional parameter along with the traditional coupling constant and cutoff frequency . This “stiffness” parameter is introduced to represent a nontrivial momentum dependence of the pairing interaction . Asymptotic behavior in the four quadrants , is explored in Ref. 19, pointing to the existence of solutions with behavior quite distinct from the familiar relation (5), in addition to a BCSanalog.
Another asymptotic formula of special interest (and of long standing) is that of Gorkov and MelikBarkhudarov ^{20} (GM),
(6) 
written for the zerotemperature gap rather than the critical temperature. This result was derived by fieldtheoretic methods in the limit of an infinitely dilute gas of interacting spin fermions, with . Here is the vacuum scattering length, assumed to be negative, is the Fermi energy, the fermion mass, and the Fermi momentum. The prefactor is an inmedium correction for a polarizationinduced interaction corresponding to exchange of virtual phonons. The same result without the GM prefactor was rederived several times in the 1990’s, basically by summing ladder diagrams for the bare interaction (see Ref. 21, where the GM prefactor is generalized to for an arbitrary number of fermion species).
In the recent work previously cited^{18}, it has been argued (cf. Eqs. (3.25) and (3.26)) that if one has corrections of the inmedium scattering length to the vacuum scattering length of the form
(7) 
it follows that
(8) 
The GM factor is just one of these corrections, which still assumes that the pairing matrix element at is the same as that at . Removing this assumption produces yet another correction of the same kind.
The above summary of BCS asymptotics is intended to provide deep background for the present work on neutron matter at densities occurring in the innercrust region of neutron stars, but their direct relevance is open to question. The neutron densities involved in this application are low compared to the saturation density of isospinsymmetric nuclear matter, which, in pure neutron matter, would correspond to a value of about . We will find that pairing in neutron matter is strongest at somewhat less than half that value, thus at a density an order of magnitude below . On the other hand, given the unusually large magnitude of the neutronneutron wave scattering length, , the diluteness condition implies , over three orders of magnitude lower in density than that of the physically relevant neutronstar environment. Naturally the dilutelimit asymptotics do apply for the extreme lowdensity tail of the roughly Gaussian shape of vs. in the neutron pairing problem considered here. The higherdensity tail is more relevant; it has been demonstrated in Ref. 9 that dies exponentially to naught as an upper critical density is approached.
1.3 Sensitivity issues in optimization
The foregoing subsections of this introduction provide a rather elaborate background and motivation for the work to be presented. Another motivation is more immediate. Recently, using updated modern NN interactions, gap calculations ^{22} for pure neutron matter have again been performed within the simpler version of correlated BCS theory in which the groundstate energy for a JastrowFeenberg trial function, estimated by a truncated cluster expansion, is minimized with respect to the parameters in an assumed analytic form for the Jastrow twobody correlation function . The correlation function so determined is used to generate a “tamed” effective pairing interaction for calculation of a corresponding superfluid gap in the state.
Ref. 22 presents results for the groundstate energy per particle and the corresponding energy gap, based on the Argonne (AV18) NN interaction^{23} and two trial correlation functions with analytic forms that have been employed in earlier JastrowFeenberg studies of nuclear and neutron matter. The optimal groundstate energies determined for these two choices show only minor quantitative differences over the low range of densities where a significant gap is to be expected (peaking at about 1/8 nuclear saturation density). The two curves obtained for the gap at the Fermi surface, plotted versus Fermi momentum , have a Gaussian appearance. In contrast to the close agreement of the results for the two correlation choices, the corresponding peak values for are found to differ by almost a factor two (with a value 1.8 MeV for the correlation function featuring an overshoot of unity versus 3.3 MeV for one that does not.
This finding could be interpreted as a reflection of the variational property that a small error of order in the wave function only entails an error of order in the energy expectation value, but of order for other observables, with in this case corresponding to the difference in the choices for . But the situation may actually be worse for two reasons: The most immediate one is that the gap itself shows an exponential amplification of errors in the coupling strength and density of states, at least for the standard BCS case. Some results of the present investigation indicate a similar strong sensitivity of gap behavior. The second, more subtle reason, involves the convergence of cluster expansions for correlated wave functions: Typically, the contribution of an body diagram in the energy is amplified by a factor of in its contribution to the effective interactions needed to calculate the coupling matrix elements.
Figure 6 of the paper of Pavlou et al. ^{22} shows plots of versus for the AV18 interaction as obtained in almost a dozen calculations by different theoretical methods, including various versions of Monte Carlo. (Actually, this is a summary figure taken from the review by Gezerlis et al. ^{24} of novel superfluidity in neutron stars, with a curve calculated by Pavlou et al. superimposed.) There is a spread of a factor of six in the peak values of , with a significant spread also in the peak densities. In view of what has been said above, this is hardly surprising, although the calculations may differ in the inclusion of inmedium effects.
At any rate, the message from the considerations of this subsection is that it is imperative within any variational approach to seek truly optimal correlations, without resorting to simple parametrizations, and that is what EulerLagrange FHNC (FHNCEL) can deliver, with minimal error.
The rest of this paper is organized as follows. Section 2 exemplifies what qualifies as a generic manybody theory, first with a brief review of the elements of the JastrowFeenberg theory of the normal ground state of a manyfermion system (Sec. 2.1), then with an introduction to the formalism associated with the method of correlated basis functions (CBF) (Sec. 2.2), concluding with the essentials of a coherent theory of fermion superfluidity within the CBF framework (Sec. 2.3), based on EulerLagrange Fermi HypernettedChain optimization (FHNCEL). In Sec. 3.1, we describe and discuss our application of two types of FHNCEL theory to the groundstate energetics of pure neutron matter. Sec. 3.2 is concerned with solution of the resulting CBF gap equation for pairing, which incorporates the effects of the optimal JastrowFeenberg correlations. Results for energetics (the equation of state) and BCS pairing in CBF framework are presented and discussed for two versions of the bare NN interaction, namely the Reid softcore potential ^{7} and the Argonne interaction^{25}. Well known from earlier microscopic studies of nuclear matter, these choices are quantitatively viable in the lowdensity regime where the pairing state is dominant. Only the central components of these potentials and their projections are needed for determination of the CBF pairing matrix elements, in contrast to the case of the full Argonne interaction ^{23}. (Note that Argonne is central by construction and we use only the part of the Reid potential, i.e., its tensor as well as spinorbit terms being omitted.) Manybody aspects of our findings unique to optimal incorporation of short and longrange correlations within the CBF/FHNC formalism are analyzed. Where meaningful, our predictions for the density dependence of the gap at the Fermi wave number are compared with those from other microscopic calculations. Sec. 4 summarizes ways in which the present numerical study may be improved and extended.
2 Generic ManyBody Theory
2.1 The normal ground state
In this section, we briefly describe the JastrowFeenberg variational method and its implementation to superfluid systems. (For comprehensive background on this manybody approach and its generalization to the method of correlated basis functions, see Refs. 13, 15. Recent descriptions and analysis of its applications to superfluid systems may be found in Refs. 18 and 26). We call this method “generic manybody theory” because the same equations can be derived by Green functions methods ^{27, 28}, from coupledcluster theory ^{29}, and by a generalization of density functional theory to pair distribution functions ^{30}, without mentioning a JastrowFeenberg wave function. We use the JastrowFeenberg point of view here because it is the simplest to implement and to generalize.
For a strongly interacting and translationally invariant normal system, the JastrowFeenberg theory assumes a nonrelativistic manybody Hamiltonian
(9) 
The method starts with an ansatz for the wave function, ^{31}
(10)  
(11) 
where is a normalizing constant. Here denotes a model state, which for normal Fermi systems is a Slaterdeterminant, and is a correlation operator written in general form, but to be truncated at the twobody term in a standard Jastrow calculation. There are basically two ways to deal with this type of wave function. In quantum Monte Carlo studies, the wave function (10) is referred to as “fixednode approximation,” and an optimal correlation function is obtained by stochastic means. Computationally far less demanding are diagrammatic methods, specifically the optimized EulerLagrange Fermihypernetted chain (FHNCEL) method, which is well suited for calculation of physically interesting quantities. These diagrammatic methods have been successfully applied to such highly correlated Fermi systems as He at ^{32}. We have shown in recent work ^{33} that even the simplest version of the FHNCEL theory is accurate within better than one percent at densities less than 25% of the saturation density of liquid He, and the same or better performance is expected for nuclear systems.
The correlations are obtained by minimizing the energy, i.e. by solving the EulerLagrange (EL) equations
(12)  
(13) 
Evaluation of the energy (12) for the variational wave function (10), (11) and analysis of the variational problem are carried out by cluster expansion and resummation methods. The procedure has been described at length in review articles ^{13, 32} and pedagogical material ^{15}. Here, we spell out the simplest version of the equations that is consistent with the variational problem (“FHNC//0EL”). These equations do not provide the quantitatively best implementation of this approach ^{32}. Instead, they provide the minimal version of the FHNCEL theory. In particular, they contain the indispensable physics, namely the correct description of both short and longranged correlations.
In the FHNC//0EL approximation, which contains both the random phase approximation (RPA) and the BetheGoldstone equation in a “collective” approximation, the Euler equation (13) can be written in the form
(14) 
where is the static structure factor of the interacting system, is the kinetic energy of a free particle, is the static structure factor of the noninteracting Fermi system, and
(15) 
is the socalled “particlehole interaction.” As usual, we define the Fourier transform with a density factor, i.e.,
(16) 
Auxiliary quantities are the “induced interaction”
(17) 
and the “directdirect correlation function,”
(18) 
a “dressed” analog of the Fourier inverse of . Eqs. (14)(18) form a closed set which can be solved by iteration. Note that the Jastrow correlation function has been eliminated entirely.
2.2 Correlated Basis Functions
Correlated Basis Functions (CBF) theory uses the correlation operator to generate a complete set of basis states through
(19) 
where the are Slater determinants of singleparticle orbitals. We review the method only very briefly, the diagrammatic construction of the relevant ingredients having been derived in Ref. 35 (see also Ref. 32 for further details).
To develop a BCS theory with correlated wave functions it is expedient to introduce a “secondquantized” formulation. The JastrowFeenberg correlation operator in (11) depends on the particle number, i.e. (whenever unambiguous, we omit the corresponding subscript). Starting from the conventional operators that create and annihilate singleparticle states, new creation and annihilation operators of correlated states are defined by their action on the correlated basis states:
(20)  
(21) 
According to these definitions, and obey the same commutation rules as the creation and annihilation operators and of uncorrelated states, but they are not Hermitian conjugates of one another.
For offdiagonal elements of an body operator , we sort the quantum numbers and such that is mapped onto by
(22) 
Then, the matrix elements depend only on the difference between the states and , and not on the states as a whole. Consequently, can be written as the matrix element of a body operator
(23) 
with the index indicating antisymmetrization. In homogeneous systems, the continuous parts of the quantum numbers are wave numbers ; we abbreviate their difference as .
The key quantities for the execution of the theory are diagonal and offdiagonal matrix elements of unity and ,
(24)  
(25) 
Eq. (25) defines a natural decomposition ^{15, 35} of the matrix elements of into the offdiagonal quantities and and diagonal quantities . These diagonal matrix elements are additive to leading order in the particle number, allowing us to define the CBF singleparticle energies that satisfy
(26) 
According to Eq. (23), and define particle operators and , e.g.
(27)  
Diagrammatic representations of and have the same topology ^{35}. In the next section, we will show that in dealing with pairing phenomena, only the twobody operators are needed.
In principle, and are nonlocal body operators. The leading, local contributions to these operators are readily expressed in terms of the diagrammatic quantities of FHNCEL theory ^{32}:
(28) 
For further reference we also display the coordinate space form of the interaction :
(29)  
which exhibits somewhat more clearly the physical meaning of the individual terms: The factor describes the shortranged correlations, the term describes the cost in kinetic energy for bending the wave functions at short distances, and the induced potential describes the corrections due to phonon exchange. In the local approximations spelled out in Eqs. (28), the CBF singleparticle energies (26) assume the simple form
(30) 
with
(31)  
(32) 
where is the degree of degeneracy of the singleparticle states, is the Slater exchange function, and the constant is determined by the condition . In the limit of a weakly interacting system, we have , and the reduce to the HartreeFock singleparticle energies (4).
2.3 BCS Theory with correlated wave functions
The BCS theory of fermion superfluidity generalizes the HartreeFock model by introducing a superposition of independentparticle wave functions corresponding to different particle numbers ^{36}, represented economically by Eq. (3) in terms of Bogoliubov amplitudes , .
The most natural way to deal with a strongly correlated system is to first project the bare BCS state on an arbitrary member of a complete set of independentparticle states with fixed particle numbers. Then apply the correlation operator to that state, normalize the result, and finally, sum over all particle numbers . Thus, the correlated BCS (CBCS) state is taken as
(33) 
The trial state (33) superposes the correlated basis states with the same amplitudes the model states have in the corresponding expansion of the original BCS vector. It is important to note that this state differs from the state proposed, analyzed, and applied computationally in Refs. 37, 38, 39, which fails to include the normalizing denominators present in Eq. (19). As shown in Ref. 17, this option leads to a meaningful gap equation only if specific diverging quantities are omitted.
Consider now the expectation value of an arbitrary twobody operator with respect to the superfluid state (33):
(34) 
For superfluid gaps that are small compared to the Fermi energy, it suffices to consider the interaction of only one Cooper pair at a time. In that case, one need retain only the terms of first order in the deviation and those of second order in the product . We refer to this as the “decoupling approximation”. The error introduced thereby is of order , where is the superfluid gap energy at the Fermi energy . Within this approximation, neither the pairing matrix elements nor the singleparticle energies entering the gap equation depend on the Bogoliubov parameters , .
The calculation of for correlated states is somewhat tedious^{16}. Details may be found in Refs. 16, 18; we only give the final result. The energy of the superfluid state may be derived from
(35)  
in terms of the “pairing interaction” specified by
(36)  
(37)  
(38) 
With the result (35), we have arrived at a formulation of the theory which is isomorphic with the BCS theory for weakly interacting systems. Closer inspection^{18} reveals that our approach corresponds to a BCS theory formulated in terms of the scattering matrix ^{40}. The correlation operator serves here to tame the shortrange dynamical correlations. The effective interaction is just an energyindependent approximation of the matrix.
We may now implement the standard procedure of determining the Bogoliubov amplitudes , , by variation of the energy expectation (35). This leads to the familiar gap equation
(39) 
The conventional (i.e. “uncorrelated” or “meanfield”) BCS gap equation ^{41} is retrieved by replacing the effective interaction by the pairing matrix of the bare interaction. The lowclusterorder approximations to the pairing interaction used by Benhar^{42} and Pavlou et. al. ^{22} are obtained by setting in Eqs. (28) and (29) and omitting the induced interaction .
3 Application to Neutron Matter
3.1 Energetics
We have carried out groundstate calculations for static properties and superfluid pairing gaps in neutron matter based on two representative NN interactions acting in the channel, namely the central parts of the Reid softcore potential ^{7} as formulated in Eqs. (A.1)(A.8) of Ref. 43, generally referred to as Reid , and the Argonne potential ^{23}. In the density range of interest, any interaction must give very close to the same for nuclear matter and the deuteron, as long as it fits the wave scattering data and the deuteron. We have carried out two types of calculations: Full FHNCEL calculations as described, for example, in Refs. 34 and 32, and FHNC//0EL calculations as described in section 2.1. Results for the equation of state for these two calculations, plotted as versus Fermi momentum , are shown in Fig. 1.
The picture is very similar to that found for LennardJones interactions ^{33}: The FHNC//0 approximation performs well up to about half nuclear saturation density. It is also noteworthy that the two potentials lead to very nearly the same equation of state in the density range considered, with the stars for the respective FHNC//0EL calculations overlapping.
For the Reid potential we have also examined the importance of optimized triplet correlations (i.e., nonvanishing in Eq. (11) and elementarydiagram cluster contributions as outlined in Ref. 32) and found their influence negligible. We have also tried the central part of the full Argonne potential in the channel. It turns out that this component of the interaction is too soft to lead to a stable solution of the Euler equation. The problem can be solved by an artificial enhancement of the repulsive regime, but the results depend sensitively on that enhancement factor and hence were considered unreliable.
3.2 BCS pairing
Once the groundstate correlations are known, the superfluid gap function can be determined by solving the gap equation (39). Since we are concerned with pairing, we have inserted the component of the chosen potential model into the effective interaction (29). In the phononexchange correction the central component of the interaction is the appropriate choice.
The gap equation was solved by the eigenvalue method with an adaptive mesh, as outlined in the appendix of Ref. 18. We have primarily adopted a free singleparticle spectrum for as it occurs in Eqs. (36) and (39). One could also use the actual spectrum of CBF singleparticle energies (30), calculated from the effective interactions ^{35}, in both the pairing interaction (36) and the denominator of Eq. (39). We have not done this for the reason outlined below.
At first glance, only the spectrum in the vicinity of the Fermi momentum is relevant. In that regime it can be approximated quite well in terms of an effective mass. Fig. 2 shows the effective mass obtained from the CBF singleparticle energies for both potential models.
Evidently, the effective mass ratio obtained for both potentials is very close to unity.
However, “first glance” may not be sufficient; there is a subtlety to consider: If the gap at the Fermi surface is small, we can replace the pairing interaction by its wave matrix element at the Fermi surface,
(40) 
Then we can write the gap equation as
(41) 
which is almost identical to Eq. (16.91) in Ref. 40. In particular, the second term, which originates from the energy numerator generated in Eq. (39) by the second term of in Eq. (36), has the function of regularizing the integral for large . This feature is lost if the bare mass is used in the relationship (28), and the integral (41) diverges unless a momentumdependent effective mass ratio is used that approaches unity in the limit of large momenta.
A second issue is that it has been known for a long time ^{46, 47} that the effective mass in nuclear systems has a peak around the Fermi surface, however, such a peak is absent in the CBF singleparticle spectrum. An effectivemass enhancement may be obtained by including complex selfenergy corrections; this can be done, for example, by going to higherorder terms in CBF perturbation theory ^{48}. We note that the enhancement effect is much stronger in He (see Refs. 49, 50, 51) due to the softness of the spinfluctuation mode.
In view of these considerations, we have deemed it more accurate to employ the free singleparticle spectrum , and to study the sensitivity of our results to changes in the effective mass. Our results for the superfluid gap for the two potentials are shown in Fig. 3. Evidently the difference of the gap between these two potential models is almost negligible and certainly within the accuracy of the FHNC approximations. To determine the importance of effectivemass corrections we have also solved the gap equation assuming effectivemass ratios between and in both the pairing interaction (36) and the energy denominator (2). The results for the gap define the gray area in Fig. 3; their spread provides a conservative estimate of the importance of a nontrivial singleparticle spectrum.
3.3 Consequences for many–body theory
To conclude this section, let us look more closely at different aspects of the convergence of cluster expansion and resummation techniques. Apart from cold atomic gases – which with rare exceptions like the “unitary limit” pose no challenges to modern manybody theory – pure neutron matter at subnuclear and nuclear densities is, apart from the complications introduced by the nucleonnucleon force, one of the most lenient manyparticle systems provided by nature. This is largely due to the low density of the system, as measured for example by the ratio of the pion Compton wavelength ^{7}, , or the radius of the “hard core,” , to the average particle spacing at the given density . Thus, at the density is , which corresponds to only 20 percent of the saturation density of He.
Evidence for the good convergence of manybody theory for neutron matter in the density regime relevant for pairing is already provided in Fig. 1, which shows that the very simple FHNC//0 approximation for the energy is quite accurate. In fact, even the very simple twobody cluster approximation
(42) 
in which is the kinetic energy of the free Fermi gas and its pair distribution function, yields results virtually identical to the FHNC//0 results plotted in Fig. 1. We have refrained from showing these results in order not to obscure the figure. Note that one can of course identify with in Eq. (42).
These findings are consistent with the fact that the optimal results for and are not very different. To demonstrate this, both functions are plotted in Fig. 4 for three representative densities, and . At the lowest density, the two functions are practically identical. As the density increases, becomes slightly steeper in the attractive regime of the interaction.
From these results, one might be led to conclude that loworder methods are also adequate for calculating the superfluid gap. We remind the reader, however, of the discussions in Sect. 1.3 on the both the sensitivity of quantities other than the energy to the correlation functions, and to the convergence rate of cluster expansions. Accordingly, we have examined the consequences of two approximations: Leaving out the energy numerator generated by the pairing interaction (36) and leaving out the induced interaction .
The most important function of the energy numerator is to regularize the integral in the gap equation for contact interactions, as witnessed in Eq. (41). The situation being discussed at that point is, of course, extreme. More generally, one would expect that the energy numerator term is important whenever the pairing interaction does not fall off sufficiently rapidly for large momenta. This is indeed the case: We show in Fig. 5 the interaction for three representative densities. Evidently, the pairing interaction does not fall off rapidly above . The effect is, of course, most pronounced for low densities. Although the gap is determined solely by the pairing matrix element in the limit of an infinitesimal gap, one must expect significant finiterange effects in the present case where the gap is of the order of 10 to 50 percent of the Fermi energy.
The second new aspect is the appearance of the induced interaction term appearing in the pairing interaction (cf. Eqs. (29) and (17). This term describes the exchange of particlehole excitations ^{27} and is one of the important effects introduced into the CBF version of BCS theory. While the gap equation includes the summation of ladder diagrams^{4, 52} and can, at least in principle, deal with bare hardcore interactions, the particlehole reducible diagrams described by introduce new physics.
Ignoring the induced interaction leads to the twobody approximation
(43) 
for the pairing interaction. We note that in this case one can again identify .
Fig. 6 demonstrates the impact on the calculated energy gap of the two approximations identified above, for the case of the Reid potential. Evidently, both simplifications have rather dramatic effects, being enhanced by the nominally exponential dependence of the gap on the pairing interaction. At this point we are not prepared to describe or affirm any systematics of the effects of these approximations. However, the close proximity of the “full CBCS” results and those for the bare interaction shown in Fig. 3 would seem to be coincidental, stemming from competing corrections.
3.4 Comparison with Previous Gap Calculations
The work we report represents the most rigorous calculation yet performed for nuclear systems within correlated BCS theory. It is therefore of special interest to compare its results with those of earlier calculations of the pairing gap for neutron matter based on microscopic manybody theories, where meaningful conclusions might be drawn.
Informative comparison of the predictions of previous gap calculations – as represented for example by the aforementioned summary figure in the review by Gezerlis et al.^{24} – is rendered problematic by the diversity of methods applied, interactions adopted, and assumptions made (e.g., whether or not polarization effects from exchange of density and/or spindensity fluctuations are included). Nevertheless, some specific and nonspecific comparisons may be useful.
Fig. 6 includes data plotted for a pureBCS calculation in which the pairing interaction is the bare potential in the channel of the Argonne interaction, used along with free singleparticle . The BCS result for Argonne , calculated by the separation method of Ref. 9, was taken from Ref. 53. (For this present purpose, the distinction between the original Argonne interaction and Argonne should be immaterial.) Corresponding bareBCS results for the Reid choice (displayed in Ref. 9 but not plotted here) are very close to those shown for Argonne , as expected. What is unexpected is that our CBF results for the Argonne case show only a modest suppression (about 15%) of the maximum, which occurs slightly above in both calculations. The approximately Gaussian shape of vs. shifts to slightly lower in the absence of JastrowFeenberg correlations. It is obvious from Fig. 6 that this near concurrence cannot be attributed to unimportance of the correlations introduced in the CBF treatment. It is possible that this feature is due in part to the presence of the inducedinteraction term in the effective pairing interaction coming from density fluctuations, which is expected to enhance the gap relative to that given by direct part of .
The most recent microscopic calculations of the gap incorporating JastrowFeenberg twobody correlations are those of Pavlou et al.^{22}, described briefly in Sec. 1.3. Their variational CBF study was carried out for each of two different parametrized forms the Jastrow factor , subject to a constraint on its “wound,” as outlined briefly in Sec. 1.3. Both of these forms have been used in earlier work: one, referred to as the “Benhar” choice, has one free parameter but allows to overshoot unity, whereas the “Davé” choice has two free parameters but no overshoot. Restricted minimization was performed on an approximation to the energy expectation that retains only the leading (zeroth) order of its cluster expansion, neglecting terms of first and higher orders in a dimensionless small parameter that grows with density. (While its value remains well below 0.05 in the relevant density regime, the implied rate of convergence of the expansion of does not extend to approximants of .)
The approach adopted in Ref. 22 may be considered the simplest implementation of correlated CBF theory. The effective pairing interaction it generates differs from its FHNCEL counterparts in two essential respects: it lacks precisely those ingredients that are the subjects of the above discussion of the “manybody consequences” of our work based on FHNCEL theory, namely the energy numerator term and the induced interaction entering the effective interaction .
This same statement applies to the variational component of the correlated BCS approach applied much earlier by Chen et al.^{8}, in which the gap in neutron matter was estimated based on the central, part of the Reid softcore interaction. In that case was found to peak at about with a maximum value close to 3.2 MeV, a result based in fact on the Davé form for . It should then be no surprise that with negligible differences, Pavlou et al. obtained to almost exactly the same result for versus , although the Argonne interaction was assumed. It should be said that all of the tests we have made support the assertion that, when the gap is calculated by the same method with different inputs for the bare NN interactions but otherwise the same assumptions (e.g., for the single particle energies ), virtually identical results will be obtained , provided the NN interaction chosen reproduces the NN scattering data up to laboratory energies relevant for below about . Indeed, this a well established property for the BCS gap^{54}.
As already pointed out, the maximum gap value obtained by Pavlou et al. with AV for their two optimized correlation functions differ by nearly a factor two (3.3 MeV for the Davé form at and 1.8 MeV at for the Benhar choice) – reflecting the extreme sensitivity of the gap to inputs for the effective interaction. Recognizing that the induced interaction and energynumerator terms are absent in these two calculations, the information provided by Fig. 6 on the relative contributions of these additional terms suggests that the results obtained in Ref. 22 for the Benhar correlation function are to be favored over those for the Davé form.
Turning to other microscopic calculations designed to provide accurate predictions for the gap in neutron matter, we first single out the study of Cao, Lombardo, and Schuck ^{55}, carried out within the framework of Brueckner theory. Meanfield theory for the superfluid state, as represented by the pureBCS treatment, was modified by replacement of the bare pairing interaction with a proper vertex part, which includes an induced interaction describing the competition between the attractive density excitations and their repulsive spindensity counterparts (i.e., screening or polarization corrections). Inmedium corrections were also introduced into the selfenergies , corresponding to both dispersion and Fermisurface depletion. The quenching of the gap due to exchange of spindensity fluctuations was found to be less extreme than indicated by some previous studies. Results based on a free singleparticle spectrum were also reported, allowing more direct comparison with our results. For the free spectrum, Cao et al.^{55} find a maximum of about 2.7 MeV occurring close to , based on the Argonne interaction. The close agreement with our CBCS results shown in Fig. 3 is remarkable, but provocative, as our treatment does not include an induced interaction term corresponding to spindensity fluctuations. On the other hand, the treatment of screening effects in the two approaches is not directly comparable. For a recent intensive computational analysis of medium polarization in asymmetric nuclear matter, see Ref. 56.
The auxiliaryfield diffusion Monte Carlo algorithm (AFDMC) purports to yield accurate results for pairing gaps in neutron matter and other manyfermion systems ^{44, 57}. For the present considerations, this algorithm has two distinctive features.

Unlike the BCS state, the trial wave function that is propagated in imaginary time describes a definite number of particles , even or odd. The part of that describes pairing is essentially the projection of the BCS state onto the particle Hilbert space, which is a Pfaffian. Correlations are otherwise introduced into by a Jastrow factor.

The pairing gap is constructed as a difference of energies obtained for different particle numbers,
(44) 
Whatever the merits and deficits of this approach, they are generally different from those of traditional manybody theory; consequently, AFDMC tends to be regarded as an essentially independent arbiter in judging the quality of such traditional method when comparison can be made. The number of data points shown in Ref. 57 for the gap in neutron matter do not allow a precise identification of the peak value of , but it lies slightly above 2 MeV, reached slightly above . The error bars shown are roughly half an MeV. We present this result only to provide a balanced perspective on the current status of the problem, but restrain from drawing any conclusions about its bearing on the quality of our calculations.
4 Summary and Prospects
In this paper we have described new calculations of the pairing gap in the partialwave channel. Our findings have been analyzed and discussed in the preceding section.
The most interesting result of previous^{18} work along these lines is the appearance of a divergence of solutions of the FHNCEL equations that occurs, as a function of potential strength, well before the divergence of the vacuum scattering length of the interaction potential. This divergence of solutions of the FHNCEL equations is analogous to the spinodal instability often discussed in earlier literature, with the principled and practical conclusion that the FHNCEL equations for the homogeneous system have no solutions if , i.e. if the system is unstable in the particlehole channel. In Ref. 18, divergence of the FHNCEL equations in the case of a diverging inmedium scattering length gave evidence that the ground state is unstable against dimerization. The appearance of such instabilities whenever the assumptions on the state of the system fail – here, assumption of a nondimerized phase; in the case of particlehole instabilities, of a uniform system – is a unique feature of theories such as FHNCEL that enjoy the topological completeness of parquet diagrams.
In the calculations being reported, we have not encountered such an instability, which could be taken as evidence that mediumdriven formation of dineutrons in lowdensity neutron matter does not occur, or, in current terminology, that a BECBCS crossover ^{58, 59, 60, 61} does not take place. This is remarkable in view of the fact that, at low densities (), the gap reaches 0.45 times the Fermi energy which is comparable to what is found in the unitary Fermi gas at the BCSBEC crossover^{62}. It should be mentioned that a recent study ^{63} of the phase diagram of spinpolarized neutron matter revealed signatures that can be interpreted ^{52} as a precursor of such a crossover.
There are four areas where the present calculation can be improved:
(a) As pointed out above, the FHNCEL method sums all ring and ladder diagrams. It does that, however, in a “collective approximation” of the particlehole and the particleparticle propagators ^{32} that treats the correlations between particles within the Fermi sea in an average way. Since pairing occurs between particles at the Fermi surface, it must be examined to what extent the average treatment of correlations is appropriate. The route to improve upon this aspect is well charted within CBF theory, and earlier studies ^{64, 65} have demonstrated that CBF corrections to the pairing matrix elements can indeed be significant.
(b) Related to (a): whereas the effect of density fluctuations (exchange of virtual phonons) has been included in the CBF pairing interaction in an averagepropagator sense, effects of spindensity fluctuations are not taken into account. Based on Landau parameters and some microscopic efforts ^{8, 66, 67, 68}, density fluctuations produce a modest enhancement of the pairing gap, whereas the spindensity channel generates a dominant suppression. Without introducing explicit spindependent correlations into the basis functions of the CBF treatment, their perturbative treatment within the CBF framework would be required.
(c) In the present work, inmedium effects on the selfenergy input to the gap equation have not been pursued quantiatively. This shortcoming warrants further attention in subsequent applications of correlated BCS theory.
(d) The most severe approximation made in this work is the use of stateindependent correlation functions, albeit the twonucleon interaction is exquisitely statedependent. Introduction of a correlation operator in Eq. (10) that contains spin, isospin, tensor, and more complicated operators in the twobody correlation vehicle figuratively opens Pandora’s box. This complexity has been largely dealt with in rather simple approximations that either completely omit commutator terms ^{69, 70} or in a “singleoperatorchain” approximation ^{71}, which only sums the ring diagrams of statedependent correlations. Unfortunately, for modern nucleonnucleon interactions, which may have different core sizes in the singlet and triplet channels, the contributions of commutator diagrams can be huge ^{72}. It remains to be seen how important these effects are in the problem considered here, but at higher densities they can be decisive.
Acknowledgements.
This work was supported, in part, by the College of Arts and Sciences, University at Buffalo SUNY, and the Austrian Science Fund project I602 (to EK). JWC acknowledges support from the McDonnell Center for the Space Sciences, and expresses gratitude to the University of Madeira and its Center for Mathematical Sciences for gracious hospitality during periods of extended residence.Footnotes
 journal: Journal of Low Temperature Physics
References
 J. Bardeen, L.N. Cooper, J.R. Schrieffer, Phys. Rev. 108, 1175 (1957)
 R. Broglia, V. Zelevensky, Fifty Years of Nuclear BCS (World Scientific, Singapore, 2013)
 A. Bohr, B.R. Mottelson, D. Pines, Phys. Rev. 110, 936 (1958)
 L.N. Cooper, R.L. Mills, A.M. Sessler, Phys. Rev. 114, 1377 (1959)
 J.R. Schrieffer, Theory Of Superconductivity (Advanced Books Classics), revised edn. (Perseus Books, 1999)
 R.L. Mills, A.M. Sessler, S.A. Moszkowski, D.G. Shankland, Phys. Rev. Lett. 3, 381 (1959)
 R.V. Reid, Jr., Ann. Phys. (NY) 50, 411 (1968)
 J.M.C. Chen, J.W. Clark, R.D. Davé, V.V. Khodel, Nucl. Phys. A 555, 59 (1993)
 V.A. Khodel, V.V. Khodel, J.W. Clark, Nucl. Phys. A 598, 390 (1996)
 J.W. Clark, C.H. Yang, Lett. Nuovo Cimento 3, 272 (1970)
 C.H. Yang, J.W. Clark, Nucl. Phys. A 174, 49 (1971)
 C.H. Yang, Theory of superfluidity in stronglyinteracting many fermion systems with applications to pure neutron matter and symmetrical nuclear matter. Ph.D. thesis, Washington University (1971)
 J.W. Clark, in Progress in Particle and Nuclear Physics, vol. 2, ed. by D.H. Wilkinson (Pergamon Press Ltd., Oxford, 1979), pp. 89–199
 N.C. Chao, J.W. Clark, C.H. Yang, Nucl. Phys. A 179, 320 (1972)
 E. Krotscheck, in Introduction to Modern Methods of Quantum Many–Body Theory and their Applications, Advances in Quantum Many–Body Theory, vol. 7, ed. by A. Fabrocini, S. Fantoni, E. Krotscheck (World Scientific, Singapore, 2002), pp. 267–330
 E. Krotscheck, R.A. Smith, A.D. Jackson, Phys. Rev. B 24, 6404 (1981)
 E. Krotscheck, J.W. Clark, Nucl. Phys. A 333, 77 (1980)
 H.H. Fan, E. Krotscheck, T. Lichtenegger, D. Mateo, R.E. Zillich, Phys. Rev. A 92, 023640 (2015)
 J.W. Clark, A. Sedrakian, M. Stein, X.G. Huang, V.A. Khodel, V.R. Shaginyan, M.V. Zverev, Journal of Physics: Conference Series 702, 012012 (2016)
 L. Gorkov, T.K. MelikBarkhudarov, Sov. Phys. JETP 13, 1018 (1961)
 H. Heiselberg, C.J. Pethick, H. Smith, L. Viverit, Phys. Rev. Lett. 85, 2418 (2000)
 G.E. Pavlou, E. Mavrommatis, C. Moustakidis, J.W. Clark, Eur. Phys. J. A 53, 96 (2017)
 R.B. Wiringa, V.G.J. Stoks, R. Schiavilla, Phys. Rev. C 51, 38 (1995)
 A. Gezerlis, C.J. Pethick, A. Schwenk, in Novel Superfluids, vol. 2, ed. by K.H. Bennemann, J.B. Ketterson (Oxford University Press, 2014), chap. 22
 R.B. Wiringa, S.C. Pieper, Phys. Rev. Lett. 89, 182501 (2002)
 J.W. Clark, in Fifty Years of Nuclear BCS, ed. by R.A. Broglia, V. Zelevinsky (World Scientific, Singapore, 2013), chap. 27, pp. 360–376
 A.D. Jackson, A. Lande, R.A. Smith, Physics Reports 86(2), 55 (1982)
 A.D. Jackson, A. Lande, R.A. Smith, Phys. Rev. Lett. 54, 1469 (1985)
 R.F. Bishop, in Condensed Matter Theories, vol. 10, ed. by M. Casas, J. Navarro, A. Polls (Nova Science Publishers, Commack, New York, 1995), vol. 10, pp. 483–508
 E. Krotscheck, Phys. Lett. A 190, 201 (1994)
 E. Feenberg, Theory of Quantum Fluids (Academic, New York, 1969)
 E. Krotscheck, J. Low Temp. Phys. 119, 103 (2000)
 J. Egger, E. Krotscheck, R.E. Zillich, J. Low Temp. Phys. 165, 275 (2011)
 E. Krotscheck, Ann. Phys. (NY) 155, 1 (1984)
 E. Krotscheck, J.W. Clark, Nucl. Phys. A 328, 73 (1979)
 S.T. Beliaev, in Lecture Notes of the 1957 Les Houches Summer School, ed. by C. DeWitt, P. Nozières (Dunod, 1959), pp. 343–374
 S. Fantoni, Nucl. Phys. A 363, 381 (1981)
 A. Fabrocini, S. Fantoni, A.Y. Illarionov, K.E. Schmidt, Phys. Rev. Lett. 95, 192501 (2005)
 A. Fabrocini, S. Fantoni, A.Y. Illarionov, K.E. Schmidt, Nucl. Phys. A 803, 137 (2008)
 C.J. Pethick, H. Smith, BoseEinstein Condensation in Dilute Gases, second edition edn. (Cambridge University Press, Cambridge, UK, 2008)
 A.L. Fetter, J.D. Walecka, Quantum Theory of ManyParticle Systems (McGrawHill, New York, 1971)
 O. Benhar, G.D. Rosi, G. Salvi, J. Low Temp. Phys. (2017). (this volume, arXiv:1305.4659)
 B.D. Day, Phys. Rev. C 24, 1203 (1981)
 S. Gandolfi, A. Illarionov, K.E. Schmidt, F. Pederiva, S. Fantoni, Phys. Rev. C 79, 054005 (2009)
 M. Baldo, A. Polls, A. Rios, H.J. Schulze, I. Vidaña, Phys. Rev. C 86, 064001 (2012)
 G. Brown, J. Gunn, P. Gould, Nucl. Phys. 46, 598 (1963)
 P. Quentin, H. Flocard, Annual Review of Nuclear and Particle Science 28, 523 (1978)
 E. Krotscheck, R.A. Smith, A.D. Jackson, Phys. Lett. B 104, 421 (1981)
 G.E. Brown, C.J. Pethick, A. Zaringhalam, J. Low Temp. Phys. 48, 349 (1982)
 B.L. Friman, E. Krotscheck, Phys. Rev. Lett. 49, 1705 (1982)
 E. Krotscheck, J. Springer, J. Low Temp. Phys. 132, 281 (2003)
 P. Noziéres, S. SchmittRink, J. Low Temp. Phys. 59, 195 (1985)
 L. Yuan, Threebody pairing interaction effect on superfluidity with applications to neutron star matter. Ph.D. thesis, Washington University (2011)
 Ø. Elgarøy, M. HjorthJensen, Phys. Rev. C 57, 1174 (1998)
 L.G. Cao, U. Lombardo, P. Schuck, Phys. Rev. C 74, 064301 (2006)
 S.S. Zhang, L.G. Cao, U. Lombardo, P. Schuck, Phys. Rev. C 93, 044329 (2016)
 S. Gandolfi, A. Illarionov, F. Pederiva, K.E. Schmidt, S. Fantoni, Phys. Rev. C 80, 045802 (2009)
 J. Margueron, H. Sagawa, K. Hagino, Phys. Rev. C 76, 064316 (2007)
 F. Isaule, H.F. Arellano, A. Rios, Phys. Rev. C 94, 034004 (2016)
 M. Stein, A. Sedrakian, X.G. Huang, J.W. Clark, Phys. Rev. C 90, 065804 (2014)
 V.A. Khodel, J.W. Clark, V.R. Shaginyan, M.V. Zverev, Physics of Atomic Nuclei 77, 1145 (2014)
 M. Randeria, E. Taylor, Annual Review of Condensed Matter Physics 5(1), 209 (2014)
 M. Stein, A. Sedrakian, X.G. Huang, J.W. Clark, Phys. Rev. C 93, 015802 (2016)
 A.D. Jackson, E. Krotscheck, D. Meltzer, R.A. Smith, Nucl. Phys. A 386, 125 (1982)
 J.M.C. Chen, J.W. Clark, E. Krotscheck, R.A. Smith, Nucl. Phys. A 451, 509 (1986)
 J.W. Clark, C. G.Källman, C.H. Yang, D.A. Chakkalakal, Phys. Lett. B 61(4), 331 (1976)
 J. Wambach, T. Ainsworth, D. Pines, Nucl. Phys. A 555, 128 (1993)
 H.J. Schulze, J. Cugnon, A. Lejeune, M. Baldo, U. Lombardo, Phys. Lett. B 375, 1 (1996)
 S. Fantoni, S. Rosati, Nuovo Cimento 43A, 413 (1977)
 R.B. Wiringa, V.R. Pandharipande, Nucl. Phys. A 299, 1 (1978)
 V.R. Pandharipande, R.B. Wiringa, Rev. Mod. Phys. 51(4), 821 (1979)
 E. Krotscheck, Nucl. Phys. A 482, 617 (1988)