Comparison between the Thomas–Fermi and Hartree–Fock–Bogoliubov Methods in the Inner Crust of a Neutron Star: The Role of Pairing Correlations

Abstract

We investigated the role of a pairing correlation in the chemical composition of the inner crust of a neutron star with the extended Thomas–Fermi method, using the Strutinsky integral correction. We compare our results with the fully self-consistent Hartree–Fock–Bogoliubov approach, showing that the resulting discrepancy, apart from the very low density region, is compatible with the typical accuracy we can achieve with standard mean-field methods.

1. Introduction

The quest to find the equation of state (EoS) that best describes the properties of neutron stars (NS) [1] is one of the major challenges in nuclear physics. NSs are extremely compact objects, and so the density and pressure through their interior spans several orders of magnitude, and consequently it is important to use a theoretical model that can accurately describe such a large density range.

The tool of choice to describe both finite nuclei and NS is the nuclear energy density functional (NEDF) [2]. By carefully adjusting the parameters of the functional [3], it is possible to obtain a unified equation of state [4] that can describe the entire NS, from the low-density outer region to the core. Thanks to the latest advances both in the way one observes them [5] and the technique used [6], it is possible to provide additional constraints to the EoS [7,8]. By combining those with more traditional constraints based on heavy-ion collision experiments [9,10] it is possible to obtain interesting information about the properties of nuclear matter at high densities. By combining all this information, together with the most recent measurements of finite nuclei [11], it is possible to construct accurate models to describe the physics of such massive objects.

Due to the extreme pressure gradient, NS matter is not homogeneous. With current models [1], one can identify two main regions: the crust, which has a crystalline structure, and the core, which is in a liquid phase whose composition is still under debate [12,13,14,15,16]. The crust spans a density range of $\approx {10}^{-8}{\rho }_0$ to $\approx 0.5{\rho }_0$, where ${\rho }_0=0.16$ fm${}^{-3}$ is the nuclear saturation density.

The crust can be further divided into two regions: the outer crust and the inner crust. From the earliest models [17] to more modern ones [18,19], it has been predicted that the outer crust has a crystalline structure of finite nuclei, surrounded by a gas of free electrons. The inner crust has a similar structure, but because of the higher density, neutrons start to drip out from the nuclei and form a gas [20,21,22].

Modeling the inner crust structure is particularly challenging, since it requires the treatment of the clusters and gas on an equal footing. Ideally one should use band theory, as typically used in solid state physics, to discus the properties of crystal with delocalized electrons [23,24].

A widely accepted simplification consists of defining a spherical Wigner–Seitz (WS) cell centered on each cluster at a given baryonic density, and to assume no interaction between cells. The size of the cell is determined by assuming charge neutrality, and that protons, neutrons and electrons are at $\beta$-equilibrium within the WS cell. We refer to reference [23] for a detailed discussion on the validity of such approximation. Having defined the main features of the system, one can then determine the proton number Z within the cell, and the cell radius $R_{WS}$, after minimizing the total energy per particle e of the system. This then yields the total particle number A, and the proton fraction $Y_p=Z/A$, a quantity of astrophysical importance.

Within the NEDF framework, one should calculate the nuclear contribution to the total energy of the WS cell using the Hartree–Fock–Bogoliubov (HFB) [25] equations, since they provide a fully quantum mechanical description of the system without making an artificial distinction between bound and unbound neutrons. The downsides of the HFB method are its computational cost and the possible numerical inaccuracies related to the boundary conditions adopted in the calculation [26,27,28]. We refer to reference [29] for a more detailed discussion. More recently, the HFB method has been used to investigate heating of the crust [30] and the crust-core transition [31], though these are not systematic studies of the crust.

To overcome the numerical difficulties that arise in the HFB method, but also to reduce the computational cost, several authors have adopted semiclassical methods based on the extended Thomas–Fermi (ETF) [32] approximation. To account for nuclear shell structure, an energy correction is added using the Strutinsky integral (SI) method [33]. Combined, these are named the ETFSI method.

