Effects of pairing gap and band gap on superfluid density in the inner crust of neutron stars

Calculations of the superfluid density in the inner crust of neutron stars by different approaches are in strong disagreement, which causes a debate on the accountability of pulsar glitches based on superfluidity. Taking a simple unified model, we study the dependence on approximation of the superfluid density in a periodic potential. In comparison with the Hartree-Fock-Bogoliubov (HF-Bogoliubov) theory which treats the effects of the band gap and the pairing gap on equal footing, we examine the HF-BCS-type approximation in which the former is incorporated in priority, and another approximation in which the latter is incorporated in priority. We find that, when the pairing gap and the band gap are comparable as in the inner crust of neutron stars, they need to be treated on equal footing, and the HF-BCS approximation can considerably underestimate the superfluid density even if the pairing gap is much smaller than the Fermi energy. Our result suggests that the validity of the HF-BCS approximation for evaluating the superfluid density in neutron star crusts is questionable.


I. INTRODUCTION
Neutron stars are very accurate natural clocks spinning at a regular period. However, it is also known that their spinning rate sometimes suddenly increases and then relaxes to that before the spin-up. This phenomenon is called "glitch", which is a long-standing problem of neutron star physics. A promising model of glitches is based on neutron superfluidity in the crust of neutron stars [1][2][3][4][5]. In the inner crust, nuclei form a crystalline lattice, and neutrons permeating the lattice forms a superfluid which possesses quantized vortices [6,7]. It is considered that, while the vortices are usually pinned to the nuclei consisting of a normal fluid, the sudden collective unpinning of the vortices releases the angular momentum of the superfluid to the normal fluid, which is observed as a glitch [8][9][10][11][12][13]. The neutron superfluid density is a fundamental quantity for understanding glitches since it is directly related to the strength of the glitches.
Although the superfluid can flow without friction, it can be still entrained by the lattice of nuclei mainly due to the band gap caused by the periodic lattice. To estimate the entrainment effect, Chamel [14] has performed a three-dimensional band calculation of neutrons in the inner crust based on the Hartree-Fock (HF) approximation. He has obtained the superfluid density that is about one tenth of the believed value. Such a low value of the superfluid density is not large enough to explain the observed strength of the glitches, and the standard glitch model based on neutron superfluidity in the crust runs into serious difficulties. Consequently, the entrainment effect on the superfluid density attracts great interests [15][16][17][18][19][20].
However, the pairing effect is neglected in the HF approximation and its importance has been pointed out in Ref. [16]. The equal-footing treatments of the pairing gap and the band gap leads to recovery of the superfluid density that is around 70 percent of the believed value. On the other hand, Chamel [17] reported an updated calculation showing that the superfluid density is still low even if the pairing effect is included by the BCS approximation. Therefore, the problem of the superfluid density in neutron star crusts is still controversial. The disagreement between these studies suggests that the superfluid density depends on the approximation method whether or not the pairing gap and the band gap are included on equal footing.
There is a recent development in the study of the entrainment effect on the slab phase of neutron star matter [18,20]. The slab phase has one-dimensional periodic layers of nuclei like lasagna pasta [21,22] and provides a simple setup for calculations. A self-consistent band calculation without pairing in Ref. [18] yields the reduction of the effective mass of neutrons in the slab phase by a factor of 0.65-0.75 from their bare mass. In addition, a time-dependent self-consistent band calculation in Ref. [20] has also reported the reduction of the collective mass of neutrons. In the latter work, the reduction of the mass is interpreted to be caused by a counter flow of the dripped neutrons along the opposite direction to the motion of the slab. These studies claim that the entrainment effect enhances the mobility of the dripped neutrons and the conventional scenario of the glitch phenomena based on neutron superfluidity is still tenable for the slab phase. However, the slab phase is a thin layer of the inner crust and the interpretation of the glitch is still unclear for the entire inner crust.
Another recent direction is on the impurity disorder of the lattice [19,23]. An effect of structural disorder on neutron superfluidity in the crust has been studied using a model of an amorphous metallic alloy [19]. The resulting superfluid density is sufficient to explain glitch phenomena in contrast to that of the regular lattice.
In the present paper, we will shed light on the dependence of the superfluid density on the approximation methods and the importance of the equal-footing treatment of the pairing and the band mixing. The previous results have been obtained by large-scale self-consistent band calculations, and thus it is not clear what causes the differences among these results. Therefore, in the present work, we take a simple minimal model of a fermionic superfluid under a periodic potential with minimal three bands including the pairing ∆ and the band mixing V . With this minimal model, we compare the following three methods: (i) We directly solve the Bogoliubov-de Gennes (BdG) equations with V and ∆ on equal footing, which is so-called the HF-Bogoliubov theory. (ii) We first incorporate only the band mixing by V , and then include the particle-hole mixing by ∆. This method can be regarded as the HF-BCS approximation [24]. (iii) We first incorporate only the particle-hole mixing by ∆, and then include the band mixing by V . Methods (ii) and (iii) are approximations of Method (i). We shall demonstrate that the resulting superfluid density shows a significant difference among the three methods when ∆ ∼ V even if ∆ is much smaller than the Fermi energy E F as in the case of superfluid neutrons in the inner crust. Especially, the superfluid density by Method (ii) can be considerably smaller than that by Method (i) in such a situation. We also discuss the implication of our result to the superfluid density in the inner crust of neutron stars. This paper is organized as follows. We present our minimal model in Sec. II, and explain the details of the three methods to calculate the superfluid density in Sec. III. In Sec. IV, we present the results showing that there is a considerable difference in the superfluid density among the three methods when ∆ ∼ V even if ∆ ≪ E F . Implications of our result to the superfluid density of the inner crust of neutron stars are discussed in Sec. V. Finally, concluding remarks are given in Sec. VI.

