C Ducoin

Ph Chomaz

F Gulminelli

Isospin-dependent clusterization of Neutron-Star Matter

Keywords:

Because of the presence of a liquid-gas phase transition in nuclear matter, compactstar matter can present a region of instability against the formation of clusters. We investigate this phase separation in a matter composed of neutrons, protons and electrons, within a Skyrme-Lyon mean-field approach. Matter instability and phase properties are characterized through the study of the free-energy curvature. The effect of β-equilibrium is also analyzed in detail, and we show that the opacity to neutrinos has an influence on the presence of clusterized matter in finite-temperature proto-neutron stars.

I. INTRODUCTION

Nuclear matter is known to present a phase transition of the liquid-gas type, which occurs at sub-saturation density up to a critical temperature [1,2,3]. In laboratory experiments, this phase diagram can be explored only through heavy-ion collisions; in particular the multifragmentation phenomenon, observed for collisions around the Fermi energy and above, has been interpreted as the occurrence of the nuclear liquid-gas phase transition rounded by finite-size effects [4,5,6,7,8]. From the theoretical point of view, the isospin properties of the phase diagram have been studied extensively [9,10,11] and we know that the phase transition concerns isospin-symmetric (ρ n = ρ p ) as well as asymmetric (ρ n = ρ p ) matter.

Then the transition phenomenology may also affect the physics of compact stars, which are made of dense, neutron-rich matter [12,13]. In particular neutron-star properties may be influenced by the presence of exotic "pasta phases" in their inner crust [14,15,16,17].

Density fluctuations due to the liquid-gas instability may also play a role in the dynamics of supernova explosions, by affecting neutrino propagation [18,19].

In this paper, we study the instability of star matter against the liquid-gas phase transition. At the densities of interest for our study, it can be considered as a medium composed of neutrons, protons and electrons, under the constraint of strict electroneutrality. We will describe nuclear matter (neutrons and protons, denoted "NM") in a mean-field approach, using different Skyrme effective interactions. The Skyrme-Lyon parameterization Sly230a, already employed for astrophysical applications in ref. [20], is taken as a reference, since it has been especially fitted to describe neutron-rich matter [21,22]. It is compared to two other Skyrme forces : the earlier SGII [23], whose parameters are optimized to describe nuclear collective motion, and SIII, one of the first and more simple Skyrme parameterizations [24]. Fig. 1 shows the density behavior of these different forces in the isoscalar (upper part) as well as isovector (lower part) channel.

We can see that realistic interactions, fitted on measured nuclear properties as Sly230a and SGII, are perfectly superposed for symmetric nuclear matter, while SIII presents an abnormally large incompressibility, given by the curvature at saturation :

K ∞ = 9ρ 2 0 ∂ 2 E/A ∂ρ 2 ρ=ρ 0 .
The density dependence of the symmetry energy is less well constrained by experimental data [10], and shows very different behaviors with the different interactions.

Electrons are described as a Fermi gas, whose density is set equal to proton density.

Nuclear matter with electrons (representing star matter) will be denoted "NMe".

In section II, we present nuclear-matter thermodynamics and compare the spinodalinstability region for NM and NMe. In section III, we consider the instability of homogeneous NMe against finite-size density fluctuations, leading to matter clusterization. Some qualitative features of the clusters are deduced from the free-energy curvature properties of homogeneous matter in the framework of a spinodal-decomposition scenario, where cluster characteristics are determined by the dominant fluctuation mode at the onset of phase separation. Finally, in section IV we investigate the effect of neutrinos on star-matter instability when β-equilibrium is assumed. , with I = ρ 3 /ρ.

E/A (MeV)

II. GLOBAL INSTABILITY OF HOMOGENEOUS MATTER

In this section, we establish the thermodynamic instability properties of infinite and homogeneous matter. For a fixed temperature, the free-energy curvature in the density plane (ρ n , ρ p ) determines the spinodal region. Inside the spinodal, a thermodynamic fluctuation of the homogeneous system can lead to a reduction of the free energy : this fluctuation is then amplified, leading to phase separation. Since the resulting phases are both homogeneous and infinite, any even infinitesimal charge density would lead to a divergent Coulomb energy:

the thermodynamic limit thus imposes that each phase is exactly neutral [25].

In the actual cooling of a proto-neutron star, the characteristic time scales may be long enough respect to the typical relaxation times associated with the strong interaction, for the system to be in thermodynamic equilibrium at all steps of its evolution. If this is the case, the spinodal zone is not directly pertinent for star-matter thermodynamics. However the spinodal region contains points of locally inversed curvature, and is always by construction included inside the region of globally inversed curvature, namely the coexistence zone [9,11].

Therefore our results on the thermodynamic instabilities will keep being relevant for the equilibrium mixed phase.

We will first evaluate the spinodal region in nuclear matter, in the absence of Coulomb interaction. Then, electroneutrality will be insured by the presence of electrons, which results in a strong quenching of instability.

A. Nuclear-matter spinodal region

Nuclear matter involves an infinite quantity of nucleons interacting through the nuclear force. The Coulomb interaction is not taken into account in this model, but neutrons and protons are distinguished by their isospin.

