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

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

S. Burrello INFN, Laboratori Nazionali del Sud and Dip. di Fisica e Astronomia, Università di Catania, 95123 Catania, Italy    F. Gulminelli CNRS and ENSICAEN, UMR6534, LPC, 14050 Caen cédex, France    F. Aymard CNRS and ENSICAEN, UMR6534, LPC, 14050 Caen cédex, France    M.Colonna INFN, Laboratori Nazionali del Sud, 95123 Catania, Italy    Ad.R.Raduta NIPNE, Bucharest-Magurele, POB-MG6, Romania
July 13, 2019

Superfluidity in the crust is a key ingredient for the cooling properties of proto-neutron stars. Present theoretical calculations employ the quasi-particle mean-field Hartree-Fock-Bogoliubov theory with temperature dependent occupation numbers for the quasi-particle states.


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.


Following a recent work, the Wigner-Seitz cell is mapped into a model with cluster degrees of freedom. The finite temperature distribution is then given by a statistical collection of Wigner-Seitz 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.


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.


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.


I 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 X-ray transients Pag12 (), as well as in the understanding of glitches Pie14 (). Moreover, it is well-known 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 non-negligible 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 Hartree-Fock-Bogoliubov (HFB) equations in the so-called Wigner-Seitz approximation Neg73 (), meaning that the assumption is done 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 conditions. Moreover, at the extremely low proton fractions associated to the inner crust, deformed nuclear structures and beyond drip-line light nuclear resonances can participate to the statistical equilibrium, and might be too exotic to be well described through standard mean-field calculations. Non-spherical pasta structures, which are not accessible to the spherical mean-field, 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 beyond-mean 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 self-consistently under conditions of statistical equilibrium. In some of these models both the gas-cluster interaction and the self-interaction of the gas are included, though within semi-classical approximations Gul15 (). In particular in Ref.Gul15 () it is shown that a proper definition of the cluster self-energies in the NSE cluster distribution allows recovering the zero temperature limit of a single Wigner-Seitz approximation in the (Extended) Thomas-Fermi 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 non-homogeneity 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 non-negligible effects of the cluster distribution are seen close to the drip point, and close to the crust-core 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 neutrino-less 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 in-medium 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 Wigner-Seitz 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 self-interacting 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 ():


where is the symmetry energy at saturation, is the surface stiffness coefficient extracted from a semi-infinite 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


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:


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

Figure 1: (Color online) pairing gap as a function of density for homogeneous neutron matter at zero temperature as obtained from Brueckner-Hartree Fock calculations (full line from  Cao06 ()). The figure also shows the energy gap deduced solving the BCS gap equations at finite temperature (symbols from Bur14 ()).

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):




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 spin-saturated matter and is the particle occupation number:


with and . and denote, respectively, chemical potential and reduced chemical potential, and


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 () :


where the parameters , , are fixed imposing to reproduce the pairing gap of pure neutron matter as obtained in Brueckner-Hartree-Fock 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


The reduced chemical potential is moreover obtained by fixing the particle number density:


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 mean-field approximation:


where the entropy density is given by :


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 ():


In this expression, is the electron free-energy density, is the total proton density and denotes the Wigner Seitz volume. is the free energy of the cluster immersed in the nucleon gas:


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:


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 in-medium effects.

The modification of the nuclear free energy consists of a bulk term


due to the presence of the gas in the same spatial volume, , occupied by the cluster, and a surface term which accounts for the isospin-dependent 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 in-medium effect.

The screening effect of the electron density which neutralizes the Wigner-Seitz cell leads to a modification of the cluster free energy according to:


with the Coulomb screening function in the Wigner-Seitz approximation:


The volume of the Wigner-Seitz cell associated to each nuclear species is univocally defined by the charge conservation constraint:


leading to


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 ():


The result is a NSE-like expression for the cluster multiplicities Gul15 ():


where the chemical potentials can be expressed as a function of the gas densities only:


The numerical solution of Eqs.(21),(22) for the two unknown , closes the model.

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 ():


with the asymmetry energy coefficient:


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 even-even (odd-odd) nuclei.

It is important to observe that this formula,Eq.(26), similar to any other mean-field 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 , mean-field 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 back-shifted 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 Wigner-Seitz cell and write .

In our model, the full distribution of clusters is accounted for and the volume fraction accessible to the gas results:


where the total volume can be written as , being the average size of the Wigner-Seitz volume and the total cluster multiplicity. Indicating with the normalized cluster multiplicity, the contribution of the gas to the energy density is therefore:


and the total baryonic energy density of star matter is


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 ():


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 Wigner-Seitz 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 non-superfluid, 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 in-medium 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 in-medium correction is given by the Coulomb screening by the electron gas, and by the Pauli-blocking 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 in-medium modification of the cluster energy can be computed by subtracting to the total energy in each Wigner-Seitz cell the contribution of the gas alone and of the nucleus alone, following Aym14 ():


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 beta-equilibrated 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):


A similar decomposition can be applied to the total LDA energy :


where is the radius of the Wigner-Seitz cell, is the hard-sphere 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 in-medium modification simply as


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 in-medium corrections were evaluated in Ref.Aym14 () adding to the LDA also higher orders in in the semiclassical Thomas-Fermi 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 spin-orbit terms are therefore neglected in the surface correction.

In order to evaluate the in-medium surface correction through Eq.(34), a model for the density profiles , has to be assumed. We use the simple Wood-Saxon analytical profiles proposed in Ref.Pan13 () and successfully compared to full Hartree-Fock calculations in spherical symmetry in Refs.Pan13 (); Aym14 ():


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


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 in-medium surface correction can be finally expressed as


