The $^1$S$_0$ pairing gap in neutron matter

We report ab initio calculations of the S wave pairing gap in neutron matter calculated using realistic nuclear Hamiltonians that include two- and three-body interactions. We use a trial state, properly optimized to capture the essential pairing correlations, from which we extract ground state properties by means of auxiliary field diffusion Monte Carlo simulations. We extrapolate our results to the thermodynamic limit by studying the finite-size effects in the symmetry-restored projected Bardeen-Cooper-Schrieffer (PBCS) theory and compare our results to other ab initio studies done in the past. Our quantum Monte Carlo results for the pairing gap show a modest suppression with respect to the mean-field BCS values. These results can be connected to cold atom experiments, via the unitarity regime where fermionic superfluidity assumes a unified description, and they are important in the prediction of thermal properties and the cooling of neutron stars.

Strongly paired Fermi systems offer a unique regime for quantum many-body physics as their relevance spans many physical settings of various scales: From the structure of neutron stars (NSs) and the physics of neutron-rich nuclei to cold-atom experiments. Neutron matter (NM), one of the most strongly interacting Fermi systems found in nature, is an important ingredient of NSs, playing an essential role in their structure [1] while strongly interacting fermionic atoms are now routinely used in experiments shedding light on the properties of strongly interacting superfluids [2]. While it initially appears different, the description of these systems can be unified via their proximity to the unitary Fermi gas, connecting atomic experiments on Earth to the NS matter.
Neutrons found in the inner crust of quiescent NSs are known to form 1 S 0 pairs, turning low-density NM to a S wave superfluid [1,3,4]. The correct description of such neutron fluids is integral to the understanding of NS physics. Properties of low-density NM can explain observations such as irregularities in the periods of NSs and their cooling [5][6][7][8], while the equation of state (EoS) of high-density NM impacts the mass-radius relations of NSs [1,9] and the hydrodynamic description of their inner crust [10]. Neutron matter of the same densities is also found on the exterior of neutron-rich nuclei [1,4]. Therefore, a correct description of low-density NM is crucial to our understanding of nuclear systems of various sizes.
Strongly interacting cold Fermi atoms have been the subject of many theoretical investigations [11][12][13][14][15][16][17][18]. Experimentally, they have been studied extensively since the beginning of the century, owing in part, to the simplicity of these experiments compared to those for their bosonic counterparts [2]. In these cold atom experiments, the strength of the interaction can be tuned through Feshbach resonances to yield a specific scattering length. Many experimental studies of strongly interacting Fermi gases utilize a 6 Li gas, which exhibits a very broad Feshbach resonance, with a vanishing effective range r e [19]. This allows one to perform studies of atomic superfluids from close to the Bardeen-Cooper-Schrieffer (BCS) limit (small arXiv:2201.01308v2 [nucl-th] 4 Feb 2022 −k F a) to unitarity (−k F a 1), where k F is the Fermi momentum and a is the scattering length of the inter-particle interaction. For a comprehensive review of experimental techniques of cold Fermi atoms at unitarity, see Ref. [2].
The neutron-neutron (NN) interaction is in principle a very complicated one. At large distances, it is described by the exchange of a pion, and at intermediate distances it is spindependent and attractive, dominated by a two-pion exchange, and it turns repulsive at short distances. However, at the very low densities found in the inner crust of NSs or the exterior of neutron-rich nuclei, the long inter-particle distance allows us to ignore most of these details, and capture the physics of the system with the scattering length and effective range of the NN interaction, while the repulsive core only guarantees that the system does not collapse to a higher-density state. This is known as the shape independence and it asserts that, at low energies, the two-body scattering phase-shift δ 0 can be described solely by a and r e , The parameters a and r e can be thought of as the strength and the range of the interaction that drives the low-energy scattering. Therefore, for low densities, we can model the NN interaction with any potential that reproduces the correct scattering length and effective range.
The unitarity limit refers to Fermi gases with an attractive interaction strong enough to create a bound state of vanishing bound energy. In terms of the scattering length and effective range, this corresponds to −k F a = ∞ and k F r e = 0 making the inter-particle distance, introduced by the density, the only length scale of the system. Therefore, Fermi gases are expected to exhibit universal behavior at this limit [20]. Past unitarity, for k F a > 1, lies the Bose-Einsein condensation (BEC) regime, where pairing has grown strong enough to create bound molecular states. At a zero temperature, pairing is present throughout the three regimes, with unitarity being a smooth crossover between the BCS and BEC regime. For a discussion on the BCS-BEC crossover and how it can be used to connect cold atoms to nuclear systems, see Refs. [12,21].
While the scattering length of cold atoms can be tuned arbitrarily close to unitarity (and beyond), NM comes with fixed a and r e . In the densities considered here, the NN interaction exhibits a very long scattering length a ≈ −18.7 fm being attractive enough to almost form bound states (dineutrons), while its effective range is finite, but much smaller than the scattering length, at r e ≈ 2.7 fm. While the scattering length is much larger than the inter-particle spacing for low-density NM, very low densities are needed to turn the inter-particle spacing to be much larger than r e . Therefore, to the extent that the effects of a finite range can be neglected, neutron matter and cold atoms exhibit 'universal' behavior; their properties depend only on the product of the Fermi momentum with the scattering length. This can be clearly seen in Figure 1, where for low densities (and large inter-particle distances) NM and cold-atoms produce identical pairing gaps. This changes as the density increases and the effects of the finite range NN interaction start becoming important.
The exponential suppression of pairing on the BCS limit allows for a mean-field description of superfluidity, namely the BCS theory (see Section 1). While the mean-field approach can be applied to stronger pairing as well, it is expected to be only qualitatively correct. For an accurate description of strongly paired systems, such as the unitary gas, we have to look past the mean-field approach. One way is paved by the various beyond-mean field theories, which include correlations neglected by the mean-field treatment (e.g., see Refs. [22][23][24][25][26][27]). Another path are ab initio approaches where the problem is tackled from first principles without any uncontrolled approximations. We have performed auxiliary field diffusion Monte Carlo (AFDMC) calculations of the ground-state properties of superfluid NM for an improved variational wavefunction. We compare our results to previous exact diffusion Monte Carlo (DMC) calculations for very low-density NM and we extend our results to higher densities. Finally, we compare our results to these predicted by the BCS theory and its particle-conserving version, for the same densities. The rest of the manuscript is organized as follows: In Section 1, we present a brief overview of the BCS theory and related techniques which provide a qualitative understanding of pairing, in Section 2 we present the overview of the AFDMC approach that was used and a comparison with the DMC in past calculations, and we conclude with Sections 3 and 4, presenting and discussing our results for the pairing gap and the EoS of low-density NM.