The results obtained using the two methods seem to qualitatively disagree. Using HFB, one typically observes clusters with a variety of different Z through the crust [29,34]. On the other hand, ETFSI points towards an inner crust comprising clusters with only $Z=40$, both with Skyrme forces [4,33,35] and with finite-range Gogny forces [36]. The main reason for this disagreement is probably related to the properties of the functional used to perform the calculations. In previous work [37], we compared the structure of some WS cells performing both ETFSI and HFB calculations using the same functional. We observed an energy discrepancy far larger than the estimated error of the HFB method [29], and we identified the cause as a lack of neutron pairing correlations [38] in the ETFSI approach [35].

In the present article, we therefore present a simple methodology to incorporate the neutron pairing energy contribution into the ETFSI method, and we compare the results obtained with those from HFB calculations.

The article is organized as follows: in Section 2 we briefly discuss the energy contributions within a WS cell, and we also illustrate the ETFSI method with the extension to include pairing correlations. We present our main results for the inner crust structure in Section 3, and we provide our conclusions in Section 4.

2. The Wigner–Seitz Cell

Adopting the Wigner–Seitz (WS) approximation, we divide the periodic lattice of nuclei of the inner crust into spherical charge-neutral non-interacting cells. For a given baryonic density ${\rho }_b$, all cells are assumed identical. Each contains a nuclear cluster at its center surrounded by a gas of superfluid neutrons, and also contains a near-homogeneous ultra-relativistic electron gas. The radius of a cell $R_{WS}$ is defined as half of the distance between neighboring clusters. Under the condition of charge neutrality, the number of protons Z in the WS cell must equal the number of electrons. For given values of ${\rho }_b$ and $R_{WS}$, the total number of particles A is fixed by the relation

$$A={\rho }_b\frac{4\pi }{3}{R}_{WS}^3.$$

(1)

Consequently, a WS cell is uniquely defined by the parameters {${\rho }_b,Z,R_{WS}$}. At zero temperature, the total energy per particle of the system is given by

$$e=e_{Sky}+e_e-Y_p{Q}_{n,\beta },$$

(2)

where $e_{Sky}$ is the contribution arising from the baryons interacting via the strong force and from the Coulomb interaction between protons, while $e_e$ is sum of the kinetic and potential energies of ultra-relativistic electrons [39] and the proton-electron interaction [34]. The last term accounts for the mass difference between neutrons and protons, $Q_{n,\beta }=0.782$ MeV. $Y_p=Z/A$ is the proton fraction of the cell.

The terms related to electrons ($e_e$) and to the Coulomb component of $e_{Sky}$ are treated on equal footing for the HFB and ETFSI methods, as described in [33]. The energy contribution per particle for a Skyrme-type nucleon-nucleon interaction [40] is expressed as functional of local densities

$$\begin{array}{c}\hfill e_{Sky}=\frac{1}{A}\int \mathcal{E}\left({\rho }_q\left(\mathbf{r}\right),{\tau }_q\left(\mathbf{r}\right),\overrightarrow{J}\left(\mathbf{r}\right)\right)d\mathbf{r},\end{array}$$

(3)

where $q=n,p$ stands for the nuclear charge. Considering only time-reversal invariant systems, the Skyrme functional depends only on a linear combination of the matter densities ${\rho }_q\left(r\right)$, kinetic densities ${\tau }_q\left(r\right)$, the spin current densities $\overrightarrow{J}\left(r\right)$ and their derivatives [41]. The procedure used to calculate the local densities differs in the HFB and semiclassical methods, and it is outlined in the Section 2.1 and Section 2.2.

2.1. Hartree–Fock–Bogoliubov

In the HFB approach, the densities (and thus the fields) are calculated directly using the quasi-particle wave-functions, which are the solutions of the HFB equations [25].

$$\begin{array}{c}\hfill {\displaystyle \sum _{n^{\prime }}}(h_{n^{\prime }nlj}^q-{\epsilon }_{F,q})U_{n^{\prime }lj}^{iq}+{\displaystyle \sum _{n^{\prime }}}{\Delta }_{n{n}^{\prime }lj}^q{V}_{n^{\prime }lj}^{iq}=E_{ilj}^q{U}_{nlj}^{iq}\end{array}$$

(4)