In a mean-field approach, each nucleon is an independent particle evolving in a potential resulting from the averaged effect of other nucleons. Let us denote Ĥ the hamiltonian of the system, and ρq the density matrix of each nucleon species, the index q representing neutrons (q = n) or protons (q = p). The mean field Ŵq is defined by the relation δ Ĥ = T r( Ŵq δ ρq ).

Within this approximation, the average energy density of homogeneous nuclear matter H h = Ĥ /V is a functional of one-body densities only: particle densities ρ q and kinetic densities τ q = p2 q / 2 . It is useful to introduce isoscalar and isovector densities :

ρ = ρ n + ρ p , τ = τ n + τ p , ρ 3 = ρ n -ρ p , τ 3 = τ n -τ p . (1)
Using a Skyrme effective interaction, the energy density of homogeneous, spin-saturated matter with no Coulomb interaction, can be written as :

H = K + H 0 + H 3 + H ef f (2)
In this expression, K is the kinetic-energy term, H 0 a density-independent two-body term, H 3 a density-dependent term, and H ef f a momentum-dependent term:

K = 2 2m τ H 0 = C 0 ρ 2 + D 0 ρ 2 3 H 3 = C 3 ρ σ+2 + D 3 ρ σ ρ 2 3 H ef f = C ef f ρτ + D ef f ρ 3 τ 3 (3)
where the coefficients C i an D i are linear combinations of the standard Skyrme parameters [11,21]. The mean field is given by the expression :

Ŵq = ∂H h ∂τ q p2 2 + ∂H h ∂ρ q = 1 2m * q p2 + U q , (4)
where m * q is the effective mass of particles of type q. The individual energy levels read

ǫ i q = p 2 i 2m * q + U q .
In the grand-canonical formalism, they are occupied according to Fermi-Dirac statistics with occupation numbers:

n q (p) = 1/ 1 + exp(β(p 2 /2m * q + U q -µ q )) , (5)
where β is linked to the temperature by β = 1/T , and µ q is the chemical potential of particles of type q. All one-body densities are then given by Fermi integrals and the meanfield grand-canonical partition sum can be expressed in the same way as for an ideal gas:

ln Z 0 V = 2 q ∞ 0 ln(1 + e -β( p 2 2m * q -µ ′ q ) ) 4πp 2 h 3 dp = q 2 3m * q βτ q , (6)
where µ ′ q = µ q -U q . For fixed values of ρ n and ρ p at a given temperature, the value of chemical potentials as well as kinetic densities are obtained by iteratively solving the self-consistent relations between ρ q , µ q and m * q . We can then calculate the entropy as a function of densities. Indeed the mean-field approach approximates the real entropy S by the corresponding mean-field entropy S 0 [26] :

s = S/V ≃ S 0 /V = ln Z 0 /V + β( Ŵ 0 /V -µ n ρ n -µ p ρ p ). ( 7
)
For a given temperature T = β -1 , statistical equilibrium corresponds to the minimization of the free energy F = V f , related to the entropy S by a Legendre transform. Per unit volume, this relation is f = H h -s/β and f is a function of nucleon densities ρ n and ρ p . If the free-energy surface presents a concavity, it can be minimized by linear interpolation (phase mixing). The equilibrium configuration corresponds to a phase coexistence determined by Gibbs contruction, which is a several-dimension generalization of the well-known Maxwell construction [9].

We now turn to spinodal instabiliy. Let us consider an homogeneous system submitted to a thermodynamic fluctuation, consisting in a separation into two portions distinguished by an infinitesimal density shift. Both portions correspond to infinite and homogeneous systems. If the fluctuation leads to a reduction of the global free energy, it will be amplified in the successive time evolution: the initial system is locally unstable. The spinodal region is then defined as the ensemble of unstable points, i.e. points for which the free-energy surface is convex in at least one direction. This is determined by the curvature matrix :

C h N M =   ∂ 2 f /∂ρ 2 n ∂ 2 f /∂ρ n ∂ρ p ∂ 2 f /∂ρ p ∂ρ n ∂ 2 f /∂ρ 2 p   =   ∂µ n /∂ρ n ∂µ n /∂ρ p ∂µ p /∂ρ n ∂µ p /∂ρ p   (8)
This matrix is defined for each point (ρ n , ρ p ), and its lower eigen-value C < gives the minimal free-energy curvature in this point. The eigen-vector u < associated with the minimal curvature, named "instability direction", gives the direction of phase separation.

We have represented on the left part of fig. 2 the minimal curvature C < as a function of total density, for different values of the proton fraction Z/A, at a fixed temperature of 10

MeV : some of the curves have a negative region. On the right part, the minimal value of C < obtained for each asymmetry is represented as a function of Z/A. This shows that, for a finite range of asymmetry, homogeneous matter presents a spinodal instability, i.e. a negative curvature of the free energy. This is observed until a critical temperature of about 14 MeV.

B. Electron quenching of instability

Let us now take the electric charge into account. Matter must be strictly neutral at a macroscopic scale, in order to avoid the divergence of Coulomb energy. In this paper, we consider star matter where the proton charge is neutralized by a Fermi gas of electrons.

Studying thermodynamic fluctuations, each portion in which the system separates is submitted to the constraint of strict electroneutrality : the electron density is then everywhere exactly equal to the proton density, ρ e = ρ p . In this case, electrons do not add any degree of freedom, and the free-energy curvature matrix remains two-dimensional :

