Statistical description of complex nuclear phases in supernovae and proto-neutron stars

Statistical description of complex nuclear phases in supernovae and proto-neutron stars


We develop a phenomenological statistical model for dilute star matter at finite temperature, in which free nucleons are treated within a mean-field approximation and nuclei are considered to form a loosely interacting cluster gas. Its domain of applicability, that is baryonic densities ranging from about g cm to normal nuclear density, temperatures between 1 and 20 MeV and proton fractions between 0.5 and 0, makes it suitable for the description of baryonic matter produced in supernovae explosions and proto-neutron stars. The first finding is that, contrary to the common belief, the crust-core transition is not first order, and for all subsaturation densities matter can be viewed as a continuous fluid mixture between free nucleons and massive nuclei. As a consequence, the equations of state and the associated observables do not present any discontinuity over the whole thermodynamic range. We further investigate the nuclear matter composition over a wide range of densities and temperatures. At high density and temperature our model accounts for a much larger mass fraction bound in medium nuclei with respect to traditional approaches as Lattimer-Swesty, with sizeable consequences on the thermodynamic quantities. The equations of state agree well with the presently used EOS only at low temperatures and in the homogeneous matter phase, while important differences are present in the crust-core transition region. The correlation among the composition of baryonic matter and neutrino opacity is finally discussed, and we show that the two problems can be effectively decoupled.

21.65.Mn, 24.10.Pa, 26.50.+x, 26.60.-c

March 21, 2018

I I. Introduction

Nuclear matter is not only a theoretical idealization providing benchmark studies of the effective nuclear interaction, but it is also believed to constitute the major baryonic component of massive objects in the universe, as exploding supernovae cores and neutron stars. The structure and properties of these astrophysical objects at baryonic densities exceeding normal nuclear matter density is still highly speculative (1); (2); (3); (4); (5); (6); (7); (8). Conversely for subsaturation densities it is well established that matter is mainly composed of neutrons, protons, electrons, positrons and photons in thermal and typically also chemical equilibrium (9); (10). Depending on the thermodynamic condition, neutrinos and anti-neutrinos can also participate to the equilibrium.

Such composite matter is subject to the contrasting couplings of the electromagnetic and the strong interaction. Because of the electron screening, the two couplings act on comparable length scales giving rise to the phenomenon of frustration (11); (12); (13); (14), well-known in condensed matter physics (15). Because of this, a specific phase diagram, different from the one of nuclear matter and including inhomogeneous components, is expected in stellar matter (16).

Many theoretical studies exist at zero temperature. In a cold neutron star, going from the dilute crust to the dense core, a transition is known to occur from a solid phase constituted of finite nuclei on a Wigner lattice immersed in a background of delocalized electrons and neutrons, through intermediate inhomogeneous phases composed of non-spherical nuclei (pasta phases), to a liquid phase composed of uniform neutrons, protons and electrons (17); (18); (19); (20); (21); (22).

At finite temperature the matter structure and properties are not as well settled. The most popular phenomenological approaches are the Lattimer-Swesty (23) (LS) and the Shen (24) equation of state, recently updated in Ref. (25). In these standard treatments currently used in most supernovae codes, the dilute stellar matter at finite temperature is described in the baryonic sector as a statistical equilibrium between protons, neutrons, alphas and a single heavy nucleus. The transition to homogeneous matter in the neutron star core is supposed to be first order in these modelizations and obtained through a Maxwell construction in the total density at fixed proton fraction.

It is clear that such single nucleus approximation (SNA) is highly schematic and improvements are possible. Concerning integrated quantities as thermodynamic functions and equations of state, such variables may be largely insensitive to the detailed matter composition (26), though we show in this paper that this is not always the case. However it is also known that the composition at relatively high density close to saturation, together with the pressure and symmetry energy, governs the electron capture rate, which in turn determines the proton fraction at bounce and the size of the homologous core, a key quantity to fix the strength of the shock-wave and the output of the supernovae explosion (27); (28); (29); (30); (31); (32). Moreover the composition may also affect the nucleosynthesis of heavy elements, which is still poorly understood (33); (34); (35); (36), as well as the neutrino scattering through the core after bounce (37); (38), and the cooling rate of neutron stars (39); (40). For these reasons, in the recent years, many efforts have been done to improve over the simplistic representation of stellar matter given by the SNA approach.

The different modelizations which consider a possible distribution of all different nuclear species inside dilute stellar matter are known under the generic name of Nuclear Statistical Equilibrium (NSE) (41); (42); (43); (44). The basic idea behind these models is the Fisher conjecture that strong interactions in dilute matter may be entirely exhausted by clusterization (45). In these approaches stellar matter in the baryonic sector is then viewed as a non-interacting ideal gas of all possible nuclear species in thermal equilibrium. The result is that thermodynamic quantities like entropies and pressure appear very similar to the ones calculated with standard approaches, while noticeable differences are seen in the matter composition. In particular an important contribution of light and intermediate mass fragments is seen at high temperature, which is neglected in standard SNA approaches.

The strongest limitation of NSE-based approaches is that they completely neglect in-medium effects, which are known to be very important in nuclear matter. Since the only nuclear interactions are given by the cluster self-energies, the homogeneous matter composing the neutron star core cannot be modelized, nor it can be the phenomenon of neutron drip in the inner crust, well described by mean-field models (46). As a consequence, these models cannot be applied at densities close to saturation and the crust-core transition cannot be described.