$$\begin{array}{c}\hfill {\displaystyle \sum _{n^{\prime }}}{\Delta }_{n{n}^{\prime }lj}^q{U}_{n^{\prime }lj}^{iq}-{\displaystyle \sum _{n^{\prime }}}(h_{n^{\prime }nlj}^q-{\epsilon }_{F,q})V_{n^{\prime }lj}^{iq}=E_{ilj}^q{V}_{nlj}^{iq}\end{array}$$

(5)

where ${\epsilon }_{F,q}$ is the Fermi energy. We used the standard notation $nlj$ for the spherical single-particle states with radial quantum number n, orbital angular momentum l and total angular momentum j. $V_{nlj}^{iq}$ and $U_{nlj}^{iq}$ are the Bogoliubov amplitudes for the i-th quasi-particle with energy $E_{ilj}^q$. The single-particle Hamiltonian h is built from the Skyrme functional, while ${\Delta }_{n{n}^{\prime }lj}^q$ are the matrix elements of the pairing gap obtained from a contact pairing interaction. In the case of vanishing pairing, these equations reduce to the Hartree–Fock (HF) one. For more details on the numerical methods used to solve these equations we refer to references [28,29,38,42]. The most relevant point for the following discussion is the choice of the boundary conditions used to solve Equation (4). In the present article, we use the Dirichlet–Neumann mixed boundary conditions: (i) even-parity wave-functions vanish at $r=R_{WS}$; (ii) the first derivative of odd-parity wave functions vanishes at $r=R_{WS}$. We call them Boundary Conditions Even (BCE), in contrast to the boundary conditions odd (BCO) where the two parity states are treated in the opposite way. We have checked that this particular choice does not affect the final results. See also reference [29] for more details.

2.2. The ETFSI Method

In this section, we briefly describe the extended Thomas–Fermi + Strutinsky integral method (ETFSI), as developed in references [33,43,44,45]. Within the semiclassical approach, the neutron and proton densities are parameterized. Assuming no proton gas, we use a generalized Fermi–Dirac distribution of the form

$$\begin{array}{c}\hfill {\rho }_q\left(r\right)=\frac{{\rho }_q^{liq}-{\delta }_{q,n}{\rho }^{gas}}{1+exp\left(\frac{r-r_q}{a_q}\right)}+{\delta }_{q,n}{\rho }^{gas}.\end{array}$$

(6)

${\rho }_q^{liq}$ are the densities of the neutrons and protons at the center of the WS cell, $r=0$, while ${\rho }^{gas}$ is the density of neutrons at the edge, $r=R_{WS}$. $r_q$ are the cluster radii of the neutrons and protons, and $a_q$ are the diffusivities of the cluster surface. These 7 adjustable parameters are determined by the minimization of the energy per particle, given in Equation (2), under the constraints of charge neutrality and $\beta$-equilibrium. See Section 2.3 for more details.

The authors of references [4,33,45] have introduced an additional damping factor in Equation (6). Although such a term may be useful for avoiding convergence problems at high ${\rho }_b$, it has a small impact on the energy per particle of the system, and for the following analysis we can safely proceed without it.

The kinetic and spin-current densities are expressed as a function of the matter density ${\rho }_q$ and its derivatives via the Wigner–Kirkwood expansion [25]. In the present work, we use the full expansion, up to 4th order in gradients, employing the explicit expressions for the 4th-order density contributions as given in the appendix of reference [44]. This differs from the approach of reference [33], explained in Section II of this reference.

We now examine the quality of our semiclassical method by comparing the densities and fields with those obtained from a fully self-consistent HFB calculation. In Figure 1, we illustrate the density profiles for a WS cell with ${\rho }_b=0.02$ fm${}^{-3}$, $Z=40$ and $R_{WS}=24.0$ fm, obtained using the SLy4 functional [46], by solving the HFB equations. In the same figure, we also illustrate the results obtained using ETFSI: the seven parameters characterizing the semiclassical matter densities (Equation (6)) were adjusted to reproduce the matter densities obtained from the HFB calculation.

We observe that the 4th order expansion works nicely, and reproduces very well the neutron kinetic and spin-current densities, as shown in Figure 1 (panel a). The proton kinetic density presents a small bump around 7 fm, as seen in panel b of Figure 1. Due to the large density gradient in the proton cluster surface compared to the neutron one, the 4th order truncation is probably not fully satisfactory. The main consequence is the poor reproduction of the proton spin current density. This is not a major issue since the spin-orbit field is not too affected by such a difference.