II. MODEL
The BdG equations in one dimension with the superfluid momentum Q are written as [16,25] where i is the shorthand index for the Bloch momentum and the band indices, u i and v i are the quasiparticle amplitudes of state i, ǫ i is the corresponding quasiparticle energy, ∆(x) is the pairing gap, and H ′ Q (x) is given by Here, m is the bare neutron mass and µ(Q) is the chemical potential. We choose the external potential V ext with the following sinusoidal form: where the magnitude V of the periodic potential is a real constant and K is the reciprocal lattice vector. We selfconsistently determine the Q-dependence of µ such that n(Q) = n(Q = 0), where n(Q) is the particle density at µ(Q). We consider the periodic solutions such that the quasiparticle amplitudes u k (x) and v k (x) for quasimomentum k are given by the Bloch states: where we restrict within three bands j = 0, ±1 for simplicity. Substituting Eqs. (4) and (5) into Eq. (6), we get Here, we have neglected the momentum dependence of ∆(k).H ′ (k, Q) in the diagonal blocks is given bỹ Equations (6), (7), and (8) constitutes our minimal model with three Bloch bands. The explicit matrix form of our model is given in Appendix A. In this simple model, both the strength V of the external periodic potential characterizing the band gap and the strength ∆ of the particle-hole mixing characterizing the pairing gap are introduced as pre-determined parameters.

III. METHODS
The superfluid density n s is defined as the second-order derivative of the energy density e with respect to Q: To obtain the superfluid density, we have to calculate the energy density. The energy density of the system described by the BdG equations is given by Here, we take a summation of i over the positive eigenvalue modes. In Method (i), we straightforwardly calculate the energy density (10) with the solutions of Eq. (6). The superfluid density by this method is denoted by n s BdG . Method (i) is so-called the HF-Bogoliubov theory which incorporates V and ∆ on equal footing, and provides the reference result to examine Methods (ii) and (iii), which are further approximations of Method (i). Below, we explain the energy density and the superfluid density by Methods (ii) and (iii).

A. Method (ii)
First, we consider only the band mixing by V and solve the eigenvalue problems forH ′ (k, Q) and −H ′ (k, −Q) in Eq. (6), respectively: where j labels the band, and the subscripts p and h in the eigenvalue represent the "particle" and "hole", respectively. This procedure corresponds to solving the HF equations for a periodic potential and ψ j (k, ±Q) corresponds to the HF basis. Next, we consider the particle-hole mixing by the pairing gap ∆, and solve the following BdG equations: We calculate the energy density (10) with the solutions of Eq. (13). The superfluid density obtained by this method is denoted by n s V ∆ . As we shall discuss in Sec. V A, Method (ii) can be identified as the HF-BCS approximation employed in [17]. This approximation incorporates the band gap in prior to the pairing gap.