To overcome this problem, different microscopic (47); (48); (49); (50) as well as phenomenological (51) approaches have been developed in the very recent past. In this paper we would like to introduce a phenomenological model that treats the nuclei component within an improved NSE, while it describes the unbound protons and neutrons in the finite temperature Hartree-Fock approximation.

The plan of the paper is as follows. The first part of the paper is devoted to the description of the model. The clusterized component, the homogeneous component, the properties of the mixture and the lepton sector are described in successive sections together with their thermodynamic properties. A particular attention is devoted to the modelization of the crust-core transition. We show that the inclusion of excluded volume is sufficient to describe the transition from the clusterized crust to the homogeneous core, and that this transition is continuous. Different generic as well as specific arguments are given against the possibility of a first-order transition. The second part of the paper gives some results relevant for the star matter phenomenology. The first section shows observables following constant chemical potential paths, in order to connect the observables with the properties of the phase diagram. Then the behavior of the different quantities for constant proton fractions is displayed in order to compare with more standard treatments of supernova matter. Finally the last section addresses the problem of neutrino trapping and the interplay between the matter opacity to neutrinos and matter composition. Conclusions and outlooks conclude the paper.

Ii II. The model

The model aims to describe the thermal and chemical properties of nuclear matter present in supernovae and (proto)-neutron stars at densities ranging from the normal nuclear density to , temperatures between 0 and 20 MeV and proton concentration between 0.5 and 0. In this regime, the star matter typically consists of a mixture of nucleons, light and heavy nuclear clusters, neutrinos (if we consider the thermodynamic stage where neutrinos are trapped), photons and a charge neutralizing background of electrons and positrons.

As there is no interaction among electrons, neutrinos, photons and nuclear matter, the different systems may be treated separately and their contributions to the global thermodynamic potential and equations of state added-up.

In the grancanonical ensemble this reads,


where the grancanonical potential,


is the Legendre transformation of the entropy with respect to the fixed intensive variables. In the previous equations is an arbitrary macroscopic volume and is the grancanonical partition sum.

The observables conjugated to the ones fixed by the reservoir and geometry can be immediately calculated as partial derivatives of . Thus, the total energy density is


the different particle densities are


where , and finally the total pressure is


ii.1 A.The baryon sector

The light and heavy nuclei are assumed to form a gas of loosely-interacting clusters which coexist in the Wigner-Seitz cell with a homogeneous background of delocalized nucleons. To avoid exceeding the normal nuclear density and naturally allow for homogeneous-unhomogeneous matter transition, nuclei and nucleons are forbidden to occupy the same volume. In the following we start describing the modelization of these two components separately, and we turn successively to the properties of the mixture obtained when the two are supposed to be simultaneously present in the Wigner-Seitz cell.

1. The homogeneous nuclear matter component

Mean-field models constitute a natural choice for approaching interacting particle systems. By introducing a mean-field potential, the physical problem is reduced to the simplified version of a system of non-interacting particles. Effective nucleon-nucleon interactions allow one to express the system average energy as a simple single-particle density functional, and to cast the nuclear matter statistics in a way which is formally very similar to an ideal Fermi gas(52).

The mean field energy density of an infinite homogeneous system is a functional of the particle densities and kinetic densities for neutrons () or protons (). At finite temperature, the mean field approximation consists in expressing the grancanonical partition function of the interacting particle system as the sum of the grancanonical partition function of the corresponding independent-particle system associated to the mean-field single-particle energies, with the temperature weighted difference between the average single-particle energy () and the mean-field energy


The one body partition sum is defined as:


and can be expressed as a function of the neutron and proton kinetic energy density :




In these equations the effective chemical potential includes the self energies according to , are the neutron (proton) effective mass, and the factor comes from the spin degeneracy.

Equation (9) establishes a self-consistent relation between the density of q-particles and their chemical potential . Introducing the single particles energies , the above densities can be written as regular Fermi integrals by shifting the chemical potential according to The Fermi-Dirac distribution indeed reads:


Eqs. (11) and (9) define a self-consistent problem since depends on the densities. For each couple a unique solution is found by iteratively solving the self-consistency between and . Then Eq. (10) is used to calculate .

At the thermodynamic limit the system volume diverges together with the particle numbers , , and the thermodynamics is completely defined as a function of the two particle densities , or equivalently the two chemical potentials .

With Skyrme based interactions, the energy density of homogeneous, spin saturated matter with no Coulomb effects is written as:


where are Skyrme parameters.

Several Skyrme potentials have been developed over the years for describing the properties of both infinite nuclear matter and atomic nuclei, and address the associated thermodynamics. It was thus in particular shown that nuclear matter manifests liquid-gas like first-order phase transitions up to a critical temperature (53); (54); (55); (56); (57).

In order to make direct quantitative comparisons between our model and the one of Lattimer and Swesty (23), through this paper nucleon-nucleon interactions are accounted for according to the SKM* parameterization (58) if not explicitely mentioned otherwise. Table 1 summarizes the force parameters while the main properties of nuclear matter are summarized in Table 2.