where and are the total and proton densities that correspond to the same cluster in the absence of the gas


and .

For the low temperatures which are of interest in the present study, the in-medium surface energy correction computed here is expected to give a small effect to the composition of the inner crust Rad14 (). The effect of the in-medium correction will therefore be estimated perturbatively. We assume that, for a given thermodynamic condition , the in-medium 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 in-medium 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 a-posteriori to the energy density.

The final expression for the total baryonic energy density at finite temperature is then given by:


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
Table 1: From left to right are given: the total baryonic density, the proton fraction at T=100 keV, the proton fraction considered in Ref.For10 (), the gas density at T=100 keV, the gas density obtained at T=0 in Ref.For10 (), the radius of the average Wigner-Seitz volume calculated at T=100 keV and the radius of the cell at T=0 shown in Ref.For10 ().

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 crust-core 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 Wigner-Seitz 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 Wigner-Seitz 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 ab-initio 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 cluster-gas 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 Wigner-Seitz 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 Wigner-Seitz 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 mean-field 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:

Figure 2: Temperature evolution of the global proton fraction obtained by imposing the neutrinoless -equilibrium condition (43) for four representative cells.

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 crust-core 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.

Figure 3: (Color online) Temperature evolution of the unbound gas component in the same representative cells as in Fig. 2. Full line: complete NSE calculation. Dashed line: the value of the global proton fraction is assumed equal to the one calculated from -equilibrium at the lowest temperature, MeV. Dotted line: as the full line, but neglecting the pairing interaction.

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 superfluid-normal 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 non-superfluid 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.

Figure 4: (Color online) Normalized cluster size distribution at different temperatures in the same four representative cells as in Fig. 2 and at -equilibrium.

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 mean-field based formalisms like HFB.

Figure 5: Normalized isotopic cluster distribution in the same four representative cells as in Fig.2, as obtained at the highest temperature considered, MeV. The different symbols indicate results at -equilibrium (triangle) or assuming a global proton fraction equal to that one corresponding to -equilibrium but at the lowest temperature, MeV.

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

Figure 6: (Color online) Temperature evolution of the baryonic energy density for the same representative cells as in Fig.2. Full line: complete NSE calculation. Dashed line: as the full line, but the value of the global proton fraction is assumed equat to that one calculated from -equilibrium at the lowest temperature, MeV. Lines with symbols: as the full line, but the NSE distribution is replaced with the most probable Wigner-Seitz cell.

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 Wigner-Seitz 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 crust-core transition, the matter is so neutron rich that standard heavy clusters are not favored any more with respect to more exotic neutron-rich 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 Wigner-Seitz 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 non-spherical 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 non-spherical 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 Wigner-Seitz cell, physically corresponds to the melting of clusters inside a hot medium. In a mean-field 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 Wigner-Seitz 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 well-known Lattimer-Swesty model Lat91 (). Indeed the baryonic density associated to the Wigner-Seitz 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.

Figure 7: (Color online) Temperature evolution of the heat capacity for the same representative cells as in Fig.2. Full line: complete NSE calculation. 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, MeV.

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 Ginzburg-Landau 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 Wigner-Seitz 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 mean-field and/or pairing model, and because of the non-negligible 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 .

Figure 8: Temperature evolution of heat capacity for Cell 6, obtained considering a pairing interaction with the same parameters of Ref.For10 (). Full line: complete NSE calculation. Dashed line: as the full line, but the value of the global proton fraction is assumed equal to that one given in Ref.For10 () (see also Table I). The inset shows a zoom at low temperature in order to facilitate the comparison with the results of Ref.For10 ().

iii.3 The effect of mass functionals

In all the calculations presented in the previous sections, we have systematically used the Skyrme-based liquid-drop 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 proto-neutron 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 non-negligible 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.

Figure 9: (Color online) Temperature evolution of heat capacity for four intermediate cells. Full lines: complete NSE calculation making use of the analytical expressions, Eq.(26), for binding energies (labeled as LDM in the figure). Line with symbols: as the full line, but experimental binding energies, from Aud12 (), are used whenever available.

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 1-2) 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 Wigner-Seitz 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.

Figure 10: (Color online) Normalized cluster size distribution at different temperatures, as obtained from complete NSE calculation using experimental binding energies, from Aud12 (), whenever available. Results are shown for two representative cells where the transition to a dominance of light clusters is observed in the range of temperatures considered.

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 in-medium free energy shift (Eq.(16)), while the surface contribution (Eq.(35)) has been added a-posteriori perturbatively, in order to consider the smearing effect of the density distribution on the pairing field. Contrary to bulk in-medium 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 in-medium 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 self-consistent 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.

Figure 11: (Color online) Temperature evolution of heat capacity for the same two cells as in Fig.10. Full line with symbols: complete NSE calculation using experimental binding energies, from Aud12 (), whenever available. Dashed line: as the full line but neglecting in-medium effects.

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 Wigner-Seitz 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 cluster-gas 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 in-medium 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 Wigner-Seitz 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 proto-neutron 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 ab-initio Brueckner-Hartree-Fock calculation. A non-relativistic Skyrme energy functional is used for the mean-field 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 crust-core interface.

The added value of the present semiclassical modelling with respect to more sophisticated HFB calculation in the Wigner-Seitz 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 neutron-rich 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 superfluid-normal 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.


This work has been partially funded by New-Compstar, COST Action MP1304.


  • (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. Schaffner-Bielich, Nucl. Phys. A837, 210 (2010); M. Hempel, T. Fischer, J. Schaffner-Bielich, 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, nucl-th/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. Roca-Maza, 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).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description