From the densities, one obtains the corresponding fields [41]: the central potentials $U_q\left(r\right)$, the effective masses $\frac{{\hslash }^2}{2m_q^{*}\left(r\right)}$ and the spin-orbit fields $W_q\left(r\right)$. They are shown in the lower panels of Figure 1, together with the corresponding ones obtained solving the HF equations. In both neutron and proton cases the agreement is very good, thereby showing the validity of the semiclassical approximation in capturing the main features of the HFB calculation. The bump in the kinetic and spin current proton densities translates to an oscillation of the central proton potential at $r\approx 7$ fm, but the impact is small.

Under the ETF approximation, shell effects are not accounted for. Consequently, the authors of reference [45] suggested including (at least for protons) a perturbative contribution to the total energy using the Strutinsky integral (SI) theorem, without acting on the densities or fields. For the sake of completeness we implemented exactly the same method.

2.3. Choice of Functionals

Before more detailed analysis on the results obtained with HFB and ETFSI, we briefly discuss the choice of the Skyrme functional used in this work. In Figure 2 we show the energy per particle e curves for pure neutron matter (PNM) as a function of the density of the system for the three functionals we investigated: SLy4 [46], BSk21 [33] and BSk24 [47]. Their parameters were all adjusted for the functionals to be applicable to both finite nuclei and higher density neutron star matter, for modeling an entire neutron star. On the same we also show the EoS calculated using ab initio methods and given in references [48,49] (APR), and the one calculated in reference [50] (LS2). Both BSk21 and BSk24 have been fit on the LS2 EoS, while SLy4 has been adjusted using APR.

The SLy4 functional was fit to doubly-magic nuclei; as a consequence, we are left with some freedom in choosing how to model pairing correlations. In this work, we use a simple density-dependent pairing interaction of the form [51]

$$\begin{array}{c}\hfill v_q^{pair}({\mathbf{r}}_1,{\mathbf{r}}_2)=v_{0q}\left[1-\eta {\left(\frac{{\rho }_q\left(r\right)}{{\rho }_0}\right)}^{\alpha }\right]\delta \left({\mathbf{r}}_1-{\mathbf{r}}_2\right).\end{array}$$

(7)

We choose the parameters $\eta =0.7$ and $\alpha =0.45$. We assume that the pairing strength is the same for neutrons and protons and we fix the pairing strength $v_{0q}$ to obtain a maximum pairing gap in PNM of $\approx 3$ MeV, hereafter named strong, or a maximum of $\approx 1$ MeV, hereafter named weak, as done in reference [34]. These choices largely cover the available range of results concerning the density evolution of the pairing gap in infinite nuclear matter [52]. To avoid the ultraviolet divergence of the interaction given in Equation (7) [53], we adopt a smooth cut-off in quasi-particle space at $E_{ijl}^q\ge 20$ MeV that is defined by an Gaussian factor $exp\left(\left(E_{ijl}^q-20\right)/100\right)$. In Figure 3, we report the density dependence of the pairing gap in PNM obtained by solving the BCS equations. The maxima of these pairing gaps are located at ${\rho }_n\approx 0.03$ fm${}^{-3}$ (SLy4), corresponding to a Fermi momentum of $k_F^n\approx 0.94$ fm${}^{-1}$.

The pairing interactions for the BSk21-24 functionals [54] have been adjusted to reproduce the ${}^1{S}_0$ pairing gaps in both symmetric nuclear matter and PNM, as obtained from Brueckner calculations using the Argonne $V18$ nucleon–nucleon potential [55]. The resulting pairing gap is also reported in Figure 3 as a function of the neutron density. We observe that in this case, the pairing gap reaches a maximum of ${\Delta }_n\approx 2.5$ MeV around ${\rho }_n\approx 0.02$ fm${}^{-3}$, corresponding to a Fermi momentum of $k_F^n\approx 0.87$ fm${}^{-1}$.

Within the literature, there is a wide consensus on the importance of pairing correlations within the inner crust of neutron stars [24,26,28,34,56,57,58,59,60,61]. While pairing correlations are naturally included within the HFB equations, Equation (4), the original ETFSI method [45] is not capable of treating such correlations. As a consequence, in reference [35] the authors have modified the ETFSI formalism to include the effects of proton pairing with the BCS approximation [25], but still with no explicit treatment of neutron pairing correlations.