B. Method (iii)
In Method (iii), we first incorporate only the particlehole mixing due to ∆ by solving the following BdG equations: with ξ k+jK,Q = ξ(k+jK, Q), the band index j = {0, ±1}, and a = + or − which labels positive and negative eigenvalues, respectively. Next, we incorporate the band mixing due to V by solving the following one-particle Schrödinger equation: with where the index l of the eigenvalue takes 0, 1, and 2 for each j band, and l = 0 labels the lowest eigenvalue.
We have determined the diagonal elements of H given by Eq. (16) by the following argument. The diagonal elements of H are the energies of the modes of k+jK (j = {0, ±1}) without the band mixing. From the expression of the energy density of the system (10) can be regarded as the energy of each mode i. In the region of ∆ ≪ µ, we can approximate that |v i | 2 ∼ 1 and . From the energy eigenvalues ǫ l ∆V obtained from Eq. (15), we calculate the energy density of the system as follows. We suppose that the N -particle wave function Ψ is given by the Slater determinant of one-particle Bloch waves ψ k with the quasimomentum k as: The energy of the lowest state is calculated as Since the Hamiltonian (16) does not have a two-body interaction, the total energy of the system is written as a sum of the single-particle energy. We approximate the summation about n in Eq. (18) by an integration of the momentum k over the Fermi sea V FS with the Fermi momentum k F as follows: . Dividing by the system size L, we get the energy density The superfluid density n s ∆V by this method is given by Opposite to Method (ii), Method (iii) incorporates the pairing gap in prior to the band gap.

IV. RESULTS
The left panel of Fig. 1 shows the superfluid density calculated by the three methods as a function of V /∆. Since the pairing gap of superfluid neutrons is typically the order of 1% of the neutron Fermi energy in the inner crust of neutron stars, here we set ∆/E F = 0.01. In addition, since we examine the three methods for the situation in which the band gap effect is nonnegligible, we set K/k F = 2 where the Fermi surface is perfectly nested [26]. Consequently, in our setup, the lowest band is almost filled and a part of the next lowest band within ∆ from the Fermi energy is slightly filled. For onedimensional systems, the particle density n is given by n = 2k F /π. The blue, orange, and green lines show the superfluid density by Methods (i), (ii), and (iii), respectively. The results show that the superfluid density by the three methods considerably differ in the region of V ∼ ∆ even though the pairing is negligibly small compared to the Fermi energy, ∆ ≪ E F . For example, at The superfluid density n s BdG by Method (i) shows an intermediate value between n s V ∆ [Method (ii)] and n s ∆V [Method (iii)]. It is close to n s ∆V in the small V /∆ region and n s V ∆ in the large V /∆ region. Furthermore, in the region of small V /∆, the behavior of the superfluid density n s V ∆ by Method (ii) is qualitatively different from the others. n s BdG and n s ∆V behave quadratically in V /∆ whereas n s V ∆ behaves linearly. The reason of this difference can be understood as follows. Let us first consider a calculation only with the band mixing, which corresponds to the first half step of Method (ii). Our setup is the nested case, K = 2k F , and the Fermi surface (i.e., two points at k = ±k F ) lies inside the band gap when V is nonzero. As a consequence, the system without pairing is an insulator and the flow density is n s /n = 0 for any nonzero V . On the other hand, for V = 0, there is no band gap and thus all the particles contribute to the flow, so that n s /n = 1 even without pairing. Therefore, n s shows step-wise behavior such that n s /n = 1 at V = 0 and n s /n = 0 otherwise, and this behavior can be regarded as the reason of the qualitative difference of the superfluid density by Method (ii) from the others. In Method (ii), while this step-wise behavior is smoothened by the pairing effect incorporated in the second half step in Method (ii), the resulting superfluid density n s V ∆ is still influenced by the step-wise behavior before including the pairing.
The right panel of Fig. 1 shows the superfluid density as a function of ∆/E F for V /E F = 0.01. The blue, orange, and green lines again show the superfluid density by Methods (i), (ii), and (iii), respectively. The superfluid density n s /n by both Methods (i) and (iii) rapidly increases and saturates to unity as ∆/E F increases. Consequently, the difference between them is negligibly small for larger ∆/E F 0.1. On the other hand, n s /n by Method (ii) slowly increases and does not saturate yet even at ∆/E F = 0.1. In the limit of ∆/E F → 0, the superfluid density by all the three methods should converge to zero while we cannot calculate for very small ∆/E F below ∼ 0.01 due to the limitation of the numerical accuracy.
In Fig. 2, we plot the ratios n s V ∆ /n s BdG (left panel) and n s ∆V /n s BdG (right panel) as functions of V /∆ for several values of ∆. The green, orange, and blue lines show these ratios for ∆/E F = 0.01, 0.05, and 0.1, respectively. For, e.g., ∆/E F = 0.01 (green lines), n s V ∆ /n s BdG takes a minimum value of 0.29 at V /∆ = 1.6 (left panel) and n s ∆V /n BdG s takes a maximum value of 8.6 at V /∆ = 4.0 (right panel). This result clearly demonstrates that the resulting superfluid density can considerably differ by the approximation when the pairing gap and the band gap are comparable. Interestingly, this difference is more prominent when the pairing gap is negligibly smaller than E F (e.g., the green lines for ∆/E F = 0.01). As ∆/E F increases, the minimum and maximum values of these ratios approach unity due to the saturation of n s /n as seen in the right panel of Fig. 1. In the region of V /∆ 1, the three lines converge to a single line and thus ∆/E Fdependence of these ratios becomes negligible.
Furthermore, the ratios n s V ∆ /n s BdG and n s ∆V /n s BdG approach to unity in the limits of V /∆ → 0 and V /∆ → ∞ in Fig. 2. In the limit of V /∆ → 0, the reduction of the superfluid density by the band gap effect vanishes. Thus, all the resulting superfluid densities by the three methods agree with the particle density, so that n s V ∆ /n s BdG = n s ∆V /n s BdG = 1. In the limit of V /∆ → ∞, on the other hand, the band mixing by V is dominant and thus it is considered that the interplay between the effects of the band gap and the pairing gap is lost. As a consequence, once the effect of V is included, the resulting superfluid density is unchanged irrespective of the order in which V and ∆ are included. In conclusion, the equal-footing treatment of V and ∆ by Method (i) is crucial for the reliable value of n s when V and ∆ are comparable. It is noted that they are indeed comparable in the inner crust of neutron stars as we will discuss in Sec. VI.