Parameter Value
(MeV fm) -2645
(MeV fm) 410
(MeV fm) -135
(MeV fm) 15595
Table 1: SKM* force parameters (58).
Observable Value
(MeV) -15.78
(MeV) 216.7
(MeV) 30.03
(fm) 0.1603
(MeV) 17.51
Table 2: Infinite nuclear matter and surface properties of SKM* force (58).

The thermodynamic properties of the system are best studied introducing the constrained entropy:


where the entropy density is nothing but the Legendre transform of the mean field partition sum:


Instabilities are then recognized as local convexities of the constrained entropy at a given temperature in the two-dimensional density plane . In the presence of an instability, the system entropy can be maximized by phase mixing, which corresponds to a linear interpolation in the density plane. This geometrical construction corresponds to the well-known Gibbs equality between all intensive parameters of two coexisting phases in equilibrium:


Phase diagrams of SKM* nuclear matter are illustrated in Figs. 1, 2 and 3 for three values of temperature 1.6, 5 and, respectively, 10 MeV in (left panels) and (right panels) representations. As it is customary in nuclear matter calculations, the system is charge neutral in the sense that the electromagnetic interaction has been artificially switched off. In the star matter case the proton electric charge has to be taken into account, and compensated by the electron charge to give a net zero charge. The resulting Coulomb interaction energy will be discussed in the next chapter. Full lines correspond to the phase coexistence region and dashed lines mark the borders of the instability region, the so-called spinodal surface.

As one may notice, a large portion of the density and chemical potential plane is characterized by the instability even at temperatures as high as 10 MeV. Inside the phase coexistence in the density representation, the mean field solution is unphysical at thermal equilibrium if matter is uniquely composed of homogeneously distributed (uncharged) protons and neutrons. In the case of stellar matter, the entropy convexity properties have to be examined only after considering all the different constituents, which can very drastically change the properties of the phase diagram. As we will show in the next sections, the introduction of electrons has the effect of making the Gibbs construction unphysical, while the introduction of clusterized matter has the effect of stabilizing the unstable mean field solution.

Figure 1: (Color online) (left) and (right) representations of the phase diagrams of the charge-neutral homogeneous infinite nuclear matter described by the SKM* interaction (lines) and net-charge neutralized clusterized nuclear matter (stars) corresponding to =1.6 MeV. In the case of the uniform system, the dashed lines mark the borders of the spinodal zone, while the solid line marks the coexistence region. These results are independent of the chosen confining volume fm except in the low density region of the phase diagram of the clusterized matter system in coordinates, which is affected by finite size effects (see text).
Figure 2: (Color online) The same as in Fig. 1 for =5 MeV. For symbol and line codes, see Fig. 1.
Figure 3: (Color online) The same as in Fig. 1 for =10 MeV.

2. The clusterized nuclear matter component

As we have discussed in the previous section, in a wide region of temperature and density, uncharged uniform nuclear matter is unstable with respect to density fluctuations. If such fluctuations occur on a macroscopic scale, a first-order liquid-gas phase transition follows. Inside the spinodal surface however, finite wavelength fluctuations may also occur, leading to spontaneous cluster formation. Spinodal decomposition and nucleation are indeed known to constitute the dynamical mechanisms leading to phase separation in macroscopic uncharged systems. As we will show in the following, in the presence of the repulsive long range Coulomb interaction among protons, which is screened by the leptons only on macroscopic scales, phase separation is quenched in stellar matter and the instability is cured by the formation of nuclear clusters, that is clusters in the femtometer scale.

One possibility to account for the thermal and phase properties of a clusterized matter sub-system is to use a statistical model with cluster degrees of freedom. Several such models have been so far proposed for the study of condensation close to the critical point (45), nuclear multifragmentation (59); (60); (61); (62) and compact stars (42); (43); (44). They basically consist in the estimation of the number of microscopic states compatible with the thermodynamic macroscopic constraints.

Considering that the center of mass of the clusters can be treated as classical degrees of freedom, the grancanonical partition function of the clustered system reads:


Here denotes a generic cluster configuration defined by the mass number , the proton number and the excitation energy of each constituent nuclear state indexed by , whose binding energy is and level density is ; represents the inter-cluster configuration dependent interaction energy, and is the volume accessible to cluster which is not excluded by the presence of the other clusters. Because of the short range of the nuclear force, this excluded volume correction can be viewed as a first-order approximation to the nuclear part of the interaction energy in the spirit of the Van der Waals gas. More sophisticated in-medium corrections can in principle be also added as modifications of the cluster self-energies (50). Because of this, we will consider that the interaction energy only contains the Coulomb electrostatic energy of the configuration.

To obtain the so-called NSE approach, one has to neglect all inter-particle interactions except the cluster self-energies, leading to :


The advantage of this simplification is that a completely factorized expression for the partition sum can be obtained, leading to analytical expressions for all thermodynamic quantities. Indeed any configuration can be ordered as


where gives the multiplicity of clusters of type . Equation (17) becomes:




Eq. (19) is not yet an analytically tractable expression, because of the infinite number of possible excited states for each nuclear isotope . The internal degrees of freedom can be however integrated over giving


where is the standard grancanonical partition sum for a single nucleus of mass and charge :


and is a temperature dependent degeneracy.