Since the neutrons in the inner crust form a gas, it is not possible to apply directly the same methodology without running into the same type of problems encountered with the HFB method, which are related to spurious shell effects in the neutron gas [27]. In this work, we have developed an additional energy correction based on the local density approximation (LDA). A similar approach was already proposed in reference [62]. In the the weak-coupling limit, the correction to the energy per particle from superfluid neutrons is [63]

$$\begin{array}{c}\hfill e_{cond}=-\frac{3{\Delta }_n^2}{8{\mu }_n}.\end{array}$$

(8)

The chemical potential ${\mu }_n$ is approximated by the corresponding Fermi energy ${\epsilon }_{F,n}=k_{F,n}^2{\hslash }^2/2m_n^{*}$, where $k_{F,n}={\left(3{\pi }^2{\rho }_n\right)}^{1/3}$. ${\Delta }_n={\Delta }_n\left({\rho }_n\left(r\right)\right)$ is the local pairing gap as extracted from PNM calculations at a given density ${\rho }_n$ as discussed in reference [64] and illustrated in Figure 3. We refer to Appendix A for more details. The energy correction in Equation (8) can be easily implemented in the ETFSI formalism without a major increase in the computational cost.

In Figure 4 we illustrate the different contributions to the energy per particle e for various WS cell at fixed baryonic density ${\rho }_b=0.02$ fm${}^{-3}$, but for the cases of weak and strong pairing. The ETFSI+pairing results have been obtained using a complete minimization of the total energy of the WS cell using the SLy4 functional. In the upper-left panel we show the nuclear contribution, referring to the ETF energy including proton-proton Coulomb interaction. The higher e with stronger pairing simply reflects the larger neutron number obtained. The higher the proportion of neutrons in the system, the more neutron pairing can occur, lowering the total energy. When we include also the electron-electron and electron-proton interactions, as shown in the lower-left panel, we see almost no difference between the weak and strong pairing cases. For strong pairing, the increase in the nuclear energy is offset by a decrease in the energy contribution from the electrons, a result of the larger WS cell.

In the upper-right panel of Figure 4, we see that the contribution to e from neutron pairing is almost flat with respect to Z, for weak and strong. Note that the weak case has been multiplied by 10 in this panel. In Figure 3, where the PNM gap is about three times higher for the strong interaction at ${\rho }_b=0.02$ fm${}^{-3}$, and Equation (8), showing the quadratic dependence of $e_{cond}$ on the PNM gap, the factor of nine increase in the neutron condensation energy was expected.

The energy per particle coming from the SI correction and proton pairing energy, as calculated in reference [35], are shown in the lower-right panel of Figure 4. The shell effects give rise to local minima at $Z=20,28,40,50,58$ for both the weak and strong interactions. The effect of increasing the proton pairing strength is to partially smooth out these shell effects.

Energy Minimization with ETFSI

A key ingredient of the ETFSI calculation is the determination of the seven parameters of the density profiles, using a minimization procedure. For this case we have used the Python library SciPy [65]. For a given ${\rho }_b$ and Z, a initial guess is made for the number of neutrons, and the corresponding $R_{WS}$ is calculated using Equation (1). The parameters of the density profiles in Equation (6) are varied, subject to constraints outlined in reference [45], to minimize the total energy of the WS cell (Equation (2)), which is calculated using a code written in Fortran 90. The SI correction and proton pairing energy is then added perturbatively, which is necessary to prevent anomalously large values for the SI correction [33]. This energy minimization is systematically repeated with different neutron numbers to find the cell configuration with the minimum energy.

After repeating this process for every even value in $20\le Z\le 60$, one finds the Z that yields the lowest energy per particle e; this is the optimum Z for a given ${\rho }_b$.