A. Correspondence to the HF-BCS approximation
We show that Method (ii) is essentially the same as the HF-BCS approximation employed in Refs. [17,27]. We first diagonalize the diagonal blocks ofH ′ (k, Q) and −H ′ (k, −Q) in Eq. (6). Here, we obtain two HF basis functions ψ(k, Q) and ψ(k, −Q) defined in Eqs. (11) and (12) for the particle Hamiltonian and the hole Hamiltonian, respectively. We also obtain the two orthogonal matricesÕ(k, Q) andÕ(k, −Q) which diagonalize each of them, and introduce the following matrix: where O is a 3 × 3 zero matrix. We rotate the BdG equations (6) by O, and obtain Here,H ′ V (k, Q) andH ′ V (k, −Q) are the diagonalized particle and hole Hamiltonians given by where ǫ i p,V and ǫ i h,V are the eigenvalues of Eqs. (11) and (12), respectively.
The rotated pairing block is given by∆ V (k, Q) = O −1 (k, Q)∆Õ(k, −Q). We note thatÕ −1 (k, Q) is not the inverse ofÕ(k, −Q) for nonzero Q, and the rotated pairing block∆ V (k, Q) has off-diagonal elements and depends on Q and k. In the HF-BCS approximation, we discard these off-diagonal elements: which could be a reasonable approximation when the system is almost uniform. Furthermore, in [17,28], the diagonal elements (∆ V (K, Q)) ii for all i are replaced by a single value ∆, so that their dependence on k and Q are neglected. This treatment is essentially the same as Eq. (13) in Method (ii).

B. Implications to three-dimensional systems
We discuss implications of our results to the superfluid density in the inner crust of neutron stars. Our model is in one dimension while the real system is in three dimensions. The reduction of the superfluid density by the band gap is the most prominent in the onedimensional, perfectly nested case of K = 2k F considered in the present work as we provide two reasons below. Therefore, our result of n s /n can be regarded as the most conservative value in which the suppression by the band gap is strongest: in the neutron star crusts, such suppression of the superfluid density should be weaker.
The first point is that the superfluid flow Q and the reciprocal lattice vector K can take various directions in three dimensions. The effect of the band gap suppressing n s should be weaker for a flow whose Q is not along K. The actual periodic lattice potential in the inner crust of neutron stars has many Fourier components with various directions and magnitudes of K. On the other hand, in one dimension, there is no orthogonal component of K to Q, and thus the reduction of the superfluid density is expected to be the most prominent.
Next, the reduction of the superfluid density is caused by the increase of the effective mass of the dripped neutrons around avoided-crossing points (ACPs) in the band structure. Curvature of the energy dispersion is the inverse of the effective mass, and the slope of the dispersion is zero at the ACP as shown at k/k F = ±1 in Fig. 3. To get large reduction of the superfluid density, almost all of the particles near the Fermi surface have to feel the effect of the ACPs. The Fermi surface in the one-dimensional momentum space is just two points at k = ±k F and lie on the ACP when K = 2k F as shown in Fig. 3. It is noted that the Fermi surface in three dimensions is a twodimensional surface while the ACPs are points. Although three-dimensional systems have many reciprocal lattice vectors, they are just a finite number of points and cannot completely cover the two-dimensional Fermi surface. On the other hand, the Fermi surface in a one-dimensional system consists of two points, and the reciprocal lattice vectors completely overlap with these two points in the nested case. Therefore, the effect of the ACPs is the most prominent in our one-dimensional setting. Figure 4 shows that the superfluid density of our model as a function of K/k F . The superfluid density rapidly recovers as the Fermi surface moves away from the ACP. This result suggests that, in three-dimensional systems, the suppression of the superfluid density by the band structure is reduced compared to the one-dimensional systems due to the large part of the Fermi surface uncovered by the ACPs.

