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

prev

next

out of 28

Published on

13-Apr-2017View

213Download

0

Transcript

PHYSICAL REVIEW C 82, 065801 (2010)Statistical description of complex nuclear phases in supernovae and proto-neutron starsAd. R. Raduta1 and F. Gulminelli21NIPNE, Bucharest-Magurele, POB-MG6, Romania2LPC (IN2P3-CNRS/Ensicaen et Universite), F-14076 Caen cedex, France(Received 10 September 2010; revised manuscript received 22 October 2010; published 9 December 2010)We develop a phenomenological statistical model for dilute star matter at finite temperature, in which freenucleons are treated within a mean-field approximation and nuclei are considered to form a loosely interactingcluster gas. Its domain of applicability, that is, baryonic densities ranging from about > 108 g/cm3 to normalnuclear density, temperatures between 1 and 20 MeV, and proton fractions between 0.5 and 0, makes it suitablefor the description of baryonic matter produced in supernovae explosions and proto-neutron stars. The firstfinding is that, contrary to the common belief, the crust-core transition is not first order, and for all subsaturationdensities matter can be viewed as a continuous fluid mixture between free nucleons and massive nuclei. Asa consequence, the equations of state and the associated observables do not present any discontinuity over thewhole thermodynamic range. We further investigate the nuclear matter composition over a wide range of densitiesand temperatures. At high density and temperature our model accounts for a much larger mass fraction boundin medium nuclei with respect to traditional approaches as Lattimer-Swesty, with sizable consequences on thethermodynamic quantities. The equations of state agree well with the presently used EOS only at low temperaturesand 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 showthat the two problems can be effectively decoupled.DOI: 10.1103/PhysRevC.82.065801 PACS number(s): 21.65.Mn, 24.10.Pa, 26.50.+x, 26.60.cI. INTRODUCTIONNuclear matter is not only a theoretical idealization provid-ing benchmark studies of the effective nuclear interaction, butit is also believed to constitute the major baryonic componentof massive objects in the universe, as exploding supernovaecores and neutron stars. The structure and properties ofthese astrophysical objects at baryonic densities exceedingnormal nuclear matter density is still highly speculative [18].Conversely, for subsaturation densities, it is well establishedthat matter is mainly composed of neutrons, protons, electrons,positrons, and photons in thermal and typically also chemicalequilibrium [9,10]. Depending on the thermodynamic condi-tion, neutrinos, and antineutrinos can also participate to theequilibrium.Such composite matter is subject to the contrasting cou-plings of the electromagnetic and the strong interaction.Because of the electron screening, the two couplings act oncomparable length scales giving rise to the phenomenon offrustration [1114], well-known in condensed matter physics[15]. Because of this, a specific phase diagram, differingfrom the one of nuclear matter and including inhomogeneouscomponents, is expected in stellar matter [16].Many theoretical studies exist at zero temperature. In a coldneutron star, going from the dilute crust to the dense core, atransition is known to occur from a solid phase constituted offinite nuclei on a Wigner lattice immersed in a backgroundof delocalized electrons and neutrons through intermediateinhomogeneous phases composed of nonspherical nuclei(pasta phases) to a liquid phase composed of uniform neutrons,protons, and electrons [1722].At finite temperature the matter structure and propertiesare not as well settled. The most popular phenomenologicalapproaches are the Lattimer-Swesty [23] (LS) and the Shen[24] equation of state, recently updated in Ref. [25]. In thesestandard treatments currently used in most supernovae codes,the dilute stellar matter at finite temperature is described inthe baryonic sector as a statistical equilibrium among protons,neutrons, s, and a single heavy nucleus. The transition tohomogeneous matter in the neutron star core is supposed tobe first order in these modelizations and obtained througha Maxwell construction in the total density at fixed protonfraction.It is clear that such single nucleus approximation (SNA) ishighly schematic and improvements are possible. Concerningintegrated quantities as thermodynamic functions and equa-tions of state, such variables may be largely insensitive tothe detailed matter composition [26], though we show in thisarticle that this is not always the case. However, it is alsoknown that the composition at relatively high density close tosaturation, together with the pressure and symmetry energy,governs the electron capture rate, which in turn determinesthe proton fraction at bounce and the size of the homologouscore, a key quantity to fix the strength of the shock-wave andthe output of the supernovae explosion [2732]. Moreoverthe composition may also affect the nucleosynthesis of heavyelements, which is still poorly understood [3336], as well asthe neutrino scattering through the core after bounce [37,38]and the cooling rate of neutron stars [39,40]. For these reasons,in recent years, many efforts have been done to improve thesimplistic representation of stellar matter given by the SNAapproach.The different modelizations which consider a possibledistribution of all different nuclear species inside dilute stellarmatter are known under the generic name of nuclear statistical0556-2813/2010/82(6)/065801(28) 065801-1 2010 The American Physical Societyhttp://dx.doi.org/10.1103/PhysRevC.82.065801AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)equilibrium (NSE) [4144]. The basic idea behind thesemodels is the Fisher conjecture that strong interactions indilute matter may be entirely exhausted by clusterization[45]. In these approaches stellar matter in the baryonicsector is then viewed as a noninteracting ideal gas of allpossible nuclear species in thermal equilibrium. The resultis that thermodynamic quantities like entropies and pressureappear very similar to the ones calculated with standardapproaches, while noticeable differences are seen in the mattercomposition. In particular, an important contribution of lightand intermediate mass fragments is seen at high temperature,which is neglected in standard SNA approaches.The strongest limitation of NSE-based approaches is thatthey completely neglect in-medium effects, which are knownto be very important in nuclear matter. Since the onlynuclear interactions are given by the cluster self-energies, thehomogeneous matter composing the neutron star core cannotbe modelized, nor can it be the phenomenon of neutron drip inthe inner crust, well described by mean-field models [46]. As aconsequence, these models cannot be applied at densities closeto saturation and the crust-core transition cannot be described.To overcome this problem, different microscopic [4750]as well as phenomenological [51] approaches have beendeveloped in the very recent past. In this article we wouldlike to introduce a phenomenological model that treats thenuclei component within an improved NSE, while it describesthe unbound protons and neutrons in the finite temperatureHartree-Fock approximation.The plan of the article is as follows. The first part ofthe article is devoted to the description of the model. Theclusterized component, the homogeneous component, theproperties of the mixture, and the lepton sector are describedin successive sections together with their thermodynamicproperties. Particular attention is devoted to the modelizationof the crust-core transition. We show that the inclusion ofexcluded volume is sufficient to describe the transition fromthe clusterized crust to the homogeneous core and that thistransition is continuous. Different generic as well as specificarguments are given against the possibility of a first-ordertransition. The second part of the article gives some resultsrelevant for the star matter phenomenology. The first sectionshows observables following constant chemical potentialpaths, in order to connect the observables with the properties ofthe phase diagram. Then the behavior of the different quantitiesfor constant proton fractions is displayed in order to comparewith more standard treatments of supernova matter. Finally,the last section addresses the problem of neutrino trappingand the interplay between the matter opacity to neutrinos andmatter composition. Conclusions and outlooks conclude thearticle.II. THE MODELThe model aims to describe the thermal and chemicalproperties of nuclear matter present in supernovae and(proto-)neutron stars at densities ranging from the normalnuclear density 0 to 1060, temperatures between 0 and20 MeV, and proton concentration between 0.5 and 0. Inthis regime, the star matter typically consists of a mixtureof nucleons, light and heavy nuclear clusters, neutrinos (ifwe consider the thermodynamic stage where neutrinos aretrapped), photons, and a charge-neutralizing background ofelectrons and positrons.As there is no interaction among electrons, neutrinos, pho-tons, and nuclear matter, the different systems may be treatedseparately and their contributions to the global thermodynamicpotential and equations of state added up.In the grand-canonical ensemble this reads,G(, n, p, e, V )= G(bar)(, n, p, V ) +G(lep)(, e, V )+G( )(, V ), (1)where the grand-canonical potential,G(, n, p, e, V ) = lnZgc(, n, p, e, V )= S[,n,p], (2)is the Legendre transformation of the entropy S with respect tothe fixed intensive variables. In the previous equations V is anarbitrary macroscopic volume and Zgc is the grand-canonicalpartition sum.The observables conjugated to the ones fixed by thereservoir and geometry can be immediately calculated aspartial derivatives of G. Thus, the total energy density ise = 1V(G)n,p,V= e(bar) + e(lep) + e( ), (3)the different particle densities arei = 1V[G(i)],j ,V, (4)where i = n, p, e, and, finally, the total pressure isp = GV= p(bar) + p(lep) + p( ). (5)A. The baryon sectorThe light and heavy nuclei are assumed to form a gas ofloosely interacting clusters which coexist in the Wigner-Seitzcell with a homogeneous background of delocalized nucleons.To avoid exceeding the normal nuclear density and naturallyallow for homogeneous-unhomogeneous matter transition,nuclei and nucleons are forbidden to occupy the same volume.In the following we start describing the modelization of thesetwo components separately, and we turn successively to theproperties of the mixture obtained when the two are supposedto be simultaneously present in the Wigner-Seitz cell.1. The homogeneous nuclear matter componentMean-field models constitute a natural choice for approach-ing interacting particle systems. By introducing a mean-fieldpotential, the physical problem is reduced to the simplifiedversion of a system of noninteracting particles. Effectivenucleon-nucleon interactions allow one to express the system065801-2STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)average energy as a simple single-particle density functionaland to cast the nuclear matter statistics in a way which isformally very similar to an ideal Fermi gas [52].The mean-field energy density of an infinite homogeneoussystem e(HM) = H 0/V is a functional of the particle densitiesq and kinetic densities q for neutrons (q = n) or protons(q = p). At finite temperature, the mean-field approximationconsists in expressing the grand-canonical partition functionof the interacting particle system as the sum of the grand-canonical partition function of the corresponding independent-particle system associated to the mean-field single-particleenergies, with the temperature weighted difference betweenthe average single-particle energy (W 0 = lnZ0) and themean-field energylnZ (HM) lnZ (HM)0 + (W 0 H 0). (6)The one-body partition sum is defined asZ (HM)0 = Tr[e(W0nNnpNp)] = Zn0Zp0 (7)and can be expressed as a function of the neutron and protonkinetic energy densitylnZq0V= 2 0ln[1 + e(p22mq q)]4p2h3dp = h23mqq,(8)withq = 2 0nq(p)4p2h3dp, (9)q = 2 0p2h2nq(p)4p2h3dp. (10)In these equations the effective chemical potential qincludes the self-energies according to q = q q e(HM),mn,p = (eHM/n,p)1/2 are the neutron (proton) effectivemass, and the factor 2 comes from the spin degeneracy.Equation (9) establishes a self-consistent relation betweenthe density of q particles q and their chemical potential q .Introducing the single particles energies iq = p2i2mq+ q e(HM),the above densities can be written as regular Fermi integralsby shifting the chemical potential according to q = q q eHM. The Fermi-Dirac distribution indeed reads:nq(p) = 11 + exp[(p2/2mq q)]. (11)Equations (11) and (9) define a self-consistent problemsince mq depends on the densities. For each couple (n, p)a unique solution (n, p) is found by iteratively solving theself-consistency between n,p and mn,p. Then Eq. (10) is usedto calculate n,p.At the thermodynamic limit the system volume V divergestogether with the particle numbers Nn, Np, and thethermodynamics is completely defined as a function of the twoparticle densities (n, p) or, equivalently, the two chemicalpotentials (n,p).With Skyrme-based interactions, the energy density ofhomogeneous, spin-saturated matter with no Coulomb effectsTABLE I. SKM* force parameters [58].Parameter Valuet0 (MeV fm3) 2645t1 (MeV fm5) 410t2 (MeV fm5) 135t3 (MeV fm3+3 ) 15 595x0 0.09x1 0.0x2 0.0x3 0.0 1/6is written as:e(HM) = h22m(n + p)+ t0(x0 + 2)(n + p)2/4 t0(2x0 + 1)(2n + 2p)/4+ t3(x3 + 2)(n + p)+2/24 t3(2x3 + 1) (n + p)(2n + 2p)/24 + [t1(x1 + 2)+ t2(x2 + 2)](n + p)(n + p)/8 + [t2(2x2 + 1) t1(2x1 + 1)](nn + pp)/8, (12)where t0, t1, t2, t3, x0, x1, x2, x3, and are Skyrme parameters.Several Skyrme potentials have been developed over theyears for describing the properties of both infinite nuclearmatter and atomic nuclei and address the associated thermo-dynamics. It was thus in particular shown that nuclear mattermanifests liquid-gas like first-order phase transitions up to acritical temperature [5357].In order to make direct quantitative comparisons betweenour model and the one of Lattimer and Swesty [23], throughthis article nucleon-nucleon interactions are accounted foraccording to the SKM* parametrization [58] if not explicitelymentioned otherwise. Table I summarizes the force parameterswhile the main properties of nuclear matter are summarized inTable II.The thermodynamic properties of the system are beststudied introducing the constrained entropy:s(HM)c = s(HM) eHM, (13)where the entropy density is nothing but the Legendretransform of the mean-field partition sum:s(HM) = lnZ (HM)/V + [e(HM) nn pp]. (14)TABLE II. Infinite nuclear matter and surfaceproperties of SKM* force [58].Observable ValueE/A (MeV) 15.78K (MeV) 216.7m/m 0.79Csym (MeV) 30.030 (fm3) 0.1603as (MeV) 17.51ks 3.74065801-3AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)-70-60-50-40-30-20-1001020-70 -60 -50 -40 -30 -20 -10 0 10 20HM coex.HM spinodalCM coex.n (MeV) p (MeV)10-1410-1210-1010-810-610-410-210-1410-1210-1010-810-610-410-2n (fm-3) p (fm-3)FIG. 1. (Color online) p n (left) and p n (right) representations of the phase diagrams of the charge-neutral homogeneousinfinite nuclear matter described by the SKM* interaction (lines) and net-charge-neutralized clusterized nuclear matter (stars) correspondingto T = 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 thecoexistence region. These results are independent of the chosen confining volume V = 2.9 104 fm3 except in the low-density region of thephase diagram of the clusterized matter system in p n coordinates, which is affected by finite size effects (see text).Instabilities are then recognized as local convexities ofthe constrained entropy at a given temperature T = 1/ inthe two-dimensional density plane (n, p). In the presenceof an instability, the system entropy can be maximized byphase mixing, which corresponds to a linear interpolation inthe density plane. This geometrical construction correspondsto the well-known Gibbs equality between all intensiveparameters of two coexisting phases in equilibrium:(1) = (2);(1)q = (2)q , q = n, p; (15)P (1) = P (2).Phase diagrams of SKM* nuclear matter are illustrated inFigs. 1, 2, and 3 for three values of temperature T = 1.6,5, and, respectively, 10 MeV in n p (left panels) andn p (right panels) representations. As it is customary innuclear matter calculations, the system is charge neutral in thesense that the electromagnetic interaction has been artificiallyswitched off. In the star matter case the proton electric chargehas to be taken into account, and compensated by the electroncharge to give a net zero charge. The resulting Coulombinteraction energy will be discussed in the next chapter. Fulllines correspond to the phase coexistence region and dashedlines mark the borders of the instability region, the so-calledspinodal surface.As one may notice, a large portion of the density andchemical potential plane is characterized by the instabilityeven at temperatures as high as 10 MeV. Inside the phasecoexistence in the density representation, the mean-field solu-tion is unphysical at thermal equilibrium if matter is uniquelycomposed of homogeneously distributed (uncharged) protonsand neutrons. In the case of stellar matter, the entropy convexityproperties have to be examined only after considering all thedifferent constituents, which can very drastically change theproperties of the phase diagram. As we will show in the nextsections, the introduction of electrons has the effect of makingthe Gibbs construction unphysical, while the introduction ofclusterized matter has the effect of stabilizing the unstablemean-field solution.2. The clusterized nuclear matter componentAs we have discussed in the previous section, in a wideregion of temperature and density, uncharged uniform nuclearmatter is unstable with respect to density fluctuations. Ifsuch fluctuations occur on a macroscopic scale, a first-order-70-60-50-40-30-20-10010-70 -60 -50 -40 -30 -20 -10 0 10n (MeV) p (MeV)10-510-410-310-210-110-510-410-310-210-1n (fm-3) p (fm-3)FIG. 2. (Color online) The same as in Fig. 1 for T = 5 MeV. For symbol and line codes, see Fig. 1.065801-4STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)-70-60-50-40-30-20-10010-70 -60 -50 -40 -30 -20 -10 0 10HM coex.HM spinodaln (MeV) p (MeV)10-310-210-110-310-210-1n (fm-3) p (fm-3)FIG. 3. (Color online) The same as in Fig. 1 for T = 10 MeV.liquid-gas phase transition follows. Inside the spinodal surface,however, finite wavelength fluctuations may also occur, lead-ing to spontaneous cluster formation. Spinodal decompositionand nucleation are indeed known to constitute the dynamicalmechanisms leading to phase separation in macroscopicuncharged systems. As we will show in the following, inthe presence of the repulsive long-range Coulomb interactionamong protons, which is screened by the leptons only onmacroscopic scales, phase separation is quenched in stellarmatter, and the instability is cured by the formation of nuclearclusters, that is, clusters in the femtometer scale.One possibility to account for the thermal and phaseproperties of a clusterized matter subsystem is to use astatistical model with cluster degrees of freedom. Several suchmodels have been so far proposed for the study of condensationclose to the critical point [45], nuclear multifragmentation[5962], and compact stars [4244]. They basically consist inthe estimation of the number of microscopic states compatiblewith the thermodynamic macroscopic constraints.Considering that the center of mass of the clusters can betreated as classical degrees of freedom, the grand-canonicalpartition function of the clustered system reads:Z (cl),n,p=CWC =C1NC!d3r1 d3rNCd3p1 d3pNCh3NC exp[(EC nNC pZC)]=C1NC!exp[Eint(C)] NCi=1[V freei (C)(mAiT2h2)3/2(Ai, Zi, i) exp(Bi i + nNi + pZi)]. (16)Here C denotes a generic cluster configuration C ={(A1, Z1, 1), . . . , (ANC , ZNC , NC )} defined by the mass num-ber Ai , the proton number Zi , and the excitation energyi of each constituent nuclear state indexed by i, whosebinding energy is Bi and level density is (Ai, Zi, i); Eint(C)represents the intercluster configuration-dependent interactionenergy, and V freei (C) < V is the volume accessible to clusteri which is not excluded by the presence of the other clusters.Because of the short range of the nuclear force, this excludedvolume correction can be viewed as a first-order approximationto the nuclear part of the interaction energy in the spirit of theVan der Waals gas. More sophisticated in-medium correctionscan in principle be also added as modifications of the clusterself-energies Bi [50]. Because of this, we will consider thatthe interaction energy only contains the Coulomb electrostaticenergy of the configuration.To obtain the so-called NSE approach, one has to neglectall interparticle interactions except the cluster self-energies,leading toZ (cl),n,p Z(cl)0 =CV NCNC!NCi=1[(mAiT2h2)3/2(Ai, Zi, i) exp(Bi i + nNi + pZi)]. (17)The advantage of this simplification is that a completelyfactorized expression for the partition sum can be obtained,leading to analytical expressions for all thermodynamic quan-tities. Indeed any configuration (C) can be ordered as(C) = {T1, T1, . . . , T1 n1(C), T2, T2, . . . , T2 n2(C), . . .}, (18)where ni(C) = nAi,Zi ,i (C) gives the multiplicity of clusters oftype Ti = {Ai, Zi, i}. Equation (17) becomesZ (cl)0 =i=1ni=01ni![z,n,p (Ai, Zi, i)]ni , (19)withz,n,p (Ai,Zi,i) =(mAiT2h2)3/2(Ai, Zi, i) exp(Bi i +nNi +pZi).(20)Equation (19) is not yet an analytically tractable expression,because of the infinite number of possible excited states foreach nuclear isotope (Ai, Zi). The internal degrees of freedom065801-5AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)can be, however, integratedZ (NSE),n,p =A,Z=1nA,Z=01nA,Z![z,n,p (A,Z)]nA,Z , (21)where z,n,p (A,Z) is the standard grand-canonical parti-tion sum for a single nucleus of mass A and charge Z:z,n,p (A,Z) =(mAT2h2)3/2 d(Ai, Zi, ) exp(Bi + nN + pZ)=(mAT2h2)3/2g(A,Z) exp(Bi + nN + pZ), (22)and g(A,Z) is a temperature-dependent degeneracy.While elegant and analytically tractable, Eq. (21) representsa correct estimation of Eq. (16) if and only if the interactionsare fully exhausted by clusterization. For this reason, inthis work we will consider only Eq. (16), which we willestimate numerically with a precise, efficient, and well-testedMetropolisMonte Carlo technique [62], and calculate relevantthermodynamic observables other than the ones fixed by theexterior as ensemble averages. For the configuration-definedquantities, XC = Atot(C), Ztot(C), Etot(C):Atot(C) =NCi=1AiZtot(C) =NCi=1Zi (23)Etot(C) =NCi=1(Bi + i),the corresponding ensemble averaged density reads:x = 1VC XCWCC WC, (24)and stability is systematically checked against an increase ofthe considered volume V .For each sampled configuration (C), we evaluate theCoulomb interaction energy in the Wigner-Seitz approxima-tion [17],VCoulomb(C) =NCi=135c()e2Z2ir0A1/3i, (25)withc() = 1 32(e0p)1/3+ 12(e0p), (26)accounting for the screening effect of electrons. 0p = Z/A0denotes the proton density inside the nuclei and e isthe electron density. Net charge neutrality of the systemimposes the electron density e to be equal with the protondensity p.The total pressure results:p =(G(cl)V)|,n,p= NCV 352V(1n1/3 1n) ie2Z2ir0A1/3i= [p(cl) + p(lattice)]. (27)Here p(cl) originates from clusters translational motion andp(lattice) is the so-called lattice Coulomb pressure.Concerning the excluded volume, the result of the NCintegrals over the coordinate space corrected for the volumeexcluded by each cluster, Vi , can be recasted as a globalconfiguration-dependent factor (C)V NC as,V (V V1)(V V1 V2) (V V1 V2 VNC1) = V NC . (28)For the clusters we assume the following properties:(i) cluster binding energies are described by a liquid-dropparameterization,B(A,Z) = (avA as(1 Tf (T ))A2/3) (1 aI (A 2Z)2/A2) (29)with av = 15.4941 MeV, as = 17.9439 MeV, and aI =1.7826 [63], where the surface term is additionallymade to vanish at the critical temperatureTC = 12 MeV[64]; the Coulomb energy contribution is accountedfor in VCoulomb(C) as discussed above; structure effectsthoroughly considered by other authors [51] are dis-regarded because of the relatively high temperatureswe are focusing on, but they can be straightforwardlyincluded;(ii) in order to avoid multiple counting of the free nucleons,already considered in the uniform background, theclusters should have Z 2; this arbitrary limit maybe replaced by A > 1, such as to allow also for 2H and3H, but this choice is not expected to be important forthe presently considered quantities;(iii) no limit is imposed for the largest cluster which, inprinciple, may reach the total system size; given thefact that the fragment mass partition enters in Eq. (16)not only in the translational energy factor and bindingenergy but also in the calculations of excluded volumeand Coulomb energy, we expect it to play a dramaticrole in the limit of low temperatures and densities closeto 0;(iv) to allow creation of exotic species, there is no restrictionfor the neutron/proton composition; a cluster is allowedto exist as soon as B(A,Z) > 0;(v) cluster internal excitation energy is upper limited bythe cluster binding energy and the corresponding leveldensity is of Fermi-gas type with cut-off correction [65],() =12a1/45/4exp(2a) exp(/ ), (30)with a = 0.114A+ 0.098A2/3 MeV1 and = 9 MeV.065801-6STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)The phase diagram of clusterized stellar matter can beinferred from the bimodal behavior of the total particle numberin the grand-canonical ensemble [64]. This model presents afirst-order liquid-gas phase transition similar to the case ofuniform nuclear matter, but, contrary to the latter, it is neverunstable. The two phase diagrams are compared for T = 1.6and 5 MeV in Figs. 1 and, respectively, Fig. 2. Figure 3presents exclusively the phase diagram of neutral uniformmatter at T = 10 MeV, as this temperature is supracriticalfor the clusterized matter.We can see that the thermodynamics of a clusterized systemdiffers substantially from the one of neutral nuclear matterin the mean-field approximation, even though the employedeffective interaction parameters are typically tuned on the sameexperimentally accessible observables of cold nuclei close tosaturation density, which are implemented in the cluster energyfunctional. As a consequence, the two phase diagrams arecompatible only at low temperatures and at densities close tosaturation.The extension of the coexistence region in density andtemperature is noticeably reduced for clusterized matter. Thiseffect is partially not physical. Indeed the numerical procedureadopted for estimating Eq. (16) implies that clusterized matterresults depend on the system size for the lowest densities:the lowest accessible density is determined by the chosenfinite volume and the minimum allowed cluster size (Amin)eventually multiplied by the minimum allowed number ofclusters (Nmin). Clusterized matter phase diagrams consideredin this section have been obtained working within a cell ofvolume V = 2.9 104 fm3. If Nmin = 1 and Amin = 3, thelowest accessible density min = NminAmin/V is of the order1.1 104 fm3, which implies that no point below this limitis present. Actually this value is even higher as for numericalefficiency purposes the minimum allowed cluster number islarger than 1.By consequence, the low density limits of the coexistencezone have to be considered as a numerical artifact. To havetrustable predictions at lower densities, the simulation volumewill have to be increased in future calculations. Conversely,the shrinking of the coexistence zone with increasing temper-ature and asymmetry is a physical effect.For a given baryonic density and temperature, we can alsosee that coexisting phases of clusterized matter are muchmore isospin symmetric with respect to the ones of theuniform matter, even if the energetics of the two systemscorresponds to compatible values for the symmetry energy.This can be physically understood from the fact that higherdensities are locally explored if a system is clusterized,and the symmetry energy is an increasing function ofdensity.Concerning the reduction of the limiting temperature in theclusterized system, this is understood as a cumulated effectof Coulomb and cluster center of mass translational degreesof freedom, which are not accounted for in the homogeneousnuclear matter calculation [64].These general observations already show the importance ofproperly accounting for the cluster properties of matter, evenif one is only interested in global integrated thermodynamicquantities as equations of state.From the p n representation it comes out that thecoexistence line of the clusterized matter is embedded in thespinodal region of the uniform matter and shifted toward highervalues of respective to the coexistence of nuclear matter; thatis, it corresponds to the low-density part of the spinodal regionof homogeneous matter. This fact, together with the absence ofinstability of the clusterized system, is already an indicationthat the homogeneous matter instability is effectively curedby accounting for the possibility of cluster formation, as wewill develop in the next section. In the end, we mention thatthe phase diagram of clusterized matter might slightly changewhile modifying cluster properties.3. Mixture propertiesIn the previous section we have examined the thermo-dynamics of purely uniform and purely clusterized nuclearsystems. In the physical situation of supernovae and proto-neutron star matter it is clear that these two components haveto be simultaneously present, and the question arises of howdescribing their mutual equilibrium conditions. In many pastworks [17,19,20,51,66], the evolution of clusterized mattertoward homogeneous matter at high density was described asa first-order phase transition, governed by Gibbs or Maxwellequilibrium rules. The simultaneous presence in the Wigner-Seitz cell of localized clusters and homogeneously distributedfree protons and neutrons in the nucleonic drip region isthen also depicted as a manifestation of the same first-ordertransition, though corrected by surface and Coulomb effectsarising from the finite size of the cell.However, as we discuss now, the proper way of treating den-sity inhomogeneities occurring at a microscopic (femtometerscale) level is the concept of a gas mixture, if the coexistingcomponents (here nucleons and fragments) are noninteracting.If interaction effects arising from the excluded volume areincluded, we show that such a mixture naturally producesa continuous (second order) transition to uniform matter atdensities close to saturation.Restricting for the moment to the baryonic sector, thethermodynamics of a system with two conserved charges nand p is a two-dimensional problem. The other charges donot affect the thermodynamics at the baryonic level becausethey are not coupled to the baryons [see Eq. (1)]. The netelectron density is coupled to the proton density throughthe electromagnetic interaction, but this does not imply anadditional degree of freedom as charge neutrality imposesp = e. This two-dimensional problem can be easily taggedif we switch from the grand-canonical ensemble, where allthe intensive parameters are constrained and the densities canfreely vary, to a statistical ensemble where one single extensiveobservable (i.e., a single density) is free to fluctuate [53].In our two-dimensional problem, two equivalent possibilitiesexist, namely the proton-canonical, neutron-grand-canonicalensemble where the thermodynamic potential is Legendretransformed to its conjugated extensive variable p,spc(, n, p) = 1VS[,n]= 1VlnZgc(, n, p) pp, (31)065801-7AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)-45-40-35-30-25-2010-510-410-310-210-1p (fm-3) p (MeV)1010 210 3-0.15 -0.1 -0.05 0 0.05HMclustersP (MeV/fm3)1/ p (fm3 )-43-42-410 0.02FIG. 4. (Color online) Properties of the uniform (solid lines) andnet-charge-neutralized clusterized (solid circles connected by lines)matter at T = 5 MeV and n = 0 MeV. Homogeneous matter isdescribed by SKM*.and the neutron-canonical, proton-grand-canonical ensembledescribed by the thermodynamic potential,snc(, n, p) = 1VS[,p]= 1VlnZgc(, n, p) nn. (32)We now come to illustrate the way of mixing clustersand homogeneous matter and the differences between suchmixture and a possible phase coexistence which would beassociated to a first-order phase transition. Figure 4 showsthe equations of state p(p) and P (1/p) for the proton-canonical ensemble, as well as the functional relation betweenthe two intensive parameters p and P for the uniform (solidlines) and clusterized (solid circles connected by lines) matterat T = 5 MeV and a representative neutron chemical potentialn = 0 MeV. The clusterized system is confined within avolume V = 2.9 104 fm3.The back-bendings of p(p)|(HM) and P (1/p)|(HM) andthe tilting shape of p(P )|(HM) indicate that in the consideredp range uniform matter, if it would exist alone, would beunstable against phase separation. If we consider that, atthe microscopic scale of the Wigner-Seitz cell, clusters and freenucleons may be simultaneously present, as expected from thephenomenon of neutron drip, then the statistical equilibriumconditions read:(HM)p = (cl)p ; (HM)n = (cl)n ;(33)T (HM) = T (cl); p(bar) = p(HM)wHM + p(cl)wcl,where wHM and wcl measure the volume fractions associatedto the homogeneous matter and clusters,wHM = (V vcl)/V,(34)wcl ={V v0[(HM)n + (HM)p](V vcl)}/V,with v0 standing for the intrinsic volume of a nucleon calcu-lated such as v0 = 1/0 and vcl = v0NCi=1 Ai representing thesummed-up clusters intrinsic volumes.Let us denote by (HM)p(H ) and (HM)p(L) the proton chemicalpotentials corresponding to the limits of the spinodal regionof homogeneous matter. For p < (HM)p(L) 46 MeV, therewill be one single low-density solution of uniform matterwhich can constitute a mixture at thermal equilibrium withclusterized matter according to Eqs. (33) and which, forsufficiently small values of p, may differ by orders ofmagnitude from the density of clusters. For (HM)p(L) < p (HM)p(H ) 24 MeV the thermodynamic condition is inprinciple the same as for high and negative chemical potentials.Since, however, the sum of the two densities corresponding tohomogeneous and clusterized matter exceeds 0, the situationgoes beyond the applicability domain of the present model andis, thus, devoid of interest.Given that one of the purposes of our article is to addressthe unhomogeneous (crust)-homogeneous (core) transition, itis possible now to anticipate that for these two representativetemperatures T = 1.6 and 5 MeV this will take place when,following a constant proton fraction trajectory, the uniformmatter state will not correspond anymore to the low-densitybut rather to the high-density solution.At higher temperatures the situation is slightly differentas shown by Fig. 5, which shows the same quantities plottedin Fig. 4, for a representative neutron chemical potential atT = 10 MeV. In this case, in the relatively narrow chemicalpotential region where three different homogeneous mattermacrostates are found, the corresponding cluster states arealways at high density close to saturation. Then, the onlymixture solution not overcoming 0 is the one given by thelow-density solution for uniform matter. At this temperaturethen, the transition to homogeneous matter does not occuras an increasing proportion of homogeneous matter thatbecomes abruptly dominant in the mixture as we havejust seen at lower temperatures but rather as an increasingsize of the largest cluster which abruptly starts to fill thewhole available volume. The two transition mechanisms arephysically very close to each other: an infinitely large clusterhas the same energetic and entropic properties of homogeneousmatter at saturation density, provided the parameters of theeffective interaction are compatible with the cluster energyparameters.Since the thermodynamic potential of the mixture canbe expressed as the sum of the thermodynamic potentials065801-8STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)-55-50-45-40-35-30-25-20-15-1010-410-310-210-1HMclustersp (fm-3) p (MeV)1010 210 3-0.2 -0.1 0 0.1 0.2P (MeV/fm3)1/ p (fm3 )FIG. 5. (Color online) The same as in Fig. 4 for T = 10 MeV andn = 20 MeV.corresponding to each subsystem:G(bar)(, n, p, V ) = G(cl)(, n, p, V )+G(HM)(, n, p, V ), (35)the total neutron and proton density of the mixture is obtainedby summing up the corresponding particle numbers,Nn,p =[G(cl)(n,p)+ G(HM)(n,p)],p,n,V= N (cl)n,p +N (HM)n,p , (36)with the extra condition that nucleons cannot be found in thespace occupied by the clusters and vice versa.The densities then result:n,p = Nn,pV= N(cl)n,pV v0[(HM)n + (HM)p](V vcl) V v0[(HM)n + (HM)p](V vcl)V+ N(HM)n,pV vcl V vclV= N(cl)n,pV v0[(HM)n + (HM)p](V vcl)wcl +N (HM)n,pV vcl wHM. (37)The excluded volume correction plays a sizable roleonly in the vicinity of 0, and is in particular responsiblefor the unhomogeneous-homogeneous matter transition athigh temperature. More precisely, pure homogeneous nuclearmatter phase occurs if [(HM)n + (HM)p ] > 0 and, respectively,(V vcl) 0.The equation of state and constrained entropy of themixture for (T = 5 MeV, n = 0 MeV) and (T = 10 MeV,10-310-210-110-510-410-310-210-1Spc/V (fm-3)-50-45-40-35-30-25-2010-510-410-310-210-1HMclustersmixture (pmix< 0)coex.mixture (pmix> 0)p (fm-3) p (MeV)00.20.40 0.02 0.04FIG. 6. (Color online)p versus p (bottom) and S[,n]/Vversus p (top) for the homogeneous-clusterized matter coexistence(open diamonds) and mixture (open circles if pmix > 0 and starsif pmix < 0) at T = 5 MeV along the constant n = 0 MeV path.The homogeneous matter and net-charge-neutralized cluster mattercurves are plotted with solid line and, respectively, small solidsymbols.n = 20 MeV) are indicated by circles in Figs. 6 and 7. Wecan see that the simultaneous presence of the two componentsconstitutes always an entropy gain respect to both purely10-310-210-110-410-310-210-1Spc/V (fm-3)-55-50-45-40-35-30-25-20-1510-410-310-210-1HM (SKM*)clustersmixturecoex.p (fm-3) p (MeV)00.10 0.05FIG. 7. (Color online) The same as in Fig. 6 for T = 10 MeV andn = 20 MeV.065801-9AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)homogeneous and purely clusterized matter and correspondstherefore always to the equilibrium solution.At the lowest temperature the p versus p equationof state is bivalued in a given density domain. This canbe understood as follows. At temperatures low enough thatthe cluster component shows a liquid-gas transition, whichcan be recognized from the plateau in the p(p) relationat constant n, two very different cluster-matter mixturescan be obtained from Eqs. (33) and (37) at the same totaldensity but different chemical potentials. These two possiblemixtures correspond to a homogeneous matter solution lyinginside the spinodal combined respectively to a low-densitylow-chemical-potential cluster solution (stars in Fig. 6) orto a higher-densityhigher-chemical-potential solution. Sincethese two possibilities correspond to the same total density, theequilibrium solution will be the one maximizing the associatedentropy. It comes out that the lowest chemical potentialsolution always correspond to a negative total pressure for themixture, while the pressure associated to the highest chemicalpotential solution is always positive, meaning that the entropyis higher (see the inset in Fig. 6).The states denoted by stars in Fig. 6 are therefore unstable ormetastable and will not be further considered. The presence ofsuch multiple solutions in this density and temperature domainmay be an indication that some degrees of freedom are missingin the model and cluster deformations should be added in thisregion where pasta structures are expected [1114].To summarize, we have modelized the simultaneous pres-ence of free nucleons and light and heavy nuclei as two inde-pendent components which are allowed to exchange particlesand energy following the laws of statistical equilibrium andwhose proportions are further determined by the geometricalconstraint given by the excluded volume. This means that thedensities will have different values in each subsystem and thatthe partial pressures will be also different, following Eqs. (33)and (37).The chemical potential equality between the two compo-nents insures the global maximization of the total entropy of themixture. As a result, a smooth decreasing weight of the clustercomponent is first observed with increasing density, followedby a successive smooth but steeper decrease as densityapproaches 0. This last result can be intuitively understoodfrom the simple physical fact that when density increases, thespace accessible to a clusterized system decreases and so doesthe effective number of degrees of freedom. As we will discussin greater detail in the following, this smooth behavior of themixture does in no way imply that the crust-core transitionis smooth. Indeed, because of the characteristic back-bendingbehavior of the homogeneous matter equation of state, whichreflects the well-known instability respect to spontaneousdensity fluctuations, the density content of the homogeneouscomponent can discontinuously change as a function of thechemical potential.4. Mixture versus phase coexistenceThe crust-core transition of a cold neutron star is often[1720] depicted as a first-order transition from a solid phaseof nuclei on a Wigner lattice to a liquid phase of homogeneousnuclear matter. By continuity, a similar phenomenology isthought to occur at the finite temperature and more isospinsymmetric conditions involved in supernova matter [23,51].This hypothesis arises from essentially two considerationsthat we now turn to examine: (i) The intermediate configurationbetween the solid crust and the liquid core, namely finitenuclei immersed in a homogeneous neutron background,intuitively resembles to a coexistence between a clustered anda homogeneous baryonic phase and (ii) normal nuclear matteris known to present an instability respect to cluster formationand phase separation. This spinodal instability is associated toa first-order transition. By continuity the same is expected inthe baryonic part of stellar matter.Let us first concentrate on the first statement. As wehave already mentioned, the simultaneous presence of twodifferent states characterized by a different average value foran order parameter (here: baryonic density) constitutes a phasecoexistence if and only if the order parameter fluctuations(here: density inhomogeneities) occur on a length scale whichis 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 becausethe Gibbs construction neglects the interface energy andentropy, which is essential to describe thermodynamics offinite systems [67]. Since this is clearly not the case in stellarmatter, 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 ofa microscopic mixture between clusters and homogeneousmatter is also a way to describe the intermediate configurationsbetween core and crust. Therefore we may ask if the twomodelizations of phase mixture and phase coexistence areequivalent or at least lead to equivalent results for theobservables of interest.To clarify this point, let us consider homogeneous andclusterized matter as two different possible phases for stellarmatter which may coexist at a first-order transition point. Theadvantage of working with a hybrid ensemble as the proton-canonical ensemble defined above is that the two-dimensionaltangent construction implied by Gibbs conditions in a multi-component first-order phase transition [68,69] is reduced to asimple one-dimensional Maxwell construction on the uniquefree extensive variable [53]. It is important to stress that thisprocedure is just a technical simplification and the obtainedequilibrium is exactly identical to the usual two-dimensionalGibbs construction in the grand-canonical ensemble. Thisconstruction in a Legendre-transformed potential should not beconfused with a Maxwell construction on the p() equationof state performed on trajectories characterized by constantproton fractions Yp = p/, which is currently used in theliterature [23] and has no physical justification. Indeed, ina first-order phase transition in a two-component system thechemical composition of the two coexisting phases is neverthe same and the total pressure is a monotonically increasingfunction of the total density [68]. This means that a Maxwellconstruction on p() at constant Yp does not give the correctGibbs equilibrium.The inset in the right upper panel of Fig. 4 showsthat for the representative (T ,n) choice of the figure the065801-10STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)two components are found at the same value of protonchemical potential and pressure in a single point locatedat (P t , tp) (0.005 MeV/fm3 ,41.8 MeV). Only at thispoint the Gibbs equilibrium rule for coexisting phases isfulfilled,(HM)p = (cl)p ; (HM)n = (cl)n ; T (HM) = T (cl); p(HM) = p(cl).(38)As it can be seen in Fig. 4, at this point the density of thetwo components differs of more than two orders of magnitude.If clusters and homogeneous matter could be viewed as twodifferent phases for stellar matter, the global system would becomposed of homogeneous matter only for p < tp, and ofclusters only for p > tp, with a discontinuous (first-order)transition at the transition point [P t (T ,n), tp(T ,n)]. Theeffect of this construction on the equation of state is shown onFig. 6 by the open diamonds. Within the picture of a first-orderphase transition, stellar matter at the chosen representative(T ,n) point would be homogeneous at low density upto p = mp , with mp 1.2 105/fm3; for intermediatedensities mp p Mp , with Mp 2 103/fm3, it wouldbe a combination of matter and clusters varying in linearproportions; and at high density p Mp it would be solelycomposed of clusters, which would extend over the wholevolume at the highest densities. As we can see, this equationof state differs completely from the physically meaningful onecorresponding to a continuous mixture of matter and clustersover the whole density domain. However, the two models leadto an almost equivalent entropy, as it is shown by the upperpart of Fig. 6.Similar considerations apply at other temperatures andchemical potentials. As a second example, Fig. 7 comparesthe properties of the mixture with the ones of a hypotheticalfirst-order transition for T = 10 MeV andn = 20 MeV, thesame thermodynamic conditions as in Fig. 5. Only within themixture picture free protons in weak proportions are presentin the matter, as it is physically expected. An inspection of thebehavior of the constrained entropy at this high temperature(higher part of Fig. 7) clearly shows that the continuousmixture of nucleons and clusters not only correctly recoversthe homogeneous matter properties at high density but is alsovery effective in compensating the mean-field instability ofhomogeneous matter at lower densities, leading to a concaveentropy over the whole density range (see the inset in the upperpanel). Thus the usual argument that a phase coexistence wouldbe needed to maximize the entropy and cure the convexityanomaly of the entropy cannot be applied: the same entropicgain is obtained considering a continuous mixture, withoutapplying an artificial construction which neglects all finite-sizeeffects.We note by passing that at high temperature there are two(p,p) points satisfying the Gibbs conditions Eq. (38) for agiven value of n (see Fig. 5). However, at the second pointat high pressure the two components correspond to almost thesame density. Making a mixture or a phase coexistence doesnot therefore produce any sensitive difference in the equationof state (lower part of Fig. 7) in this density region.-60-40-200mixturecoex. p (MeV)02468Acl00.050.10.1510-210-1tot (fm-3)Ncl/AbFIG. 8. (Color online) Mixture versus coexistence for T =10 MeV and n = 20 MeV. (Top) p versus tot; (middle) averagemass of a cluster versus tot; (bottom) average cluster number pertotal baryon number versus tot.The matter composition obtained with a continuous mixtureis compared to the one corresponding to a discontinuoustransition obtained with the Gibbs construction in Fig. 8.We can see, together with the different behavior of theequation of state, that the predictions for the average clustersize and cluster multiplicity show quantitative differences.It is difficult to know a priori the consequence of suchaltogether slight differences on the macroscopic modelizationof supernova matter, and we cannot exclude that the useof a mixture equation of state would not be really distin-guishable from the currently used approximations employingfirst-order phase transitions. Even if this would be the case,the discontinuities in all thermodynamic quantities impliedby these artificial first-order transitions lead to numericalinstabilities in supernova codes [70] which are difficult tohandle and which could be avoided using continuous equationsof state.To defend the possibility of a first-order crust-core tran-sition, one may argue that the behavior of the entropy as afunction of the density shown in Figs. 6 and 7 is certainlymodel dependent and that in general residual convexities couldstill be present after considering the mixture of nucleonsand clusters, especially at the lowest temperatures. Suchconvexities would represent residual instabilities and shouldthen be cured with a Gibbs construction, leading again to a first-order phase transition from a low-density, cluster-dominated065801-11AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)-1.6-1.4-1.2-1-0.8-0.6-0.4-0.200.20 0.01 0.02 0.03 0.04 0.05 0.06HM (SKM*)mixturep (fm-3)Spc/V (fm-3)-1.6-1.4-1.2-1-0.8-0.6-0.4-0.200 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09HM (SKM*)mixturep (fm-3)Spc/V (fm-3)FIG. 9. (Color online) Constrained entropy for the ensemble baryons plus electrons at T = 5 MeV along the constant n = 0 MeV path(left panel) and at T = 10 MeV and n = 20 MeV (right panels) for the homogeneous matter (full line) and the mixture of matter and clusters(open circles).mixture to a high-density, matter-dominated mixture as inthe usual representation of a neutron star. Such an instabilityis known to be present in homogeneous nuclear matter, asshown by the back-bendings of Figs. 4 and 5, and it maywell survive if clusterization was allowed in a microscopicway [50] going beyond the mean-field approximation ina microscopically consistent way. The argument exposedabove is the second argument which can be invoked tojustify a discontinuous transition, which was mentioned atthe beginning of this section and that we now come todiscuss.It is certainly correct that our model depicting the simul-taneous presence of nucleons and clusters as a noninteractingmixture of two different components associated to completelyindependent degrees of freedom and self-energies is highlyschematic, and improvements are in order. In particular, theway the excluded volume is parametrized has an influence onthe quantitative balance between nucleons and clusters, with apossible consequence on the precise location of the transitionwhich is not completely under control.Fully self-consistent treatments of the quantal many-bodyproblem at finite temperature as in Refs. [4750], includingthe possibility of self-consistent clustering at low densities, areextremely promising and will in the long run provide a muchbetter description of stellar matter than the phenomenologicalmodel proposed here. Concerning the issue of the order of thetransition, however, we believe that the availability of suchsophisticated calculations would not modify our conclusions.Indeed, as we have already stated, the thermodynamic stabilityhas to be checked on the global thermodynamic potentialEq. (1), after all the different constituents are added up. If s and neutrinos are completely decoupled from the baryonicsector, this is not true concerning electrons and positronsbecause of the neutrality constraint p = e. Because of this,the properties of the electron free energy will influence theconvexity properties of the baryonic entropy as a function ofthe density.Figure 9 shows the thermodynamic potential of the proton-canonical, neutrongrand-canonical ensemble in the samethermodynamic conditions as in Figs. 6 and 7, including thecontributions of electrons (detailed in the next section) andexplicitly including the neutrality constraints totpc (, n, p) =1VlnZ (cl)gc (, n, p)+ 1VlnZ (HM)gc (, n, p)+ 1VlnZ (el)(,e) (p + e)p. (39)We can see that taking into account the presence of clustershas virtually no effect on the convexity of the global entropy,which is largely dominated by the convexity properties of theelectron thermodynamic potential. Since this latter is highlyconcave, whatever the treatment of the baryonic sector, theglobal entropy will stay concave in the whole density domain.The effect of the electrons is thus to compensate the instabilityof nuclear matter and quench the possibility of a first-orderphase transition.This result can be physically understood [16,71] as follows:Because of the low electron mass, the electron gas is muchmore incompressible than the nuclear component. A fluctua-tion in the baryonic density creates a fluctuation in the chargedensity, which because of the neutrality constraint impliesa fluctuation in the net electron density. Such fluctuationsare naturally present in the length scale of the Wigner-Seitzcell and are at the origin of the electrostatic interactionenergy which is nonzero even if the system is neutral; seeEq. (25). However, these fluctuations cannot extend overmacroscopic lengths because they would imply macroscopicelectron inhomogeneities, which are contrasted by the highelectron incompressibility.To conclude, phases of different global baryon densitycannot exist in stellar matter if baryonic matter is constituted ofheavy particles as nucleons and clusters, as it is the case at thedensities of interest in the present article. Since this statementis valid for all chemical potentials and temperatures, includingT = 0, we conclude that the crust-core phase transition cannotbe of first order at any temperature.065801-12STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)B. Leptons and photonsIn the temperature and density domains relevant for thisarticle, T > 1 MeV and > 106 g/cm3, leptons (electronsand neutrinos) are relativistic, in particle-antiparticle pairequilibrium and in thermal equilibrium with nuclear matter[17,23].Given that the net charge of the electron gas (that isthe number of electrons minus the number of positrons)has to neutralize the positive charge of protons in uniformand clusterized matter, the electron chemical potential eis determined by the equality among the electron density,calculated in the relativistic Fermi gas model, and the protondensity,Yp = ge62(ehc)3 [1 + 12e(2T 2 32m2ec4)], (40)where represents the baryon density, ge = 2 is the spindegeneracy, and me is the rest mass.The electron gas pressure, entropy, and energy densities arefunctions of T and e and their expressions are the following,p(el) = gee242(ehc)3 [1 + 12e(22T 2 3m2ec4)+ 2T 24e(7152T 2 12m2ec4)], (41)s(el) = geT 2e6(hc)3[1 + 12e(7152T 2 12m2ec4)], (42)and, respectively,e(el) = gee82(ehc)3 [1 + 12e(22T 2 m2ec4)+ 2T 24e(7152T 2 12m2ec4)]. (43)Neutrinos also form a relativistic Fermi gas in pair equilib-rium and thermal equilibrium with nuclear matter. This meansthat Eqs. (40), (41), (42), and (43) with m = 0 and g = 1instead of me and ge describe their density, pressure, entropyand energy densities. Under the -equilibrium hypothesis,neutrinos chemical potential is dictated by the chemicalpotential of neutrons, protons, and electrons,(p +mpc2) + (e +mec2) = (n +mnc2) + . (44)If not explicitly mentioned otherwise (Sec. III C), all overthis article we shall work out of equilibrium. Neutrinoscontribution to the total energy, entropy, and pressure willbe disregarded as well. The motivation of this choice isthat equilibration with respect to weak interaction is oftennot achieved over the time scales of many astrophysicalphenomena, but we will discuss this point further in Sec. III C.Photons are assumed to be in thermal equilibrium with theother star ingredients and their pressure, entropy and energydensities arep( ) = 2T 445(hc)3, (45)s( ) = 4p( )T, (46)and, respectively,e( ) = 3p( ). (47)III. STAR MATTER RELEVANT RESULTSA. Constant chemical potential pathsIn the evolution of supernovae explosions, nuclear com-position plays an extremely important role. First, isotopicabundances in the light and medium mass region influencenucleosynthesis via r, rp, and photodisintegration processes.Then, nuclei with masses larger than 4060 determine the rateof electron capture [33,34], the most important weak nuclearinteraction in the dynamics of stellar core collapse and bounce,responsible for the neutron enrichment of baryonic matter andneutrino emission. On the other hand, nuclei absorb energyvia internal excitation making the increase in temperature andpressure during star collapse be less strong than in the case of auniform nucleon gas [27]. Nuclei influence also the maximumdensity reached during the collapse. In the absence of nuclei thecollapse stops before reaching 0, while the limiting value islarger than 0 if nuclei are present [28]. Finally, by dissociationunder the shock wave, nuclear composition affects the shockpropagation [30].Under this perspective, it is obvious that realistic starmatter equations of state require accurate treatment of thenuclear statistical equilibrium. Working within the SNA, thefirst generation models [17,23] do not offer a satisfactorydescription of nuclear miscellanea and may be to someextent responsible for the failure of present supernovaesimulations to produce explosions, even if realistic neutrinotransport [72] and convective instabilities in multidimensionalhydrodynamics [73] are probably the key ingredient of asuccessful explosion of highly massive cores. Most recentstatistical models [42,44,51] fix this inconvenience by workingin a grand-canonical approximation (GCA) which allows forcontinuous distributions in the nuclear cluster mass and charge.This means that practically everywhere in the temperature-density domain, the cluster gas may contain a variety ofnuclei of very different sizes. In our model these, additionally,may change from one configuration to the other, meaning thatthe composition fluctuation is also accounted for.To understand the correlation between the thermodynamicproperties and the matter composition obtained with ourmodel, Fig. 10 presents the evolution with density of theaverage cluster size and isotopic composition at the threetemperatures T = 1.6, 5, and 10 MeV analyzed before, and forsome chosen characteristic values of the chemical potentials,in comparison to the Lattimer-Swesty (LS) model employingthe SNA [74]. As a technical detail, we mention that atT = 1.6 and 5 MeV the clusterized matter was treated withina canonical-constant Z ensemble, while a grand-canonical onewas employed at T = 10 MeV. In all considered cases, the cellvolume spanned the domain (1.9 104) (1.4 108) fm3.The left part of the figure shows the evolution withtemperature and density of the average cluster size, comparedto the single nucleus size of LS, corresponding to the sameeffective nucleon-nucleon potential, SKM*. The general trend,065801-13AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)10 210 310-310-210-1Acl0.20.30.40.510-310-210-1Zcl/Acl1010 210-310-210-1Acl0.20.30.40.510-310-210-1Zcl/Acl1010-310-210-1tot (fm-3)Acl0.20.30.40.510-310-210-1tot (fm-3)Zcl/AclFIG. 10. (Color online) Evolution with totalbaryonic density of the cluster size (left) andisospin composition (right) along trajectories ofconstant p for T = 1.6, 5, and 10 MeV. (Opencircles, squares, and triangles) Predictions of thepresent model corresponding to a positive totalpressure of the mixture pmix > 0; (open stars)the same as above for a negative total pressureof the mixturepmix < 0; (dashed lines) Lattimer-Swesty data from Ref. [74].which is independent of the chosen values of the chemicalpotential, is an increase with density and a decrease withtemperature of the cluster size. This behavior is smooth, eventhough from the thermodynamic viewpoint at T = 1.6 MeVand T = 5 MeV the calculation corresponds to the coexistencezone for the cluster component (stars in the left part of Figs. 1and 2), while T = 10 MeV is supercritical for clusters. Thissmoothness is due to the presence of the homogeneous mattercomponent as a continuous mixture, as we have discussedin Sec. II A3. The only indication of the underlying clusterphase transition can be seen from the fluctuations of the sizedistribution, indicated by the vertical bars in Fig. 10, which aremaximal at T = 5 MeV in the middle of the coexistence regionas it is well known from multifragmentation studies in finitenuclei [7577]. The open stars in Fig. 10 denote the unstable ormetastable mixture solutions discussed in Sec. II A3, obtainedwhen the homogeneous matter component is characterizedby a density very similar to the density of the clusters. Aswe have already mentioned, the presence of these multiplesolutions with similar entropy content may be an indicationthat allowing for pasta-like clusters may provide a higherentropy solution.Coming to the comparison to LS, the two models nicelyagree at low density and temperature, but important differencesare seen elsewhere. In the LS model the crust-core transition isa discontinuous (first-order) disappearance of the single heavycluster which merges at high density and temperature withhomogeneous matter: this is why Acl decreases with densityat T = 5 MeV in LS, and clusters abruptly disappear at thetemperature-dependent transition points. This is at variancewith our calculation, where both components exist, thoughin very different proportions, in the whole temperature andchemical potentials space.Considering particles together with the heavy nucleuswhen discussing the cluster component in the LS case anddefining Acl = (Aheavy + 4N)/(N + 1) would not changethe situation. The explanation relies on the fact that, for thesetrajectories, when a heavy cluster exists the number of particles in a cell is negligible.The right part of Fig. 10 gives the average cluster chemicalcomposition, while the fluctuations of this quantity are givenby the vertical bars. We can see that at low temperature (belowthe cluster limiting point) the isospin content of clusters isfairly constant both in our calculation and in the LS model. Thisis expected, since by definition a trajectory moving inside acoexistence region implies a transformation p = transp = cteand n = transn = cte, that is Z/A = cte. Since a hugepart of the density plane at low temperature belongs to thecoexistence region if the cluster component is consideredalone (see Figs. 1 and 2), the chosen trajectory can span thewhole relevant density domain staying inside the coexistence,that is, at constant chemical composition. The little deviationsobserved in our calculation from this constant behavior athigh density have to be considered as numerical accidentsdue to the nonperfect achievement of the thermodynamiclimit.In the case of T = 10 MeV, the chemical composition isentirely determined by the chosen transformation and hasno specific physical meaning. It is interesting to remarkthat, except at the lowest temperature, this quantity presentsvery huge fluctuations that cannot be addressed within theSNA. Changing the chosen chemical potential value obviouslymodifies the behavior of the chemical composition withdensity, but does not change our qualitative conclusions.Figure 11 displays the characteristics of the Wigner-Seitzcell in terms of baryonic number, proton fraction, and linear065801-14STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)050010001500200010-310-210-1AWS00.20.40.610-310-210-1YWS020406010-310-210-1L WS (fm)020040060080010-310-210-1AWS00.20.40.610-310-210-1YWS020406010-310-210-1L WS (fm)11010 210 310-310-210-1tot (fm-3)AWS00.20.40.610-310-210-1tot (fm-3)YWS05010010-310-210-1tot (fm-3)L WS (fm)FIG. 11. (Color online) Evolution with total baryonic density of the total baryonic number (left), proton fraction (middle), and lineardimension (right) of the Wigner-Seitz cell for T = 1.6, 5, and 10 MeV along trajectories of constant p . (Open symbols) Predictions of thepresent model; (dashed lines) Lattimer-Swesty data from Ref. [74].size in the same thermodynamic conditions of Fig. 10. TheWigner-Seitz size is defined as the total baryonic massassociated to each fragment,AWS = i AiNcl, (48)where the sum runs over clusters and free nucleons and Ncl =NC is the average cluster number obtained in the simulation.It is important to stress that the concept of a Wigner-Seitz celldiffers substantially in our approach and in SNA. While inSNA the Wigner-Seitz cell is a primitive cell, i.e., the smallestvolume element which contains complete information on theinfinite system and from which the last one can be built, in GCAwe work with an ensemble of tens or hundreds of elementarycells. None of these cells is primitive, as their cluster partitionsmay differ. Only within fluctuations they may be seen as replicaof an average elementary cell, the counterpart of the Wigner-Seitz cell in SNA.The decrease of AWS with increasing density, as observedboth in LS and in our calculations at high temperature,can be understood as a trivial compression effect, whichimplies that reducing the volume the matter composition doesnot drastically change. A different behavior is seen at lowtemperature, due to the rapid increase of the heaviest clustersize toward the transition to the core.As before, open stars denote unstable mixture solutionsand we cannot exclude that the crust-core transition wouldbe more abrupt if additional degrees of freedom were addedin our model at high density. Indeed, at the peak a smallcluster coexists with a huge quantity of slightly diluted matter.This situation is energetically similar to a very big bubblein a dense = 0 medium, but the absence of deformationdegrees of freedom may suppress the entropy of such exoticconfigurations.Concerning the comparison with the LS model, the twocalculations fairly agree only at low temperature. As a generaltrend, taking into account the whole equilibrium distributionof clusters (GCA) produces smaller fragments (see Fig. 10, leftpanel) with higher multiplicity respect to the SNA. This is inqualitative agreement with other authors findings [42,44,51].The absence of points at intermediate density at T = 5 MeVis not an effect of the discontinuous transitions in LS butis simply due to the fact that LS calculations are performedin discrete Yp steps, and therefore some points are missingwhen looking for = cte trajectories. The much higher AWSseen in LS calculations at high temperature is due to the factthat no fragments other than s are present there in the SNA(see Fig. 10). In this situation (tot < 3 103 fm3 for T =5 MeV and the whole considered total density domain at T =10 MeV), the cell is built around an particle such that thetotal baryonic mass of a cell becomes,AWS = (Nn +Np + 4N)/N= 12 + 4(n + p)/, (49)065801-15AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)where Nn, Np, N , and n, p, denote respectively thenumber of free neutrons, free protons, and particles perunit volume and the corresponding mass fractions. This makesthe Wigner-Seitz mass increase because of our definitionAtot/V = AWSNcl/V .The global proton fraction of the Wigner-Seitz cell is shownin the middle panels of Fig. 11 as a function of temperature andtotal density. At the lowest temperature fragments dominatein the global composition, which implies that the behavior ofthe cell closely follows the behavior of the proton fraction ofthe clusters shown in Fig. 10. At higher temperature, clustersstart to be immersed in a nuclear gas which is neutron rich.As density increases, proton drip becomes allowed, slightlyincreasing the proton fraction of the free nucleons. Since matterdominates over fragments at high density and temperature,YWSreflects this behavior.Finally the right part of Fig. 11 shows the average distanceamong fragments, which can be physically interpreted as thelinear size of the Wigner-Seitz cell:LWS =(VNcl)1/3. (50)In the case of SNA, at temperatures where the heavy nucleusis no longer present, this quantity is replaced with the averagedistance among alphas. LWS decreases with density andtemperature as expected. Again, only in the thermodynamicconditions where a single massive cluster dominates in theequilibrium distribution (that is, at the lowest temperature) theresults agree with the SNA which assumes a single cluster inthe whole phase diagram.B. Results at constant proton fractionThe construction of the phase diagram and mixture prop-erties requires a grand-canonical formulation, that is, to fixthe intensive observables T , n, and p, or, alternatively, acanonical formulation working in terms of T , n, and p,where n and p stand for the total neutron and protondensities, i.e., the summed up contributions of uniform andclusterized matter. For the purpose of studying the stabilityproperties of the different phases, hybrid ensembles whereonly one charge is let free to fluctuate have also been employed.These variables are, however, not convenient for expressing theequation of state of stellar matter, as the chemical potentialsare not measurable and the total baryonic density is of outmostimportance for astrophysical applications.The mostly spread choice is then to work in terms ofT , , and Yp, where = n + p is the total baryonicdensity and Yp = p/ is the proton fraction. The net chargeneutrality requires the proton and electron fractions to be equal,Yp = Ye.Let us adopt these coordinates and investigate which arethe domains of n and p spanned by the total system withtwo arbitrarily chosen values of Yp = 0.2 and 0.3 if the totaldensity ranges from 0 to 0/1000. The results are plottedin Fig. 12 for the same values of temperature consideredbefore, T = 1.6, 5, and 10 MeV. Present model predictions-80-70-60-50-40-30-20-100102010-410-310-210-1T=1.6 MeVT=5 MeVT=10 MeV n (MeV)-80-70-60-50-40-30-20-100102010-410-310-210-1tot (fm-3) n (MeV)-100-90-80-70-60-50-40-30-20-10010-410-310-210-1 p (MeV)-80-70-60-50-40-30-20-10010-410-310-210-1tot (fm-3) p (MeV)FIG. 12. (Color online) n (left panels) and p (right panels) versus total baryonic density for different values of the constant protonfraction Ye = 0.2, 0.3 and temperature T = 1.6, 5, and 10 MeV. Present model predictions are illustrated with open circles (T = 1.6 MeV),squares (T = 5 MeV), and triangles (T = 10 MeV); dashed lines correspond to Lattimer-Swesty results [74]. In both cases the uniform mattercomponent is described by SKM*. The thin dotted line corresponds to the case in which, at T = 10 MeV, nuclear matter would exclusivelyconsist out of free nucleons.065801-16STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)-60-50-40-30-20-1001020HMmixture n (MeV)00.10.20.30.40.50.60.70.80.9110-310-210-1tot (fm-3)Yp-70-60-50-40-30-20-10010 p (MeV)00.10.20.30.40.50.60.70.80.9110-310-210-1tot (fm-3)YpFIG. 13. (Color online) n,p()|p,n (top) and Yp()|p,n (bottom) curves at T = 10 MeV. The values of the constant chemical potentialsare expressed in MeV. Solid lines correspond to the homogeneous matter; open circles correspond to the mixture between homogeneous andclusterized matter.(open circles, squares, and triangles) are systematically com-pared to LS predictions [74] (dashed lines). While for n()our predictions agree with LS, for p() this is true onlyin the limit of low and medium temperatures. At the high-est considered temperature (T = 10 MeV) and > 0/100present model predictions differ from the LS ones by few tensof MeV.This difference between the models can be understood as aneffect of the presence of clusters at high temperature implied byour approach. The thin dotted line corresponds to a calculationwhere the cluster contribution in the mixture is artificiallyswitched off, and the nuclear matter is made out exclusively offree nucleons. Quite remarkably, in this case our calculationagrees with LS predictions, which can be understood knowingthat the amount of mass contained in clusters in LS at hightemperature is negligible, as we have seen in the previoussection. We want to stress that the total entropy of the mixtureis always higher than the total entropy of homogeneous matter(see Fig. 7), meaning that our mixture solution is preferred atequilibrium.Surprisingly, this difference between the models is not seenfor the neutron chemical potential (left side of Fig. 12). Avery similar behavior was already observed in Ref. [51]. Thiscan be understood from inspection of Fig. 13, which displaysthe relation among n, p, Yp, and at T = 10 MeV for thehomogeneous matter and clusterized component.Let us look at the homogeneous matter first. Because of thesymmetry properties of the effective interaction, the behaviorof n() at constant p is the same as the behavior of p() atconstantn and shows the characteristic back-bending impliedby the homogeneous matter instability. Because of that, thebehavior of the proton fraction Yp() at constant n is alsoback-bending, and the behavior of Yp() at constant p isinverted respect to it. A constant Yp = 0.3 path then explores,as a function of the density, the back-bending region of theiso-p, while an analogous behavior for the iso-n would beseen in the (unphysical) symmetric trajectory Yp = 0.7. Asa consequence, p shows a maximum as a function of thedensity for this low value of the proton fraction while n is amonotonically increasing function, as observed in Fig. 12.Turning to the clusterized component, this latter does notshow any instability, meaning that the chemical potentials aremonotonically increasing functions of the density. Then thebehavior of Yp as a function of the density at constant (neutronor proton) chemical potential is essentially flat, and the neteffect is to translate the maximum of p in the constant protonfraction trajectories, while no effect is seen in n becauseof its monotonic behavior at this low Yp. This explains thedifferent behavior of the two chemical potentials at constantproton fraction seen in Fig. 12.We now turn to analyze some quantitative predictions ofour model concerning both matter composition and equations065801-17AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)10-410-310-210-1110-510-310-1X(Afrag)00.250.50.75110-510-310-1X(nuclei)00.050.10.1510-510-310-1Ncl/Ab10-410-310-210-1110-410-2X(Afrag)00.250.50.75110-410-2X(nuclei)00.050.10.1510-410-2Ncl/Ab10-410-310-210-1110-410-2tot (fm-3)X(Afrag)00.250.50.75110-410-2tot (fm-3)X(nuclei)00.050.10.1510-410-2tot (fm-3)Ncl/AbFIG. 14. (Color online) Evolution with total baryonic density of the average fragment mass fraction (left panels), cumulated cluster massfraction (middle panels), and cluster number per baryonic number (right panels) for T = 1.6, 5, and 10 MeV. Large and small open symbolsstand for present model predictions corresponding to Ye = 0.2 and, respectively, 0.3; thick (Ye = 0.2) and thin (Ye = 0.3) dashed lines illustrateLattimer-Swesty data from Ref. [74]; thick (Ye = 0.2) and thin (Ye = 0.3) dot-dashed lines correspond to the Horowitz-Schwenk results basedon the virial expansion of the equation of state [78].of state in thermodynamic conditions which are relevant forsupernova matter.1. Matter compositionFigure 14 illustrates the total baryonic density dependenceof the average fragment mass fraction Acl/AWS (left panels),cumulated fragment mass fractionNCi=1 Ai/AWS (middlepanels), and average number of clusters per baryonic number(right panels) for the three values of temperature alreadyconsidered, T = 1.6, 5, and 10 MeV, densities ranging from105 to 101 fm3, and Yp = 0.2 and 0.3. LS predictions(dashed lines) correspond to the unique heavy nucleus ofSNA whose mass fraction is immediately calculated out ofthe mass fractions of free neutrons, free protons, and s,heavy = 1 (n + p + ) (left panels) and, respectively,summed-up contributions of the heavy nucleus and particles(heavy + ) (middle panels). The fact that particles areconsidered for LS also when computing Ncl/AWS is respon-sible for the double-humped structure of the correspondingdistributions (right panels).The huge difference between SNA and GCA in whatregards Acl/Atot vs. is due to the different acceptations ofthe cell concept in the two frameworks. As we have alreadymentioned, contrary to SNA in GCA no primitive cell exists:because of the fluctuations implied by finite temperature, thesystem cannot be viewed as an infinite number of replicasof a single elementary volume element containing one singlecluster.Quite remarkably, the discrepancies originating from cellfluctuations are to a large extent washed-out when particlesare included (middle panels) and SNA and GCA agree muchbetter. However, at T = 5 MeV and, to a larger extent, atT = 10 MeV GCA accounts for more mass bound in massiveand light nuclei than SNA. This result may be considereda step forward in making numerical simulations of stellarevolution able to reproduce the supernovae explosions. Thesame qualitative information on cell population is revealedby the cluster number per baryon number depicted in theright panels, for which the maximum deviation is observed atthe highest temperature. The bimodal character of the SNAdistributions is an artefact of the fact that only two typesof nuclei, s and the heavy nucleus, are allowed to exist,each of them being preferentially produced in another totaldensity domain. A fully microscopic description of nuclearmatter composed of neutrons, protons and particles has beenrecently proposed by Horowitz and Schwenk [78] and relieson the virial expansion of the equation of state. The resultscorresponding to the density dependence of the cumulatedcluster mass fraction 4n/(nn + np + 4n) and average clusternumber per baryonic number n/(nn + np + 4n) are plottedwith dot-dashed lines in the middle and, respectively, rightpanels. The density domains over which this calculation ispossible was fixed by the condition of having subunitarynucleon fugacities. From a qualitative point of view, our modelagrees with HS in the sense that both of them show a smoothbehavior of Xnuclei and Ncl/Ab as a function of . From aquantitative point of view, the absence of nuclei heavier than particles makes HS systematically underestimate, in thedensity domains where such heavy cluster may exist, the per-centage of mass bound in clusters with respect to both presentmodel and LS. On the other hand, in the temperature-density065801-18STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)00.250.50.75110-610-510-410-310-210-1 HM/tot00.050.10.150.210-610-510-410-310-210-1 HM p/HM00.250.50.75110-410-310-210-1 HM/tot00.050.10.150.210-410-310-210-1 HM p/HM00.250.50.75110-410-310-210-1tot (fm-3) HM/tot00.050.10.150.210-410-310-210-1tot (fm-3) HM p/HMFIG. 15. (Color online) Evolution with totalbaryonic density of relative baryonic densityin uniform matter (left panels) and free protonfraction (right panels) for T = 1.6, 5, and 10 MeVand Ye = 0.2. Open circles, squares, and trianglescorrespond to the predictions of the present modelwhen uniform matter is calculated according toSKM*; dashed lines stand for Lattimer-Swestydata from Ref. [74]; dot-dashed lines correspondto Horowitz-Schwenk results based on the virialexpansion of the equation of state [78]; opendiamonds at T = 5 MeV correspond to presentmodel predictions in the case in which the uniformmatter is calculated according to Sly230a [79].regions where the largest cluster allowed by LS is the particle, there is no clear relation among the predictions of thetwo models. Thus, at T = 5 MeV and 3 103 fm3 HSand LS predictions agree perfectly while for T = 10 MeV and1 103 1.8 102 fm3 HS results exceed thosecorresponding to LS.2. Homogeneous-inhomogeneous transitionWe have already reported that, in the vicinity of 0, thestrong interplay between clusters and uniform matter whichis not allowed to occupy the same volume is responsiblefor a homogeneous-unhomogeneous matter transition. Thisis visible in Figs. 10, 11, and 14.A more suited representation is given in the left panels ofFigs. 15 and 16, where the relative baryonic density of uniformmatter (HM)/(w(HM)(HM) + w(cl)(cl)) is expressed as a func-tion of total baryonic density. The considered temperatures andproton fractions are T = 1.6, 5, and 10 MeV and, respectively,Yp = 0.2 (Fig. 15) and 0.3 (Fig. 16). The common pattern ofthese curves is a sudden jump toward (HM)/ = 1 takingplace at 0 and a much smoother increase toward 1 atmuch smaller values of total density. This means that closeto normal nuclear density and at low total densities, uniformnuclear matter dominates. While for the low-density domainthe SNA statement is true in the strict sense, in GCA it shouldbe rather understood as the dominance of free nucleons overclusters. As we have already seen, in GCA clusters survive, byconstruction, in the whole density-temperature domain, but atlow densities they represent a negligible percentage of matter.With the exception of T = 1.6 MeV where Monte CarloGCA convergence problems prevent accurate determinationof the transition point, SNA and GCA predictions agree well.Both models show that the transition density depends on thetemperature but, at least in the presently investigated range,only slightly on Yp.The interesting point here is that this transition is sponta-neously obtained in our model without performing an externalartificial Maxwell construction at high density, as is necessaryin standard NSE models [4244].References [80,81] show that the homogeneous-unhomogeneous transition density, as well as the neutron skinthickness, depends on the density dependence of the symmetryenergy and, more precisely, on L = 30[S2()/()]|=0 ,where S2() = [2(E/A)/22]|=0 with = (n p)/.To check to what extend this cold matter result holds atfinite temperature, we have calculated the transition regionfor T = 5 MeV also in the case in which homogeneousmatter is described by Sly230a [79] instead of SKM*. Theresult plotted with open diamonds in Figs. 15 and 16 showsno significative modification with respect to the previousones. Given that L(SKM) = 46 MeV and L(Sly230a) =43.9 MeV, we partially confirm the conclusions of Refs.[80,81]. In addition to this, within our model the transitiondensity as well as the proton enrichment of the uniform mattercould, at least in principle, depend also on how the clusterizedmatter component is implemented.The isospin composition of uniform matter as a function of is illustrated in the right panels of Figs. 15 and 16 for the samevalues of temperature and Yp as above. Outside the transitionregion, where homogeneous matter dominates, (HM)p /(HM)approaches the total system proton fraction, Ye. Inside thetransition region, uniform matter presents a proton fractioninferior to the one of the (clusters+uniform matter) system,as expected given the relative proton enrichment manifestedby clusters (right panels of Fig. 10). Moreover, pHM/HMis sensitive to Ye. In the clusterized subsystem, free protonfractions keep the memory of the total system proton fraction.065801-19AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)00.250.50.75110-610-510-410-310-210-1 HM/tot00.10.20.310-610-510-410-310-210-1 HM p/HM00.250.50.75110-410-310-210-1 HM/tot00.10.20.310-410-310-210-1 HM p/HM00.250.50.75110-410-310-210-1tot (fm-3) HM/tot00.10.20.310-410-310-210-1tot (fm-3) HM p/HMFIG. 16. (Color online) The same as in Fig. 15for Ye = 0.3.With regard to the sensitivity on how the homogeneous matteris implemented, the open diamonds show that SKM* andSly230a give comparable results.In what concerns the virial expansion of the equation ofstate as performed in Ref. [78], it is clear that the condition ofsubunitary nucleon and fugacities prevents the explorationof HM/tot and pHM/HM vs. up to the normal nucleardensity. Nevertheless, over the accessible density domainthe behavior of the two quantities clearly reproduces thealready discussed patterns, that is, (i) within a restraineddensity domain, which shrinks with increasing temperature,clusterized matter ( matter, in this case) dominates overhomogeneous matter and (ii) homogeneous matter is depletedin protons. As expected, the evolution with density of the twoquantities is smooth and the lack of nuclei heavier than leadsto less mass bound in clusters with respect to both LS and thepresently proposed model. The relative amount of matterand its isospin symmetry result in homogeneous matter that-10-5051010-510-410-310-210-1Ebaryonic/Ab (MeV)-10-5051010-510-410-310-210-1-10-5051010-410-310-210-1Ebaryonic/Ab (MeV)-10-5051010-410-310-210-10102010-410-310-210-1tot (fm-3)Ebaryonic/Ab (MeV)0102010-410-310-210-1tot (fm-3)FIG. 17. (Color online) Evolution with totalbaryonic density of the baryonic energy per baryon atdifferent temperatures T = 1.6, 5 and 10 MeV forYe = 0.2 (left panels) and Ye = 0.3 (right panels).Open symbols correspond to present model predic-tions; dashed lines stand for Lattimer-Swesty data[74]; dot-dashed lines illustrate Horowitz-Schwenkresults [78]; dotted lines correspond to the casein which only homogeneous nuclear matter wouldexist.065801-20STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)0204010-510-410-310-210-1Etot/Ab (MeV)0102030405010-410-310-210-1Etot/Ab (MeV)102030405010-410-310-210-1tot (fm-3)Etot/Ab (MeV)FIG. 18. (Color online) Evolution with total baryonic density ofthe total energy per baryon at different temperatures T = 1.6, 5, and10 MeV and Ye = 0.2 and 0.3. Open symbols illustrate present modelpredictions while dashed lines stand for Lattimer-Swesty data [74].Large (small) symbols and thick (thin) lines correspond to Ye = 0.2(0.3).has a proton fraction similar to (T = 1.6 MeV), higher than(T = 5 MeV), or between (T = 10 MeV) the values predictedby LS and the present model.3. Equations of stateIn astrophysical simulations, the equations of state, i.e., thefunctional dependences between different state observables,are extremely important ingredients.Figures 17 and 18 present the baryonic energy per baryon,E(bar)/Ab = [E(cl) + E(HM)]/Ab, and, respectively, total en-ergy per baryon, E(tot)/Ab = [E(bar) + E(el) + E( )]/Ab, as afunction of total baryonic density for T = 1.6, 5, and 10 MeVand Yp = 0.2 and 0.3.As one may note, at all considered temperatures theE(bar)/Ab() curves behave similarly: At low densities, theyare flat as dominated by the ideal gas character of dilutedhomogeneous matter, while in the transition density, theymanifest a sharp fall. LS results additionally show that for > 0 E(bar)/Ab rises steeply with . The dependence withisospin asymmetry and temperature is trivial: more symmetricsystems are more bound and average baryon energy increaseswith the temperature. Concerning the comparison between01020304010-510-410-310-210-1Esym./Ab (MeV)01020304010-510-410-310-210-1Esym./Ab (MeV)01020304010-510-410-310-210-1tot (fm-3)Esym./Ab (MeV)FIG. 19. (Color online) Evolution with total baryonic density andtemperature of the symmetry energy per baryon calculated accordingto Eq. (52). See the caption to Fig. 17 for the line and symbol codes.our model and LS, one can say that at T = 1.6 MeV theyagree perfectly. At higher temperatures the two models resultsdeviate in the total baryonic density range corresponding to thehomogeneous-unhomogeneous matter transition. Moreover,the discrepancy augments with the temperature. Taking intoaccount that uniform matter is described by the same equationof state in LS and present model, the only possible explanationof the observed behavior concerns the clusterized component.As we have already discussed that the present model accountsfor a systematically larger mass fraction bound in clusterswith respect with LS, the clusters binding energy is expectedto act in the sense of diminishing E(cl). This is confirmed byinspection of the thin lines in Fig. 17, which represent thebaryonic energy per baryon number, obtained by artificiallyswitching off the cluster component. For T = 10 MeV, thethin line perfectly agrees with LS results, showing that inSNA clusters practically do not contribute to the EOS at hightemperatures.We have already seen (middle panels of Fig. 14) that inthe high-density subregion of the transition density region HSaccounts for a smaller fraction of mass bound in clusters withrespect LS, while in the low-density subregion the oppositeholds. This relation is replicated in the density dependence ofthe baryonic energy per baryon: in the high-density subregion(Eb/Ab)|HS > (Eb/Ab)|LS while in the low-density region065801-21AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)024610-510-410-310-210-1Sbaryonic/Ab024610-510-410-310-210-10246810-410-310-210-1Sbaryonic/Ab0246810-410-310-210-10246810-410-310-210-1tot (fm-3)Sbaryonic/Ab0246810-410-310-210-1tot (fm-3)FIG. 20. (Color online) The same as in Fig. 17but for the baryonic entropy. See the caption toFig. 17 for the line and symbol codes.(Eb/Ab)|HS (Eb/Ab)|LS. Concerning the comparison of HSversus the present model, a similar conclusion may be drawn:as far as our model accounts for a higher percentage of boundmass, it predicts smaller values of Eb/Ab. In addition to this,HS shows less Ye sensitivity than the other two models.As the electron and photon contributions are trivial, it isobvious that all features of the baryonic energy density will bethe same in the representation of the total energy per baryon.Indeed, Fig. 18 shows that at low temperatures and whateverdensities, the present model results agree with the LS ones. Thesame is true for high temperatures and densities correspondingto homogeneous matter. Otherwise, GCA results deviate fromthe SNA ones.Information on the density dependence of the baryonicenergy for different system asymmetry allows one to inferthe density dependence of the symmetry energy, an extremelyimportant quantity which is so far largely unknown. Indeed,recasting the first term in the expansion of the internal energyper baryon in powers of asymmetry = (N Z)/A,Esym(,T )/A = 122(E/A)2=0, (51)in a finite difference formula, one gets,Esym(,T )/A =E(, T , Y 1p)/A E(, T , Y (2)p )/A(1 2Y 1p)2 (1 2Y 2p)2 . (52)Equation (52) is exact because of the quadratic dependenceof the baryonic density with the isospin asymmetry. The modelpredictions on the baryonic density dependence of the symme-try energy per baryon (open symbols) are illustrated in Fig. 19along with the Lattimer-Swesty (dashed lines), Horowitz-Schwenk (dot-dashed lines), and homogeneous matter(dotted lines) results. The common feature of the three modelsis that, when a significant amount of matter is bound in clusters,Esym/A deviates from expectations based on homogeneousmatter behavior. The deviation is visible in diluted matter andis washed out while approaching 0, as nuclei and homo-geneous matter are described by similar EOS. In additionto this remarkable symmetry energy increase, clusterizationinduces also a temperature dependence of Esym/A(), notpresent in homogeneous matter. At the highest temperature,where our model accounts for much more clusters withrespect to LS and HS, the agreement is mainly qualitative. AtT = 1.6 and 5 MeV the results of the present model perfectlyagree with those of LS and both exceed the HS predictions.The same conclusions have been recently pointed out inRef. [50]. There it is shown that clusterization also reducesthe sensitivity of the EOS to the parameters of the effectiveinteraction.Figures 20 and 21 illustrate the total baryonic densitydependence of the baryonic entropy per baryon, S(bar)/Ab =(S(cl) + S(HM))/Ab, and total entropy per baryon, Stot/Ab =[S(bar) + S(el) + S( )]/Ab, for T = 1.6, 5, and 10 MeV andYe = 0.2, 0.3. Both curve families present a monotonicdecrease with , a monotonic increase with T and an isospinasymmetry dependence which vanishes with increasing T .Again, at high temperatures and densities corresponding tothe transition region, the predictions of our model are typicallydifferent than LS results,while the two calculations agree welloutside the transition region and at low temperatures. Thefact that at T = 10 MeV the thin line, corresponding to thecase in which baryonic matter would exclusively consist ofuniform matter, sits perfectly on top of LS reflects the differentroles played by the clusterized component in the two models.The fact that at low temperatures HS results exceed the LS065801-22STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)024610-510-410-310-210-1Stot/Ab0246810-410-310-210-1Stot/Ab024681010-410-310-210-1tot (fm-3)Stot/AbFIG. 21. (Color online) The same as in Fig. 18 but for the totalentropy. See the caption to Fig. 18 for the line and symbol codes.ones may be explained by the much-increased cluster numberwithin the first model (see the right panels of Fig. 14). This dis-crepancy diminishes with increasing temperature such that atT = 10 MeV, HS and LS results are practically identical. As inthe case of energy, the total entropy inherits the characteristicsof the baryonic term. Finally, baryonic and total entropies showless sensitivity to clusters respect to energy.Figures 22 and 23 present the dependence on the totalbaryonic density of the baryonic pressure p(bar) = wclp(cl) +wHMp(HM) and, respectively, total pressure p = p(bar) +p(lattice) + p(el) + p( ) for the same temperatures and protonfractions considered before. For dilute matter p(bar)() is alinear function, as the system recovers the ideal gas limit. Aspectacular feature of the LS results is that, up to a certaintemperature, at a certain value of which depends on thetemperature p(bar)() manifests a sudden fall.This behavior is due to those solutions of the uniform mattercomponent which are unstable and, thus, characterized by neg-ative values ofp(HM). If baryonic matter were exclusively madefrom homogeneous matter, p(bar)(tot)|Ye would correspond tothe dotted lines in Fig. 22. Their shape is qualitatively similarto LS but, except for T = 10 MeV, quantitative discrepanciesexist over the whole density range, as one would have actuallyexpected given the cluster contribution.It is interesting that the falling pattern of HM and LSpressure curves practically does not exist within our model.The explanation relies on the mixture stability criterion wehave adopted in order to maximize the (constrained) entropy(see Sec. II A1). A reminiscence of the fall may, nevertheless,be noticed in the short plateaus or back-bendings. HS resultsplotted with dot-dashed lines in the left panels do not manifestsignatures of negative pressure either, meaning that the nuclearmatter obtained by performing a virial expansion of theequation of state is stable.The disappearance of the sudden fall at high temperatures(in our case T = 10 MeV) may be understood looking againat Fig. 3. This figure shows that, if the n p trajectory10-510-410-310-210-1110-510-410-310-210-1Pbaryonic (MeV/fm3 )10-510-410-310-210-1110-510-410-310-210-110-410-310-210-1110-410-310-210-1Pbaryonic (MeV/fm3 )10-410-310-210-1110-410-310-210-110-310-210-1110-410-310-210-1tot (fm-3)Pbaryonic (MeV/fm3 )10-310-210-1110-410-310-210-1tot (fm-3)FIG. 22. (Color online) The same as in Fig. 17but for the baryonic pressure, pb = wclp(cl) +wHMp(HM). See the caption to Fig. 17 for the lineand symbol codes.065801-23AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)10-510-410-310-210-1110-510-410-310-210-1Ptot (MeV/fm3 )10-310-210-1110-410-310-210-1Ptot (MeV/fm3 )10-210-111010-310-210-1tot (fm-3)Ptot (MeV/fm3 )FIG. 23. (Color online) The same as in Fig. 18 but for thetotal pressure,ptot = wclp(cl) + wHMp(HM) + p(lattice) + p(el) + p( ). Itis remarkable to notice that by summing up the electron and contribution, the total EOS shows almost no sensitivity to the baryonicEOS. See Fig. 18 for the line and symbol code.through the phase diagram fixed by the clusterized system issystematically situated under the high-density spinodal limit,then only one solution for the uniform nuclear matter willexist, namely the low-density stable solution.Figure 23 corresponds to ptot vs. , and shows thatafter considering the electron and photon contributions thediscrepancies originating from the baryonic component of theEOS are to a large extent washed out. The same conclusionwas reached by Refs. [42,51].C. Neutrino opacityAbundant neutrino fluxes are emitted during the collapseof a supernova core and in the first 105 years of life of thesubsequently born neutron star, thus being mainly responsiblefor the star cooling. The mean free path of the neutrinosdepends on the EOS and baryonic composition. During thecollapse phase, it decreases as the star radius shrinks from100 to 10 km, such that it eventually becomes smaller thanthe star radius. This final regime corresponds to the so-calledfull neutrino trapping, while the initial stage corresponds toneutrino transparency and is characterized by a vanishingneutrino chemical potential. Full neutrino trapping is expectedto occur also in postbounce supernovae, when the matter is stilldense. The transition between neutrino trapping and neutrinofree streaming should in principle be a continuous process,obtained by solving the Boltzmann equation for neutrinotransport coupled with the hydrodynamics of the supernovaevolution [82]. In many actual simulations, schematic flux-limiting schemes are adopted in the diffusion approximation[83] allowing regimes of partial neutrino trapping. In thesemodelizations the neutrino diffusion coefficient depends onthe neutrino mean free path, which in turn obviously dependson the matter composition through the different elastic- andinelastic-scattering processes of neutrinos on protons andnuclei [32,82].In Ref. [71] it was additionally suggested that the mattercomposition itself may depend in turn on the percentageof trapped neutrinos, thus introducing a self-consistencyproblem. To explore this issue, we study in this section thecomposition of the baryonic matter in our model as a functionof the trapped neutrino density.The total number of neutrinos produced per unit volumeprod can be calculated from the (local) proton fraction con-sidering hot star matter as produced from the deleptonizationof an initial iron core according to [71]prod = [(ZA)0(ZA)], (53)where (Z/A)0 stands for the Z/A in56Fe.The percentage of trapped neutrinos x depends on time andon the local neutrino mean free path. Determining this quantityconsistently would require a full hydrodynamic calculationwith a complete treatment of the weak processes implied,which is completely out of scope of the present article. To studythe interplay between trapping and matter composition wesimply take the neutrino density at the given thermodynamiccondition as a free parameter, linked to the opacity x byx /prod . The scope of this analysis will then be tocheck whether the matter composition depends on the neutrinodensity, which would imply an extra self-consistent couplingbetween the baryon and the lepton sector. Indeed if the clusterproperties strongly vary during the transition between neutrinotrapping and neutrino streaming, it would be necessary tocouple the neutrino propagation with the calculation of theEOS, which would be a very heavy numerical task.Equation (40) in Sec. IID shows that for a relativistic Fermigas the particle density is univocally linked to its chemicalpotential. In the case of neutrinos, one can safely use thesimpler T = 0 expression of an ultrarelativistic degenerateFermi gas: = hc(62)1/3. (54)The -equilibrium equation (44) with a neutrino chemicalpotential () fixed from Eq. (54) such as to produce thechosen neutrino density = xprod defines a trajectory inthe (n,p) or, equivalently, (n,p) plane. These trajectoriesare shown in Fig. 24 for different values of x , covering thewhole domain between full trapping x = 1 and zero trappingx = 0.065801-24STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)10-510-410-310-210-110-510-410-310-210-1x=1x=0.1x=0.01x=0.001n (fm-3) p (fm-3)10-410-310-210-110-410-310-210-1n (fm-3) p (fm-3)10-310-210-110-310-210-1n (fm-3) p (fm-3)FIG. 24. (Color online) Constant neutrinoopacity (x = 1, 101, 102, 103) paths in thetotal density plane for T = 1.6, 5, and 10 MeVunder the assumption of equilibrium. Neutrinoopacity is calculated according to Eq. (53). Thelines correspond to paths of constant Z/A, whosevalues are indicated on the figure.10 210 310-510-410-310-210-1x=1x=0.1x=0.01x=0.001tot (fm-3) A 0.20.250.30.350.40.450.510-510-410-310-210-1tot (fm-3) Z / A11010 210-410-310-210-1tot (fm-3) A 00.10.20.30.40.50.610-410-310-210-1tot (fm-3) Z / A11010-310-210-1tot (fm-3) A 00.10.20.30.40.50.610-310-210-1tot (fm-3) Z / AFIG. 25. (Color online) Average cluster size (left panels) andaverage cluster isospin (right panels) as a function of total baryonicdensity for T = 1.6, 5, and 10 MeV along constant neutrinoopacity (x = 1, 101, 102, 103) paths under the assumption of equilibrium. No selection is performed with respect to the mixturestability.As one may note, the constant x paths roughly corre-spond to constant Z/A paths (lines on the figure). More-over, the global isospin asymmetry of the matter increasesmonotonically with the percentage of trapped neutrinos andno significant dependence on temperature is observed. Forinstance, over the whole considered temperature domain thefull neutrino trapping path lies along the Z/A = 0.1 path,while x = 103 corresponds to Z/A slightly higher than0.4. These results are in good agreement with Ref. [71],where only homogeneous nuclear matter was consideredand no clusters were included. The extra information borneby the present work is that we can correlate the percent-age of trapped neutrinos with the characteristics of theclusters.Figure 25 displays the average cluster size (left panels) andisospin composition (right panels) as a function of total bary-onic density at different temperatures along constant neutrinoopacity paths. For the sake of completeness, no selection isnow made with respect to the stability of the mixture. Wecan see that the size strongly depends on the thermodynamicconditions (temperature and density) but, in the cases in whichhomogeneous matter is stable, it is virtually independent ofthe opacity to neutrinos. This result can be understood fromthe fact that changing the percentage of trapped neutrinosat equilibrium essentially amounts to changing the matterisospin asymmetry Yp for a given total density (see Fig. 24).In turn, this affects the chemical composition of clustersbut not the cluster size, as we have already observed indiscussing Fig. 10. Since the interaction probability of neu-trinos with matter essentially depends on isoscalar quantitiesas the average cluster density and size [29,32], our resultindicates that the problem of neutrino-baryon interaction canbe effectively decoupled from the problem of the mattercomposition.065801-25AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)IV. CONCLUSIONSIn this article we have presented a phenomenologicalmodel for stellar matter at finite temperature and subsatu-ration densities which describes matter inhomogeneities asa continuous mixture of a distribution of loosely interactingclusters in statistical equilibrium, with free nucleons treatedwithin the nonrelativistic mean-field approximation. The twocomponents interact through the electrostatic energy in theWigner-Seitz approximation and a configuration-dependentexcluded volume term in the spirit of the Van-der-Waals gas.Such a mixture is demonstrated to minimize (maximize) theassociated energetic (entropic) thermodynamic potential in thewhole considered range of densities and temperatures such thata stable solution can always be found. As a consequence, thetransition from the homogeneous core to the inhomogeneouscrust of a finite temperature neutron star is naturally obtainedin the model, at variance with other approaches wherediscontinuous transitions have to be invoked.In qualitative agreement with previous works based on aNSE, we find that the inclusion of a statistical distribution ofclusters in low-density stellar matter modifies in an importantway the average matter composition respect to standardtreatments as the LS equation of state. For densities close to thetransition density which are not accessible to NSE models ina thermodynamically consistent way, we show that clusteringhas also sizable consequences of thermodynamic quantitiesand equations of state. In particular, both the symmetryenergy and the baryonic energy are considerably altered bythe presence of clusters, in agreement with the microscopicresults of Ref. [50]. These results open different perspectivesfor future work which will be pursued in the near future.Our formalism where the cluster partition sum is sampledwithin an exact, though numerical, Monte Carlo procedureallows us to account for cluster correlations. We have alreadyshown in a previous article [64] that going beyond thenoninteracting NSE approximation for cluster equilibriumhas an important effect on the cluster phase diagram, andin the present work the configuration-dependent excludedvolume correction, which phenomenologically accounts forthe nuclear interaction, was shown to be responsible for thecrust-core transition. The Coulomb part of the interaction,which is presently treated in the usual one-body Wigner-Seitzapproximation, could also be calculated exactly accountingfor charge fluctuations inside the Wigner-Seitz cell. From theexisting literature [16,84] we expect that properly accountingfor Coulomb correlations induced by the charge fluctuationsat the microscopic level may have an important influence onmatter composition and thermodynamic features, particularlyat high density and temperature where the configurationfluctuations are maximal.In the regime of baryonic density close to saturation, weobserve metastable and unstable solutions corresponding tobubblelike structures where small spherical clusters coexistwith slightly diluted homogeneous matter. The possible per-sistence at high temperature of pasta-like structures, welldocumented at low temperatures [1214], is still an openquestion and it will be interesting to see if such exotic structuresappear if deformation degrees of freedom are allowed in theconfiguration geometry and cluster energy functional.Concerning the cluster energy functional, in this ex-ploratory work we have employed a very simple and simplisticliquid-dropbased expression. It will, however, be very easyto implement a composite functional giving the experimen-tal masses at zero temperature and vanishing density andevolving toward a parametric formulation at high density andtemperature and for unknown neutron-rich nuclear species.Concerning the homogeneous matter functional, we plan toinclude pairing correlations within the finite temperature HFBapproach as a natural extension of the present formalism,which will allow having a trustworthy equation of state atthe lowest temperatures.Finally, once the different physical inputs are settled, alonger-range project consists of constructing a complete EOStable for direct implementation in supernova and neutron starscodes in the framework of the COMPSTAR project [85].ACKNOWLEDGMENTSA.R.R acknowledges partial support from the RomanianNational Authority for Scientific Research under Grant No.IDEI nr. 267/2007 and PN 09370105/2009 and kind hospitalityfrom LPC-Caen within IFIN-IN2P3 agreement No. 07-44.F.G. acknowledges partial support from ANR under ProjectNExEN.[1] R. Knorren, M. Prakash, and P. J. Ellis, Phys. Rev. C 52, 3470(1995).[2] T. Takatsuka, S. Nishizaki, Y. Yamamoto, and R. Tamagaki,Prog. Theor. Phys. 115, 355 (2006).[3] D. Page, M. Prakash, J. M. Lattimer, and A. W. Steiner, Phys.Rev. Lett. 85, 2048 (2000).[4] D. Blaschke, T. Klahn, and D. N. Voskresensky, Astrophys. J.533, 406 (2000).[5] P. K. Panda, D. P. Menezes, and C. Providencia, Phys. Rev. C69, 058801 (2004).[6] O. E. Nicotra, M. Baldo, G. F. Burgio, and H. J. Schulze, Phys.Rev. D 74, 123001 (2006).[7] K. Nakazato, K. Sumiyoshi, and S. Yamada, Phys. Rev. D 77,103006 (2008).[8] N. Yasutake, T. Maruyama, and T. Tatsumi, Phys. Rev. D 80,123009 (2009).[9] J. M. Lattimer and M. Prakash, Science 304, 536(2004).[10] P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, Neutron Stars:Equation of State and Structure (Springer, Berlin, 2007).[11] C. J. Horowitz, M. A. Perez-Garcia, D. K. Berry, andJ. Piekarewicz, Phys. Rev. C 72, 035801 (2005).[12] H. Sonoda, G. Watanabe, K. Sato, K. Yasuoka, and T. Ebisuzaki,Phys. Rev. C 77, 035806 (2008); G. Watanabe, H. Sonoda,065801-26http://dx.doi.org/10.1103/PhysRevC.52.3470http://dx.doi.org/10.1103/PhysRevC.52.3470http://dx.doi.org/10.1143/PTP.115.355http://dx.doi.org/10.1103/PhysRevLett.85.2048http://dx.doi.org/10.1103/PhysRevLett.85.2048http://dx.doi.org/10.1086/308664http://dx.doi.org/10.1086/308664http://dx.doi.org/10.1103/PhysRevC.69.058801http://dx.doi.org/10.1103/PhysRevC.69.058801http://dx.doi.org/10.1103/PhysRevD.74.123001http://dx.doi.org/10.1103/PhysRevD.74.123001http://dx.doi.org/10.1103/PhysRevD.77.103006http://dx.doi.org/10.1103/PhysRevD.77.103006http://dx.doi.org/10.1103/PhysRevD.80.123009http://dx.doi.org/10.1103/PhysRevD.80.123009http://dx.doi.org/10.1126/science.1090720http://dx.doi.org/10.1126/science.1090720http://dx.doi.org/10.1103/PhysRevC.72.035801http://dx.doi.org/10.1103/PhysRevC.77.035806STATISTICAL DESCRIPTION OF COMPLEX NUCLEAR . . . PHYSICAL REVIEW C 82, 065801 (2010)T. Maruyama, K. Sato, K. Yasuoka, and T. Ebisuzaki, Phys.Rev. Lett. 103, 121101 (2009).[13] W. G. Newton and J. R. Stone, Phys. Rev. C 79, 055801(2009).[14] F. Sebille, S. Figerou, and V. de la Mota, Nucl. Phys. A 822, 51(2009).[15] M. Grousson, G. Tarjus, and P. Viot, Phys. Rev. E 62, 7781(2000); 64, 036109 (2001).[16] P. Napolitani, Ph. Chomaz, F. Gulminelli, and K. H. O. Hasnaoui,Phys. Rev. Lett. 98, 131102 (2007); C. Ducoin, K. H. O.Hasnaoui, P. Napolitani, Ph. Chomaz, and F. Gulminelli, Phys.Rev. C 75, 065805 (2007).[17] J. M. Lattimer, C. J. Pethick, D. G. Ravenhall, and D. Q. Lamb,Nucl. Phys. A 432, 646 (1985).[18] M. Lassaut, H. Flocard, P. Bonche, P. H. Heenen, and E. Suraud,Astron. Astrophys. 183, L3 (1987).[19] C. J. Pethick, D. G. Ravenhall, and C. R. Lorenz, Nucl. Phys. A584, 675 (1995).[20] F. Douchin and P. Haensel, Phys. Lett. B 485, 107 (2000).[21] G. Watanabe, K. Iida, and K. Sato, Nucl. Phys. A 726, 357(2003).[22] T. Maruyama, T. Tatsumi, D. N. Voskresensky, T. Tanigawa, andS. Chiba, Phys. Rev. C 72, 015802 (2005).[23] J. M. Lattimer and F. Douglas Swesty, Nucl. Phys. A 535, 331(1991).[24] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, Nucl.Phys. A 637, 435 (1998); Prog. Theor. Phys. 100, 1013(1998).[25] G. Shen, C. J. Horowitz, and S. Teige, Phys. Rev. C 82, 015806(2010).[26] A. Burrows and J. M. Lattimer, Astrophys. J. 307, 178(1986).[27] H. A. Bethe, G. E. Brown, J. Applegate, and J. M. Lattimer,Nucl. Phys. A 324, 487 (1979).[28] Ya. B. Zeldovich and I. D. Novikov, Relativistic Astrophysics(University of Chicago Press, Chicago, 1971).[29] C. J. Horowitz, M. A. Perez-Garcia, J. Carriere, D. K. Berry,and J. Piekarewicz, Phys. Rev. C 70, 065806 (2004).[30] G. Martnez-Pinedo, M. Liebendorfer, and D. Frekers, Nucl.Phys. A 777, 395 (2006).[31] H. T. Janka, K. Langanke, A. Marek, G. Martnez-Pinedo, andB. Muller, Phys. Rep. 442, 38 (2007).[32] H. Sonoda, G. Watanabe, K. Sato, T. Takiwaki, K. Yasuoka, andT. Ebisuzaki, Phys. Rev. C 75, 042801(R) (2007).[33] W. R. Hix, O. E. B. Messer, A. Mezzacappa, M. Liebendorfer,J. Sampaio, K. Langanke, D. J. Dean, and G. Martnez-Pinedo,Phys. Rev. Lett. 91, 201102 (2003).[34] K. Langanke and G. Martnez-Pinedo, Nucl. Phys. A 731, 365(2004).[35] J. J. Cowan, F. K. Thieleman, and J. W. Truran, Phys. Rep. 208,267 (1991).[36] Y. Z. Qian, Prog. Part. Nucl. Phys. 50, 153 (2003).[37] J. Margueron, J. Navarro, and P. Blottiau, Phys. Rev. C 70,028801 (2004).[38] O. L. Caballero, C. J. Horowitz, and D. K. Berry, Phys. Rev. C74, 065801 (2006).[39] D. Page, U. Geppert, and F. Weber, Nucl. Phys. A 777, 497(2006).[40] D. G. Yakovlev, O. Y. Gnedin, M. E. Gusakov, A. D. Kaminker,K. P. Levenfish, and A. Y. Potekhin, Nucl. Phys. A 752, 590(2005).[41] A. C. Phillips, The Physics of Stars (John Wiley & Sons,Chichester, 1994).[42] A. S. Botvina and I. N. Mishustin, Nucl. Phys. A 843, 98 (2010).[43] S. I. Blinnikov, I. V. Panov, M. A. Rudzsky, and K. Sumiyoshi,arXiv:0904.3849 [astro-ph].[44] S. R. Souza, A. W. Steiner, W. G. Lynch, R. Donangelo, andM. A. Famiano, Astrophys. J. 707, 1495 (2009).[45] M. E. Fisher, Physics (NY) 3, 255 (1967).[46] J. Negele and D. Vautherin, Nucl. Phys. A 207, 298 (1973).[47] S. K. Samaddar, J. N. De, X. Vinas, and M. Centelles, Phys.Rev. C 80, 035803 (2009).[48] A. Arcones, G. Martnez-Pinedo, E. OConnor, A. Schwenk,H. Th. Janka, C. J. Horowitz, and K. Langanke, Phys. Rev. C78, 015806 (2008).[49] S. Heckel, P. P. Schneider, and A. Sedrakian, Phys. Rev. C 80,015805 (2009).[50] S. Typel, G. Ropke, T. Klahn, D. Blaschke, and H. H. Wolter,Phys. Rev. C 81, 015803 (2010).[51] M. Hempel and J. Schaffner-Bielich, Nucl. Phys. A 837, 210(2010).[52] D. Vautherin, Adv. Nucl. Phys. 22, 123 (1996).[53] C. Ducoin, Ph. Chomaz, and F. Gulminelli, Nucl. Phys. A 771,68 (2006).[54] M. Barranco and J. R. Buchler, Phys. Rev. C 24, 1191 (1981).[55] S. Shlomo and V. M. Kolomietz, Rep. Prog. Phys. 68, 1 (2005).[56] F. Douchin, P. Haensel, and J. Meyer, Nucl. Phys. A 665, 419(2000).[57] A. Rios, Nucl. Phys. A 845, 58 (2010).[58] J. Bartel, P. Quentin, M. Brack, C. Guet, and H. B. Hakansson,Nucl. Phys. A 386, 79 (1982).[59] J. P. Bondorf, A. S. Botvina, A. S. Iljinov, I. N. Mishustin, andK. Sneppen, Phys. Rep. 257, 133 (1995).[60] D. H. E. Gross, Rep. Progr. Phys. 53, 605 (1990).[61] S. E. Koonin and J. Randrup, Nucl. Phys. A 474, 173 (1987).[62] A. H. Raduta and A. R. Raduta, Phys. Rev. C 55, 1344 (1997);65, 054610 (2002).[63] G. Audi and A. H. Wapstra, Nucl. Phys. A 595, 409 (1995).[64] Ad. R. Raduta and F. Gulminelli, Phys. Rev. C 80, 024606(2009).[65] A. S. Iljinov, M. V. Mebel, N. Bianchi, E. de Sanctis,C. Guaraldo, V. Lucherini, V. Muccifora, E. Polli, A. R. Reolon,and P. Rossi, Nucl. Phys. A 543, 517 (1992).[66] C. Ishizuka, A. Ohnishi, and K. Sumiyoshi, Prog. Theor. Phys.Suppl. 146, 373 (2002).[67] D. H. E. Gross, Microcanonical Thermodynamics: Phase Transi-tions in Finite Systems, Lecture Notes in Physics, vol. 66 (WorldScientific, Singapore, 2001).[68] N. K. Glendenning, Phys. Rep. 342, 394 (2001).[69] M. Hempel, G. Pagliara, and J. Schaffner-Bielich, Phys. Rev. D80, 125014 (2009).[70] J. Margueron and A. Fantina (private communication).[71] C. Ducoin, Ph. Chomaz, and F. Gulminelli, Nucl. Phys. A 789,403 (2007).[72] F. S. Kitaura, H. Th. Janka, and W. Hillebrandt, Astron.Astrophys. 450, 345 (2006).[73] A. Marek and H. Th. Janka, Astrophys. J. 694, 664 (2009).[74] J. M. Lattimer [http://www.astro.sunysb.edu/lattimer/EOS/key.txt].[75] M. DAgostino et al., Phys. Lett. B 473, 219 (2000).[76] F. Gulminelli and M. DAgostino, Eur. Phys. J. A 30, 253 (2006).[77] E. Bonnet et al., Phys. Rev. Lett. 103, 072701 (2009).065801-27http://dx.doi.org/10.1103/PhysRevLett.103.121101http://dx.doi.org/10.1103/PhysRevLett.103.121101http://dx.doi.org/10.1103/PhysRevC.79.055801http://dx.doi.org/10.1103/PhysRevC.79.055801http://dx.doi.org/10.1016/j.nuclphysa.2009.02.013http://dx.doi.org/10.1016/j.nuclphysa.2009.02.013http://dx.doi.org/10.1103/PhysRevE.62.7781http://dx.doi.org/10.1103/PhysRevE.62.7781http://dx.doi.org/10.1103/PhysRevE.64.036109http://dx.doi.org/10.1103/PhysRevLett.98.131102http://dx.doi.org/10.1103/PhysRevC.75.065805http://dx.doi.org/10.1103/PhysRevC.75.065805http://dx.doi.org/10.1016/0375-9474(85)90006-5http://dx.doi.org/10.1016/0375-9474(94)00506-Ihttp://dx.doi.org/10.1016/0375-9474(94)00506-Ihttp://dx.doi.org/10.1016/S0370-2693(00)00672-9http://dx.doi.org/10.1016/S0375-9474(03)01601-4http://dx.doi.org/10.1016/S0375-9474(03)01601-4http://dx.doi.org/10.1103/PhysRevC.72.015802http://dx.doi.org/10.1016/0375-9474(91)90452-Chttp://dx.doi.org/10.1016/0375-9474(91)90452-Chttp://dx.doi.org/10.1016/S0375-9474(98)00236-Xhttp://dx.doi.org/10.1016/S0375-9474(98)00236-Xhttp://dx.doi.org/10.1143/PTP.100.1013http://dx.doi.org/10.1143/PTP.100.1013http://dx.doi.org/10.1103/PhysRevC.82.015806http://dx.doi.org/10.1103/PhysRevC.82.015806http://dx.doi.org/10.1086/164405http://dx.doi.org/10.1086/164405http://dx.doi.org/10.1016/0375-9474(79)90596-7http://dx.doi.org/10.1103/PhysRevC.70.065806http://dx.doi.org/10.1016/j.nuclphysa.2006.02.014http://dx.doi.org/10.1016/j.nuclphysa.2006.02.014http://dx.doi.org/10.1016/j.physrep.2007.02.002http://dx.doi.org/10.1103/PhysRevC.75.042801http://dx.doi.org/10.1103/PhysRevLett.91.201102http://dx.doi.org/10.1016/j.nuclphysa.2003.11.049http://dx.doi.org/10.1016/j.nuclphysa.2003.11.049http://dx.doi.org/10.1016/0370-1573(91)90070-3http://dx.doi.org/10.1016/0370-1573(91)90070-3http://dx.doi.org/10.1016/S0146-6410(02)00178-3http://dx.doi.org/10.1103/PhysRevC.70.028801http://dx.doi.org/10.1103/PhysRevC.70.028801http://dx.doi.org/10.1103/PhysRevC.74.065801http://dx.doi.org/10.1103/PhysRevC.74.065801http://dx.doi.org/10.1016/j.nuclphysa.2005.09.019http://dx.doi.org/10.1016/j.nuclphysa.2005.09.019http://dx.doi.org/10.1016/j.nuclphysa.2005.02.061http://dx.doi.org/10.1016/j.nuclphysa.2005.02.061http://dx.doi.org/10.1016/j.nuclphysa.2010.05.052http://arXiv.org/abs/arXiv:0904.3849http://dx.doi.org/10.1088/0004-637X/707/2/1495http://dx.doi.org/10.1016/0375-9474(73)90349-7http://dx.doi.org/10.1103/PhysRevC.80.035803http://dx.doi.org/10.1103/PhysRevC.80.035803http://dx.doi.org/10.1103/PhysRevC.78.015806http://dx.doi.org/10.1103/PhysRevC.78.015806http://dx.doi.org/10.1103/PhysRevC.80.015805http://dx.doi.org/10.1103/PhysRevC.80.015805http://dx.doi.org/10.1103/PhysRevC.81.015803http://dx.doi.org/10.1016/j.nuclphysa.2010.02.010http://dx.doi.org/10.1016/j.nuclphysa.2010.02.010http://dx.doi.org/10.1007/0-306-47065-9_4http://dx.doi.org/10.1016/j.nuclphysa.2006.03.005http://dx.doi.org/10.1016/j.nuclphysa.2006.03.005http://dx.doi.org/10.1103/PhysRevC.24.1191http://dx.doi.org/10.1088/0034-4885/68/1/R01http://dx.doi.org/10.1016/S0375-9474(99)00397-8http://dx.doi.org/10.1016/S0375-9474(99)00397-8http://dx.doi.org/10.1016/j.nuclphysa.2010.05.057http://dx.doi.org/10.1016/0375-9474(82)90403-1http://dx.doi.org/10.1016/0370-1573(94)00097-Mhttp://dx.doi.org/10.1088/0034-4885/53/5/003http://dx.doi.org/10.1016/0375-9474(87)90199-0http://dx.doi.org/10.1103/PhysRevC.55.1344http://dx.doi.org/10.1103/PhysRevC.65.054610http://dx.doi.org/10.1016/0375-9474(95)00445-9http://dx.doi.org/10.1103/PhysRevC.80.024606http://dx.doi.org/10.1103/PhysRevC.80.024606http://dx.doi.org/10.1016/0375-9474(92)90278-Rhttp://dx.doi.org/10.1143/PTPS.146.373http://dx.doi.org/10.1143/PTPS.146.373http://dx.doi.org/10.1016/S0370-1573(00)00080-6http://dx.doi.org/10.1103/PhysRevD.80.125014http://dx.doi.org/10.1103/PhysRevD.80.125014http://dx.doi.org/10.1016/j.nuclphysa.2007.03.006http://dx.doi.org/10.1016/j.nuclphysa.2007.03.006http://dx.doi.org/10.1051/0004-6361:20054703http://dx.doi.org/10.1051/0004-6361:20054703http://dx.doi.org/10.1088/0004-637X/694/1/664http://www.astro.sunysb.edu/lattimer/EOS/key.txthttp://www.astro.sunysb.edu/lattimer/EOS/key.txthttp://dx.doi.org/10.1016/S0370-2693(99)01486-0http://dx.doi.org/10.1140/epja/i2006-10121-xhttp://dx.doi.org/10.1103/PhysRevLett.103.072701AD. R. RADUTA AND F. GULMINELLI PHYSICAL REVIEW C 82, 065801 (2010)[78] C. J. Horowitz and A. Schwenk, Nucl. Phys. A 776, 55(2006).[79] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer,Nucl. Phys. A 627, 710 (1997).[80] C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647(2001).[81] I. Vidana, C. Providencia, A. Polls, and A. Rios, Phys. Rev. C80, 045806 (2009).[82] M. Liebendorfer, M. Rampp, H. Th. Janka, and A. Mezzacappa,Astrophys. J. 620, 840 (2005).[83] R. B. Bowers and J. R. Wilson, Astrophys. J. 263, 366(1982).[84] G. Watanabe and K. Iida, Phys. Rev. C 68, 045801 (2003).[85] COMPSTAR, European program 06-RNP-106 of the EuropeanScience Foundation, EOS construction task coordinated byM. Baldo, Y. Illarionov, and M. Oertel.065801-28http://dx.doi.org/10.1016/j.nuclphysa.2006.05.009http://dx.doi.org/10.1016/j.nuclphysa.2006.05.009http://dx.doi.org/10.1016/S0375-9474(97)00596-4http://dx.doi.org/10.1103/PhysRevLett.86.5647http://dx.doi.org/10.1103/PhysRevLett.86.5647http://dx.doi.org/10.1103/PhysRevC.80.045806http://dx.doi.org/10.1103/PhysRevC.80.045806http://dx.doi.org/10.1086/427203http://dx.doi.org/10.1086/160510http://dx.doi.org/10.1086/160510http://dx.doi.org/10.1103/PhysRevC.68.045801