Superfluid neutron matter with a twist

Superfluid neutron matter is a key ingredient in the composition of neutron stars. The physics of the inner crust is largely dependent on that of its $S$-wave neutron superfluid which has made its presence known through pulsar glitches and modifications on the neutron star cooling. Moreover, with recent gravitational-wave observations of neutron star mergers, the need for an equation of state for the matter of these compact stars is further accentuated and a model-independent treatment of neutron superfluidity is important. \textit{Ab initio} techniques developed for finite systems can be guided to perform extrapolations to the thermodynamic limit and attain this model-independent extraction of various quantities of infinite superfluid neutron matter. To inform such an extrapolation scheme, we performed calculations of the neutron $^1S_0$ pairing gap using the model-independent odd-even staggering in the context of the particle-conserving, projected BCS theory under twisted boundary conditions. While the practice of twisted boundary conditions is standard in solid state physics and has been used repeatedly in the past to reduce finite-size effects, this is the first time it is employed in the context of pairing. We find that a twist-averaging approach results in a substantial reduction of the finite-size effects, bringing systems with $N\gtrapprox 50$ within a $2\%$ error margin from the infinite system. This can significantly reduce extrapolation-related errors in the extraction of superfluid neutron matter quantities.


Introduction
Historically, the idea of neutron stars (NSs), i.e., ultra-compact objects where "the atomic nuclei come in close contact, forming a gigantic nucleus", was first proposed by L. Landau [1] shortly before [6] the discovery of the neutron [2]. This idea was further explored by W. Baade and F. Zwicky in two seminal publications where they identified the birth of NSs with supernova explosions, a term also coined therein. These theoretical ideas were substantiated by the observation of the first pulsar by J. Bell [5] which, interpreted by T. Gold [8] shortly after, marked the first observation of a NS. While the interest in compact stars only grew after that, a parallel thread involved nuclear superfluidity which was first proposed by A. Bohr, B. R. Mottelson, and D. Pines [9] a year after the microscopic theory of superfluidity was introduced (see Sec. 2.2) and a year before A. B. Migdal [10] remarked that similar mechanisms could take place in the interior of NSs. Today, the presence of supefluidity in NSs is a widely accepted fact (for a more detailed account on the history of NSs see Refs. [6] & [7], as well as other contributions to this special issue [11][12][13]).
While the idea of compact stars was first discussed during the beginning of the modern nuclear theory, the details of the composition of these objects came with later advancements in the physics of nuclei and the further understanding of the nuclear interaction [7]. The current understanding of NSs separates the compact object in layers of different density which also correspond to regions of different physical phenomena owing to the density-dependent nature of the nuclear interaction. The outer crust of a NS is composed of a lattice of nuclei bathed in a sea of electrons and neutrons at densities ρ ≈ 10 6 g cm −3 . Starting from 52 Fe, as one goes deeper the neutron fraction in the nuclei of the lattice increases as higher density forces an increasing number of neutrons in each nucleus. The free neutrons permeating the lattice create a state of matter similar to terrestrial binary alloys while also providing an effective attraction between nuclei resulting in clumps similar to those found in metallic alloys [14]. Reaching the innermost region of the outer crust, the composition of matter becomes uncertain largely due to the uncertainties in measurements of nuclear masses or the lack of relevant arXiv:2012.04663v1 [nucl-th] 8 Dec 2020 data [15]. The inner crust of the NS starts at densities where the neutrons start dripping out of the nuclei (ρ ≈ 4 × 10 11 g cm −3 ) resulting in a state of neutron-rich nuclei in a background of ultrarelativistic electrons and a dilute fluid of neutrons [16]. The inner crust extends to densities up to half the nuclear saturation density, ρ 0 = 2.8 g cm −3 . At these densities the 1 S 0 channel of the neutron-neutron (NN) interaction becomes attractive [17] ensuring superfluidity for the neutron fluid. With the protons mainly confined in the nuclei, in this depth one can approximate the neutron superfluid as pure neutron matter (NM) coupled to the crustal phonons or band structure [18][19][20]. At this point superfluidity exists both inside [21] and outside the neutron-rich nuclei but this boundary becomes increasingly vague as the density of the neutron fluid increases, i.e., as one reaches the bottom of the inner crust where the deformed neutron-rich nuclei mark the onset of the nuclear pasta. These extended clusters of neutrons, protons, and electrons can be shown to be the result of a balance between Coulomb and surface effects, in the context of a liquid drop model [22]. The first order transition leading to a continuous pasta phase can be considered as the beginning of the outer core of the star at a density close to nuclear saturation. Here the neutron density, while ensuring an S-wave repulsion between the neutrons, allows for attractions through other channels of the NN interaction creating PF-superfluidity while the proton density reaches that of the neutrons in the inner crust turning the proton fluid into an S-superconductor [23]. With densities ranging from 2ρ 0 to 10ρ 0 , the consistency of the star in the core is largely uncertain. A popular conjecture is the appearance of hyperons such as Σ ± , Λ, and Ξ ± , when the Fermi energy of the neutron and electron gas surpasses their corresponding rest masses. This can also lead to hyperon superfluidity if their interaction is attractive. Another popular conjecture involves quark deconfinement which, due to the uncertainty in the density marking its onset, might come before or after the formation of hyperons. As the quark degrees of freedom start becoming important one might also expect the appearance of quark-gluon superconductivity [24]. Finally, Bose-Einstein condensation of pions and kaons have also been theorized to exist in the core [25,26]. For a more detailed review on the current consensus on the composition of NS see Ref. [18] & [15].
Nuclear superfluids (and superconductors) created by pairing of nucleons in various channels, permeate most layers of a NS and are responsible for a variety of astrophysical phenomena, such as impacts on the cooling of the star [27][28][29][30][31][32] and the anomalous glitches observed in the rotation of pulsars [33,34]. The physics of the superfluids found in a NS can also affect its seismic properties and, hence, the neutrino and gravitational radiation emitted from it [18,[35][36][37]. Furthermore, pairing has direct consequences on the Equation of State (EoS) of NM [38] which in turn determines the mass-radius relation of NSs and their maximum mass. Finally, the radiation from various superfluid phases in NSs can be used for constraints on the coupling of exotic particles, proposed as extensions to the Standard Model (e.g., axions [39]), with Standard Model matter. Thus, the correct description of superfluid NM is an important step in the understanding of the physics of NSs and their connection to the cosmos.
Many models have been proposed for a concise description of superfluid NM generating a polyphony of results and a landscape of pairing gaps [18,40]. This only underlines the need for a model-independent extraction of the properties of superfluid NM. Promising candidates for such an extraction are ab initio approaches which calculate quantum many-body properties by attacking the problem from first principles. With a sizable part of these being techniques developed for tackling finite systems, one is faced with the task of creating a well-informed extrapolation scheme able to map the properties of the finite system to the infinite one, where NS matter lies. The rest of this paper is organized as follows. In section 2.1 we present a brief overview of superfluidity in NM, in section 2.2 we discuss the fundamentals of the BCS theory of superconductivity and its particle-conserving variation PBCS, and how these can form an extrapolation scheme to the TL for ab initio approaches. Finally, in section 3, we present a brief overview of the finite-size effects (FSE) in superfluid NM and we apply techniques of manipulating the periodic boundary conditions (PBC) of a finite system to decrease these FSE and further improve the extrapolation schemes already in use for NM. We present results demonstrating that approaches such as twisting the boundary conditions (BC) or averaging properties calculated with different twisted boundary conditions (TBC) can significantly improve the extrapolation to the TL.

The variety of approaches in superfluid neutron matter
Neutron matter is a strongly interacting Fermi gas making for both intriguing physics as well as systems where known weak-interaction approximations break down. Historically, the strongly interacting nature of neutron matter gave rise to the study of the unitarity regime which was conceptualized as a model for the dilute neutron gas by Bertsch [41] and Baker [42]. This is the regime where a Fermi gas, while still dilute (the interparticle distance is much larger than the range of the interaction) is very strongly interacting (the scattering length of the interaction diverges, k F a → ∞) [43]. This ensures that all length scales associated with the interaction disappear making the two-body problem scale invariant and parameter-free. Today, unitarity is known to lie in the crossover separating the BCS regime of a Fermi gas from the Bose-Einstein Condensate (BEC) regime, aptly named the BCS-BEC crossover. For vanishingly small effective range r e , the BCS-BEC crossover can be parametrized by the inverse of the coupling constant, 1/k F a. Starting on the BCS side with 1/k F a < 0, one finds particles arranged in Cooper pairs bound with energy ∆, also known as the pairing gap. Their binding increases with the coupling, reaching unitarity at 1/k F a = 0. Past the unitarity lies the BEC regime, with 1/k F a > 0, where the tightly bound Cooper pairs have become bosonic diatomic molecules and are condensed in the familiar BEC state. Neutron matter found in the inner crust of quiescent neutron stars, with its large S-wave scattering length a s ≈ −18.5 fm and relatively small effective range r e ≈ 2.7 fm, is situated close to unitarity, on the BCS side -a proximity that has lead to connections between neutron matter and cold Fermi atoms [42,[44][45][46][47][48][49][50] (for a generalization of the crossover for finite r e see [51]). In the latter, the scattering length of the interatomic interaction and its sign can be changed by tuning an external magnetic field (Feshbach resonances [52][53][54][55]) allowing one to navigate the BCS-BEC crossover and perform experiments at unitarity, where the details of the interaction are irrelevant and universal conclusions can be drawn. At the same time, the proximity of neutron matter to the unitary Fermi gas makes it the most strongly paired fermion superfluid system known in Nature with calculations and experiments suggesting pairing gaps peaking at 30 per cent of the Fermi energy [56]. The physics of these strongly-interacting superfluids is not understood nearly as well as its weakly-interacting counterparts.
The superfluid neutron matter of the inner crust does not exhibit the universal behavior expected from unitary gases entirely, despite its proximity to the unitarity regime. The reason behind this is its finite effective range [57]. The low-energy phenomena of an interacting gas depend only on the large scale characteristics of the interaction since higher energy is needed to probe finer details in the potential. This is the essence of the effective range expansion, where the physics of the 2-body phase shift, namely δ 0 are captured by the scattering length a and the effective range r e in The pairing responsible for superfluidity happens predominantly in momenta lying in a narrow band centered at the Fermi level k F and so the relevant scattering amplitudes are those involving momenta of that order. That allows us to investigate the range of validity of the effective range expansion by considering only momenta k ≈ k F . Using this, one can see that, in neutron matter interacting through the S-channel, at k F a ≈ −1 the momentum dependent terms become of the order of 1/a already setting it apart from a unitary gas. As k F a increases further, the importance of the momentum-dependent terms becomes higher reaching k F a ≈ −10 where the presence of these terms becomes vital for the correct description of neutron matter. This is not far from the interpretation of  higher-energy phenomena as probes for the finer-scale characteristics of the interaction. As the density increases, k F is pushed to higher momenta and so is the the band of momentum states relevant to pairing. With scattering in higher momentum states present, finer details of the potential become relevant distinguishing neutron matter from unitary gases. This can be seen in QMC calculations of the pairing gap in neutron matter and cold atoms [48]. Another feature of the effective range is the initiation of an effective repulsion at higher densities. The positive second term in Eq. (1) decreases in magnitude the preceding negative term whose magnitude and sign are related to the "attractiveness" of the interaction. This reduction of the attractive interaction leaves its footprint at higher densities where the gap of the superfluid state decreases after k F a -5, (see Fig. 1) instead of saturating as expected from unitary gases. Finally the attractive S-channel of the NN interaction turns repulsive at k F ∼ 1.5 fm −1 [58] instigating the closure of the pairing gap in that channel at the corresponding densities [59].
Having identified neutron matter's place in the BCS side of the BCS-BEC crossover, what follows will focus on that region. There, on the weak-coupling limit, 1/k F a → −∞, the (mean-field) BCS theory of superfluidity yields a correct qualitative and quantitative description of pairing [60] with an analytic expression for the pairing gap: When including polarization effects, the above expression acquires an extra k F a-dependent factor which tends to (1/4e) 1/3 ≈ 0.45 at the limit of 1/k F a → −∞ [61,62]. As one moves to strongly-paired superfluids, such mean-field approaches, while still valid qualitatively, fail to provide a quantitative description. The BCS description, while accurately predicting the two-body bound states that lie on the BEC side of the crossover, it does not provide precise results in that regime. For these systems one needs to turn to beyond-mean-field approaches. These approaches are built by considering correlations neglected by mean-field theories. The inclusion of short-and long-range correlations by a summation of ladder diagrams in the context of self-consistent Green's function method results in a ∼ 0.75 MeV reduction of the pairing gap compared to the BCS value, and in a closure of the pairing gap at lower densities [63]. Similar results were found by replacing the bare interaction with the G-matrix in the Bruckner theory to include screening effects [64]. In-medium corrections to the effective interaction can also be included "adiabatically" by solving the Renormalization Group (RG) flow equations in a RG approach [65] where polarization effects of the particle-hole kind result in a smaller gap compared to BCS. Similarly, reduced gaps are found when the particle-particle and particle-hole polarization effects are added by phenomenologically modifying the short-distance behavior of the bare fermion interaction [66]. A significant quenching of the gap is also found when treating short-range correlations by means of Correlated Basis Functions (CBF) [67]. The addition of three-body forces has also been explored in the context of chiral effective interactions where a reduced 1 S 0 pairing gap and an increased 3 PF 2 pairing gap were reported [68]. Finally the complete quantum many-body problem can be solved using the Quantum Monte Carlo (QMC) family of stochastic approaches where the ground-state of a many-body Hamiltonian is identified given a suitable trial wave-function [38,48]. Calculations of the 1 S 0 pairing gap by the means of QMC simulations have also resulted in a reduction of the gap as compared to the BCS value. An important distinction of the QMC methods is that one is obliged to work with finite systems and, therefore, calculations of quantities of infinite matter require the extra step of extrapolation, i.e., dealing with the FSE (see section 3). These techniques use the Rayleigh-Ritz principle to find a state with minimum energy combined with an imaginary-time propagation to purify that state of any excited-state contributions. When fermionic, the imaginary-time evolution of such states contains the additional complication of the fermion-sign problem arising from anti-symmetric states losing their anti-symmetry due to the statistical sampling inherent in the stochastic method (for a review of QMC methods see Ref [69]). With the advancements of the QMC approach over the past decade having rendered the fermion-sign problem under control, QMC results can be treated as exact yielding an error margin of typically ∼ 1% [57]. For a more detailed review on various approaches for the calculation of the pairing gap in neutron matter see Ref. [57]. The polyphony of theoretical results for the pairing gap in NM comes in the era of multi-messenger astronomy in which NSs have come in the spotlight with observations of mergers via graviational waves [70,71]. The interpretation of the data gathered from these events is in need of an EoS for NM which depends on a correct description of superfluidity in the crust. In this light, and with the strongly interacting nature of NM prohibiting a perturbative approach, a model-independent extraction of the pairing gap is a consequential issue in nuclear theory.

The BCS and PBCS Theories for Neutron Matter
The BCS theory introduced in 1957 by Bardeen, Cooper, and Schrieffer, describes superfluid (or superconducting) states as a result of an instability of the normal state in fermionic gases, with an attractive interaction, at low temperatures. In BCS, this instability, present below a transition temperature, instigates a formation of Cooper pairs. The attractive interaction responsible for the instability of the normal state and the eventual formation of the condensate of pairs, as the temperature is further decreased, can come from various sources. In the original formulation of the theory, in the context of supeconductivity in metals, the origin of the attractive interaction was the presence of an ionic lattice which allows for an effective interaction between electrons via the exchange of a phonon. In NM the nature of the attractive interaction is different, as it originates from the S-wave component of the NN interaction which is attractive at densities encountered in the inner crust of cold NSs. As mentioned in subsection 2.1, the strong nature of this interaction which is portrayed by its large scattering length a s ≈ −18.5 fm, makes NM one of the most strongly correlated superfluids encountered in Nature.
The BCS theory describes superfluidity on a mean-field level and so it is not expected to give accurate quantitative results for strongly interacting superfluids such as NM. However, one can use the BCS theory and its straightforward treatment of pairing to produce qualitative results and offer guidance to more accurate methods. Specifically, one can use results from a BCS treatment to help ab initio techniques, applied to finite systems, extrapolate to the thermodynamic limit (TL). For more details on the assistance that BCS theory can provide to these extrapolations see section 3. Here we present a brief introduction on the BCS theory, and its particle-conserving variation, namely PBCS theory, applied to NM.
One can use the BCS theory to describe a finite part of the bulk of a pure NM superfluid by enclosing a system of N 0 neutrons in a finite volume, e.g. a cubic box of length L, under PBC. The PBC require from the modulus of the wave-function to be the same on opposite sides of the box allowing for a complex phase difference between the values of the wave-function on the same points: Following the literature, we name this phase difference "twist" to avoid the polysemous word "phase". The choice of PBC is mandated by the translational symmetry that characterizes all quantities of infinite matter, and it gives rise to the familiar single-particle spectrum The twist angle θ is not an observable of the infinite system: at the limit of L → ∞ it drops out of the problem. Therefore, it is typical to assume untwisted PBC,i.e., θ = 0. However, ways to manipulate this extra degree of freedom, and bring finite systems closer to TL, have been developed. In section 3 we will use special twist-angles as well as a twist-averaging approach, where we average the quantities of a finite system over the twist angle, to make the finite system a better approximation of infinite matter. For the rest of this section we restrict our BC to the untwisted θ = 0 case,i.e., PBC. When studying the effects of pairing in NM, it is customary to isolate only the interaction responsible for the pairing and ignore normal state interactions, hence describing the system using the so-called pairing Hamiltonian: whereĉ † kσ ,ĉ kσ are fermionic creator and annihilation operators, respectively, creating or annihilating free single-particle states with momentum k and spin σ in a cubic box under PBC. Many state-of-the-art potentials have been used to model the NN interaction such as the CD-Bonn potential [72], the Nijmegen I and II potentials [73], the Argonne family of potentials (AV4, AV8, etc) [38,60], chiral potentials (NLO, N 2 LO, and N 3 LO) [74], etc. At low energies, where the details of the interaction do not matter, any potential that reproduces the 1 S 0 scattering length and effective range of NM should produce identical results, as discussed already in section 2.1. In this light, we choose the simplest path and we model the NN interaction using the purely attractive Pöschl-Teller potential, where the parameters λ and β are adjusted to reproduce the scattering length and effective range of the S-wave of the NN interaction. The matrix element in Eq. (7) is that of the 1 S 0 channel of this potential, It is worth noting that assuming a purely attractive potential, such as the PT potential, we ignore the repulsive core of the NN interaction but, again consistent with the shape independence outlined in section 2.1, for low-density NM, the exact form of the potential is irrelevant and the results produced by all the potentials mentioned above should be identical [56,60].

Even particle-number superfluid
The BCS theory describes the ground-state of a superfluid with an even number of particles N 0 as a coherent state of N 0 /2 pairs of time-reversed states: The state |0 stands for the vacuum while the distributions v k and u k represent the probability amplitude of finding or not finding, respectively, a pair of states k ↑, −k ↓, and so they are subject to the normalization v 2 Furthermore, with v 2 k representing the probability of finding a pair characterized by momentum k, it is connected to the average number of particles, This average is to be understood in the context of a grand canonical ensemble: we enforce the particle-number conservation in the box only on average. In fact, the BCS ground state in Eq. (10) is defined in a grand canonical ensemble and so it describes a superlfuid with an average particle-number equal to N 0 . Therefore, the probability distributions v 2 k and u 2 k , which characterize the pair occupation of the condensate in k-space, should be such that they minimize the free energy of the BCS state, where ξ k = ( k − µ). This minimization yields the gap equation: where every momentum-dependent quantity is replaced by an angle-averaged version. This is because only the first term of the corresponding Laplace series of each quantity couples with the S-wave of the potential which reflects the fact that momentum states scattered by an S-wave interaction can only change in magnitude and not in direction. Thus, momentum shells can be defined, each containing momenta corresponding to the same k-magnitude which are treated as identical by the interaction. The population of the shell corresponding to |k| = k is denoted by M(k) and it used in the discrete sums over k-space to allow for a one-dimensional sum over the discrete k-magnitudes instead of the entire three-dimensional momentum space. The populations of the different momentum shells depends on the twist of the BC as well, which we have taken to be the trivial θ = 0 in this section. See section 3 for the dependence of M(k) on θ and how it affects the FSE. Here the gap function ∆ k is the binding energy of a pair with momentum k and the energy needed to break a pair and create an excitation is the quasi-particle excitation energy E k . It is defined as and it defines the occupation probability distributions, Given these, Eq. (12) becomes It is worth noting that the normal state can be written in the form of the BCS ground state with v k = 1 for |k| ≤ k F and v k = 0 otherwise (this corresponds to the ∆ k = 0 solution). This means that when minimizing the free energy of the BCS state, i.e., Eq. (10), with respect to the distribution v k , the normal state is a viable candidate. In other words, if the solution of the BCS gap equations yields solutions other than the normal state solution, then a pair condensate yields lower free energy than the normal state. This is an equivalent statement of the pairing instability proposed in the initial BCS theory.
As seen in Fig. 2, the distributions v 2 k is smeared over k-space, compared to the Fermi distribution. This is a consequence of pairing: taking ∆ k → 0 in Eq. (16) one can retrieve the Fermi distribution, i.e., a free Fermi gas. This is demonstrated by the condensation amplitude which can also be seen in Fig. 2. The product v k u k is non-zero only when v k and u k are simultaneously non-zero. Hence, the spread of the condensation amplitude is strongly related to pairing with a wide spread F k characterizing a strongly paired superfluid and a normal state yielding an infinitely sharp F k around |k| = k F .
A key property of superfluidity is the energy gap in the quasi-particle excitation spectrum, e.g., see Fig. 2. The minimum energy required to break a pair and create an excitation is called the pairing gap and it is defined as Figure 3 shows the pairing gap as a function of N . Apart from the familiar oscillations seen in all intensive quantities of finite systems (see section 3), the pairing gap suffers from additional FSE due to its definition as a minimum of the discrete spectrum E k . In more detail, the quantization of momenta, imposed on finite systems by their BC, ensures that all functions of momentum be a discrete version of their TL-counterparts. Quantities defined as the minima of such discrete functions, in the manner of Eq. (21), will have to compromise with the lowest available point. The position of this available minimum might change to a neighboring k-state as N changes and thus make the N -dependence appear more random. The definition of the pairing gap, namely Eq. (21), makes it hard to compare the BCS approach to superfluidity with other techniques. A prime example of this difficulty is QMC simulations where one does not have access to the excitation spectrum of the system. Drawing from the odd-even mass staggering in nuclear physics, one can define a quantity similar to the pairing gap which, at the TL, is identical to the pairing gap: where N is an odd number of particles since it provides better decoupling of the result from the underlying mean-field theory [83], i.e., in this case, BCS. The odd-even staggering (OES) in neutron superfluidity has been demonstrated to reproduce the BCS pairing gap for systems far from the TL as well [82]. This provides OES with a phenomenological value: being a quantity readily available to a variety of ab initio approaches to superfluidity, it can be used to compare them and quantify differences between them. To have access to OES in the BCS theory via Eq. (22) one needs an expression for the energy of a system with N particles. As discussed above, the BCS treatment breaks the particle-conservation of the pairing Hamiltonian resulting in a state part of a grand canonical ensemble with a constant average particle-number. Hence, the BCS state can be understood as a superposition of eigenstates of the number operator, i.e., states describing superfluids with fixed particle number: The spectrum of eigenstates |λ N | 2 is sharply peaked around N = N 0 , as seen in Fig. 4, with a spread proportional to N −1/2 0 . In other words, the main contribution in Eq. (23) largely comes from the state ψ N 0 and that leads to the idea of the projected BCS (PBCS) theory where the ground-state of a superfluid with an exact (even) particle-number N 0 is described by ψ N 0 . The sharper the peak of |λ N | 2 the more particle-conserving the BCS ground state |ψ BCS is and the closer PBCS is to BCS. It follows that for strongly interacting superfluids, of which NM is a perfect specimen, the particle-conserving PBCS theory is expected to agree with the non-particle-conserving BCS theory.
In PBCS theory one projects out of the sum in Eq. (23) the state corresponding to N 0 = N , namely ψ N 0 . To execute this projection we need to allow v k and u k to take complex values and with the BCS ground state being a coherent state of pairs, the complex phase of v k and u k does not depend on k. Furthermore, having the freedom to define the BCS ground state up to an overall complex phase leads to The infinite product in Eq. (24) generates an infinite sum in which each term comes with as many factors of e iφ as pair creation operatorsĉ † k↑ĉ † −k↓ . Therefore, the term creating N 0 /2 pairs can be projected out by integrating the rest of the terms over a integer number of periods in φ: where the normalization R 0 0 is one of the residuum integrals: The residuum integrals are related to the probabilities of arrangements of pairs in k-space. For example, R 0 0 is the sum of the probabilities of all possible pair-arrangements in k-space and as such it is the probability of the realization of the state in Eq. (25) and, therefore, its normalization. A detailed derivation of the PBCS theory can be found in Ref. [82]. The energy corresponding to the state in Eq. (25) is where we have again used the angle-averaged distributions. Owing to the projection, E BCS even describes systems with a fixed particle number and, therefore, it can be used to calculate the gap using the OES.
The v k and u k distributions in Eq. (27) are the ones in Eqs. (16) & (17), respectively, since they originate from the already optimized, BCS ground-state in Eq. (24). One could re-optimize the PBCS state and find the distributions v k and u k that minimize the energy in Eq. (27). This corresponds to the FBCS approach [84] of which BCS is a saddle point approximation valid for 2 ∑ F 2 k 1. Following the preceding discussion on the condensation amplitude F k , it is clear to see that FBCS for NM should be qualitatively and quantatively very close to BCS and, a fortiori, to PBCS, allowing one to use the latter when probing quantities such as the OES.

Odd particle-number systems
The BCS ground state describes a condensation of pairs and as such it is not suitable for systems with an odd number of particles. With the energy of an odd particle-numbered system being central to the OES (cf. Eq. (22)), a variation of the BCS ground state is essential for the calculation of the gap. This is the blocked BCS state, where the extra, unpaired, particle occupies a momentum state b blocking it from the rest of the underlying condensate. The optimization of this state gives rise to its own blocked BCS equations where the quantities ∆(k), and E(k) retain their definitions as the gap distribution and quasi-particle excitation energy, respectively. The minimum of E(k) is defined as the BCS pairing gap for the odd particle-number systems. Following the same principles as the ones used for the projection in the even particle-number systems, one can define a fixed particle-number state for odd N. This is done by projecting out of the underlying condensate an N − 1 particle-conserving state, Finally, the energy of an odd number of particles in PBCS is the energy that corresponds to the state in Eq. (31): where, similarly to Eq. (27), we use the angle-averaged distributions along with the a slight abuse of notation: the blocking of the momentum state b does not block an entire shell but rather reduces its occupancy by 1. The unpaired particle's contribution in E PBCS odd is that of a free particle since all normal state interactions are excluded in the pairing Hamiltonian (cf. Eq. (7)). However, the presence of the extra particle is "felt" by the condensate via the blocking of the momentum state occupied by the particle. The energy in Eq. (32) depends on the blocked state b and so when used in an OES prescription, the right b should be chosen to minimize E PBCS odd . This is because the odd particle-number system in OES is to be understood as the lowest quasi-particle excitation of its even particle-number vacuum. For the sake of completeness, the corresponding energy of an odd particle-numbered superfluid in BCS is:

The Thermodynamic Limit
The treatment of the TL in BCS is straightforward underlying the worth of the theory as a tool for qualitative studies of the approach of finite systems to their TL. A system at the TL corresponds to the limiting case L → ∞ of the periodic box used for finite systems. In that limit, the discretization of the momenta becomes infinitely fine, turning all previously discretized distributions into continuous functions of k, and the sums in the corresponding equations into integrals. Hence the gap equations describing the system at the TL are: where n is the number density, and ∆(k) and E(k) are the continuous version of the gap distribution and quasi-particle excitation energy, respectively. Hence, a BCS pairing gap can be again defined in an identical way to the even and odd particle-number systems: From the discussion on the spread around the peak of the eigenstate spectrum λ 2 N it follows that for N → ∞, i.e., at the TL, the peak is infinitely sharp and the particle number projection loses its meaning and so does the distinction between even and odd particle numbers.
The energy of the infinite system diverges, being an extensive quantity. The relevant intensive quantities are the energy density and the energy per particle which are connected through where the energy density is calculated as the TL limit of Eqs. (19) or (33): It should be noted that, both the BCS and PBCS ground-states converge to the same state at TL [82]. The difference lies in that BCS approaches the infinite system through a grand canonical ensemble while PBCS does so in a canonical ensemble. This means that Eq. (38) could be calculated as the TL limit of Eqs. (27) or (32), equivalently.

The solution of the BCS gap equations
The gap equations of BCS come from minimizing the free energy of the BCS ground state with respect to the pair occupation distributions. As mentioned above, the normal state corresponds to ∆(k) = 0 and so when solving the gap equations, for a given interaction, one is inquiring about the existence of pairing correlations and gets an answer in the form of the gap distribution: the trivial solution represents the absence of pairing. The gap equations for a finite system with an even number of particles are Eqs.  (35). Irregardless of the nature of the superfluid system (even, odd, or infinite), the gap equations are a set of non-linear coupled equations and so they have to be solved self-consistently. We will describe only the method used for the even particle-numbers, since it is identical to the ones used for odd particle-numbers and the TL. With the solutions of the equations being ∆(k) and µ, one can solve the first gap equation, Eq. (14), iteratively for a given µ and use the resulting ∆(K) in the second gap equation, Eq. (18), and calculate N , thus reducing the problem to the solution of the equation: The root of the Eq. (39) can be found numerically concluding the solution of the BCS gap equations. While it is straightforward to apply identical schemes to the gap equations for odd particle-numbers as well as at the TL, when it comes to blocking an additional step is needed to properly describe the ground state of an odd particle-numbered system. As mentioned above, the blocked state b that appears in Eqs. (29) & (30) should be chosen so that it minimizes the energy of the system. For PBCS this energy corresponds to the energy in Eq. (32). This minimization entails the computationally expensive task of solving the blocked gap equations multiple times to acquire the v k and u k distributions needed for the blocked PBCS energy. An alternative route, suitable for strongly correlated systems, is a perturbative scheme. That is, one can attain an approximation of the pair probability distributions for an odd system by solving the gap equations that correspond to an even system with an odd number of particles. Consequently these distributions can be used repeatedly in Eq. (32) for a range of momentum states b to identify the state that yields the minimum energy. The error introduced in the distributions v k and u k through this perturbative approach is inversely proportional to the number of pairs N 0 /2 and so one can expect that it will not alter the structure of the excitations of the state [82,85]. With the gap equations solved only once, the use of this perturbative scheme allows for the identification of the proper blocked state in a computationally inexpensive way.

Finite Size Effects and Twisted Boundary Conditions
As a finite system approaches its TL, all intensive quantities of the system approach their TL values and the random trends that they follow are called the Finite-Size Effects (FSE). These effects can be seen in quantities like the energy per particle [82], or the pairing gap, in Fig. 3. One can typically attribute this systematic error to the quantization of k-space induced by the BC: the distributions describing an infinite system are defined on a continuous phase-space and so any approximation that uses distributions defined on a discrete phase-space comes with an error. The origin of one such error can be seen in the FSE observed in the pairing gap in Fig. 3, where the curve drawn by the values of E(k) over the grid had a minimum at a momentum which did not coincide with any point on the grid site. That resulted in compromising with the closest "available" minimum which can jump to the next or the previous k-magnitude with a change in N . This is not the only source of FSE in the pairing gap but it is responsible for an additional "jiggle" seen in the pairing gap compared to other quantities. From the perspective of the OES, these additional FSE have a similar origin: the energy of the odd particle-numbered system is the minimum with respect to the blocked state b which, in finite systems, increases discontinuously with the particle number N. In fact, in systems with odd particle numbers, we observe the location of the minimum of the quasi-particle energy to almost always coincide with the momentum state b that yields the minimum energy E PBCS odd (b; N) thus synchronizing the additional FSE seen in the pairing gap and the OES. Quantities like the energy per particle experience less pronounced FSE due to their nature as sums of distributions over the discrete k-space, which makes them less sensitive to small changes in the grid size as the particle number increases. The systematic errors that are the FSE have to be dealt with when a formulation of the TL is not available and one has to extrapolate from studies of finite-size systems. In the case of NM, such difficulties are faced by the QMC family of approaches. Quantum Monte Carlo deals with finite systems performing, for the case of NM, computationally expensive simulations whose time complexity scales unfavorably with the particle-number. That allows for tackling systems up to particle numbers of the order of ∼ 100 and so, when studying infinite matter, one needs to extrapolate to the TL from, oftentimes limited, data of finite systems. This has motivated studies of the FSE using theories which, even though not quantitatively accurate, are qualitatively reliable and can identify a successful extrapolation scheme. For superfluid NM such studies have been carried out both in the BCS [60] and the PBCS [82] theories and QMC calculations have been conducted using guidance from such studies [48] at great precision. These studies use systems of N = 66 particles to simulate a system at the TL, an assumption justified by results similar to Fig. 3 in the context of BCS, and extract the pairing gap using the model-independent OES. With the studies of the FSE rendering them under control for NM, one might wonder whether we can improve the extrapolations employed and actively reduce the observed FSE. The answer lies in the very source of the error: the BC. Indeed, techniques have been developed, initialy in the context of solid state physics, that allow for active reduction of the FSE by manipulating the twist θ introduced in Eq. (5). In band structure, calculations of properties of infinite periodic solids integrate over the first Brillouin zone, an integration which, in insulators, has been replaced by the choice of a "special k-point" [75] or a grid of k-points [76]. The integration of the results over the twist is called Twist Averaged Boundary Conditions (TABC) and it has been proven that it will produce exact results in the Hubbard model, in a grand canonical ensemble, for non-interacting systems [77], and that it will make certain lattice models converge faster to the TL [78]. In the context of nuclear matter, TABC has been used for the reduction of FSE in QMC calculations regarding the EoS of nucleonic matter [79] and in Skyrme-Hartree-Fock calculations of the energy in various nuclear pasta phases [80]. Such a study has not yet been done in superfluid NM.
Motivated by their success in other nuclear systems, we applied TBC and TABC in superfluid NM aiming for a reduction of the FSE and to further improve extrapolation schemes. In 3D systems under TBC, the twist is a vector with three components each restricted to a 2π circle. All identical particles are characterized by the same twist, and all quantities are triply periodic [86] in the twist, so that: −π < θ i ≤ π , i = 1, 2, 3 .
When applying the BCS or PBCS theory in the way described in Section 2.2, TBC would give rise to a different k-grid through Eq. (6), i.e., a grid translated by θ i /L in the i axis. This is depicted in Fig. 5 where the k-grid of a 2D system of 26 free neutrons has shifted after twisting the BC by θ = 2π(1, 2)/7. This does not only shift the momentum states available to the different distributions describing the condensate, it also decreases the population M(k) of each k-magnitude. The reason behind this becomes clear if we view the shift of the grid as a shift of the origin in k-space, the reference point of the momentum magnitudes. By shifting the origin one destroys the symmetry once present in k-space making it less likely for momentum states to fall on the same circle with radius k. With the momentum cut-off remaining the same, the number of points on the grid is the same and so less population on the different k-magnitudes translates to a larger and denser set of momentum values. This is demonstrated by Fig. 6 where the serial number of the k-magnitude is plotted versus the k-magnitude that it corresponds to, showing the increase in the number of distinct k-magnitudes in TBC. Not all twist vectors with components in the (−π, π] range are independent: the set of twists generated by symmetries of the cubic lattice (i.e., the k-grid) will give rise to identical k-spaces, i.e., identical sets of k-magnitudes and populations. This can be seen directly from Eq. (6) where, if we let R be a representation of the symmetry group of a cubic lattice, R maps k to a site of a differently twisted lattice: where θ = Rθ. With the translations being taken care of by restricting θ in [0, π], the only symmetry left to exploit is rotations of π/2 around each axis. Therefore, one twist can be taken to represent all eight rotated twists that would generate the same k-space. Such transformations conserve length further justifying why twists connected through them generate momentum spaces that yield the same populations M(k). These arguments can be used to significantly decrease the number of distinct θ-points needed for a twist-integration. With different TBC giving rise to differently shifted k-spaces, they also generate different FSE since a shifted k-space gives access to different k magnitudes for the discrete distributions describing the condensate. In principle, as noted above, the resulting k-magnitude structure will be denser and so the FSE for a finite system under TBC are expected to be smaller than those under PBC. Under this light, we performed calculations of various BCS and PBCS quantities under TBC for twist-angles on a grid defined as θ = 2π 7 (l, m, n) , l, m, n = 0, 1, . . . , 6 .
As noted above, twists connected through permutations of l, m, n in Eq. (43) yield k-spaces connected through a rotation of π/2 and so they can be considered identical. Thus the values of l, m, n (without an ordering) can uniquely define a set of TBC. In the spirit of Ref. [75], we have identified a special twist-angle that yields minimum FSE in the OES, seen in Fig. 8. This specific twist does not only significantly decrease the overall FSE observed but it also decreases the minimum of FSE seen in the low-N region. These particle-numbers can be used to simulate infinite systems. It is noteworthy that the altered k-magnitude structure created by twisting the BC can give rise to tricky edge cases. For instance, the twists that come closest to the edges of the first Brillouin zone create k-magnitudes with large discontinuities at low-k. This proves problematic for the iterative solution of the gap equations, outlined in section 2.2. The large discontinuities at low momenta in the gap equations mean that a large part of the non-trivial region of the potential stays unexplored, as seen in Fig. 7 where the kernel of the gap equation in Eq. (14) is plotted for an edge-case together with that of a normal one. Even though this does not impact the solution of the fully-paired gap equation, i.e., Eq. (14) when used for even systems or for odd systems (in the context of a perturbative scheme), it can hinder the convergence of the iterative scheme for the blocked gap equation, Eq. (29). There, blocking the k-state that yields the lowest blocked energy results in decreasing the k-magnitude's population, which is very likely to be 1 for the reasons discussed above, thus removing one more probing point from the non-zero region of the potential, making the situation worse. Effectively, this translates to a non-convergence of the iterative scheme for the blocked gap equation for a range of values of the chemical potential wherein the function N (µ) in Eq. (39) becomes ill-defined. In these cases one is forced to use the solution of the perturbative scheme described in section 2.2. Even though this introduces a small systematic error, this error can be quantified by comparing perturbative solutions to exact ones in normal cases which indicate a maximum error of ∼ 2.2% in the OES, encountered at low-N. With these edge cases making up about 2% of the total number of twist angles considered in the grid of Eq. (43), this error is irrelevant for approaches like the TABC applied below. While twisting the BC can significantly reduce FSE, the averaging of the various twist-dependent quantities over a grid of twist-angles, such as the one in Eq. (43), has been shown to be even more promising in reducing the amplitude of the FSE and flattening the otherwise semi-random trend towards the TL [81]. These are the Twist-Averaged BC where one creates twist-averages of properties of the system: This is done by studying the system under TBC and calculating twist-dependent properties, which are then integrated over the twist angle. Therefore, one can create a twist-averaged OES by calculating the OES under TBC and subsequently integrating Eq. (22) over the twist space. It has been shown [81] that a grand canonical ensemble under TABC produces exact single particle properties for non-interacting systems which is because each momentum state of the infinite system, over the range of the twist angles, occurs precisely once. This asserts that the non-interacting kinetic energy in the grand canonical ensemble under TABC, is exactly equal to its TL value. With Fermi liquid theory guaranteeing a one-to-one mapping between the low-lying excited states of an interacting system and those of a non-interacting one, we can expect a substantial reduction of FSE for interacting systems in a grand canonical ensemble; argued before in the case of non-superfluid systems [77,78,81]. Furthermore, since the projection described in section 2.2 does commute with the twist-averaging, we anticipate the nice qualities of TABC in a grand canonical ensemble to be present here as well: one can imagine a process where the integration of the twist-averaging, seen in Eq. (44), is carried out before the integration of the projection, seen in Eq. (26), while still employing the grand canonical ensemble.
The averaging of the OES over the twist-angles achieves substantially reduced FSE, as seen in Fig. 8. This was calculated by numerically integrating over the twist-angle the OES in Eq. (22) which, under TBC, has a θ dependence inherited by the energy. The TABC does not only reduce the amplitude of the oscillations observed in the untwisted case, it also brings the gap to its TL value faster. We observe that for N 50 the OES is in a 2% error margin. To achieve similar accuracy over a range of N values under PBC, one has to look at systems well beyond N = 200. This result is specifically useful to ab initio approaches where the study of arbitrarily high particle numbers is not an option. For instance, in QMC approaches to NM it is customary to use systems of N = 66 to extract the TL value of the pairing gap, as per Fig. 3. However, under TABC, a better approximation of the TL is provided by N = 80 neutrons; a system well within reach of modern QMC simulations. This delineates a program of accurately extracting the pairing gap of pure NM in a model-independent fashion via a twist-averaging of the model-independent OES, a quantity accessible by a wide range of ab initio techniques.

Summary & Conclusions
We have performed the first study of superfluid NM under TABC, a technique used widely in the field of solid state physics, in an effort to prescribe a well-informed extrapolation scheme for ab initio approaches that do not have direct access to the infinite system. We applied TABC by averaging calculations of the pairing gap under TBC, in a range of twist-angles. Our results demonstrate that a twist-averaging approach, compared to PBC, substantially reduces the FSE of the pairing gap both in amplitude and variation. With the pairing gap being a quantity distinctively sensitive to finite-size, this reduction of the FSE can significantly improve the existent extrapolation schemes in ab initio methods, like the QMC family of approaches, where similar studies have dictated systems of N = 66 neutrons optimal for extrapolation to the TL, under PBC. Given our results, these systems can be replaced by ones with N = 80 neutrons under TABC which would result in extractions of higher precision.