While elegant and analytically tractable, Eq. (21) represents a correct estimation of Eq. (16) if and only if the interactions are fully exhausted by clusterization.

For this reason, in this work we will consider only Eq. (16), which we will estimate numerically with a precise, efficient and well-tested Metropolis-Monte-Carlo technique (62), and calculate relevant thermodynamic observables other than the ones fixed by the exterior as ensemble averages. For the configuration-defined quantities, =, , :

the corresponding ensemble averaged density reads:


and stability is systematically checked against an increase of the considered volume .

For each sampled configuration , we evaluate the Coulomb interaction energy in the Wigner-Seitz approximation (17),




accounting for the screening effect of electrons. denotes the proton density inside the nuclei and is the electron density. Net charge neutrality of the system imposes the electron density to be equal to the proton density .

The total pressure results:


Here originates from clusters translational motion and is the so-called lattice Coulomb pressure.

Concerning the excluded volume, the result of the integrals over the coordinate space corrected for the volume excluded by each cluster, , can be re-casted as a global configuration-dependent factor as,


For the clusters we assume the following properties:

  • clusters binding energies are described by a liquid-drop parameterization,


    with =15.4941 MeV, =17.9439 MeV, (64) where the surface term is additionally made to vanish at the critical temperature =12 MeV (63); the Coulomb energy contribution is accounted for in as discussed above; structure effects thoroughly considered by other authors (51) are disregarded because of the relatively high temperatures we are focusing on, but they can be straight-forwardly included;

  • in order to avoid multiple counting of the free nucleons, already considered in the uniform background, the clusters should have ; this arbitrary limit may be replaced by , such as to allow also for H and H, but this choice is not expected to be important for the presently considered quantities;

  • no limit is imposed for the largest cluster which, in principle, may reach the total system size; given the fact that the fragment mass partition enters in Eq. (16) not only in the translational energy factor and binding energy, but also in the calculations of excluded volume and Coulomb energy, we expect it to play a dramatic role in the limit of low temperatures and densities close to ;

  • to allow creation of exotic species, there is no restriction for the neutron/proton composition; a cluster is allowed to exist as soon as ;

  • cluster internal excitation energy is upper limited by the cluster binding energy and the corresponding level density is of Fermi-gas type with cut-off correction (65),


    with MeV and =9 MeV.

The phase diagram of clusterized stellar matter can be inferred from the bimodal behavior of the total particle number in the grancanonical ensemble (63). This model presents a first-order liquid-gas phase transition similar to the case of uniform nuclear matter, but contrary to the latter it is never unstable. The two phase diagrams are compared for =1.6 and 5 MeV in Figs. 1 and, respectively, 2. Fig. 3 presents exclusively the phase diagram of neutral uniform matter at =10 MeV, as this temperature is supra-critical for the clusterized matter.

We can see that the thermodynamics of a clusterized system is very different from the one of neutral nuclear matter in the mean field approximation, even though the employed effective interaction parameters are typically tuned on the same experimentally accessible observables of cold nuclei close to saturation density, which are implemented in the cluster energy functional. As a consequence, the two phase diagrams are only compatible at low temperatures and at densities close to saturation.

The extension of the coexistence region in density and temperature is noticeably reduced for clusterized matter. This effect is partially not physical. Indeed the numerical procedure adopted for estimating Eq. (16) implies that clusterized matter results depend on the system size for the lowest densities: the lowest accessible density is determined by the chosen finite volume and the minimum allowed cluster size () eventually multiplied by the minimum allowed number of clusters (). Clusterized matter phase diagrams considered in this section have been obtained working within a cell of volume fm. If and , the lowest accessible density is of the order fm, which implies that no point below this limit is present. Actually this value is even higher as for numerical efficiency purposes the minimum allowed cluster number is larger than 1.

By consequence, the low density limits of the coexistence zone have to be considered as a numerical artifact. To have trustable predictions at lower densities, the simulation volume will have to be increased in future calculations. Conversely the shrinking of the coexistence zone with increasing temperature and asymmetry is a physical effect.

For a given baryonic density and temperature, we can also see that coexisting phases of clusterized matter are much more isospin symmetric with respect to the ones of the uniform matter, even if the energetics of the two systems corresponds to compatible values for the symmetry energy. This can be physically understood from the fact that higher densities are locally explored if a system is clusterized, and the symmetry energy is an increasing function of density.

Concerning the reduction of the limiting temperature in the clusterized system, this is understood as a cumulated effect of Coulomb and cluster center of mass translational degrees of freedom, which are not accounted for in the homogeneous nuclear matter calculation (63).

These general observations already show the importance of properly accounting for the cluster properties of matter, even if one is only interested in global integrated thermodynamic quantities as equations of state.

From the representation it comes out that the coexistence line of the clusterized matter is embedded in the spinodal region of the uniform matter, and shifted towards higher values of respect to the coexistence of nuclear matter, that is it corresponds to the low density part of the spinodal region of homogeneous matter. This fact, together of the absence of instability of the clusterized system, is already an indication that the homogeneous matter instability is effectively cured by accounting for the possibility of cluster formation, as we will develop in the next section. In the end, we mention that the phase diagram of clusterized matter might slightly change while modifying cluster properties.

3. The mixture properties

