Statistical description of complex nuclear phases in supernovae and protoneutron stars
Abstract
We develop a phenomenological statistical model for dilute star matter at finite temperature, in which free nucleons are treated within a meanfield 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 protoneutron stars. The first finding is that, contrary to the common belief, the crustcore 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 LattimerSwesty, 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 crustcore 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.
pacs:
21.65.Mn, 24.10.Pa, 26.50.+x, 26.60.cMarch 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 antineutrinos 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), wellknown 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 nonspherical 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 LattimerSwesty (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 shockwave 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 noninteracting 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 NSEbased approaches is that they completely neglect inmedium effects, which are known to be very important in nuclear matter. Since the only nuclear interactions are given by the cluster selfenergies, 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 meanfield models (46). As a consequence, these models cannot be applied at densities close to saturation and the crustcore 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 HartreeFock 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 crustcore 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 firstorder 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 addedup.
In the grancanonical ensemble this reads,
(1) 
where the grancanonical potential,
(2)  
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
(3)  
the different particle densities are
(4) 
where , and finally the total pressure is
(5) 
ii.1 A.The baryon sector
The light and heavy nuclei are assumed to form a gas of looselyinteracting clusters which coexist in the WignerSeitz cell with a homogeneous background of delocalized nucleons. To avoid exceeding the normal nuclear density and naturally allow for homogeneousunhomogeneous 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 WignerSeitz cell.
1. The homogeneous nuclear matter component
Meanfield models constitute a natural choice for approaching interacting particle systems. By introducing a meanfield potential, the physical problem is reduced to the simplified version of a system of noninteracting particles. Effective nucleonnucleon interactions allow one to express the system average energy as a simple singleparticle 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 independentparticle system associated to the meanfield singleparticle energies, with the temperature weighted difference between the average singleparticle energy () and the meanfield energy
(6) 
The one body partition sum is defined as:
(7) 
and can be expressed as a function of the neutron and proton kinetic energy density :
(8) 
with
(9)  
(10) 
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 selfconsistent relation between the density of qparticles 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 FermiDirac distribution indeed reads:
(11) 
Eqs. (11) and (9) define a selfconsistent problem since depends on the densities. For each couple a unique solution is found by iteratively solving the selfconsistency 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:
(12)  
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 liquidgas like firstorder 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 nucleonnucleon 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 
0.09  
0.  
0.  
0.  
1/6 
Observable  Value 

(MeV)  15.78 
(MeV)  216.7 
0.79  
(MeV)  30.03 
(fm)  0.1603 
(MeV)  17.51 
3.74 
The thermodynamic properties of the system are best studied introducing the constrained entropy:
(13) 
where the entropy density is nothing but the Legendre transform of the mean field partition sum:
(14) 
Instabilities are then recognized as local convexities of the constrained entropy at a given temperature in the twodimensional 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 wellknown Gibbs equality between all intensive parameters of two coexisting phases in equilibrium:
(15) 
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 socalled 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.
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 firstorder liquidgas 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 subsystem 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:
(16)  
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 intercluster 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 firstorder approximation to the nuclear part of the interaction energy in the spirit of the Van der Waals gas. More sophisticated inmedium corrections can in principle be also added as modifications of the cluster selfenergies (50). Because of this, we will consider that the interaction energy only contains the Coulomb electrostatic energy of the configuration.
To obtain the socalled NSE approach, one has to neglect all interparticle interactions except the cluster selfenergies, leading to :
(17) 
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
(18) 
where gives the multiplicity of clusters of type . Equation (17) becomes:
(19) 
with
(20) 
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
(21) 
where is the standard grancanonical partition sum for a single nucleus of mass and charge :
(22)  
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 welltested MetropolisMonteCarlo technique (62), and calculate relevant thermodynamic observables other than the ones fixed by the exterior as ensemble averages. For the configurationdefined quantities, =, , :
the corresponding ensemble averaged density reads:
(24) 
and stability is systematically checked against an increase of the considered volume .
For each sampled configuration , we evaluate the Coulomb interaction energy in the WignerSeitz approximation (17),
(25) 
with
(26) 
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:
(27)  
Here originates from clusters translational motion and is the socalled 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 recasted as a global configurationdependent factor as,
(28) 
For the clusters we assume the following properties:

clusters binding energies are described by a liquiddrop parameterization,
(29) 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 straightforwardly 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 Fermigas type with cutoff correction (65),
(30) 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 firstorder liquidgas 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 supracritical 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 protoneutron 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 firstorder phase transition, governed by Gibbs or Maxwell equilibrium rules. The simultaneous presence in the WignerSeitz 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 noninteracting. 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 twodimensional 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 twodimensional 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 twodimensional problem, two equivalent possibilities exist, namely the protoncanonical, neutrongrancanonical ensemble where the thermodynamic potential is Legendre transformed to its conjugated extensive variable ,
(31)  
and the neutroncanonical, protongrancanonical ensemble described by the thermodynamic potential,
(32)  
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 firstorder phase transition. Fig. 4 shows the equations of state and for the protoncanonical 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 backbendings 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 WignerSeitz cell, clusters and free nucleons may be simultaneously present, as expected from the phenomenon of neutron drip, then the statistical equilibrium conditions read:
(33) 
where and measure the volume fractions associated to the homogeneous matter and clusters,
(34) 
with standing for the intrinsic volume of a nucleon calculated such as and representing the summedup 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 subsystem:
(35) 
the total neutron and proton density of the mixture is obtained by summing up the corresponding particle numbers,
(36)  
with the extra condition that nucleons cannot be found in the space occupied by the clusters and viceversa.
The densities then result:
(37)  
The excluded volume correction plays a sizeable role only in the vicinity of , and is in particular responsible for the unhomogeneous  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 liquidgas transition, which can be recognized from the plateau in the relation at constant , two very different clustermatter 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 densitylow chemical potential cluster solution (stars in Fig. 6), or to a higher densityhigher 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 subsystem 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 crustcore transition is smooth. Indeed because of the characteristic backbending behavior of the homogeneous matter equation of state, which reflects the wellknown instability respect to spontaneous density fluctuations, the density content of the homogeneous component can discontinuously change as a function of the chemical potential.
4. Mixture versus phase coexistence
The crustcore transition of a cold neutron star is often (17); (18); (19); (20) depicted as a firstorder 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 firstorder 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 finitesize 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 firstorder transition point. The advantage of working with an hybrid ensemble as the protoncanonical ensemble defined above, is that the twodimensional tangent construction implied by Gibbs conditions in a multicomponent firstorder phase transition (67); (68) is reduced to a simple onedimensional 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 twodimensional Gibbs construction in the grancanonical ensemble. This construction in a Legendretransformed 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 firstorder phase transition in a twocomponent 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,
(38) 
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 (firstorder) 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 firstorder 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 firstorder 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 meanfield 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 finitesize 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.
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 firstorder phase transitions. Even if this would be the case, the discontinuities in all thermodynamic quantities implied by these artificial firstorder 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 firstorder crustcore 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 firstorder phase transition from a lowdensity, 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 backbendings of Figs. 4 and 5, and it may well survive if clusterization was allowed in a microscopic way (50) going beyond the meanfield 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 noninteracting mixture of two different components associated to completely independent degrees of freedom and selfenergies 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 selfconsistent treatments of the quantal manybody problem at finite temperature as in Refs. (47); (48); (49); (50), including the possibility of selfconsistent 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.
Fig. 9 shows the thermodynamic potential of the protoncanonical, neutrongrancanonical 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
(39) 
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 firstorder 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 WignerSeitz cell, and are at the origin of the electrostatic interaction energy which is nonzero 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 crustcore phase transition cannot be of firstorder 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 particleantiparticle 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,
(40) 
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,
(41) 
(42) 
and, respectively,
(43) 
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,
(44) 
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,
(45) 
(46) 
and, respectively,
(47) 
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 4060 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 multidimensional