Overview of BCS Theory
The BCS theory describes superfluid fermionic gases as condensates of Cooper pairs. The ground state of these systems is the fully paired condensate that minimizes the free energy corresponding to the pairing Hamiltonian, i.e., the one neglecting all normal state interactions: Here, k signifies the single-particle energy of a free particle with momentum k, in a box of length L under Periodic Boundary Conditions (PBC), and the matrix elements V kl are those of the pairing interaction, the attractive interaction responsible for the creation of the Cooper pairs. The BCS ground state for systems with even and odd particle numbers, can be written as: where, the distribution v k (u k ) is the probability amplitude for (not) finding a pair with momenta and spin k ↑, −k ↓ and so it is normalized to unity, u 2 k + v 2 k = 1. Hence, the distribution v k completely defines the BCS ground state and it is determined so that the state in Equation (3) minimizes the free energy. This condition leads to the BCS gap equations which determine the gap function ∆ k , corresponding to half the binding energy of a k-pair, and the quasiparticle excitation energy, where µ is the chemical potential. These in turn define the A key signature of pairing correlations is the existence of a minimum non-zero energy required to create an excitation by breaking apart a pair. This is called the pairing gap and, in the BCS theory, it corresponds to the minimum of the quasi-particle excitation energy: Here the subscript "MF" refers to the mean-field nature of the BCS theory. This can be seen by performing a Bogolyubov transformation and describing the theory in terms of non-interacting quasi-particles moving in an average field [28].