In the previous section we have examined the thermodynamics of purely uniform and purely clusterized nuclear systems. In the physical situation of supernovae and proto-neutron star matter it is clear that these two components have to be simultaneously present, and the question arises of how describing their mutual equilibrium conditions. In many past works (17); (19); (20); (51); (66), the evolution of clusterized matter towards homogeneous matter at high density was described as a first-order phase transition, governed by Gibbs or Maxwell equilibrium rules. The simultaneous presence in the Wigner-Seitz cell of localized clusters and homogeneously distributed free protons and neutrons in the nucleonic drip region is then also depicted as a manifestation of the same first order transition, though corrected by surface and Coulomb effects arising from the finite size of the cell.

However, as we discuss now, the proper way of treating density inhomogeneities occurring at a microscopic (femtometer scale) level is the concept of a gas mixture, if the coexisting components (here: nucleons and fragments) are non-interacting. If interaction effects arising from the excluded volume are included, we show that such a mixture naturally produces a continuous (second order) transition to uniform matter at densities close to saturation.

Restricting for the moment to the baryonic sector, the thermodynamics of a system with two conserved charges and is a two-dimensional problem. The other charges do not affect the thermodynamics at the baryonic level because they are not coupled to the baryons (see Eq. (1)). The net electron density is coupled to the proton density through the electromagnetic interaction, but this does not imply an additional degree of freedom as charge neutrality imposes . This two-dimensional problem can be easily taggled if we switch from the grancanonical ensemble, where all the intensive parameters are constrained and the densities can freely vary, to a statistical ensemble where one single extensive observable (i.e. a single density) is free to fluctuate (53). In our two-dimensional problem, two equivalent possibilities exist, namely the proton-canonical, neutron-grancanonical ensemble where the thermodynamic potential is Legendre transformed to its conjugated extensive variable ,


and the neutron-canonical, proton-grancanonical ensemble described by the thermodynamic potential,

Figure 4: (Color online) Properties of the uniform (solid lines) and net-charge neutralized clusterized (solid circles connected by lines) matter at =5 MeV and =0 MeV. Homogeneous matter is described by SKM*.
Figure 5: (Color online) The same as in Fig. 4 for =10 MeV and MeV.

We now come to illustrate the way of mixing clusters and homogeneous matter and the differences between such mixture and a possible phase coexistence which would be associated to a first-order phase transition. Fig. 4 shows the equations of state and for the proton-canonical ensemble, as well as the functional relation between the two intensive parameters and for the uniform (solid lines) and clusterized (solid circles connected by lines) matter at =5 MeV and a representative neutron chemical potential =0 MeV. The clusterized system is confined within a volume fm.

The back-bendings of and and the tilting shape of indicate that in the considered range uniform matter, if it would exist alone, would be unstable against phase separation.

If we consider that, at the microscopic scale of the Wigner-Seitz cell, clusters and free nucleons may be simultaneously present, as expected from the phenomenon of neutron drip, then the statistical equilibrium conditions read:


where and measure the volume fractions associated to the homogeneous matter and clusters,


with standing for the intrinsic volume of a nucleon calculated such as and representing the summed-up clusters intrinsic volumes.

Let us denote by and the proton chemical potentials corresponding to the limits of the spinodal region of homogeneous matter. For MeV, there will be one single low density solution of uniform matter which can constitute a mixture at thermal equilibrium with clusterized matter according to Eqs. (33) and which, for sufficiently small values of , may differ by orders of magnitude from the density of clusters. For , three different density states of uniform matter can exist in equilibrium with the single solution of clusterized matter, leading to three different mixtures corresponding to the same chemical potentials but different total density, chemical composition and pressure. It is important to stress that these three uniform matter solutions may all lead to thermodynamically stable mixtures, as the stability properties of the entropy has to be examined after adding all the different matter components. This is at variance with the case of normal nuclear matter, where only one solution out of the three corresponds to a thermodynamically stable state, while the other two are respectively metastable and unstable. For the thermodynamic condition is in principle the same as for high and negative chemical potentials. Since however the sum of the two densities corresponding to homogeneous and clusterized matter exceeds , the situation goes beyond the applicability domain of the present model and is, thus, devoid of interest.

Given that one of the purposes of our paper is to address the unhomogeneous (crust) - homogeneous (core) transition, it is possible now to anticipate that for these two representative temperatures =1.6 and 5 MeV this will take place when, following a constant proton fraction trajectory, the uniform matter state will not correspond anymore to the low density, but rather to the high density solution.

At higher temperatures the situation is slightly different as shown by Fig. 5, which shows the same quantities plotted in Fig. 4, for a representative neutron chemical potential at MeV. In this case, in the relatively narrow chemical potential region where three different homogeneous matter macrostates are found, the corresponding cluster states are always at high density close to saturation. Then, the only mixture solution not overcoming is the one given by the low density solution for uniform matter. At this temperature then, the transition to homogeneous matter does not occur as an increasing proportion of homogeneous matter that becomes abruptly dominant in the mixture as we have just seen at lower temperatures, but rather as an increasing size of the largest cluster which abruptly starts to fill the whole available volume. The two transition mechanisms are physically very close to each other: an infinitely large cluster has the same energetic and entropic properties of homogeneous matter at saturation density, provided the parameters of the effective interaction are compatible with the cluster energy parameters.