VI. SUMMARY AND CONCLUDING REMARKS
There is a serious disagreement among previous works [14,16,17] in the superfluid density in neutron star crusts, and this problem is still controversial. To get a better understanding of this issue, we have studied the dependence on approximation of the superfluid density in a periodic potential using a minimal model including the key competing effects of the band gap V and the pairing gap ∆. We have demonstrated three calculation methods: (i) We directly solve the BdG equations with V and ∆ on equal footing. (ii) We first incorporate only the band mixing by V , and then include the particle-hole mixing by ∆. (iii) We first incorporate only the particlehole mixing by ∆, and then include the band mixing by V . Methods (ii) and (iii) are not equivalent to Method (i), but they are different approximations of the latter such that Method (ii) incorporates V in prior to ∆ and Method (iii) incorporates ∆ in prior to V . Method (i) is called the HF-Bogoliubov theory employed in [16] and Method (ii) corresponds to the HF-BCS approximation employed in [17]. We have found that the resulting superfluid density is considerably different among the three methods when V is comparable to ∆ even if ∆ is negligibly smaller than E F . For example, in the perfectly nested case of K/k F = 2 in our one-dimensional minimal model, the superfluid density by Method (ii) is about one third of that by Method (i) for ∆ = V = 0.01E F [Eqs. (22) and (23)].
Our results imply a possibility that the HF-BCS approximation considerably underestimates the superfluid density in the inner crust of neutron stars. In the inner crust, the pairing gap is estimated to be comparable to the magnitude of the lattice potential [16]. Namely, the lattice potential has many Fourier components V K with wave vectors K corresponding to the reciprocal lattice vectors, and the magnitude |V K | for the primitive reciprocal lattice vector is largest whose value is around 1-2 MeV [16,28,29]. On the other hand, 1 S 0 pairing gap in the density region of the inner crust is typically ∼ 1 MeV [30]. Therefore, we conclude that the equalfooting treatment of the pairing and the band mixing by the HF-Bogoliubov theory is essential to calculate the superfluid density of the inner crust in neutron stars. The HF-BCS approximation tends to underestimate the superfluid density, especially under the condition of the neutron star inner crust.
In the density region of the inner crust, where the su-perfluid density is drastically reduced in the calculation by Chamel [14], typical values of the Fermi energy of dripped neutrons and the distance between the neighboring nuclei are 15 MeV and 40-50 fm, respectively. Therefore, the typical magnitude of the ratio between K and k F is K/k F ∼ 10 −1 . This value is very different from the value K/k F = 2 taken in our calculations. The real system allows an infinite number of bands and many ACPs exist around the Fermi surface, so that the band gap effect is nonnegligible even at K/k F ∼ 10 −1 as in Refs. [14,16]. On the other hand, since our minimal model has only three bands and the band gap is much smaller than the Fermi energy in the current problem, the effect of the band gap is negligible in our model unless K/k F = 2. To discuss the band gap effect as well as the pairing gap effect by this minimal model, we thus have set K/k F = 2.
In our results, the superfluid density is always smaller than the particle density whereas some previous studies on the slab phase reported that the effective mass of the neutrons is reduced from their bare mass [18,20]. In our work, we consider permeating neutrons under the external periodic potential while these works have studied the system consisting of free neutrons, lattice neutrons (i.e., neutrons bound in nuclei), and protons with the Bloch boundary conditions. In the latter case, there is an ambiguity in the definition of the free neutrons [18]. The reduction of the effective mass suggests that some of the lattice neutrons turn to the conducting neutrons by the flow. Our setup does not have such an ambiguity, and the superfluid density is necessarily smaller than the particle density.
To understand the discrepancies in the superfluid density among previous calculations, we have employed a simple toy model. Once the difference among the approximations has been clarified within the minimal model in the present work, it is now interesting to study the approximation dependence using a model taking account of more realistic setup of the inner crust such as the lattice structure in three dimensions [31], defects and disorder of the lattice [32,33], realistic nuclear forces [34] and so on. We leave these issues for future works.