BCS Theory for NM and Cold Atoms
The BCS theory can provide a phenomelogical description of the 1 S 0 pairing in lowdensity NM, where neutrons form spin-singlet (S = 0) pairs. In such a description, we can model the neutron-neutron (NN) interaction with a simple potential that is tuned to reproduce the 1 S 0 scattering length (a ≈ −18.5 fm) and effective range (r e ≈ 2.7 fm) of the NN interaction: Such a potential should produce indistinguishable results for low-density NM where the details of the functional form of the interaction are irrelevant and the physics is captured by a and r e . We choose the form of a modified Pöschl-Teller potential: where the parameters λ and β are tuned to reproduce the 1 S 0 scattering length and effective range of NM. This potential has been used successfully before in phenomenological studies of pairing in NM [29][30][31].
For cold atoms, the interaction's scattering length can be tuned using Feshbach resonances [2]. For low densities, the functional form of the interaction is irrelevant and the inter-atomic potential can be modeled by Equation (7). To study atoms close to unitarity, we tune the parameters β and λ to yield a small effective range (smaller than the inter-particle distance) and the appropriate scattering length. Therefore we can treat NM and cold atoms, at low densities, on equal footing when formulating a description within the BCS theory.
Since we aim for a description of the 1 S 0 pairing, all BCS equations need to be expanded in partial waves where only the S wave is to be kept, leading to the angle-averaged BCS gap equations, for even systems, and for odd systems, and the angle-averaged version of the ground state energies, In this angle-averaged expression of BCS, the population function M(k) counts the number of k states in the momentum shell with |k| = k, and all functions of k have been replaced by their angle-averaged counterparts which are functions of just the momentum's magnitude k. It should be noted that, in the expressions for odd systems, only a single state with momentum b is excluded from the sum rather than the entire |k| = b shell. The potential V 0 is that of the 1 S 0 channel of the potential in Equation (7), i.e., the S wave term in its partial wave expansion, V 0 (k, p) = ∞ 0 drr 2 j 0 (kr)V(r)j 0 (pr), where j 0 is the zero-th order spherical Bessel function. For a complete derivation of the angle-averaged version of BCS, see Ref. [29].
While the formulation of BCS presented until this point describes a finite number of neutrons or atoms in a cubic box of length L, a formulation for a superfluid system at the Thermodynamic Limit (TL), such as NM found in the inner crust of NSs, can be retrieved by taking the limit of L → ∞ while keeping the particle density constant, i.e., n = N /L 3 = const. The distinction between even and odd systems is irrelevant at the TL therefore Equations (8), (9), and (12), and Equations (10), (11), and (13) have the same large-N limit.

PBCS: Particle-Number Projected BCS
The BCS ground state in Equations (3) and (4) does not conserve the particle number and, therefore, it describes a linear combination of states with an average particle number of N . Quantum Monte Carlo calculations deal with states of a well-defined particle number. To connect such calculations to values at the TL, we need to study the superlfuid FSE in a particleconserving theory. Since, particle-number conservation is a symmetry of the Hamiltonian (cf. Equation (2)), we can use symmetry-restoration techniques to restore this symmetry in the BCS ground state, which amounts to projecting out of the ground state the component that corresponds to the proper eigenstate of the number operator [32,33]: Here, C and C(b) are the normalization of the projected BCS state for even and odd systems, respectively. In principle, one must treat the states in Equations (14) and (15) as variational wavefunctions and determine the distributions u k and v k that minimize the energy of each state, an approach called variation after projection (VAP) or full BCS (FBCS). An alternative route is to use the distributions u k and v k that solve the BCS gap equations, which amounts to the projection after variation (PAV) or projected BCS (PBCS). The PBCS ground states are an approximation of the FBCS ones, in the sense that the distributions u k and v k are not optimized to yield minimum energy. However, the error introduced is small for strongly paired systems [33], such as NM. Dealing with S wave superfluidity, we must again isolate the S wave terms from a partial wave expansion, and so the energy of the PBCS ground states for even and odd systems are: where we defined the residuum integrals, The energy per particle calculated with Equations (16) and (17), or Equations (12) and (13), for a given density is not equal to the TL value of the energy, as is the case with any intensive quantity of a finite system. This is known as the finite-size effects (FSE) and they can be seen in the left panel of Figure 2. By studying the FSE of intensive quantities, we can prescribe extrapolation schemes to the TL and estimate their accuracy. In the case of superfluidity, pairing tends to create smooth FSE for the energy with no abrupt changes, compared to the FSE of the free Fermi gas. This can be most clearly seen in a comparison of the superfluid kinetic energy with that of the free Fermi gas, plotted in Figure 3. This can be attributed to the pairing's smearing of the Fermi surface (see Figure 4). When studying the FSE, the particle-number projection of PBCS should be seen as separating the contribution of a system with N neutrons from the linear combination of systems with average particle-number N , that is the BCS state. As such, the PBCS curves in Figure 2 and 3 represent more well-defined FSE curves compared to the BCS ones. A detailed study of the FSE in BCS and PBCS for NM was carried out in Ref. [29]. The prescription of the pairing gap defined in the BCS theory in Equation (6) cannot be applied in the PBCS theory. Alternatively, one can define the odd-even staggering (OES), inspired by the odd-even mass staggering of nuclei. It has been demonstrated that for the 1 S 0 pairing gap in NM, ∆ MF and ∆ are probing the same physical quantity and Equations (6) and (19) can be used interchangeably [29], as shown in the right panel of Figure 2. Since the pairing gap is an intensive quantity, it generally suffers from larger FSE than the energy (see the right panel of Figure 2).