Since the thermodynamic potential of the mixture can be expressed as the sum of the thermodynamic potentials corresponding to each sub-system:


the total neutron and proton density of the mixture is obtained by summing up the corresponding particle numbers,


with the extra condition that nucleons cannot be found in the space occupied by the clusters and vice-versa.

The densities then result:


The excluded volume correction plays a sizeable role only in the vicinity of , and is in particular responsible for the un-homogeneous - homogeneous matter transition at high temperature. More precisely, pure homogeneous nuclear matter phase occurs if and, respectively, .

The equation of state and constrained entropy of the mixture for (=5 MeV, MeV) and (=10 MeV, =-20 MeV) are indicated by circles in Figs. 6 and 7. We can see that the simultaneous presence of the two components constitutes always an entropy gain respect to both purely homogeneous and purely clusterized matter, and corresponds therefore always to the equilibrium solution.

At the lowest temperature the versus equation of state is bivalued in a given density domain. This can be understood as follows. At temperatures low enough that the cluster component shows a liquid-gas transition, which can be recognized from the plateau in the relation at constant , two very different cluster-matter mixtures can be obtained from Eqs. (33) and (37) at the same total density but different chemical potentials. These two possible mixtures correspond to a homogeneous matter solution lying inside the spinodal combined respectively to a low density-low chemical potential cluster solution (stars in Fig. 6), or to a higher density-higher chemical potential solution. Since these two possibilities correspond to the same total density, the equilibrium solution will be the one maximizing the associated entropy. It comes out that the lowest chemical potential solution always correspond to a negative total pressure for the mixture, while the pressure associated to the highest chemical potential solution is always positive, meaning that the entropy is higher (see the inset in Fig. 6).

The states denoted by stars in Fig. 6 are therefore unstable or metastable and will not be further considered. The presence of such multiple solutions in this density and temperature domain may be an indication that some degrees of freedom are missing in the model and cluster deformations should be added in this region where pasta structures are expected (11); (12); (13); (14).

To summarize, we have modelised the simultaneous presence of free nucleons and light and heavy nuclei as two independent components which are allowed to exchange particles and energy following the laws of statistical equilibrium, and whose proportions are further determined by the geometrical constraint given by the excluded volume.

This means that the densities will have different values in each sub-system and that the partial pressures will be also different, following Eqs. (33) and (37).

The chemical potential equality between the two components insures the global maximization of the total entropy of the mixture. As a result, a smooth decreasing weight of the cluster component is first observed with increasing density, followed by a successive smooth but steeper decrease as density approaches . This last result can be intuitively understood from the simple physical fact that when density increases, the space accessible to a clusterized system decreases, and so does the effective number of degrees of freedom. As we will discuss in greater detail in the following, this smooth behavior of the mixture does in no way imply that the crust-core transition is smooth. Indeed because of the characteristic backbending behavior of the homogeneous matter equation of state, which reflects the well-known instability respect to spontaneous density fluctuations, the density content of the homogeneous component can discontinuously change as a function of the chemical potential.

Figure 6: (Color online) versus (bottom) and versus (top) for the homogeneous-clusterized matter coexistence (open diamonds) and mixture (open circles if and stars if ) at =5 MeV along the constant MeV path. The homogeneous matter and net-charge neutralized cluster matter curves are plotted with solid line and, respectively, small solid symbols.
Figure 7: (Color online) The same as in Fig. 6 for =10 MeV and MeV.

4. Mixture versus phase coexistence

The crust-core transition of a cold neutron star is often (17); (18); (19); (20) depicted as a first-order transition from a solid phase of nuclei on a Wigner lattice to a liquid phase of homogeneous nuclear matter. By continuity, a similar phenomenology is thought to occur at the finite temperature and more isospin symmetric conditions involved in supernova matter (23); (51).

This hypothesis arises from essentially two considerations that we now turn to examine:

  • The intermediate configuration between the solid crust and the liquid core, namely finite nuclei immersed in an homogeneous neutron background, intuitively resembles to a coexistence between a clustered and a homogeneous baryonic phase.

  • Normal nuclear matter is known to present an instability respect to cluster formation and phase separation. This spinodal instability is associated to a first-order transition. By continuity the same is expected in the baryonic part of stellar matter.

Let us first concentrate on the first statement. As we have already mentioned, the simultaneous presence of two different states characterized by a different average value for an order parameter (here: baryonic density) constitutes a phase coexistence if and only if the order parameter fluctuations (here: density inhomogeneities) occur on a length scale which is macrosocopic respect to the range of the interaction (here: much larger than the femtometer scale of typical nuclei size). This condition on macroscopic scales is necessary because the Gibbs construction neglects the interface energy and entropy, which is essential to describe thermodynamics of finite systems (69). Since this is clearly not the case in stellar matter, surface and Coulomb terms are sometimes added (23) to approximately account for finite-size corrections. However we have shown in the last section that the modelization of a microscopic mixture between clusters and homogeneous matter is also a way to describe the intermediate configurations between core and crust. Therefore we may ask if the two modelizations of phase mixture and phase coexistence are equivalent or at least lead to equivalent results for the observables of interest.

