The heat capacity of the neutron star inner crust within an extended NSE model
Abstract
 Background

Superfluidity in the crust is a key ingredient for the cooling properties of protoneutron stars. Present theoretical calculations employ the quasiparticle meanfield HartreeFockBogoliubov theory with temperature dependent occupation numbers for the quasiparticle states.
 Purpose

Finite temperature stellar matter is characterized by a whole distribution of different nuclear species. We want to assess the importance of this distribution on the calculation of heat capacity in the inner crust.
 Method

Following a recent work, the WignerSeitz cell is mapped into a model with cluster degrees of freedom. The finite temperature distribution is then given by a statistical collection of WignerSeitz cells. We additionally introduce pairing correlations in the local density BCS approximation both in the homogeneous unbound neutron component, and in the interface region between clusters and neutrons.
 Results

The heat capacity is calculated in the different baryonic density conditions corresponding to the inner crust, and in a temperature range varying from 100 KeV to 2 MeV. We show that accounting for the cluster distribution has a small effect at intermediate densities, but it considerably affects the heat capacity both close to the outer crust and close to the core. We additionally show that it is very important to consider the temperature evolution of the proton fraction for a quantitatively reliable estimation of the heat capacity.
 Conclusions

We present the first modelization of stellar matter containing at the same time a statistical distribution of clusters at finite temperature, and pairing correlations in the unbound neutron component. The effect of the nuclear distribution on the superfluid properties can be easily added in future calculations of the neutron star cooling curves. A strong influence of resonance population on the heat capacity at high temperature is observed, which deserves to be further studied within more microscopic calculations.
pacs:
26.60.kp,26.60.Gj,64.10.+h,74.20.FgI Introduction
Superfluidity in the crust is a key ingredient in the understanding of many different phenomena in compact star physics, from the cooling of young neutron stars Lat94 (); Gne01 (), to the afterburst relaxation in Xray transients Pag12 (), as well as in the understanding of glitches Pie14 (). Moreover, it is wellknown that pairing correlations reduce the crust thermalization time by a large fraction Gne01 (); For10 (). The specificity of the inner crust is the simultaneous presence of clusters and homogeneous matter, which are both influenced by pairing interactions. Indeed the occurrence of dishomogeneities has a nonnegligible influence on the pairing properties of the inner crust For10 (); Bar98 (); San04 (); Cha10 (); Pas15 (), and consequently on the time evolution of the surface temperature of the neutron star.
Present studies of crust superfluidity at finite temperature are typically done solving HartreeFockBogoliubov (HFB) equations in the socalled WignerSeitz approximation Neg73 (), meaning that the assumption is done that the cluster component is given by a single representative quasiparticle 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 conditions. Moreover, at the extremely low proton fractions associated to the inner crust, deformed nuclear structures and beyond dripline light nuclear resonances can participate to the statistical equilibrium, and might be too exotic to be well described through standard meanfield calculations. Nonspherical pasta structures, which are not accessible to the spherical meanfield, have been reported to be only marginally populated in equilibrium Ava08 (). However, at sufficiently high temperatures, light particles can appear and even become dominant in the composition of matter Roe13 (); Ava12 () and can modify the local distribution of neutron density, and the associated pairing field.
A way to include these beyondmean field effects is given by finite temperature Nuclear Statistical Equilibrium (NSE) models. In the most recent NSE implementations Hec09 (); Hem10 (); Rad10 (); Rad14 (); Fur13 (); Buy13 () the distribution of clusters is taken into account and obtained selfconsistently under conditions of statistical equilibrium. In some of these models both the gascluster interaction and the selfinteraction of the gas are included, though within semiclassical approximations Gul15 (). In particular in Ref.Gul15 () it is shown that a proper definition of the cluster selfenergies in the NSE cluster distribution allows recovering the zero temperature limit of a single WignerSeitz approximation in the (Extended) ThomasFermi limit.
This kind of approaches are however not adequate to describe the heat capacity of the crust because they do not consider the presence of pairing correlations. The aim of this paper is to analyze how the nonhomogeneity of crust matter and the associated wide distribution of nuclear species, affects the superfluid properties of the crust. Specifically, we introduce pairing correlations both in the cluster and homogeneous matter component of the NSE model in the local BCS approximation and study the effect of the cluster distribution on the heat capacity of the inner crust. We will show that the single nucleus approximation is perfectly adequate in some regions of the inner crust, but nonnegligible effects of the cluster distribution are seen close to the drip point, and close to the crustcore transition.
In most HFB calculations for the cooling problem For10 (); Bar98 (); San04 (); Cha10 (); Pas15 (), the approximation is made 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 of reference calculations Neg73 (). Even with the inclusion of pairing, the NSE model is still much less numerically demanding than a full HFB calculation at finite temperature. For this reason, we have released this approximation and imposed equilibrium at each finite temperature. This condition is justified by the fact that the time scale of cooling is sufficiently slow to insure the chemical equilibrium of weak processes at all times Gne01 (). The temperature dependence of the proton fraction is shown to have considerable effects on the heat capacity.
The paper is organized as follows. The improved NSE model with inclusion of pairing correlations in the neutron gas, is presented in section II. In this section are detailed the superfluid neutron gas (section II.1) and cluster distribution (section II.2) modelling, the calculation of the total energy (section II.3) and the calculation of the inmedium modification of the cluster surface and pairing properties due to the presence of the gas (section II.4). Section III is devoted to the presentation of the results. The composition of the inner crust in terms of cluster distribution and unbound neutrons, the temperature evolution of the energy and the heat capacity are given in the dedicated subsections (III A and III B). Section III.3 discusses the importance of a highly predictive model for the binding energies of the different nuclear species. Finally section IV gives a summary and conclusions.
Ii The improved Nuclear Statistical Equilibrum model
The complete formalism that we use can be found in Refs.Gul15 (); Rad14 (). Here, we recall the main equations and detail the inclusion of pairing both in the bulk and in the surface region inside the WignerSeitz cells, which was not considered in Ref.Gul15 ().
The model is based on a statistical distribution of compressible nuclear clusters immersed in a homogeneous background of selfinteracting nucleons and electrons. We label each nuclear species composed by neutrons and protons by their mass number and bulk asymmetry . Even below drip, the asymmetry in the bulk for a nucleus in the vacuum, differs from the global asymmetry of the nucleus, , because of the presence of a neutron skin and Coulomb effects. The relation between and is given by Mye80 (); Cen98 (); Cen09 ():
(1) 
where is the symmetry energy at saturation, is the surface stiffness coefficient extracted from a semiinfinite nuclear matter calculation and is the Coulomb parameter taken equal to MeV. The continuum states leading to the existence of a free nucleon gas can in first approximation be modelled as leading to a constant density contribution. As a consequence, the bulk asymmetry inside the clusters can be decomposed into the asymmetry of the gas weighted by the gas fraction inside the cluster, plus the asymmetry of the cluster in the vacuum weighted by the complementary mass fraction , namely
(2) 
In the previous equation denotes the bulk density. Variational arguments lead to the conclusion that, independent of the presence of an external gas, the equilibrium bulk density corresponds to the saturation density at the corresponding bulk asymmetryPan13 (). This means that the following expression, at the second order in asymmetry, can be used for the bulk density:
(3) 
In this equation, is the saturation density of symmetric nuclear matter, and are the slope and curvature of the symmetry energy at saturation. Then, solving the coupled Eqs.(2) and (3), it is possible to extract the bulk density and asymmetry.
ii.1 The free energy of the superfluid gas
The energy density of a nuclear gas, of density and asymmetry , at finite temperature T, in the mean field approximation reads Chamel (); Bur14 () (q=n,p):
(4) 
with
(5) 
Hereafter the acronym HM stands for Homogeneous Matter. We will use the Sly4 parametrization Cha98 () of the Skyrme energy functional for the local energy density and the effective nucleon mass , for the numerical applications of this paper. In Eq.(4), is the spin degeneracy in spinsaturated matter and is the particle occupation number:
(6) 
with and . and denote, respectively, chemical potential and reduced chemical potential, and
(7) 
is the single particle energy. is the temperature dependent pairing gap and denotes the anomalous density. In Eq.7, the derivative with respect to is taken at constant .
The pairing interaction is given by Chamel (); Bur14 () :
(8) 
where the parameters , , are fixed imposing to reproduce the pairing gap of pure neutron matter as obtained in BruecknerHartreeFock calculations Cao06 (). The resulting gap is displayed in Fig.1. Then the density dependence of the pairing strength can be calculated exactly in the BCS approximation by inverting the gap equation
(9) 
The reduced chemical potential is moreover obtained by fixing the particle number density:
(10) 
It should be noticed that, owing to the zero range of the pairing interaction, a cutoff has to be introduced in the gap equations to avoid divergences. Following Refs.Chamel (); Bur14 (), we adopt the energy cutoff MeV.
The free energy density is obtained adding the entropy term in the meanfield approximation:
(11) 
where the entropy density is given by :
(12)  
with .
ii.2 The cluster distribution
A given thermodynamic condition in terms of temperature, baryonic density and proton fraction is characterized by a mixture of configurations defined by with a free energy given by Gul15 ():
(13)  
In this expression, is the electron freeenergy density, is the total proton density and denotes the Wigner Seitz volume. is the free energy of the cluster immersed in the nucleon gas:
(14)  
where the total volume has been introduced and we have defined the bound fraction of the cluster by , , with .
The temperature dependent degeneracy factor includes the sum over the cluster excited states as:
(15) 
where is the density of states of the cluster, is the average particle separation energy, is the nucleon mass. See Ref.Gul15 () for details.
We can observe that the cluster energy is modified with respect to the corresponding vacuum energy both because of nuclear and Coulomb inmedium effects.
The modification of the nuclear free energy consists of a bulk term
(16) 
due to the presence of the gas in the same spatial volume, , occupied by the cluster, and a surface term which accounts for the isospindependent modification of the surface tension due to the presence of the gas at the surface of the cluster. The calculation of this last term will be detailed in section II.4. Moreover, to insure additivity of the cluster and the gas component, only the bound part of the cluster appears in the translational entropy term of Eq.(14), which also can be considered as an inmedium effect.
The screening effect of the electron density which neutralizes the WignerSeitz cell leads to a modification of the cluster free energy according to:
(17) 
with the Coulomb screening function in the WignerSeitz approximation:
(18) 
The volume of the WignerSeitz cell associated to each nuclear species is univocally defined by the charge conservation constraint:
(19) 
leading to
(20) 
The equilibrium distribution is obtained by minimizing the total free energy corresponding to an arbitrary collection of different cells , subject to the constraint of total baryonic and charge density conservation Gul15 ():
(21)  
(22) 
The result is a NSElike expression for the cluster multiplicities Gul15 ():
(23) 
where the chemical potentials can be expressed as a function of the gas densities only:
(24)  
(25) 
Concerning the cluster binding energies , theoretical coherence with the treatment of the gas demands that they are evaluated with the same Skyrme energy functional employed for the gas component. We will use for the binding energy the analytical expressions proposed in Ref.Dan09 ():
(26) 
with the asymmetry energy coefficient:
(27) 
where the different parameters are fitted from numerical Skyrme calculations in slab geometry Dan09 (). The pairing contribution to the cluster energy is evaluated according to the phenomenological expression: , where the +() sign refers to eveneven (oddodd) nuclei.
It is important to observe that this formula,Eq.(26), similar to any other meanfield model, systematically underbinds light particles, which will then tend to be underestimated in the calculations. We will discuss the effect of this limitation in section III.3. Concerning the density of states , meanfield models are known to be far off in the reproduction of these observables, and empirical adjustements have to be done. For this reason we use a backshifted Fermi gas model with parameters fitted from experimental data Buc05 ().
ii.3 Computation of the total energy
For the calculation of the heat capacity, the simplest approximation consists in considering the contribution of the cluster and the gas as simply additive, which corresponds to neglecting the term in Eq.(14).
In some early studies Lat94 (); Piz02 (), the cluster contribution was completely ignored, and the nonuniform distribution was replaced with a uniform gas formed by the total number of neutrons in the cell or by taking only the number of the unbound neutrons. It was shown in Ref.For10 () that these approximations very poorely reproduce the total heat capacity of a complete HFB calculations: the shape of the peak is too sharp, and the transition temperature is underestimated.
The displacement of the transition temperature can be simply understood as an effect of the gas density, which for a given particle number and cell size is obviously modified in the presence of the cluster, because this latter occupies a finite volume. This effect cannot be simply accounted if the same boundary conditions of the complete HFB calculations are applied to the uniform neutron gas configuration For10 ().
To compute the contribution of the gas to the energy density in the simplified hypothesis , we have to consider the total volume, , accessible to the gas, i.e. the volume left after excluding the volume of the clusters, and evaluate the corresponding gas volume fraction . Then the gas energy density can be simply written as: . In the Single Nucleus Approximation (SNA) employed in HFB calculations, we can consider a single representative WignerSeitz cell and write .
In our model, the full distribution of clusters is accounted for and the volume fraction accessible to the gas results:
(28) 
where the total volume can be written as , being the average size of the WignerSeitz volume and the total cluster multiplicity. Indicating with the normalized cluster multiplicity, the contribution of the gas to the energy density is therefore:
(29) 
and the total baryonic energy density of star matter is
(30) 
The energy entering Eq.(30) is given by the vacuum energy, Eq.(26), shifted by the electron screening effect, , and augmented of the average translational energy, , and average excitation energy corresponding to the considered temperature and cluster density of states Gul15 ():
(31)  
Here, (Eq.(17)) because the Coulomb shift is determined by the electrons. These latter being independent fermions with respect to the nucleons, they have no effect on the baryonic state counting.
As shown by Eq.(14), this simple excluded volume effect can be also formulated as the additivity of the gas with the bound part of the clusters.
This decomposition of the WignerSeitz cell between cluster and gas, accounting for the excluded volume effect, was tried in the recent HFB analysis of Ref.Pas15 (). It was shown that the transition temperature of the superfluid gas is correctly recovered, but the peak is still sharper than in the full HFB calculation. Moreover, in the full calculation a second peak at higher temperature can appear in the outer part of the inner crust, depending on the energy functional For10 (); Pas15 (). This peak corresponds to the temperature at which the whole system becomes nonsuperfluid, due to the pairing effect in the surface of the cluster. It is clear that this peak cannot be reproduced in the hypothesis of energy additivity.
These observations show the importance of accounting for the inmedium pairing corrections of the interface between the cluster and the gas, that we examine in the next section.
ii.4 In medium effects
In Eq.(30) we have assumed that the bound part of the cluster and the gas contribution are additive. This hypothesis is based on the approximation that : (i) the most important part of the inmedium correction is given by the Coulomb screening by the electron gas, and by the Pauliblocking effect of high energy cluster single particle states due to the gas Typ10 (); (ii) this latter effect can be approximately accounted for by subtracting from the local energy density the contribution of the unbound gas states. The Coulomb screening is indeed considered in the functional , and the bulk part of the in medium correction is accounted for by considering the cluster excluded volume, Eq.(28).
The approximation of Eq.(30) neglects the modification of the cluster surface tension due to the presence of an external neutron gas.
This residual inmedium modification of the cluster energy can be computed by subtracting to the total energy in each WignerSeitz cell the contribution of the gas alone and of the nucleus alone, following Aym14 ():
(32)  
Considering that this correction is expected to be a surface effect, it appears reasonable to compute it in the local density approximation (LDA), .
Since the proton contribution to the nucleon gas is very small for betaequilibrated matter and in the temperature regime concerned by our study, we can safely neglect any Coulomb effect to , meaning that we can consider solely the nuclear part of the energy in Eq.(32). Then the cluster energy in the LDA can be decomposed in an isospin dependent bulk part and residual terms varying with slower than linear (surface, curvature and higher order):
(33) 
A similar decomposition can be applied to the total LDA energy :
(34)  
where is the radius of the WignerSeitz cell, is the hardsphere radius, associated with the cluster volume , and represents a surface term since the bulk parts have been highlighted.
Using Eqs. (32), (33) and (34), we can express the residual inmedium modification simply as
(35) 
Since is related to the two surface terms deduced from Eqs. (33) and (34), we can expect the following relation to hold: , where the temperature dependent parameter should have a weak dependence on , revealing the small effect of the curvature terms.
These inmedium corrections were evaluated in Ref.Aym14 () adding to the LDA also higher orders in in the semiclassical ThomasFermi developement of the energy functional, but neglecting the pairing interaction and the temperature dependence. It was shown that is indeed a surface term , but it displays a very complex behavior with the cluster bulk asymmetry , the gas density and the gas asymmetry .
In this work we include the temperature effect and the pairing interaction according to Eq.(4), but we limit ourselves to the simple LDA. Gradient and spinorbit terms are therefore neglected in the surface correction.
In order to evaluate the inmedium surface correction through Eq.(34), a model for the density profiles , has to be assumed. We use the simple WoodSaxon analytical profiles proposed in Ref.Pan13 () and successfully compared to full HartreeFock calculations in spherical symmetry in Refs.Pan13 (); Aym14 ():
(36)  
(37) 
such that the local asymmetry is given by . The radius parameters entering the density profile (37) are related to the equivalent hard sphere radii by
(38) 
and a similar relation holds for . The diffuseness parameters of the total density profile are assumed to depend quadratically on the bulk asymmetry , , where and were fitted from HF calculations in Ref.Pan13 ().
Using Eq.(37) the inmedium surface correction can be finally expressed as
(39)  
where and are the total and proton densities that correspond to the same cluster in the absence of the gas
(40)  
(41) 
and .
For the low temperatures which are of interest in the present study, the inmedium surface energy correction computed here is expected to give a small effect to the composition of the inner crust Rad14 (). The effect of the inmedium correction will therefore be estimated perturbatively. We assume that, for a given thermodynamic condition , the inmedium surface correction affects only slightly the gas density and composition, and consequently the chemical potentials. This correction will then be taken using the values for obtained from a NSE calculation where the inmedium effect is not considered. With this assumption, the modified binding energies solely depend on the cluster and on the thermodynamic condition and can therefore be simply added aposteriori to the energy density.
The final expression for the total baryonic energy density at finite temperature is then given by:
(42)  
where all terms depend on the temperature, and on the gas density and composition.
Cell  [fm]  [fm]  [fm]  [fm]  [fm]  