Ab Initio: DMC and AFDMC
While the mean-field description, given by BCS, provides a qualitative understanding of strongly paired systems, accuracy demands that we treat pairing correlations from first principles. For NM, this can be done by numerically solving Schrödinger's equation for the nuclear Hamiltonian, to find the ground state of a finite number of neutrons, and then extrapolating the results to the TL. We employ the non-relativistic nuclear Hamiltonian: where m is the mass of the neutron, and v ij and V ijk are two-and three-body potentials. All the results presented in this paper have been obtained using the Argonne AV8' and the Urbana-IX (UIX) [34,35]. The AV8' belongs to the Argonne family of realistic two-nucleon potentials, which are generated by high-precision fitting of experimental scattering data. The functional form of the AV8', contains eight two-particle operators: The four central components 1, τ i · τ j , σ i · σ j , (τ i · τ j )(σ i · σ j ), the tensor S ij and the tensor-τ S ij (τ i · τ j ) components, the spin-orbit L ij · S ij , and the spin-orbit-τ L ij · S ij (τ i · τ j ) components (where S ij = 3(σ i ·r ij )(σ j ·r ij ) − (σ i · σ j ) and L ij , S ij are the relative angular momentum and the total spin of the particle ij). The UIX is a three-body potential, which describes the exchange of two pions between three nucleons via a spin-isospin dependent term [36] and it is fit to reproduce the correct triton energy in Green's function Monte Carlo calculations and the expected saturation energy of nuclear matter in the Fermi hypernetted-chain approximation [35]. The remaining term is a phenomenological part that sums other neglected terms. We have also considered two-and three-body local interactions constructed within chiral effective field theory [37,38], but the results are very similar. This is expected as this study is dedicated to low-density NM. The ground-state is calculated using the AFDMC method, and then compared to earlier ab initio results obtained by the DMC method. Both of these methods are members of the QMC family of ab initio approaches, which solve the many-body Schrödinger's equation stochastically. The DMC approach is a projector method, which uses imaginary-time propagation to extract the ground-state of a Hamiltonian from a trial state. The method relies on the fact that, in imaginary time, Schrödinger's equation turns into a diffusion equation where the large-time limit of any initial state (not orthogonal to the ground state) is the ground state of the system: The speed of the convergence depends on the choice of the trial state, which should capture the qualitative features of the problem at hand. In practice, the DMC method is applied by distributing particles according to the trial wavefunction and then propagating them in space by sampling the short-time propagator. The spin of each particle is considered 'frozen', i.e., different spin-projections define different species of particles. This means that this method is most suitable for systems with spin-independent interactions. The antisymmetry of the fermionic wavefunctions creates regions where the trial wavefunction is negative, which complicates the interpretation of Schrödinger's equation as a diffusion equation. This is known as the fermion sign problem and in DMC, it is typically addressed by fixing the nodal surface of the wavefunction (Ψ = 0) to be the same as that of the trial wavefunction. This corresponds to separating the simulation space in regions where the trial wavefunction has a definite sign (nodal pockets) and evolving each one independently. Therefore, the large-time limit of the initial configuration of particles corresponds to the state with the lowest energy and the same nodal surface as the trial wavefunction. The energy of this state provides an upper-bound to the true ground state energy. The DMC method has been used in the past to calculate the 1 S 0 pairing gap of low-density NM and cold atoms [39]. For a review of DMC, see Ref. [40].
The AFDMC approach is built on the same underlying principle as the DMC method: A ground state is projected out of a trial wavefunction via imaginary-time propagation. In contrast to DMC, the AFDMC time propagation can alter the spin-projection of the particles. When done naively, the inclusion of spin degrees of freedom has an unfavorable scaling behavior. This is mediated in AFDMC by applying a Hubbard-Stratonovich transformation to the short-time propagator, expressing the action of an operator exp (−λÔ 2 /2) as: thus improving the scaling behavior at the cost of additional integrations over auxiliary fields. This makes AFDMC the method of choice for systems with spin-dependent interactions, like nuclear systems (for a more detailed description of AFDMC, see Refs. [41,42]). This method has been applied in the past to calculate the energy and pairing gap of low density neutron matter [43][44][45].
The trial wavefunction Ψ T of N neutrons in AFDMC is typically of the form: where R and S represent the 3N spatial coordinates and 3N up-and down-spin components of the neutrons, and f (r) is a two-body spin-independent correlation (see Ref. [41] for details). The antisymmetric part Φ of the trial wave function is usually given by the ground state of noninteracting Fermions (Fermi gas), which is written as a Slater determinant of single particle functions. In the case of superfluid systems, the function Φ must include pairing correlations. For spin-independent interactions, like ultra-cold Fermi gases, pairing correlations can be included by using a Slater determinant [46]. However, in the case of a realistic nuclear interaction that can exchange the spin of neutrons, a pfaffian wave function must be used instead: The details on how to calculate the pfaffian above are given in Ref. [45]. The pairing functions φ(r) are: where sum over α indicates the k-space shells of the cube with k values: for integer n x , n y , and n z , and L is the simulation box size. The function χ is the spin-singlet wave function for two neutrons: Note that if the pairing coefficients c α are zero for all |k α | > k F , the pfaffian in Equation ( 26) turns into a Slater determinant of spin-up and spin-down neutrons filling the Fermi sea, and the pfaffian form goes over to the normal liquid state.
We calculate the superfluid pairing gap and the EoS for the ground state of NM. The pairing gap is evaluated by taking the difference of the total energy of systems with even and odd particle numbers, as per Equation (19), which makes the calculation very sensitive to the error-bars of the energy. The AFDMC method involves the constrained-path approximation to control the sign problem, and the results depend upon the quality of the trial wave function used to model the system. In several cases, an unconstrained path evolution is possible, giving exact results for the ground-state energy within (usually large) statistical error bars [42]. Thus, unconstrained-path calculations of the pairing gap would give very large error bars. In addition, it is reasonable to assume that the quality of the constrained-path approximation is similar for systems with N − 1, N, and N + 1 neutrons, so the systematic uncertainties would cancel. However, the results for the energy, and in particular for the pairing gap, ultimately depend on the trial wave function.
In the old calculations of Refs. [43][44][45], the parameters c α were chosen by performing a correlated basis function calculation [43,47]. Old results showed a pairing gap very close to the one predicted by the BCS theory, and predicted a higher EoS with respect to other calculations.
In this study, we improved the trial wave function by performing a variational search of the optimal c α parameters using the stochastic reconfiguration method [48]. As will be discussed in Sections 3 and 4, this yields a state that captures the appropriate correlations of superfluid NM and thus yields a ground state with lower energy and results consistent with the accurate DMC calculations at very low-densities performed by Gezerlis and Carlson [31,39]. The pairing functions in Equation (27), once properly normalized, can be connected to a BCS state, like the one in Equation (3), through the c α parameters, This allows us to compare the pairing functions and pair probability distributions corresponding to the pairing function of Equation (27) for 40 neutrons to those of the BCS ground state in Figure 4 and 5. We find that the BCS state defined by the c α parameters yields a pair probability distribution with a slightly higher spread around the Fermi surface, describing higher pairing correlations. This is reflected in the optimized state's pair function which decreases slowly with the particle separation (see Figure 5). This suggests that a pair of particles remains correlated at larger separations. Additionally, we find that the optimized state's pair function is less isotropic compared to the BCS ground state, owing to the noncentral components of the nuclear force described by the richer operator-structure of the AV8' potential.