To clarify this point, let us consider homogeneous and clusterized matter as two different possible phases for stellar matter which may coexist at a first-order transition point. The advantage of working with an hybrid ensemble as the proton-canonical ensemble defined above, is that the two-dimensional tangent construction implied by Gibbs conditions in a multi-component first-order phase transition (67); (68) is reduced to a simple one-dimensional Maxwell construction on the unique free extensive variable (53). It is important to stress that this procedure is just a technical simplification and the obtained equilibrium is exactly identical to the usual two-dimensional Gibbs construction in the grancanonical ensemble. This construction in a Legendre-transformed potential should not be confused with a Maxwell construction on the equation of state performed on trajectories characterized by constant proton fractions , which is currently used in the literature (23) and has no physical justification. Indeed in a first-order phase transition in a two-component system the chemical composition of the two coexisting phases is never the same and the total pressure is a monotonically increasing function of the total density (67). This means that a Maxwell construction on at constant does not give the correct Gibbs equilibrium.

The inset in the right upper panel of Fig. 4 shows that for the representative choice of the figure the two components are found at the same value of proton chemical potential and pressure in a single point located at MeV/fm MeV. Only at this point the Gibbs equilibrium rule for coexisting phases is fulfilled,


As it can be seen in Fig. 4, at this point the density of the two components differs of more than two orders of magnitude. If clusters and homogeneous matter could be viewed as two different phases for stellar matter, the global system would be composed of homogeneous matter only for , and of clusters only for , with a discontinuous (first-order) transition at the transition point . The effect of this construction on the equation of state is shown on Fig. 6 by the open diamonds. Within the picture of a first-order phase transition, stellar matter at the chosen representative point would be homogeneous at low density up to , with fm; for intermediate densities , with fm, it would be a combination of matter and clusters varying in linear proportions; and at high density it would be solely composed of clusters, which would extend over the whole volume at the highest densities. As we can see, this equation of state is very different from the physically meaningful one corresponding to a continuous mixture of matter and clusters over the whole density domain. However, the two models lead to an almost equivalent entropy, as it is shown by the upper part of Fig. 6.

Similar considerations apply at other temperatures and chemical potentials. As a second example, Fig. 7 compares the properties of the mixture with the ones of an hypothetical first-order transition for MeV and MeV, the same thermodynamic conditions as in Fig. 5. Only within the mixture picture free protons in weak proportions are present in the matter, as it is physically expected. An inspection of the behavior of the constrained entropy at this high temperature (higher part of Fig. 7) clearly shows that the continuous mixture of nucleons and clusters not only correctly recovers the homogeneous matter properties at high density, but is also very effective in compensating the mean-field instability of homogeneous matter at lower densities, leading to a concave entropy over the whole density range (see the inset in the upper panel). Thus the usual argument that a phase coexistence would be needed to maximize the entropy and cure the convexity anomaly of the entropy cannot be applied: the same entropic gain is obtained considering a continuous mixture, without applying an artificial construction which neglects all finite-size effects.

We note by passing that at high temperature there are two points satisfying the Gibbs conditions Eq. (38) for a given value of (see Fig. 5). However at the second point at high pressure the two components correspond to almost the same density. Making a mixture or a phase coexistence does not therefore produce any sensitive difference in the equation of state (lower part of Fig. 7) in this density region.

Figure 8: (Color online) Mixture versus coexistence for =10 MeV and MeV. Up: versus ; Middle: average mass of a cluster versus ; Bottom: average cluster number per total baryon number versus .

The matter composition obtained with a continuous mixture is compared to the one corresponding to a discontinuous transition obtained with the Gibbs construction in Fig. 8. We can see, together with the different behavior of the equation of state, that the predictions for the average cluster size and cluster multiplicity show quantitative differences. It is difficult to know a priori the consequence of such altogether slight differences on the macroscopic modelization of supernova matter, and we cannot exclude that the use of a mixture equation of state would not be really distinguishable from the currently used approximations employing first-order phase transitions. Even if this would be the case, the discontinuities in all thermodynamic quantities implied by these artificial first-order transitions lead to numerical instabilities in supernova codes (70) which are difficult to handle, and which could be avoided using continuous equations of state.

To defend the possibility of a first-order crust-core transition, one may argue that the behavior of the entropy as a function of the density shown in Figs. 6 and 7 is certainly model dependent, and that in general residual convexities could still be present after considering the mixture of nucleons and clusters, especially at the lowest temperatures. Such convexities would represent residual instabilities and should then be cured with a Gibbs construction, leading again to a first-order phase transition from a low-density, cluster dominated mixture to a high density, matter dominated mixture as in the usual representation of a neutron star. Such an instability is known to be present in homogeneous nuclear matter, as shown by the back-bendings of Figs. 4 and 5, and it may well survive if clusterization was allowed in a microscopic way (50) going beyond the mean-field approximation in a microscopically consistent way.

The argument exposed above is the second argument which can be invoked to justify a discontinuous transition, which was mentioned at the beginning of this section, and that we now come to discuss.

It is certainly correct that our model depicting the simultaneous presence of nucleons and clusters as a non-interacting mixture of two different components associated to completely independent degrees of freedom and self-energies is highly schematic, and improvements are in order.

In particular the way the excluded volume is parametrized has an influence on the quantitative balance between nucleons and clusters, with a possible consequence on the precise location of the transition which is not completely under control.