Attempts to use the full 4th-order expressions for densities, along with the profiles using an extra damping factor as used in [4,33,45], made the minimization much more difficult. Even without the damping factor, at ${\rho }_b\approx 0.05$ fm${}^{-3}$ and above, our minimization procedure begins to fail increasingly often. This is the same difficulty reported by the authors of reference [33]. Furthermore, the presence of non-spherical clusters is expected above ${\rho }_b\approx 0.05$ fm${}^{-3}$ [66]. Since the goal of the present article is not to provide a complete EoS for calculations of an entire NS, we have limited our investigation to the range of baryonic densities $0.00025$ fm${}^{-3}\le {\rho }_b\le 0.048$ fm${}^{-3}$.

2.4. HFB vs. ETFSI+Pairing

In this section we compare the results obtained in WS cells from solving the HFB equations with those from the ETFSI method, with and without the additional correction for neutron pairing correlations.

In Figure 5, we show the results obtained with ETFSI, ETFSI+pairing and HFB, for the WS cells at a few selected ${\rho }_b$, using the functionals SLy4 and BSk24. For the SLy4 functional, the strong pairing was used. The results obtained with BSk21 were almost identical to BSk24, so we do not show them in the figure. For each WS cell, the energy minimization is performed using the full ETFSI method, with or without pairing. The parameters of the resulting WS cells (${\rho }_b,Z,R_{WS}$) were used to perform HFB calculations using the same functionals. We have previously performed fully self-consistent HFB calculations [29] using the SLy4 functional and the strong pairing interaction given in Equation (7), not using the ETFSI cell parameters. We found that, while this leads to slightly different minima, they all still lie within the error bars estimated in that work.

The energy dependence as a function of Z for a full ETF calculation would be a smooth parabola, but due to the Strutinsky integral correction, we clearly observe a modification to the total energy resulting from shell structure. By comparing the ETFSI calculations with and without pairing, we observe that the positions of the minima do not change, but we observe a general reduction in the relative energy difference between the shell or sub-shell closure values and between the other WS cells. As discussed in reference [35], this is mainly the effect of proton pairing. The neutron pairing acts to globally shift the total energy as shown in Figure 4.

The HFB results are remarkably close to the ETFSI+pairing ones: near the drip density, ${\rho }_b=0.00041$ fm${}^{-3}$, there is a larger discrepancy between HFB and ETFSI+pairing. This energy difference is of the order of ≈100 keV per particle. At these low densities, the energy from the cluster is greater than that from the neutron gas, and so the different density profiles of the ETFSI and HFB methods, seen in Figure 1, are more important. The correction for neutron pairing is also less accurate at these very low densities.

For both SLy4 and BSk24, we find the two local minima at $Z=40$ and $Z=50$ for ETFSI+pairing, as also found by other works [4,33,35]. However, the minima found with HFB take several values between $Z=36$ and $Z=50$. The small discrepancy could be related to the role of neutron shell effects of the cluster that are not taken into account within the ETFSI method. From ${\rho }_b=0.013$ fm${}^{-3}$, when the neutron gas contribution starts to become more important, there is a remarkable agreement between the full HFB calculation and the ETFSI+pairing. This result confirms our previous hypothesis concerning the discrepancy observed between the two methods in reference [37]. At this density, the total energy difference between the two methods is less than 10 keV per particle. The energy minimum obtained with HFB is at $Z=50$, while the one with ETFSI+pairing is at $Z=40$. As discussed in reference [29], the accuracy of our HFB code is ≈4–5 keV per particle. Inspecting the figure, we notice that the relative energy difference between the HFB configuration at $Z=40$ and $Z=50$ is ≈2 keV per particle, clearly within the error of our calculations.

At higher baryonic density ${\rho }_b=0.005$ fm${}^{-3}$ the agreement between the HFB and ETFSI+pairing remains, although the energy minima obtained by the two calculations do not agree. However they are still compatible due to the error we estimate for the HFB calculations. From this figure, we conclude that the inclusion of neutron pairing correlations in the ETFSI method leads to a very good agreement in the total energy per particle with the more involved HFB calculation.

3. Results

Having demonstrated the compatibility between the ETFSI+pairing and HFB methods, it is now possible to use ETFSI+pairing for systematic calculations of the inner crust. We first tested our method, without neutron pairing, against previous work. We used BSk21 and SLy4, as in reference [33], and BSk24, as in reference [4]. For all of these functionals, we found that the energy minimum occurs at $Z=40$ for all densities in the range $0.00025$ fm${}^{-3}\le {\rho }_b\le 0.048$ fm${}^{-3}$ (as in the previous work), except for at very low densities, less than $\approx 0.01$ fm${}^{-3}$, where we find $Z=50$. As noted in [33], the two local minima corresponding to $Z=40$ and 50 are very close in energy at low densities. This result is independent on the use or not of the energy correction given by neutron pairing correlations.