C h N M e =  
The presence of electrons contributes to this matrix via the electron susceptibility χ -1 e = ∂µ e /∂ρ e . This term is readily evaluated in the approximation of a degenerated, ultrarelativistic Fermi gas. Because of the low electron mass, these approximations are valid until very low densities. For a degenerated Fermi gas, the Fermi momentum p Fe is given by the relation :

ρ e = ρ p = (p Fe c) 3 3π 2 ( c) 3 . (10)
In the ultra-relativistic approximation, the chemical potential is equal to the Fermi momentum µ e = p Fe c, so it is related to the density by :

µ e = c(3π 2 ρ e ) 1/3 , (11)
which gives the electron contribution to the curvature matrix :

∂µ e /∂ρ e = χ -1 e = c( π 2 9ρ 2 e ) 1/3. (12)
In the density region of interest for our study, the numerical value of χ -1 e is large in comparison with the other terms of the free-energy curvature matrix [START_REF]For example for symmetric matter at saturation density, the susceptibility χ -1 e is 1052.45M eV , while the nuclear free-energy second derivatives (∂µ q /∂ρ q and ∂µ[END_REF]. As a result, adding an electron gas leads to a drastic effect on the eigen-modes of the matrix. Because of its positive sign, it leads to a quenching of the instability. This is illustrated on fig. 2 The value of C min < depends on the details of the nuclear interaction used for the calculation of neutron and proton chemical potentials. For one given point, the sign of C < is given by

det C h N M e = det C h N M + ∂ ρn µ n * ∂ ρe µ e .
It can be negative if det C h N M has a sufficiently strong negative value, ie nuclear-matter thermodynamic instability is strong enough to balance the positive term due to the electron gas. This is does not happen for any of the three Skyrme forces we have considered, as shown by figure 3 : the instability is totally suppressed by electron incompressibility. However, this result is model dependent. For example, in reference [27] the spinodal instability is studied in a relativistic-mean-field approach, using the NL3 set of parameters : the authors obtain a small thermodynamic instability region in the case of nuclear matter with electrons. This region, located at low density, is very reduced with respect to nuclear-matter spinodal, as expected from electron quenching.

III. CLUSTERIZATION OF N M e

In the previous section we have shown that the electron gas in stellar matter, because of its high incompressibility, acts against the formation of macroscopic charge dishomogeneities, thus completely quenching the liquid-gas phase transition, at least in the mean-field approximation adopted here. This does not mean that we should expect stellar matter at subsaturation density to be homogeneous. Indeed it has been recognized since a long time [28] that the ground state of neutron-rich stellar matter at low density should be a Coulomb lattice of exotic nuclei embedded in a neutron gas, meaning that we should expect a solid-liquid transition between the core and the inner crust of a neutron star [36]. Numerical simulations performed with (semi)-classical models [16,17] indicate that such dishomogeneous phases should persist at finite temperature and in a large range of asymmetry. Clusterisation of hot stellar matter is also considered in a statistical approach in ref. [29,30]. A way to address the possible existence of a phase of clusterized matter in mean-field approaches, is to study the instability of nuclear matter with electrons against finite-size density fluctuations [31]. If such fluctuations lead to a reduction of the free-energy, they will be amplified through the process of spinodal decomposition [32] until the formation of structures of alternatively high and low density. Homogeneous matter is then unstable against clusterization. Clusters correspond to the high-density regions, which can be seen as liquid nuclear drops immersed in a nuclear gas. Considering the short time scales associated with spinodal decomposition, the cluster characteristics are determined by the free-energy curvature properties at the onset of instability, when density fluctuations have small amplitudes [32].

A. Free-energy variation with finite-size fluctuations

The instability region of homogeneous NMe against clusterization is determined by studying how the system free energy changes, when a finite-size density fluctuation is introduced.

If fluctuations occur on a finite microscopic scale, electron and proton density can fluctuate independently : only their mean values are constrained to be equal, insuring macroscopic electroneutrality. Fluctuations thus affect independently the three species of the medium (neutrons, protons and electrons), whose densities become :

ρ q = ρ 0 q + δρ q , (13)
with q ∈ {n, p, e}. Each density variation can be expressed by a Fourier transform

δρ q = d ka q ( k)e i k• r , (14)
with a q ( k) = a * q (-k) to insure that δρ q is real. Since the different wave vectors k are decoupled in the global free-energy variation, the problem reduces to the study of planewave density fluctuations:

δρ q = A q e i k• r + c.c., (15)
where each species is associated with a different amplitude. Thermodynamic fluctuations, studied in the previous section, can be obtained as the limit k → 0 of this model.

To evaluate the free-energy variation, we consider a Thomas-Fermi approximation, i.e.

the density variation is supposed smooth enough to allow in each point the definition of a Fermi sphere corresponding to the local density. Then, in each point of density ρ q (r) = ρ 0 q + δρ q , the local bulk term of the free energy is equal to the free energy f h of an infinite homogeneous system at the same density. The global bulk free energy of the system is the space average of this local term :

f b = 1 V f h ({ρ q (r)})dr. ( 16
)
In the small-amplitude limit, the integration leads to f b = f h ({ρ 0 q }) + δf b , with :

δf b = i,j A i A * j + A * i A j 2 ∂ 2 f h ∂ρ i ∂ρ j {ρ 0 q } = i,j A i A * j + A * i A j 2 ∂µ j ∂ρ j {ρ 0 q } , ( 17
)
where i, j ∈ {n, p, e}. First-order terms have vanished in the integration because the average density variation is zero.

The variation of the entropy is contained in the bulk term, since entropy depends only on the local density [26]. However, in the case of a finite wave number k, the energy density is modified by two additional terms [31], arising from the density-gradient dependence of the nuclear force, and from the Coulomb interaction. Denoting these two contributions δE ∇ and δE c respectively, the free-energy variation is :

δf = δf b + δE ∇ + δE c . ( 18
)
Let us now explicit the gradient and Coulomb contributions. In the presence of density gradients the nuclear energy density has the form H = H h + H ∇ , where H h is given by eq.( 2) and the density-gradient term H ∇ is expressed as:

H ∇ = C ∇ nn (∇ρ n ) 2 + C ∇ pp (∇ρ p ) 2 + 2C ∇ np ∇ρ n • ∇ρ p . ( 19
)
Coefficients C ∇ ij are combinations of the Skyrme parameters given for example in ref. [21]:

C ∇ nn = C ∇ pp = 3 16 [t 1 (1 -x 1 ) -t 2 (1 + x 2 )] , C ∇ np = C ∇ pn = 1 16 [3t 1 (2 + x 1 ) -t 2 (2 + x 2 )] . (20)
The global contribution of this term to the energy density is given by the space average:

δE ∇ = 1 V H ∇ ( r)d r = k 2 i,j A i A * j + A * i A j C ∇ ij , (21)
where i, j ∈ {n, p}.

The Coulomb contribution δH c is due to the independent density fluctuations of charged particles {e, p}. Denoting q i = √ 4πǫ 0 e i the electric charge of a particle of type i, this creates the charge distribution :

ρ c ( r) = i e i A i e i k• r + c.c. ( 22
)
The consequent Coulomb energy per unit volume is :

δE c = 1 2V ρ c (r 1 )ρ c (r 2 ) |r 2 -r 1 | dr 1 dr 2 = 4π k 2 i,j A i A * j + A * i A j 2 , ( 23
)
where i, j ∈ {e, p}.

Summing the contributions ( 17), ( 21) and ( 23), we have the total free-energy variation, expressed to the second order in A q :

δf = ij A i A * j + A * i A j 2 ∂µ i ∂ρ j {ρ 0 q } + 2C ∇ ij k 2 + 4πe i e j k 2 , ( 24
)
which can be written in a matrix form in the three-dimensional space of density fluctuations à = (A n , A p , A e ) :

δf = Ã * C f Ã, (25)
where

C f =      ∂µ n /∂ρ n ∂µ n /∂ρ p 0 ∂µ p /∂ρ n ∂µ p /∂ρ p 0 0 0 ∂µ e /∂ρ e      + k 2      2C f nn 2C f np 0 2C f pn 2C f pp 0 0 0 0      + 4π 2 e 2 k 2      0 0 0 0 1 -1 0 -1 1      (26)
is the free-energy curvature matrix. The bulk term of the curvature is conserved in this new expression, with the difference that the electron term now belongs to an independent space. The density-gradient part of the nuclear interaction adds a term proportional to k 2 , while the Coulomb interaction induced by the plane-wave charge distribution adds a term inversely proportional to k 2 .

B. k-dependent instability

An instability against matter clusterization corresponds to a negative free-energy curvature in the space of density fluctuations. It is studied through the k-dependent curvature matrix C f . The eigen-mode of minimal curvature is represented in fig. 4 for one given point of average density (ρ 0 n , ρ 0 p ), as a function of k, by three quantities : its eigen-value, its direction in the (ρ p , ρ n )-plane and in the (ρ p , ρ e )-plane.

The upper part of fig. 4 gives the minimal eigen-value C < of the curvature matrix. We observe that, for the chosen average-density point, there is an interval of k allowing a negative free-energy curvature in the space of density fluctuations : the corresponding plane wave will then be amplified leading to matter clusterization. The lower limit of the instability interval is due to the Coulomb term (∝ 1/k 2 ), which tends to suppress long wavelenght instabilities which would lead to the formation of a macroscopic charge. The upper limit is determined by the nuclear density-gradient term (∝ k 2 ), which acts against the surface formation associated with short wavelenght fluctuations. Thus, both effects arising from finite-size fluctuations lead to a reduction of the instability, with respect to the case of thermodynamic phase separation in uncharged nuclear-matter.

Let us now consider the phase-separation direction, which is given by the eigen-vector associated with the minimal curvature. The middle part of fig. 4 represents this direction projected in the (ρ n , ρ p )-plane. Comparing this direction to the line of constant proton fraction, we can see that the amplification of fluctuations will lead to a dense phase (clusters)

which is more symmetric than the surrounding gas. This is the well-known phenomenon of isospin distillation [9], resulting from the minimization of the nuclear symmetry energy.

The lower part of the figure represents the eigen-vector direction in the (ρ p , ρ e ) plane, i.e. it represents the mobility of electrons. This behavior is essentially determined by the relative weight of coulomb interaction and electron incompressibility, as shown by the fact that the complete neglect of nuclear interaction (dashed line) changes the results only slightly. At low k-values, this quantity is dominated by Coulomb energy, which forces the electrons to closely follow the protons in order to lower the electric charge. Conversely, for high values of k, electron incompressibility dominates over Coulomb interaction, and the electron gas tends to stay a uniform background.

In order to better quantify the effect of electrons, the phase-separation direction is compared to two extreme idealized cases represented on the same figure. In both cases, the proton charge is neutralized by a gas of negative charges : the incompressibility of this gas is set to either zero or infinity. In the first case, the Coulomb energy is exactly cancelled : the resulting curves show the effect of nuclear density-gradient term, the point at k = 0 corresponding to homogeneous nuclear matter as discussed in section II A. In the second case, the negative gas forms a strictly homogeneous background, which suppresses the monopolar term of Coulomb interaction. In reality the electron gas has a finite incompressibility, related to the susceptibility χ -1 e . As noticed in section II B, the high value of χ -1 e strongly suppresses the electron degree of freedom. As a result, the electron dimension in the curvature matrix (26) has a only a perturbative effect with respect to the case of a uniform negative background. Since the (small) electron mobility allows to lower the Coulomb energy, it makes the fluctuation more favorable and leads to a slight increase of the instability interval.

Fixing different values of k, we have represented in fig. 5 the corresponding instability regions (k-spinodals) in the density plane (ρ n , ρ p ) at a given temperature. Each point of this plane gives the initial densities of the homogeneous system, from which fluctuations start. The effect of the finite electron incompressibility is more visible at lower k, where the Coulomb term dominates. For a fixed k, the perturbation increases as χ -1 e decreases, i.e. towards higher values of ρ 0 e = ρ 0 p . Starting from k ≈ 80 MeV/c the electron gas can safely be represented as a fixed background. The curves presented on fig. 5 can be qualitatively compared with the results of ref. [27] (obtained in a RMF approach at zero temperature).

Nuclear-matter spinodal region are very similar in both cases. When finite-size instabilities are considered, a similar behaviour is obtain for the higher values of k (k 100). For the lower values, Coulomb contribution becomes dominant and is found to affect differently the shape of the k-spinodals : fig. 5 shows a stronger reduction of instability in the proton-rich part of the plan, while the curves remain essentially symmetric in the RMF calculation. This signs a different evolution of the derivative ∂ ρn µ n through the density plane (indeed, when Coulomb interaction dominates, the direction of minimal free-energy curvature gets closer to the neutron-axis direction).

The k-interval, for which NMe presents an instability in at least one density point, is represented on fig. 6 as a function of the temperature T . It is maximal at zero temperature, and goes to zero at a limiting temperature which is about 12MeV for the realistic forces (almost 16MeV for SIII).

The global instability region (at a given temperature) is defined as the envelope of the different k-spinodals : it contains all the points presenting an instability against clusterization for at least one value of k. It is represented in fig. 7 for three different Skyrme parameterizations : Sly230a, SGII and SIII. We can see that the two realistic forces give similar envelopes. SIII presents an atypical behaviour, with a much wider instability region. This is essentilly due to the too high value of K ∞ (nuclear incompressibility at saturation density) (see Fig. 1).

C. Influence on cluster properties

In a spinodal-decomposition scenario, phase properties are mostly determined by the conditions at the onset of instability [32]. In particular matter clusterization is dominated by the development of the most unstable mode, i.e. e , compared to two limiting cases (zero and infinite incompressibility of the negatively charged gas) for Sly230a. The curve correponding to χ -1 = 0 goes to zero at the critical temperature of nuclear matter T u . For the physical N M e case, the instability interval disappears for a limit temperature T lim < T u . This temperature is slightly higher for N M e. Bottom : N M e curves for three Skyrme parameterizations.

is given by:

δρ q = A q e i( k• r-ωt) + c.c., (27)
where unstable modes are characterized by an imaginary frequence. Because of the exponential amplification of such modes, clusterization is dominated by the wave number k 0 corresponding to the highest imaginary frequence |ω (k 0 ) | = max. The dispersion relation ω(k) can be established in a hydrodynamic description. This is easily derived [32] in the case of one single fluid, from the Euler and continuity equations:

d(mρv)/dt = -∇P, ( 28
)
∂ρ/∂t + ∇(ρv) = 0, ( 29
)
where P is the system's pressure, related to the free-energy density f by P = ρdf /dρ. These equations in the limit of small amplitude fluctuations ρ = ρ 0 + δρ can be expressed as a diffusion equation

∂ 2 ρ ∂t 2 = ρ 0 C∇ 2 ρ, (30)
where C = d 2 f /dρ 2 | ρ 0 is the curvature of the homogeneous-system free-energy density. A Fourier representation of this equation leads to the dispersion relation ω 2 = ρ 0 Ck 2 . When this curvature is negative, the frequency is imaginary : it can be expressed as ω = i/τ , where τ is the time constant characterizing the development of fluctuations.

In our study, the instability direction for a given wave-number k in each point of the density plane, is associated with the eigen-value C < (k) of the curvature matrix C. Since the temporal evolution of this minimal-curvature mode is given by τ 2 ∝ 1/(-C < k 2 ), the most unstable mode for each temperature and density can be determined by the condition -Ck 2 0 = max. In the diabatic limit that the successive dynamics can, as a first approximation, be neglected [32], we can link the free-entropy curvature matrix of homogeneous matter, to final cluster properties as issued from the subsequent spinodal decomposition. In particular, an order of magnitude for cluster sizes is given by the half wavelength λ 0 /2 = π/k 0 ; cluster isotopic composition is related to the direction of the eigen-vector u < (k 0 ) associated with

C < (k 0 ).
The evolution of the most unstable mode through the density plane is given by fig. 8, for a fixed temperature. The upper part presents quantities linked to cluster size : the most unstable wave-number k 0 , and its associated half wavelength λ 0 /2. For a wide range of asymmetries, the typical cluster size depends essentially on the total density, rising towards the instability boundary. Level lines circle a minimum at low density, quasi-symmetric matter. At these low densities the very high frequency of the most unstable mode gives typical sizes of a few f m, corresponding to clusters of ≈ 10 particles. This result should be taken with care: the Thomas-Fermi approximation inherent in our treatment breaks down for such small clusters, and we may be facing the limit of validity of the mean field approach.

The lower part of the figure presents properties affecting the isotopic composition of clusters: the direction of u < (k 0 ) in the nuclear-density plane is given by the rapport δρ 3 /δρ.

It is compared to the direction of constant proton fraction by the difference ρ 3 /ρ -δρ 3 /δρ.

We can see that clusterization almost systematically involves portions of matter which are more symmetric than the global system, which is again a manifestation of isospin distillation.

A representation of the most unstable mode is given by fig. 9 and10. They show some characteristic evolutions for the half wavelentgh and the phase separation direction as a function of density and proton fraction, at different temperatures and using different Skyrme parametrizations. We can see that we can expect clusterized stellar matter with similar characteristics in a very wide range of temperatures, densities and asymmetries, and these qualitative results appear robust respect to the chosen Skyrme parametrization. In particular realistic forces, irrespective of their different symmetry energy at subsaturation densities, lead to very close results for all observables. Increasing the temperature has the effect of Bottom-right : difference respect to the constant Z/A direction, ρ 3 /ρ -δρ 3 /δρ(k 0 ), levels regularly spaced of 0.05.

ρ ρ) δ / ρ -(δ (k ) 0 ρ / ρ 0.2 0 -0.2 3 3 ρ δ / ρ δ 06 ) -3 (fm p ρ ) -3 (fm n ρ ) -3 (fm n
increasing the cluster size and decreasing the distillation. The electron participation to phase separation is given by δρ e /δρ p : this ratio is low, and correlated with λ 0 (the larger clusters lead to a higher polarization of the electron gas).

In the absence of Coulomb interaction (as for the study of normal nuclear matter), isospin symmetry imposes along the ρ n = ρ p axis the relation δρ 3 /δρ = ρ 3 /ρ = 0. In the present case, the presence of a proton charge hinders ρ p fluctuations, leading to an enhancement of δρ 3 /δρ. Level-lines δρ 3 /δρ = 0 and ρ 3 /ρ -δρ 3 /δρ = 0 are respectively shifted to the protonrich and neutron-rich part of the density plane, delimiting a region of abnormal isospin distillation , shown in fig. 11 for a temperature T = 5 MeV. On the neutron-rich side of this region, phase separation leads to a dense phase more asymmetric than the global system, and a proton-rich gas. On the proton-rich side, for a large enough final density the dense phase is neutron rich. However, the region concerned by this anomaly is thin : for the major part of the instability region, the isospin distillation phenomenon of normal nuclear matter is not qualitatively modified by the presence of electric charge.

IV. β-EQUILIBRIUM IN N M e

We now study the case of a β-equilibrated star matter, composed of neutrons, protons, electrons and neutrinos ν. For simplicity, leptons will be treated in the degenerated-gas approximation. If this equilibrium is reached at a given temperature and density, it imposes a given isotopic composition, determined by the rate equality between reactions of electronic capture and neutron disintegration :

p + e -→ n + ν n → p + e -+ ν. ( 31
)
When the core of a massive star collapses, the increase of density displaces β-equilibrium towards a more neutron-rich composition.

β-equilibrium corresponds to free-energy minimization. In the presence of different species i of mass m i , the total free energy density is given by :

f t ({ρ i }) = f ({ρ i }) + i ρ i m i . ( 32
)
f t derivatives define the chemical potentials µ t i : The grey area gives the region of abnormal distillation due to the Coulomb interaction, situated between the two lines δρ 3 /δρ = 0 and δρ 3 /δρ = ρ 3 /ρ. (δρ 3 /δρ is the phase-separation direction given by the most unstable mode). so that β-equilibrium corresponds to the relation

µ t i = ∂ ρ i f t = ∂ ρ i f + m i = µ i + m i , (33)
µ t p + µ t e = µ t n + µ t ν µ p + m p + µ e + m e = µ n + m n + µ ν . (34)
Eq. (34) supposes that neutrinos participate in β-equilibrium. This is true in the initial post-bounce supernova dynamics, when matter is dense enough for neutrinos to be completely trapped in the star. In the later phase of the cooling process, matter becomes transparent to neutrinos [33,34], which then do not equilibrate their chemical potential with the one of the other species. The transition between these two regimes is smooth, and in the following we will study post-bounce supernova matter considering various rates of neutrino trapping. The detailed relation between density, temperature and neutrino trapping in supernova matter is not a settled issue, especially concerning the crust region which is of interest for the liquid-gas phase transition. For this reason, for our study we treat the trapping fraction as a free parameter, denoted x ν : it is defined as the total fraction of neutrinos produced by electron-capture reactions which is retained in the star. We suppose an initial proton fraction Z/A 0 corresponding to 56 F e (which is the main constituent of the star core before collapse). For each point (ρ, Z/A) of the nuclear density plane, x ν and Z/A 0 determine the neutrino density inside the matter according to :

ρ ν = x ν ρ [(Z/A) 0 -Z/A] . (35)
Leaving x ν as a free parameter, we can estimate the effect of neutrino trapping on the thermodynamic properties of the system. The neutrino density is linked to the chemical potential by the ultra-relativistic degenerate Fermi gas expression :

µ ν = c(6π 2 ρ ν ) 1/3 . (36)
Eqs.( 34),( 35),(36) describe a trajectory in the plane of chemical potentials, which fixes the isotopic composition of matter at each temperature as a function of the x ν parameter, linked to neutrino opacity. The β-equilibrium compositions corresponding to different values of

x ν are represented in fig. 12 in the space of chemical potentials (µ t n , µ t pe = µ t p + µ t e ). The diagonal in this space gives the β-equilibrium curve for a matter completely transparent to neutrinos. Any point of the plane is β-equilibrated for a value of µ ν (i.e. a number of trapped neutrinos) corresponding to the vertical shift respect to this line. The lines corresponding to x ν = 0 and 1 give the region which is in principle accessible to β-equilibrium : this crosses the instability envelope, which is reported as a thick line in the chemical-potential plane. We can see that the matter composition at β-equilibrium is very sensitive to the trapping fraction x ν : only a small fraction of trapped neutrinos is sufficient to sensibly increase the stellar proton fraction. At zero temperature, the instability envelope is crossed for any value of x ν . At higher temperatures of the order of T = 10MeV , relevant for the initial dynamical steps of the proto-neutron-star, matter at β-equilibrium is too asymmetric to be unstable against clusterization if all neutrinos escape from the cooling star, while a fraction as small as x ν = 0.01 is sufficient to lead the system to explore the instability region. This result depends very slightly on the different Skyrme parametrizations, meaning that we can reasonably expect that matter will be clusterized in the whole evolution of the proto-neutron-star.

A representation of β-equilibrium composition in the plane (ρ, Z/A) is given by fig. 13 for different temperatures and Skyrme interactions. Each path corresponding to a fixed value of x ν shows a quasi-constant proton fraction. The proton fraction shows very little evolution until low densities : then it first slightly decreases before going abruptly to 1 in the limit of zero density. Indeed, at vanishing density the mass balance m n -m p -m e > 0 dominates β-equilibrium, imposing a full decay of neutrons into protons. However, since this rising happens at very low density, fig. 13 shows that for a fixed x ν β-equilibrium essentially corresponds to a constant-Z/A line.

We should stress that in this calculation nuclear matter is always considered as homogeneous. When clusters are formed, β-equilibrium does not correspond any more to the simple mean-field trajectories represented in figs.12 and 13. Thus the β-equilibrium curves displayed by these figures represent the initial condition for matter clusterization. The instability envelope we have represented corresponds to the region where clusters can be formed by spinodal decomposition. This gives a lower-limit estimate for the density interval corresponding to clusterization. Indeed the region for which completely equilibrated matter is clusterized has to be wider than the spinodal-instability envelope. However, outside this envelope, matter can in principle be homogeneous and β-equilibrated since in a situation of meta-stability the time needed for cluster formation can be longer than β-equilibrium time scale. Let us also consider the case the case of a thermally equilibrated initially homogenous dense matter. If it is submitted to a fast dilatation, Z/A remains constant while density decreases. Then, the proportion of neutrinos trapped in the initial state, which fixes the constant composition Z/A, determines wether the instability envelope is crossed or not.

Figure 14 gives a quantitative view of the neutrino fraction needed for β-equilibrated matter to reach the region of instability against clusterization, as a function of temperature.

The value x min ν corresponds to the path at constant x ν which is tangential to the instability envelope. The curves obtained show a slow evolution with very small values of neutrino fraction at low temperature, and a sharp enhancement as the limiting temperature for instability is approached. This profile is due to the typical evolution of the instability region, since, as we have seen, β-equilibrium paths depend only weakly on the temperature.

V. SUMMARY AND CONCLUSION

The nuclear-matter liquid-gas phase transition results in a region of instability against cluster formation in low-density finite-temperature star matter composed of neutrons, protons and electrons. This has been studied analyzing homogeneous-matter instabilities against a plane-wave density fluctuation, in the small-amplitude limit. In a spinodal decomposition scenario, cluster properties as sizes and chemical compositions can be predicted from the characteristics of the most unstable mode at the onset of fluctuation amplification : they are determined studying the free-energy curvature matrix in the space of density fluctuations, as a function of the wave number k.

For a cold neutron star, thermodynamic equilibrium is reached. The liquid-gas transition concerns the inner crust, where pasta phases are expected. Our approach does not include the study of such structures, since it is based on homogeneous-matter properties : however, the instability envelope which we obtain gives the minimal region where equilibrated matter is clusterized. The out-of-equilibrium phase separation may concern proto-neutron stars issuing from core-collapse in supernova explosions, since the shock-wave causes abrupt changes in the density. At this time, star matter is hot and contains trapped neutrinos, whose propagation is a key issue for the understanding of supernova dynamics [18]. It has been pointed out [17] that neutrino in-medium propagation is affected by the eventual presence of clusters. On the other hand, we have found that presence of neutrinos affects the star-matter composition in a significant way. Indeed neutrino trapping rules the crossing of the instability region, which in turn leads to the formation of clusters. This interplay is sensitive to the different time scales implied in the explosion : temperature and density evolution, as well as the realization of β and thermodynamic equilibrium.

T (MeV)

FIG. 1 : 2 ∂ 2

 122 FIG. 1: Energy proporties at zero temperature for the different Skyrme forces used in this paper. Saturation points are indicated by dots. Top : energy per nucleon of symmetric matter. Bottom : density dependence of the asymmetry coefficient a s = 1 2 ∂ 2 E/A ∂I 2 I=0

FIG. 2 :

 2 FIG.2: Minimal free-energy curvature obtained with the Sly230a interaction. Ordinary nuclear matter (N M )presents a spinodal region, which is suppressed when an electron gas is added to neutralize the proton charge (N M e). Left : lower eigen-value as a function of the total density ρ, for different Z/A from Z/A = 0.5 to Z/A = 0, equally spaced of 0.1. Right : minimal curvature as a function of Z/A. Dotted line : M N e at T = 0.

  : the middle part presents the minimal curvature as a function of total density for MNe with different values of Z/A, at T = 10MeV . The lowest minimal curvatures are reported on the right part as a function of Z/A. No negative value of the curvature is found, and the same is true at zero temperature. Since the extension of the spinodal region monotonically increases with decreasing temperatures, this result implies that NMe does not present any thermodynamic instability.

FIG. 3 :

 3 FIG.3: Suppression of spinodal instability in nuclear matter with electrons, for different Skyrme forces at T = 0. Top : minimal free-energy curvature obtained for a given proton fraction, as a function of Z/A. Bottom : density at which this minimum is obtained.

FIG. 4 :

 4 FIG.4: Eigen-mode of the minimal free-energy curvature for ρ = 0.05f m -3 , Z/A = 0.3, T = 5MeV, as a function of the wave number k, calculated with the Sly230a Skyrme interaction. Top : eigen-value. Middle : associated eigen-vector in the nuclear-density plane. The N M e-curves are compared to two limiting cases corresponding to a zero and infinite incompressibility of the negatively charged gas. When the curvature is negative (full lines) this eigen-vector gives the phase-separation direction. The dots give the points of zero curvature. The dashed line gives ρ 3 /ρ (direction of constant Z/A). Bottom : same as the middle part, in the (ρ p , ρ e )-plane. The dashed line gives δρ e /δρ p when the nuclear force is zero.

FIG. 5 :FIG. 6 :

 56 FIG.5: k-spinodals for different values of k at T = 5M eV . The external dotted curve is the thermodynamic spinodal region of nuclear matter, it is the same in all panels. The k-spinodals of N M e are compared to two limiting cases (zero and infinite incompressibility of the negatively charged gas). The thermodynamic spinodal (dotted curve) corresponds to the case χ -1 = 0 (dashed curve) at k = 0.

FIG. 7 :

 7 FIG. 7: Envelope of instability against finite-size density fluctuations at T = 5M eV . For every point inside this region, at least one fluctuation mode can be amplified. Left : some N M e k-spinodals (k values given in M eV /c) and the instability envelope (thick curve). N M thermodynamic spinodal is also given (external thin curve). Right : N M e instability envelope for three different Skyrme forces.

FIG. 8 :

 8 FIG. 8: Most unstable mode properties through the density plane, represented in level lines. Top-left : most unstable wave number k 0 (M eV /c), levels regularly spaced of 20M eV /c. Topright : corresponding half wavelength λ 0 /2(f m), levels regularly spaced of 1f m. Bottom-left : most-unstable mode direction in nuclear-density plane δρ 3 /δρ(k 0 ), levels regularly spaced of 0.1.

FIG. 9 :FIG. 10 :FIG. 11 :

 91011 FIG. 9: Cluster properties as indicated by the most unstable mode, along two axis : one at constant density (left side), the other at constant Z/A (right side). Three temperatures are considered with Sly230a.

FIG. 12 :

 12 FIG. 12: β-equilibrium in chemical-potential plane (µ t n , µ t pe = µ t p + µ t e ), for different temperatures and neutrino opacities x ν = 0, 0.01, 0.1 and 1. The lines x ν = 0 and 1 delimit the region which is in principle accessible to β equilibrium. The N M e instability envelope is also represented.

FIG. 13 :

 13 FIG.13: β-equilibrium in the plane (ρ, Z/A) for different temperatures and neutrino opacities x ν = 0, 0.01, 0.1 and 1. M N e instability envelope is also represented.

FIG. 14 :

 14 FIG. 14: Temperature dependance of the minimum neutrino opacity x min ν leading to the instability

  ∂µ n /∂ρ n ∂µ n /∂ρ p ∂µ p /∂ρ n ∂µ p /∂ρ p + ∂µ e /∂ρ e

 