Fully self-consistent treatments of the quantal many-body problem at finite temperature as in Refs. (47); (48); (49); (50), including the possibility of self-consistent clustering at low densities are extremely promising, and will in the long run provide a much better description of stellar matter than the phenomenological model proposed here. Concerning the issue of the order of the transition however, we believe that the availability of such sophisticated calculations would not modify our conclusions. Indeed, as we have already stated, the thermodynamic stability has to be checked on the global thermodynamic potential Eq. (1), after all the different constituents are added up. If gammas and neutrinos are completely decoupled from the baryonic sector, this is not true concerning electrons and positrons because of the neutrality constraint . Because of this, the properties of the electron free energy will influence the convexity properties of the baryonic entropy as a function of the density.

Figure 9: (Color online) Constrained entropy for the ensemble baryons plus electrons at =5 MeV along the constant MeV path (left panel) and at =10 MeV and MeV (right panels) for the homogeneous matter (full line) and the mixture of matter and clusters (open circles).

Fig. 9 shows the thermodynamic potential of the proton-canonical, neutron-grancanonical ensemble in the same thermodynamic conditions as in Figs. 6 and 7, including the contributions of electrons (detailed in the next section) and explicitly including the neutrality constraint


We can see that taking into account the presence of clusters has virtually no effect on the convexity of the global entropy, which is largely dominated by the convexity properties of the electron thermodynamic potential. Since this latter is highly concave, whatever the treatment of the baryonic sector, the global entropy will stay concave in the whole density domain. The effect of the electrons is thus to compensate the instability of nuclear matter and quench the possibility of a first-order phase transition.

This result can be physically understood (16); (71): because of the low electron mass, the electron gas is much more incompressible than the nuclear component. A fluctuation in the baryonic density creates a fluctuation in the charge density, which because of the neutrality constraint implies a fluctuation in the net electron density. Such fluctuations are naturally present in the length scale of the Wigner-Seitz cell, and are at the origin of the electrostatic interaction energy which is non-zero even if the system is neutral, see Eq. (25). However these fluctuations cannot extend over macroscopic lengths because they would imply macroscopic electron inhomogeneities, which are contrasted by the high electron incompressibility.

To conclude, phases of different global baryon density cannot exist in stellar matter if baryonic matter is constituted of heavy particles as nucleons and clusters, as it is the case at the densities of interest in the present paper. Since this statement is valid for all chemical potentials and temperatures, including , we conclude that the crust-core phase transition cannot be of first-order at any temperature.

ii.2 B. Leptons and photons

In the temperature and density domains relevant for this paper, MeV and g cm, leptons (electrons and neutrinos) are relativistic, in particle-antiparticle pair equilibrium and in thermal equilibrium with nuclear matter (17); (23).

Given that the net charge of the electron gas (that is the number of electrons minus the number of positrons) has to neutralize the positive charge of protons in uniform and clusterized matter, the electron chemical potential is determined by the equality among the electron density, calculated in the relativistic Fermi gas model, and the proton density,


where represents the baryon density, is the spin degeneracy and is the rest mass.

The electron gas pressure, entropy and energy densities are functions of and and their expressions are the following,


and, respectively,


Neutrinos also form a relativistic Fermi gas in pair equilibrium and thermal equilibrium with nuclear matter. This means that Eqs. (40, 41, 42 and 43) with and instead of and describe their density, pressure, entropy and energy densities. Under the -equilibrium hypothesis, neutrinos chemical potential is dictated by the chemical potential of neutrons, protons and electrons,


If not explicitely mentioned otherwise (section III.C), all over this paper we shall work out of -equilibrium. Neutrinos contribution to the total energy, entropy and pressure will be disregarded, as well. The motivation of this choice is that equilibration with respect to weak interaction is often not achieved over the time scales of many astrophysical phenomena, but we will discuss this point further in section III.C.

Photons are assumed to be in thermal equilibrium with the other star ingredients and their pressure, entropy and energy densities are,


and, respectively,


Iii III. Star matter relevant results

iii.1 A. Constant chemical potential paths

In the evolution of supernovae explosions, nuclear composition plays an extremely important role. First of all, isotopic abundances in the light and medium mass region influence nucleosynthesis via r-, rp- and photodisintegration processes. Then, nuclei with masses larger than 40-60 determine the rate of electron capture (33); (34), the most important weak nuclear interaction in the dynamics of stellar core collapse and bounce, responsible for the neutron enrichment of baryonic matter and neutrino emission. On the other hand, nuclei absorb energy via internal excitation making the increase in temperature and pressure during star collapse be less strong than in the case of a uniform nucleon gas (27). Nuclei influence also the maximum density reached during the collapse. In the absence of nuclei the collapse stops before reaching , while the limiting value is larger than if nuclei are present (28). Finally, by dissociation under the shock wave, nuclear composition affects the shock propagation (30).

Under this perspective, it is obvious that realistic star matter equations of state require accurate treatment of the nuclear statistical equilibrium. Working within the single nucleus approximation (SNA), the first generation models (17); (23) do not offer a satisfactory description of nuclear miscellanea and may be to some extent responsible for the failure of present supernovae simulations to produce explosions, even if realistic neutrino transport(72) and convective instabilities in multi-dimensional