Equation of State
We first examine the EoS of low-density NM. In Figure 6, we present our AFDMC calculations with the better optimized trial wavefunction and we compare with the old AFDMC calculations from Refs. [44,45], where correlated basis function calculations prescribed the pairing functions, and the results from Ref. [39], where the DMC method was used instead and the pairing functions were determined by minimizing the trial state's energy. The results display the correct behavior at low densities, i.e., E/E FG → 1, where pairing becomes less prevalent. Evidently, the better optimization of the pairing functions in the trial wave function of AFDMC provides a better description of superfluid NM yielding a EoS consistent with the DMC calculations of Ref. [39] for low and intermediate densities. Our work then extends these results to higher densities. Our energy results have been obtained by simulating a system using N = 40 neutrons under periodic boundary conditions. By choosing this specific particle number, we avoid closed-shell configurations e.g., N = 38 or 66. This functions to test that the pfaffian wave function with a properly optimized pairing functions that can reproduce open-shell configurations. It should be noted that, as demonstrated in Figure 4, superfluid systems do not yield a well-defined Fermi surface, and so all mentions to closed-and open-shell configurations refer to free systems with the same particle number. Certain closed-shell configurations are seen to produce minimum superfluid FSE (see Figure 3) making them a common choice for QMC calculations [39,44]. Studying the TL with open-shell configurations would otherwise require the use of twist-averaged boundary conditions, as has been demonstrated for normal-state NM [49] and superfluid NM [30]. For some density, we changed N from 40 to 66, without finding any dependency to N. This is because for such small densities, the box size is large enough to avoid FSE due to the truncation of the interaction. This has also been verified for ultra-cold Fermi gases with strong interactions [13]. The choice of N = 40 is additionally motivated by overall efficiency: The optimization of the variational parameters becomes much harder for large N. As seen in Figure 3, the FSE of the kinetic energy are smeared, compared to that of the free Fermi gas, due to the pairing correlations, providing smoother FSE to the total energy as well, as seen in the left panel of Figure 2. From these calculations, we expect that the energy of a system with N = 40 neutrons is within ∼2.5% from its TL value.

