Impact of pairing on thermodynamical properties of stellar matter

Superfluidity in the crust is a key ingredient for the cooling properties of proto-neutron stars. Investigations on crust superfluidity carried out so far typically assumed that the cluster component was given by a single representative nucleus and did not consider the fact that at finite temperature a wide distribution of nuclei is expected to be populated at a given crust pressure condition. We want to assess the importance of this distribution on the calculation of the heat capacity in the inner crust, in the framework of an extended NSE model. We additionally show that it is very important to consider the temperature evolution of the proton fraction, imposed by the β-equilibrium condition, for a quantitatively reliable estimation of the heat capacity.


Introduction
Superfluidity in the crust is essential to understand many different phenomena in compact star physics, such as the cooling of neutron stars. It is well-known indeed that pairing correlations reduce the crust thermalization time by a large fraction [1]. The specificity of the inner crust is the simultaneous presence of clusters and homogeneous matter. The studies of crust superfluidity at finite temperature have always been done solving Hartree-Fock-Bogoliubov (HFB) equations in the so-called Wigner-Seitz approximation, meaning that the cluster component is given by a single representative quasi-particle configuration, corresponding to a single representative nucleus immersed in a neutron gas. These works do not consider the fact that at finite temperature a wide distribution of nuclei is expected to be populated at a given crust pressure and temperature condition.
Finite temperature Nuclear Statistical Equilibrium (NSE) models [2,3] take into account the whole distribution of clusters, which is obtained selfconsistently under conditions of statistical equilibrium. Our aim is to analyze how the non-homogeneity of crust matter and the associated wide distribution of nuclear species affects the superfluid properties of the crust.
In most HFB calculations for the cooling problem, moreover, it is assumed that the proton fraction does not evolve with the temperature and can be estimated by the value imposed, at each baryonic density, by the condition of neutrinoless chemical equilibrium at zero temperature. As the NSE model is much less numerically demanding than a full HFB calculation, we have released this approximation and imposed the β-equilibrium condition at each finite temperature.

The extended NSE model
The complete formalism that we use can be found in [4]. We detail here our extension of the NSE model, that is the inclusion of pairing in the bulk region of the Wigner-Seitz cell.
The energy density of a superfluid and spin saturated nuclear gas at finite temperature T, in the mean field approximation, reads [5] (q=n,p): We use the Sly4 parametrization [6] of the Skyrme energy functional for the local energy density E sky and the effective nucleon mass m * q . In Eq.
(1), f q is the particle occupation number including pairing,ρ gq = 2Δ/v π (ρ gq ) denotes the anomalous density while Δ is the temperature dependent pairing gap. The pairing interaction v πq (ρ gq ), as a function of the gas density ρ gq , is given by [7] :  1 S 0 pairing gap as a function of density for homogeneous neutron matter at zero temperature, as obtained from Brueckner-Hartree Fock calculations (full line from [8]). The figure also shows the energy gap deduced solving the BCS gap equations at finite temperature (symbols from [5]).
where ρ sat is the saturation density and the parameters V π , η, α are fixed imposing to reproduce the 1 S 0 pairing gap of pure neutron matter, as obtained in Brueckner-Hartree-Fock calculations [8]. The resulting gap is displayed in Fig.1. The binding energy of the clusters is evaluated by the analytical liquid drop model (LDM) expression reported in [9], where the different parameters are fitted from numerical calculations which make use of the same Skyrme energy functional employed for the unbound component.
The equilibrium distribution is obtained by minimizing the total free energy corresponding to an arbitrary collection of different Wigner-Seitz cells, subject to the constraint of total baryonic and charge density conservation [4].

Main results
In order to facilitate a quantitative comparison with the previous literature, we have chosen ten representative values for the baryonic density which have been proposed in [10]. These values cover the inner crust of the neutron star, approximately from the emergence of the neutron gas close to the drip point (denoted as cell 10) to a density close to the crust-core transition (cell 1).

Clusters distribution
The distribution of the cluster size obtained by our NSE calculation is displayed in Fig. 2. We can see that at the lowest densities and temperatures the distribution is strongly peaked and can be approximated by a unique nucleus, but increasing the temperature and/or moving towards the inner part of the crust, many different nuclear species can appear with comparable probability. Moreover, at sufficiently high temperatures, light resonances can appear and even become dominant in the composition of matter. Such configuration cannot be addressed in mean-field based formalisms like HFB. Looking at the isotopic distribution of clusters, especially for the lowest density cells (cells 7 and 10), one observes that the β-equilibrium condition leads to a significantly larger cluster asymmetry. This happens especially for light resonances, so at the highest temperatures, where the contribution of the most unbound clusters is increased.

Energy and heat capacity
The effects of the temperature dependence of the β-equilibrium condition and of the whole distribution of cells can also be appreciated looking at the behavior of the total energy density of the system. We have seen that these effects are very important not only at the highest densities, but also at the lowest ones when the temperature gets higher, whereas the emergence of light clusters occurs. The transition temperature from the superfluid to the normal phase is signalled by a kink in the behavior of the energy density, which will lead to a peak in its derivative with respect to the temperature, that is the associated heat capacity shown in Fig. 3. The temperature derivative was performed numerically following the trajectory of β-equilibrium: this means that only the total baryonic density is constant, but the proton fraction is not. As we can see, the temperature dependence Dashed line: as the full line, but the value of the global proton fraction is assumed equal to that one calculated from β-equilibrium at the lowest temperature, that is T=100 keV. of the β-equilibrium condition has a dramatic effect on the heat capacity. In particular the peak due to the phase transition is strongly smeared out in the outer region of the inner crust, from cell 7 to 10, while on the contrary, at the highest densities (cells 1 to 3) the consideration of the temperature variation of the proton fraction increases the size of the peak. In this case, indeed, it is possible to show that β-equilibrium path favors a discontinuous trend of all thermodynamic quantities at the transition point.
Actually, the quantitative value of the heat capacity, as shown in Fig. 3, cannot be directly compared to the results of previous HFB works [1] because of the different mean-field and/or pairing model. However, we have verified that the temperature location of the heat capacity peak, its height and width are almost identical, if we take the same parameters for the pairing interaction employed in that work.

The effect of mass functional
In all the calculations presented in the previous sections, we have used the Skyrme-based liquid-drop formula given in [9]. This choice allows a consistent treatment of the bound and unbound matter component within the same energy functional. However, as it happens also in other mean-field model, light clusters are systematically underbound, with respect to heavier ones, by this formula and this effect is even more dramatic for the most neutron rich light resonances, at the limit of nuclear binding, which can appear in our distribution. We have therefore repeated the same calculations, making use of the experimental mass, whenever this value is known [11]. complete NSE calculation making use of the analytical expressions reported in [9] for binding energies (labeled as LDM in the figure). Line with symbols: as the full line, but experimental binding energies from [11] are used whenever available (EXP + LDM in the figure).
As we can see in Fig. 4, in the range of temperatures considered, the results are similar to the ones presented in Fig. 3 both for the highest (cells 1-2) and lowest (cells 7 to 10) densities. However, in the cells from 3 to 6 the consideration of the experimental masses has dramatic consequences because it modifies the height of the heat capacity peak and also its position in temperature. Moreover, a second peak appears in the case of cells 3 and 4 whose location, it is possible to show, coincides with the temperature where a sharp transition to the dominance of light clusters occurs. The same features could also appear in cells at lower densities, but at temperature values that are beyond the range considered in the present study.