9  20  
22  28  
30  33  
33  36  
36  39  
41  42  
44  44  
46  46  
47  49  
48  54 
Iii 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 the seminal paper by Negele and Vautherin Neg73 (). These values cover the inner crust of the neutron star, approximately from the emergence of the neutron gas close to the drip point (cell 10) to a density close to the crustcore transition (cell 1), where bubbles and possibly other exotic nuclear shapes start to be formed. We recall that such structures are not included in our model. The corresponding values of the baryonic density, as well as the gas density, the proton fraction and the radius of the average WignerSeitz cell volume we obtain imposing the equilibrium condition, at the lowest temperature (T=100 KeV) considered in this study, are given in Table 1. We notice that the proton fraction increases, whereas the gas density decreases moving from cell 1 to cell 10. For comparison, proton fraction, gas density and radius of WignerSeitz cell obtained at T=0 in the full HFB calculation of Ref.For10 () are also given in the table. In Ref.For10 (), the same Sly4 parametrization was used in the calculations. As far as the gas density is concerned, the difference between the HFB values and our results, at the lowest temperature considered, are of the order of or less, except for the lowest densities; in that case however the gas contribution is negligible. It should also be noticed that the HFB calculations of Ref.For10 () adopt the same proton fraction of the representative calculations of Neg73 (), i.e. the equilibrium condition is not consistently implemented. The residual variation can be partly due to the different energetic description of the clusters. Our simplified mass model from Ref.Dan09 () is augmented of a phenomenological pairing term Gul15 () but does not contain shell effects. Neutron shell effects do not play any role above drip, but proton shell closures are known to be still effective at zero temperature in the inner crust Pea12 (), which can slightly affect the neutron gas density close to the drip condition. More important, the pairing interaction of this work is not the same as in Ref.For10 (). As explained in section II.1, we have fitted the parameters of the pairing interaction from abinitio BHF calculations of infinite neutron matter at zero temperature. This choice, also employed in Refs.Bur14 (); Pas15 (), is justified by the fact that the dominant pairing contribution comes from the unbound neutrons, which constitute, at the thermodynamic limit of the neutron star, an homogeneous neutron matter system. The clustergas interface, which is treated in the present work in the local density BCS approximation (see section II.4) , gives only a correction to this dominant term. Conversely, in the finite WignerSeitz calculation of Ref.For10 (), these parameters were fitted from finite nuclear properties San04 (). The resulting maximum pairing gap MeV is very close to the one displayed in Fig.1, but the density dependence of the pairing gap (see Ref.For10 ()) is different with respect to our calculation. Concerning the radius of the average volume of the WignerSeitz cell, again our results are in good agreement with HFB, except at the highest density (Cell1). As it will be shown further in the paper, this difference is due to the dominance of light resonances in our calculation, which are not included in a meanfield approach.
iii.1 Composition of the inner crust
Most thermodynamic calculations of the inner crust For10 (); Bar98 (); San04 (); Cha10 (); Pas15 () neglect the temperature variation of the proton fraction due to the temperature dependence of the chemical potentials entering the neutrinoless equilibrium condition:
(43) 
This variation, as obtained in our calculations, is shown in Fig.2, in four representative cells spanning the density and temperature interval concerned by this study. We can see that the change of the proton fraction is indeed very small close to the crustcore transition (up to cell 4), but it cannot be neglected at lower densities (cells 5 to 10). The density corresponding to the unbound neutron component is shown in Fig.3 for the same baryonic density conditions as in Fig.2.
The impact of the equilibrium condition on the gas density can be appreciated from the difference between the full and the dashed curve in Fig.3. It is clear from this result that the equilibrium condition has to be consistently implemented at each temperature. However, the most striking feature of Fig.3 is the clear discontinuity observed at the highest densities (up to cell 4 in the present calculation), corresponding to the transition point from superfluid to normal matter.
At first sight it is surprising to observe a density discontinuity, which is characteristic of first order phase transitions, at the superfluidnormal fluid transition, which is second order. This behavior is due to the fact that we are not observing an equation of state, that is at constant chemical potential, but a specific thermodynamic transformation implied by the minimization of the system total free energy. Specifically, one should consider that the pairing gap jumps, continously but suddenly, to zero at the critical temperature. This behavior influences the energetics of the system and may create discontinuities in the solution obtained for the gas density. This is particularly evident in the cells where, as in Cell 1, the gas density is larger than the value associated with the maximum gap (see Fig.1). In this case, the gas density solution corresponding to zero temperature in the full NSE calculation is lower than the gas density obtained neglecting the pairing interaction (dotted line in Fig.3), because it corresponds to a larger gap energy. As the temperature increases, in the regime where pairing is still active, the gas density decreases because the (negative) pairing contribution to the gas energy reduces. At the critical temperature the nonsuperfluid solution is recovered as it should. This corresponds to a higher density value, leading to a discontinuity.
The distribution of the cluster size as a function of the temperature is displayed in Fig.4.
We can see that at the lowest densities and temperatures the distribution is strongly peaked and can be safely 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, light particles systematically dominate at the highest temperatures. Such configuration cannot be addressed in meanfield based formalisms like HFB.
Fig.5 shows the cluster isotopic distribution, for the same four cells, at the temperature T=2 MeV. Results obtained neglecting the proton fraction variation imposed by the equilibrium condition are also shown. One can observe that, especially for the lowest density cells (Cells 7 and 10), the cluster asymmetry is significantly larger when the equilibrium condition is imposed. It is also interesting to notice that, at the high limits of the N/Z distribution, the yield is higher than the corresponding value obtained in absence of equilibrium. These extreme N/Z values are obtained from the lightest clusters, which dominate at the temperature considered (see Fig.4). Properly accounting for the equilibrium thus increases the contribution of the most unbound clusters.
From these results we can already anticipate that neglecting the temperature evolution of equilibrium will lead to a strong underestimation of the energy density, and the associated heat capacity, at high temperature.
iii.2 Energy and heat capacity
The variation with temperature of the energy density is displayed, for the same density conditions as in the previous figures, in Fig. 6.
The effect of the temperature dependence of the equilibrium condition can be appreciated comparing the full thin lines with the dashed lines. As expected, we can see that the temperature evolution of the proton fraction has a strong effect on the energy density, especially at the lowest densities.
Finally, the lines with symbols give the energy density of the most probable WignerSeitz cell, to be compared to the complete result (full lines) where the whole distribution of cells is taken into account. We can see that the effect of properly accounting for the cluster distribution is very important at the highest densities, but also at the lowest ones when the temperature gets higher. Indeed these situations are dominated by the emergence of light clusters. Close to the crustcore transition, the matter is so neutron rich that standard heavy clusters are not favored any more with respect to more exotic neutronrich forms of matter. As it can be seen from Fig.4, in this thermodynamic conditions the mass distribution extends up to but is dominated by light resonances at the limit of the nuclear binding (heavy hydrogen, helium, or lithium). The energy density associated with the full distribution is thus very different from the one associated with the most probable cluster. Specifically, the discontinuities observed in the most probable WignerSeitz cell in Fig.6 (upper left) appear at the temperatures where a transition occurs between the different elements.
It is important to remark that in this density region in principle nonspherical pasta phases, which are not included in the present work, could dominate over the light resonances. This is certainly true for low temperatures and matter close to isospin symmetry since the breaking of spherical symmetry leads to an important gain in binding energy Wat00 (). However finite temperature calculations in equilibrium Ava12 () tend to show that nonspherical pasta phases are only marginal, meaning that the energy behavior displayed in Fig.6 might be physical.
A similar transition, from heavy cluster dominated to light resonance dominated configurations, is observed at all densities. We recall that starting from cell 2 the density is too low for pasta phases to be present. This transition, leading to a sharp discontinuity in the energy density of the most probable WignerSeitz cell, physically corresponds to the melting of clusters inside a hot medium. In a meanfield treatment, cluster disappearence can only lead to a homogeneous medium, because small wavelength fluctuations cannot be treated in these approaches. However such fluctuations are entropically favored and naturally appear in the NSE treatment at high temperature.
The transition temperature from the superfluid to the normal fluid phase is signalled by a kink in the behavior of the energy density, which will lead to a peak in the associated heat capacity. This transition occurs at the same point in the full NSE calculation and considering only the most probable WignerSeitz cell. This can be understood from the fact that the electron and nucleon gases are uniform along the different cells, meaning that by construction the density and isospin caracteristics of the gas are the same in the two calculations. On the other hand, it is interesting to observe that a temperature shift could be observed if a standard calculation considering a single representative cell (SNA) was performed Gul15 (), as in the wellknown LattimerSwesty model Lat91 (). Indeed the baryonic density associated to the WignerSeitz cell of the most probable cluster is not the same as the total baryonic density of the distribution. This is a consequence of the fact that, especially at high temperature, the most probable cluster can be very different from the average cluster, thus it is very important to consider the full cluster distribution, as in the NSE calculations.
This also means that the consideration of the cluster distribution could modify the transition temperature as predicted by finite temperature HFB, though the effect is expected to be small.
Fig.7 shows the temperature behavior of the total baryonic energy derivative with respect to temperature, in four different cells. 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 have anticipated observing the energy density behavior of Fig.6, the temperature dependence of the equilibrium condition is seen to have 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, due to the rapid variation of the unbound component with temperature implied by the equilibrium condition (see Fig.3). 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. Indeed, in this case the equilibrium path favors a discontinuous trend of all thermodynamic quantities at the transition point (see Fig.3).
The LDA approximation was compared to HFB calculations in the case of trapped fermionic atoms in Gra03 (). It was shown that this approximation nicely works even in small systems at zero temperature, but it rapidly deteriorates at finite temperature. This is expected from the GinzburgLandau theory in cases where the critical temperature is much higher than the harmonic level spacing. In particular the LDA pairing field is seen to show a sudden drop at the surface, which is not apparent in the full HFB.
We however expect this limitation of LDA to be less severe in our case, because contrary to Ref. Gra03 () we do not use the LDA to solve the variational problem, but only to calculate the energy correction. Moreover our physical system is obviously not the same as in Gra03 (). We have verified that in our WignerSeitz cells the radial profile of the pairing field does not drop off but presents a decreasing tail, similar to HFB results.
Concerning the heat capacity, as shown in Fig.7, its quantitative value cannot be directly compared to the results of previous HFB works For10 (); Bar98 (); San04 (); Cha10 (); Pas15 () because of the different meanfield and/or pairing model, and because of the nonnegligible effect of the cluster distribution that we have observed in Fig.6. However, we have verified that the temperature location of the heat capacity peak, its height and width are almost identical to the results of Ref.For10 (), if we take the same parameters for the pairing interaction employed in that work. This is illustrated in Fig.8, where we represent the corresponding results for the heat capacity, obtained imposing the equilibrium condition or neglecting it (as in the HFB calculations). It is observed that the full curve compares rather well with the results of Ref.For10 (). It is also interesting to notice that, as already pointed out, the consideration of equilibrium induces non negligible effects on the .
iii.3 The effect of mass functionals
In all the calculations presented in the previous sections, we have systematically used the Skyrmebased liquiddrop formula, Eq.(26). This choice allows a consistent treatment of the bound and unbound matter component within the same energy functional. However, light clusters are systematically underbound with respect to heavier ones. To give an example, employing the parameters extracted in Dan09 () for Sly4, the binding energy of a particles is underestimated of while it is overestimated of for . This effect is even more dramatic for the most neutron rich light resonances, at the limit of nuclear binding, which can in principle be excited in the extremely neutron rich equilibrated matter of protoneutron stars at finite temperature: the last bound hydrogen isotope is according to the simplistic formula Eq.(26), while controlled extrapolations from experimental mass measurements predict that should be bound by MeV Aud12 ().
In Figs. 4 and 6 we have seen that at sufficiently high temperature, the last bound isotopes of light elements can become dominant in the composition of matter. It is therefore interesting to see how much these results depend on the poor energy description of light clusters of our mass formula. We have therefore repeated the same calculations, replacing Eq.(26) with the experimental value of the binding energy, whenever this value is known Aud12 (). By the very definition of the inner crust, all the nuclei populated with nonnegligible probability in the different density and temperature conditions explored in this work are beyond the dripline. This means that their experimental binding energy is typically not known, and Eq.(26) is still used for those nuclei in the new calculation. However experimental or extrapolated mass values exist for all bound isotopes of the lightest elements , and in that case the experimental value is used.
We find that, in the range of temperatures considered, the results are similar to the ones presented in Fig. 7 both for the highest (cells 12) and lowest (cells 7 to 10) densities. This means that the underbinding of light clusters does not influence the heat capacity calculation.
However, as shown in Fig.9, in the cells from 3 to 6 the situation is very different and the effect of accounting for the experimental binding energy of light clusters has a dramatic consequence. Indeed we see that accounting for the whole distribution of WignerSeitz cells, including the contribution of light clusters and resonances, modifies the height of the heat capacity peak and also its position in temperature. Moreover, an extra peak appears, which was not present in the calculations of Fig. 7. This peak corresponds to the “critical” temperature of the light clusters, depending on the thermodynamical conditions of the cell. This effect is easy to understand: the temperature at which the nuclei melt into a gas of free particles and resonances depends on the energy of these latter. If resonances correspond to bound states, they will dominate over the standard nuclei component at much lower temperatures than if they lie high in the continuum. The dominance of light clusters and resonances induces a change in the temperature dependence of the energy density, leading to an additional peak in the . To better illustrate this point, Fig.10 shows the cluster distribution obtained considering the experimental binding energies, whenever available, in the case of cells 3 and 4, where the second peak in the heat capacity is observed. One can appreciate that the cluster distribution is quite different with respect to the results shown in Fig.4. Moreover, we observe that the location of the second peak of the heat capacity, shown in Fig.9, coincides with the temperature where the cluster distribution starts to be dominated by light clusters. It should also be noticed that 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.
However, a few words of caution, about employing the experimental masses, are in order. All the calculations presented in this chapter have been obtained including in the statistical weight of the clusters the bulk part of the inmedium free energy shift (Eq.(16)), while the surface contribution (Eq.(35)) has been added aposteriori perturbatively, in order to consider the smearing effect of the density distribution on the pairing field. Contrary to bulk inmedium effects which increase the binding energy of the cluster, surface interaction with the surrounding gas is strongly dependent on the cluster asymmetry, as well as on the density and proton fraction of the gas. Surface inmedium shifts for very neutron rich species immersed in a neutron gas tend in particular to decrease the binding energy of the cluster Rad14 (); it is therefore possible that a selfconsistent inclusion of this energy term in the statistical calculation will reduce the contribution of the light resonances. Moreover, our local density BCS approximation to evaluate the pairing contribution of an inhomogeneous density distribution, including the population of light resonances, is certainly a quite crude approximation.
The effects of the surface corrections, as described in Section II.4, are evidenced in Fig.11. The comparison between the full and dashed lines allows appreciating the importance of the density fluctuations inside the WignerSeitz cells. In the calculation illustrated by dashed lines, the total energy is simply given by the sum of the cluster and uniform gas component according to Eq.(30), as suggested in early papers Lat94 (); Piz02 (). The energy contribution of the clustergas interface, according to Eq.(42), is considered in the full NSE results, shown by full lines.
It should be noticed that the effect of density fluctuations is never negligible, but still represents a correction of the total energy density, thus globally justifying the perturbative treatment developed in section II.4. As mentioned above, further corrections could be necessary in the case of very neutron rich species immersed in a neutron gas Rad14 (). This point is currently under study.
In the present calculations we observe a small, though appreciable, effect of the inmedium corrections on the heat capacity. In particular, we notice that, especially in the calculations neglecting the surface effects (dashed line), the transition temperature from superfluid to normal matter is very sharply defined. This is due to the fact that the gas density is characterized by a single value at each temperature point. This artificial feature comes from the neglect of the neutron density distribution inside the WignerSeitz cells. Our procedure to introduce surface corrections can partially cure this problem and leads to the results represented by the full lines. We can see that the transition temperature is smeared, as expected, and as it is observed in HFB calculations For10 (). Moreover a third small peak appears, in cell 3, at 1.8 MeV, due to the disappearance of pairing effects on the surface of the clusters.
Iv Conclusions
In this paper we have presented a calculation of the heat capacity in the inner crust of protoneutron star, within an approach based on cluster degrees of freedom that considers the complete distribution of different nuclear species in thermal and equilibrium. Superfluidity is taken into account including the pairing contribution of the homogeneous unbound neutron component in the BCS approximation. A standard pairing interaction is employed, with parameters fitted such as to reproduce the pairing gap of infinite neutron matter, as calculated from abinitio BruecknerHartreeFock calculation. A nonrelativistic Skyrme energy functional is used for the meanfield part of the unbound nucleon energy, as well as for the energy functional of the clusters. In this modelization, interparticle interactions are explicitly accounted for the unbound particles, and implicitly for the bound ones through the cluster energy functional. Interactions between bound and unboud particles are taken into account by considering the modification of the cluster surface due to the presence of the nucleon gas, which in turn affects the pairing properties. To this purpose, we have introduced an interface correction calculated in the local density BCS approximation. The resulting heat capacity appears compatible with complete HFB calculations as far as the location and the width of the peak associated with the transition from superfluid to normal matter is concerned, if the same pairing interaction is employed and the same treatment of the equilibrium condition is performed. Indeed we show that an accurate treatment of equilibrium is important for a quantitative determination of the heat capacity, and consequently the neutron star cooling curve. Specifically, accounting for the temperature dependence of proton fraction is seen to modify the energy density in a sizeable way, and to sharpen the phase transition peak close to the crustcore interface.
The added value of the present semiclassical modelling with respect to more sophisticated HFB calculation in the WignerSeitz cell is the consideration of the full distribution of different nuclear species at finite temperature with their proper statistical weight. We show that this feature considerably affects the heat capacity. In particular, the cluster disappearence at high temperature does not lead in this model to a uniform gas of nucleons, but correlations are still present in the form of exotic neutronrich resonant states at the limit of nuclear binding. This feature may lead to the appearance of an extra peak in the heat capacity, corresponding to the cluster “critical” temperature, at which the nuclei melt into a gas of free particles and resonances. Moreover, treating the pairing properties of this inhomogeneous matter in the local density approximation leads to the prediction of a second peak in the heat capacity associated to the superfluidnormal fluid transition of this clusterized matter. A similar feature was already reported in the literature For10 (); Pas15 () with HFB: a second peak is observed when the critical temperature is attained for the cluster surface. Further calculations within a more microscopic treatment are needed to confirm this finding.
Acknowledgments
This work has been partially funded by NewCompstar, COST Action MP1304.
References
 (1) J.M. Lattimer, K.A. Van Riper, M. Prakash, M. Prakash, The Astrophysical Journal 425, 802 (1994).
 (2) O.Y. Gnedin, D.G. Yakovlev, A.Y. Potekhin, MNRAS 324, 725 (2001).
 (3) D. Page, S. Reddy, Neutron Star crust. edited by C. Bertulani and J. Piekarewicz, Nova Science Publishers, 281 (2012).
 (4) J. Piekarewicz, F.J. Fattoyev, C.J. Horowitz, Phys. Rev. C 90, 015803 (2014).
 (5) M. Fortin, F. Grill, J. Margueron, D. Page, N. Sandulescu, Phys. Rev. C 82, 065804 (2010).
 (6) F. Barranco, R.A. Broglia, H. Esbensen, E. Vigezzi, Phys. Rev. C 58, 1257 (1998).
 (7) N. Sandulescu, N.V. Giai, R.J. Liotta, Phys. Rev. C 69, 045802 (2004).
 (8) N. Chamel, S. Goriely, J.M. Pearson, M. Onsi, Phys. Rev. C 81, 045804 (2010).
 (9) A. Pastore, Phys. Rev. C 91, 015809 (2015); A. Pastore, N. Chamel, J. Margueron, Mon. Not. Roy.As.Soc. 448, 1887 (2015).
 (10) J.W. Negele, D. Vautherin, Nucl. Phys. A207, 298 (1973).
 (11) S.S. Avancini, D.P. Menezes, M.D. Alloy, J.R. Marinelli, M.M.W. Moraes, C. Providência, Phys. Rev. C 78, 015802 (2008).
 (12) G. Ropke, N.U. Bastian, D. Blaschke, T. Klahn, S. Typel, H.H. Wolter, Nucl. Phys. A897, 70 (2013).
 (13) S.S. Avancini, C.C. Jr. Barros, L. Brito, S. Chiacchiera, D.P. Menezes, C. Providência, Phys. Rev. C 85, 035806 (2012); S.S. Avancini, S. Chiacchiera, D.P. Menezes, C. Providência, Phys. Rev. C 85, 059904 (2012).
 (14) S. Heckel, P.P. Schneider, A. Sedrakian, Phys. Rev. C 80, 015805 (2009).
 (15) M. Hempel and J. SchaffnerBielich, Nucl. Phys. A837, 210 (2010); M. Hempel, T. Fischer, J. SchaffnerBielich, M. Liebendrfer, Astrophys. J. 748, 70 (2012).
 (16) Ad.R. Raduta, F. Gulminelli, Phys. Rev. C 82, 065801 (2010).
 (17) Ad.R. Raduta, F. Gulminelli, F. Aymard, Eur. Phys. J. A50, 24 (2014).
 (18) S. Furusawa, K. Sumiyoshi, S. Yamada, H. Suzuki, Astrophys. J. 772, 95 (2013).
 (19) N. Buyukcizmeci et al., Nucl. Phys. A907, 13 (2013).
 (20) F. Gulminelli, Ad.R. Raduta, nuclth/ArXiv:1504.04493.
 (21) W.D. Myers, W.J. Swiatecki, Nucl. Phys. A336, 267 (1980).
 (22) M. Centelles, M.D. Estal, Vinas X., Nucl. Phys. A635, 193 (1998).
 (23) M. Warda, X. Vinas, X. RocaMaza, M. Centelles, Phys. Rev. C 80, 024316 (2009).
 (24) P. Papakonstantinou, J. Margueron, F. Gulminelli, Ad. R. Raduta, Phys. Rev. C88, 045805 (2013).
 (25) N. Chamel, Phys. Rev. C 82, 014313 (2010).
 (26) S. Burrello, M. Colonna, F. Matera, Phys. Rev. C 89, 057604 (2014).
 (27) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, R. Schaeffer, Nucl. Phys. A635, 236 (1998).
 (28) L.G. Cao, U. Lombardo, P. Schuck, Phys. Rev. C 74, 064301 (2006).
 (29) P. Danielewicz, J. Lee, Nucl. Phys. A818, 36 (2009).
 (30) T. Von Egidy, D. Bucurescu, Phys. Rev. C 72, 044311 (2005); ibid., Phys. Rev. C 73, 049901 (2006).
 (31) P.M. Pizzochero, F. Barranco, E. Vigezzi and R.A. Broglia, Astrophys. J. 569, 381 (2002).
 (32) S. Typel, G. Roepke, T. Klahn, D. Blaschke, H.H. Wolter, Phys. Rev. C 81, 015803 (2010).
 (33) F. Aymard, F. Gulminelli, J. Margueron, Phys. Rev. C 89, 065807 (2014).
 (34) S. Typel, Phys. Rev. C89, 064321 (2014).
 (35) J.M. Pearson, N. Chamel, S. Goriely, C. Ducoin, Phys. Rev. C 85, 065803, (2012).
 (36) J.M. Lattimer, F.D. Swesty, Nucl. Phys. A535, 331 (1991).
 (37) G. Watanabe, K. Iida, K. Sato, Nucl. Phys. A676, 455 (2000).
 (38) G. Audi, M. Wang, A.H. Wapstra, F.G. Kondev, M. MacCormick, X. Xu and B. Pfeiffer, Chinese Physics C36, 1287 and Chinese Physics C36, 1603 (2012).
 (39) M. Grasso, M. Urban, Phys. Rev. A 68, 033610 (2003).