Pairing Gap
We have also performed calculations of the pairing gap using the better optimized trial state. The pairing gap is calculated using Equation (19) for an odd number of neutrons N: The results are shown in Figure 7. Since the pairing gap is calculated as differences between the total energy of the system, the result is that the uncertainty associated grows with N making the use of large systems impractical. For most of the results presented here, we chose N = 47 to get a reliable pairing gap and to minimize any possible bias due to closed shells. This is particularly important for lower densities where the pairing is suppressed. We tested the consistency of our results by performing simulations at different numbers of neutrons up to N = 68. Additionally, based on PBCS calculations of the pairing gap, we expect the pairing gap at N = 47 to be close to the TL. In fact, the FSE, which for low densities and at small N, appear to be insensitive to density variations, suggest that the pairing gap of a system with N = 47 neutrons is within ∼ 5% from the TL (see the right panel of Figure 2).
The results of the pairing gap are presented in Figure 7 where we also show the pairing gaps of Ref. [39] calculated using DMC and the old AFDMC results of Refs. [44,45]. At low densities, the AFDMC calculations are in reasonable agreement with DMC calculations. For higher densities, the AFDMC results depart from the BCS prediction as the mean-field approximation becomes less accurate. In this work, we find a sizable pairing gap up to k F ∼ 1.3 fm −1 , in contrast to the previous AFDMC calculations that seem to drop around k F = 0.8 fm −1 because of the much simpler variatinal wave function employed. This again is a consequence of the better optimization of the pairing functions that we adopted in this work; we are capturing more of the essential correlations found in the ground state of superfluid NM. Overall, we find a moderate suppression of the pairing gap predicted by the mean-field BCS.  The pairing gap calculated at constant density with (black circles) and without (red triangles) the correction described in the text. Our calculations are compared to the old AFDMC results of Refs. [44,45], and to those of Ref. [39]. The brown dashed curve shows the pairing gap predicted by BCS.
The calculations of the energies in Equation (31) are typically performed by simulating the system at a constant density in order to minimize the FSE due to the truncation of the potential energy in the periodic box. However, for the free Fermi gas, this procedure would give a non-zero pairing gap. We corrected our results of the pairing gap by subtracting that of the free Fermi gas obtained at constant densities. The correction is negligible at low densities, as can be seen in Figure 7, and within error bars at large densities. In a few cases, we also performed simulations at a constant volume instead of constant density and the results are very similar to corrected constant-density ones.
In summary, we performed a detailed ab initio study of the S wave pairing in low-density NM found in the inner crust of cold NSs. We calculated the EoS and the pairing gap by means of AFDMC simulations for a range of densities, using a variationally optimized trial state benchmarked against previous calculations at low densities. A study of the FSE within a symmetry-restored mean-field treatment shows that the AFDMC energies and pairing gaps are within ∼2.5% and ∼5% from the TL, respectively. Our AFDMC pairing gaps show a modest suppression with respect to the mean-field BCS values. These results can be used in calculations of the thermal properties of NSs and be indirectly tested in cold atom experiments utilizing the universality of the unitary Fermi gas.
Award DBI-1565180. The work of G.P. and A.G. has been supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Canada Foundation for Innovation (CFI), and the Early Researcher Award (ERA) program of the Ontario Ministry of Research, Innovation, and Science. Computational resources were provided by SHARCNET and NERSC. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under contract no. 89233218CNA000001.

Conflicts of Interest:
The authors declare no conflict of interest.