We attribute the discrepancy in our results to our treatment of the 4th-order ETF contributions to the energy (explained in Section 2.2). The smaller clusters found at these lower densities have less diffuse surfaces. Therefore, using the full 4th-order expressions with its higher-order derivatives is more vulnerable to numerical inaccuracies. A possible source of error is the way the total energy is calculated. These small numerical discrepancies are related to the method of numerical integration and we have seen that they are sufficient to explain the different minima observed at very low density, as discussed in reference [33].

In Figure 6, we present the equation of state obtained with a full ETFSI minimization, using the SLy4. ETFSI results are shown with blue dotted lines, and ETFSI+pairing results using the strong pairing interaction with red dashed lines. As expected, the inclusion of neutron pairing effects in the system decreases the energy per particle of the system, as shown in Figure 6, but it does not affect the global trend.

In Figure 7, we compare other properties of the WS cells obtained using SLy4 with and without pairing. In panel a of Figure 7, we observe a larger total nucleon number A is the case of ETFSI+pairing. This reflects in a small reduction of the proton fraction $Y_p=Z/A$ as shown in panel b, and a small increase of the WS cell radius $R_{WS}$, shown in panel c. This adds to the reliability of the ETFSI method, whose semiclassical Wigner–Kirkwood expansion is exact in the limit of PNM ($Y_p\to 0$).

The pressure, defined as

$$\begin{array}{c}\hfill P=-{\left(\frac{\partial E}{\partial V}\right)}_{T,N,Z},\end{array}$$

(9)

is calculated following closely the approach described in Appendix B of [33]. In the presence of neutron pairing, the pressure of the cell decreases, but as panel d shows this difference is negligible, certainly far smaller than the difference that would arise from using different functionals.

All of these observed differences are at their greatest around ${\rho }_b=0.015$ fm${}^{-3}$. This is near the densities where the PNM pairing gaps, shown in Figure 3, are at their maximum. The small bumps for the ETFSI+pairing case in Figure 7, most notably in panel a, can be attributed to the inclusion of the SI correction and pairing in the second minimization step, when the total particle number A is determined (see Section 2.3).

We find that, for all functionals and densities tested, the inclusion of neutron pairing does not change the optimum Z. This is easily understood by looking at the upper-right panel of Figure 4, where the contribution to e from neutron pairing is almost constant with respect to Z.

4. Conclusions

We have presented a systematic comparison between solving the HFB equations and using the ETFSI method, using exactly the same numerical conditions. We have observed that the inclusion of neutron pairing correlation in ETFSI using a simple LDA approximation leads to a remarkable reduction in the discrepancy between the two methods, thereby confirming the success of using ETFSI+pairing for the determination of the properties of the WS in the inner crust.

After neutron pairing is included in the modeling of the inner crust, we find no change in the prediction of the optimum proton number Z at a given baryonic density ${\rho }_b$, for any of the functionals SLy4, BSk21 and BSk24. The location of the minimum of the nuclear contribution, dictated by the choice of functional, has a much bigger influence on the optimum Z.

The additional neutron pairing energy contribution leads to a general energy reduction in the Wigner–Seitz cell across the range $20\le Z\le 60$ investigated. The most interesting effect of neutron pairing is the increase in the radius of the WS cells and the number of neutrons per cell, which means a small decrease in the proton fraction is observed. The pressure is also slightly decreased, but the change is not significant.

We conclude that ETFSI+pairing method is capable of giving a very accurate description of the structure of the inner crust. The results obtained are of the same quality as more advanced HFB calculations done under the same numerical conditions. We recall that given the proximity in energy of the minima, small numerical inaccuracies may lead to different minima, thereby explaining the apparent discrepancy of results within the scientific literature.

The reliability of the ETFSI+pairing method calls for a hybrid approach to inner crust calculations: HFB would be still be used at lower densities, where it has a superior quality, and at higher densities up to the crust-core transition, ETFSI+pairing would be used.