C Ducoin

J Margueron

Ph Chomaz

Cluster formation in asymmetric nuclear matter: semi-classical and quantal approaches

Keywords: numbers: asymmetric nuclear matter, compact stars, liquid-gas spinodal instabilities, Skyrme forces, semi-classical theory, random-phase approximation

Cluster formation in asymmetric nuclear matter:

semi-classical and quantal approaches

I. INTRODUCTION

In infinite nuclear matter, a phase transition of the liquid-gas type [1,2] is predicted to occur at sub-saturation density with a critical temperature of about 15 MeV, depending on the theoretical models. This phenomenon corresponds to the separation of the system into two macroscopic (infinite) phases of different densities. This thermodynamic property is at the origin of finite-size spinodal instabilities, leading to the clusterization of the homogeneous medium. In this process, in contrast with the thermodynamic case, the Coulomb interaction can be included and density-gradient terms come into play: the liquid-gas instabilities can then be related to the formation of fragments occurring in finite nuclear systems. In heavy-ion collisions performed in nuclear facilities, around the Fermi energy, several signals of the liquid-gas phase transition have been analyzed [3]. In particular, the spinodal instabilities could be at the origin of the multifragmentation process [4,5,6]. Several models have been built in order to describe the dynamics leading to the formation of clusters [7,8]. Clustering is also expected to play an important role in compact stars [9,10]. It is predicted to lead to the presence of complex structures like the "pasta phases" occurring in the inner crust of neutron stars [11,12,13,14]. At lower densities and finite temperature, star matter may clusterize into light nuclei [15,16]. Around the core of supernovae, density fluctuations could have important effects on neutrino propagation [17], which is crucial for the description of supernova dynamics [18].

This paper addresses the clusterization properties of an infinite nuclear medium: the onset of the liquid-gas instabilities is characterized by studying the incidence of a small density fluctuation introduced in the homogeneous medium. Such approach is based on a spinodal-decomposition scenario, which assumes that fragments are formed in an outof-equilibrium process. Indeed, if a homogeneous system is quenched into the region of finite-size instabilities, it will decompose into clusters: in the spinodal-decomposition scenario, the early dynamics governing this evolution consists in the rapid amplification of the unstable fluctuations, whose properties can be understood in the small-amplitude limit. For a fast reaction, fragments could be separated before equilibrium is realized, thus keeping the memory of the bulk instability. Identifying fragments produced through spinodal decomposition is then of great importance to obtain informations about the homogeneous-matter energetics at sub-saturation density. Experimental indications of this mechanism in heavy-ion collisions have been obtained for instance through the analysis of size correlations [19,20]. We should note that the study of finite-size instabilities in homogeneous matter is complementary to the description of an equilibrated, clusterized system, as can be given at very low density by the virial equation of state [15]. With I: Saturation properties of the three Skyrme forces SLy230a, SGII and RATP. Are given the values for the saturation density ρ0: the binding energy E/A; the incompressibility modulus K∞; the symmetry energy as,0; the symmetry energy slope L = 3ρ0(∂as/∂ρ); the symmetry energy curvature Ksym = 9ρ 2 0 (∂ 2 as/∂ρ 2 ); the critical temperature of the liquid-gas phase transition Tc.

this respect, the finite-size instability region we discuss can be viewed as the minimal region where the equilibrated system is formed of clusters.

The model of infinite nuclear matter provides a useful insight into the behavior of nuclear systems, revealing robust qualitative features whose origin can be clearly analysed. The spinodal decomposition studied in this framework can be related to the nuclear multifragmentation process by neglecting the finite size of the homogeneous medium formed in the early stage of the collision. It also gives a picture of the instability against clusterization in compact-star matter, where the degree of freedom associated with the electrons can be neglected [21].

Our calculations are based on Skyrme density functionals. Specific attention is paid to the isospin degree of freedom, relevant for the physics of compact stars and exotic nuclei. The study of the instabilities is treated in three different approaches. Starting from the thermodynamic definition of the liquid-gas phase transition, we first concentrate on the free-energy curvature properties. This leads us to study the finite-size density fluctuations in a Thomas-Fermi approach [21,22], defining a static criterion for the existence of instabilities based on the free-energy variations. Afterward, consistently with this local equilibrium framework, we treat the time evolution by introducing hydrodynamical equations. This defines a first semi-classical approach of the clusterization dynamics. The second approach is fully quantal and dynamical: it is based on the linearization of the Time-Dependent Hartree-Fock theory (TDHF) [START_REF] Ring | The Nuclear Many-Body Problem[END_REF]. This corresponds to the random-phase approximation (RPA). In this case, the evolution of the instability is described in the collisionless limit. The differences expected between hydrodynamics and RPA results have then two origins: the quantal effects and the collisionless description. In order to disentangle these two effects, we will also consider a third approach, the Vlasov method which is the semi-classical limit of the RPA. It then includes the collisionless description but not the quantal effects.

Different Skyrme forces have been used in order to check the dependence of the results on the choice of a specific set of parameters. We take as a reference the Skyrme-Lyon parameterization SLy230a [START_REF] Chabanat | [END_REF], a modern force designed to treat neutron-rich matter, already employed for astrophysical applications [25,26]. Two other Skyrme forces are used for comparison: RATP [27] -the first interaction which includes neutron matter into the fitting procedureand SGII [28], a typical widely used parameterization. The interaction RATP is constrained by the neutron-matter equation of state calculated by Friedman and Pandharipande in the 80ies [29], while the interaction SLy230a uses a more recent neutron-matter equation of state based on an improved realistic interaction [30]. Let us remark that these interactions have quadratic velocity-dependent terms which modify the single-particle energies and give rise to effective masses m * i (i=n, p). The effects of an isovector momentum-dependent mean field on nuclear-matter properties have been studied in Ref. [31].

For these three interactions, Fig. 1 shows as a function of the density the binding energy E/A in symmetric nuclear matter and in pure neutron matter, as well as the symmetry energy a s . The different saturation properties are reported in Tab. I. One can see that the predictions of scalar properties, saturation density, energy and incompressibility are very close for these realistic forces. The values of the symmetry energy at saturation are also quite close, the main difference being its density dependence, as shown on Fig. 1. The symmetry energy slope L and curvature K sym are also indicated in Tab. I. It should be noticed that the selected Skyrme forces present an "asy-soft" behaviour, namely values of L lower than 60 MeV, while recent constraints extracted from different experimental analysis indicate that L shall be larger than about 60 MeV [32,33,34]. This however has no drastic effects on the results we present, since the spinodal contour in asymmetric nuclear matter is non linearly related to the symmetry energy and its first and second derivatives, cf Eq. (6). A systematic analysis of the spinodal properties for various relativistic and non-relativistic interactions, "asy-soft" and "asy-stiff", is in preparation [35].

The paper is organized as follows. Before discussing the finite-size instabilities, we address the thermodynamic spinodal properties in Sect. II. The hydrodynamical description for the clusterization process is treated in Sect. III. The collisionless dynamics based on the TDHF equation is treated in Sect. IV, where the quantal RPA and Vlasov methods are addressed. We discuss the different results in Sect. V, and conclusions are drawn in Sect. VI.

II. SPINODAL REGION FROM A THERMODYNAMICAL APPROACH

In this section, we address the thermodynamic liquid-gas phase transition of nuclear matter, where each phase is infinite and homogeneous. This allows a geometrical representation of the phase equilibrium, obtained by a Gibbs construction on the free-energy surface: thus can be emphasized the fundamental role of the free-energy curvature in the liquid-gas instabilities studied in the following sections.

The thermodynamic equilibrium of a system is given by maximizing the entropy in the space of observables. For homogeneous nuclear matter, this surface presents an anomalous convexity in the phase-transition region. In the thermodynamic limit where interface contribution is neglected, this convexity can be corrected maximizing the entropy by linear interpolations which define pairs of phases in coexistence. In other words, the thermodynamic equilibrium is given by the Gibbs construction, which delimits the phase-transition boundaries. This coexistence region contains the spinodal region, where homogeneous matter is locally unstable (the entropy surface being locally convex). For a given value of the temperature, this is equivalent to a negative value of the free-energy curvature matrix, which we now determine.

In the Skyrme density-functional theory (DFT), the average energy density of homogeneous nuclear matter H = Ĥ /V is a function of the one-body density matrix. The full expression of H obtained from a Skyrme-like interaction is given for example in Ref. [START_REF] Chabanat | [END_REF]. Considering homogeneous, spin-saturated matter with no Coulomb interaction, the energy-density functional reduces to:

H = H h = 2 2m τ + C eff ρτ + D eff ρ 3 τ 3 +C 0 ρ 2 + D 0 ρ 2 3 + C 3 ρ σ+2 + D 3 ρ σ ρ 2 3 , (1)
where the superscript "h", introduced for later use, means "homogeneous". Coefficients C an D are linear combinations of the standard Skyrme parameters [START_REF] Chabanat | [END_REF]37]. We have introduced the isoscalar and isovector particle densities, ρ and ρ 3 , as well as kinetic densities, τ and τ 3 :

ρ = ρ n + ρ p ; τ = τ n + τ p ρ 3 = ρ n -ρ p ; τ 3 = τ n -τ p (2)
where, denoting i the third component of the isospin (n for neutrons and p for protons), the kinetic densities are defined by τ i = p2 i / 2 . The energy density H is then a function of neutron and proton densities {ρ i , τ i }. From the Skyrme functional is derived the mean field for each type of particle i = (n, p). In uniform matter, the mean field reduces to a potential U i and an effective mass m * i such that the single-particle energies read:

ǫ i (p) = p 2 2m * i + U i ; U i = ∂H ∂ρ i ; 2 2m * i = ∂H ∂τ i , ( 3
)
where p is the momentum quantum number of the free particles. In the grand-canonical formalism, individual levels are occupied according to the Fermi-Dirac statistics with occupation numbers n

i (p) = [1 + exp(β(ǫ i (p) -µ i ))] -1
, where β = 1/T is the inverse temperature, and µ i the chemical potential of the corresponding nucleon type i. The mean-field grand-canonical partition function Z 0 is formally identical to the free-Fermi-gas one, replacing µ i by ν i = µ i -U i and m i by m * i . At given values of T , ρ n and ρ p , chemical potentials and kinetic densities are fixed, determining the energy density H and the mean-field entropy density s, directly related to Z 0 . The free-energy density f(T, ρ n , ρ p ) is deduced from the Legendre transform f = H -T s. We denote f h the free-energy density of a homogeneous system:

f h (T, ρ n , ρ p ) = H h -T s . (4)
In the following, working at constant temperature, the dependence of the functions on T will be implicit for a more concise notation. For a given temperature, the thermodynamic spinodal region corresponds to a negative curvature of the surface f h (ρ n , ρ p ), i.e. a negative eigen-value of the curvature matrix C h = {C h ij }, with:

C h ij (ρ n , ρ p ) = ∂ 2 f h ∂ρ j ∂ρ i = ∂µ i ∂ρ j , (5)
whose lower eigen-value is denoted C h < . Figure 2 shows the spinodal contours obtained for two values of temperature (T = 0 and 10 MeV) and three Skyrmelike interactions: SLy230a, SGII and RATP. Small quantitative differences can be noticed between these interactions, especially at large asymmetry and density. These differences can be analyzed in terms of isovector properties of the density functional at (ρ n = ρ p ), which we will now discuss. We introduce the variables (ρ, y) where ρ is the total density and y = ρ 3 /ρ the asymmetry. Let us focus on the point (ρ = ρ s , y = 0), defined as the high-density border of the spinodal for symmetric matter: the extension of the spinodal region in the asymmetry direction depends on the shape (convexity) of the spinodal contour at this point. This contour is defined by C h < (ρ, y) = 0, where C h < is the lower eigen-value of the free-energy curvature matrix. The isospin symmetry imposes ∂C h < /∂y| ρ (ρ, 0) = 0. The convexity of the spinodal contour at the point (ρ s , 0) can then be characterized by the quantity Cs = ∂ 2 C h < /∂y 2 | ρ (ρ s , 0). Indeed, if Cs is positive, for a small asymmetry y we have C h < (ρ s , y) > 0, meaning that the point (ρ s , y) is outside the spinodal contour. Thus, a positive (negative) value of Cs corresponds to a concave (convexe) shape of the spinodal contour at point (ρ = ρ s , y = 0). The value of Cs depends on the symmetry energy and its density derivatives:

Cs = 4 ρ a s a s (ρ 2 a ′′ s + 2ρa ′ s ) -2(ρa ′ s ) 2 (6)
where a ′ s = [∂a s /∂ρ| y ] (ρ = ρ s , y = 0) and a ′′ s = ∂ 2 a s /∂ρ 2 | y (ρ = ρ s , y = 0). The respective contributions entering in Eq. ( 6) are shown in Tab. II for the three Skyrme-like interactions SLy230a, SGII and RATP. It is clear that none of these terms dominates the sum. Information on the density dependence of the symmetry energy (slope and curvature with respect to ρ) is then determinant to predict the qualitative behaviour of the spinodal. Conversely, information on the evolution of the spinodal extension with the asymmetry would bring constraints on the relation between these terms according to Eq. ( 6).

III. LOCAL EQUILIBRIUM DESCRIPTION OF CLUSTER FORMATION

Let us now turn to finite-size density fluctuations, relevant for the study of matter clusterization in neutron-star crust and supernova core, as well as nuclear multifragmentation. Coulomb and surface energy now contribute to reduce the instability induced by the thermodynamic spinodal. However, a region of finite-size instabilities can still be found, causing the instability of homogeneous matter against cluster formation. This region is determined considering planewave density fluctuations, varying the corresponding wave-number q: indeed, any density fluctuation can be expressed as a coherent sum of plane waves and the instability with respect to at least one q value signs the presence of finite-size instabilities.

This section presents a first approach for the study of the unstable region: we assume local equilibrium in a semiclassical framework, the local densities defining local Fermi momenta. First, an instability criterion is defined from a Thomas-Fermi calculation of the free-energy variation; then hydrodynamical equations are introduced in order to describe the temporal evolution. We detail here how the free-energy variation induced by a small density fluctuation is derived from the Thomas-Fermi approximation. Let us consider a plane-wave density fluctuation of wave number q, affecting both the neutron and proton densities (i = n, p): ρ i = ρ 0,i + δρ i , with δρ i = A i e iq•r + cc. The size of a density fluctuation, λ, is characterized by the momentum q: λ = 2π/q (note that the thermodynamic limit is obtained back with the infinitewavelength density fluctuation, namely q → 0). The Thomas-Fermi approximation supposes that density variations δρ i are smooth enough to allow the definition of a local Fermi momentum. This approximation is valid if the number of states ρ i Ω -where Ω is the volume -is very large for a typical scale of λ. Taking for the volume Ω = λ 3 , this condition leads to the following criterium: 4πg(k Fi /q) 3 /3 ≫ 1, where g is the spin degeneracy g = 2. It is then clear that this approach is valid if q is smaller than the Fermi momentum k Fi [7,36].

In the description of nuclear matter with density fluctuations, finite-size effects appear because of the momentum dependence of the nuclear force, which is also related to its range. In the Skyrme framework, the surface energy can be written as:

H ∇ = C ∇ nn (∇ρ n ) 2 + C ∇ pp (∇ρ p ) 2 + 2C ∇ np ∇ρ n • ∇ρ p , (7)
where coefficients C ∇ ij are combinations of the standard Skyrme parameters:

C ∇ nn = C ∇ pp = 3 32 [t 1 (1 -x 1 ) -t 2 (1 + x 2 )] , C ∇ np = C ∇ pn = 1 32 [3t 1 (2 + x 1 ) -t 2 (2 + x 2 )] . (8)
To get a deeper insight, we can also express H ∇ in terms of the isoscalar and isovector density gradients, namely:

H ∇ = C ∇ 11 (∇ρ) 2 + C ∇ 33 (∇ρ 3 ) 2 (9)
where the isoscalar and isovector coefficients are related to (7) as:

C ∇ nn = C ∇ pp = C ∇ 11 + C ∇ 33 , C ∇ np = C ∇ pn = C ∇ 11 -C ∇ 33 . (10)
The values of the coefficients C ∇ associated with the density-gradient dependence of the Skyrme functional are given in Tab. III for several interactions. We can notice that the main part of the surface energy is due to the isoscalar density-gradient term, associated with the coefficient C ∇ 11 , which is much smaller for SGII than for SLy230a and RATP: consequently, large density gradients should cost less energy with SGII.

The local-Fermi-momentum approximation generates a coordinate-dependent Fermi momentum k F (r). Together with the temperature T , it provides expressions for the local particle, kinetic-energy and entropy densities. Using the Skyrme energy functional and the entropy, we can directly evaluate the free-energy density, which can be decomposed into a bulk and a surface term.

The bulk term contains the bulk part of the energy functional (Eq. ( 1)) and the entropy term. It is obtained from the space average of the free energy density f h (defined in Eq. ( 4)) calculated with the local densities [21]. Performing a development up to the second order in the amplitudes A i associated with the density fluctuations -the first order being null due to particle conservation -the bulk variation of the free-energy is given by:

δf b = i,j A i A * j + A * i A j 2 ∂µ i (ρ 0,n , ρ 0,p ) ∂ρ j , (11)
which depends on the curvature-matrix elements of the homogeneous free-energy C h ij = ∂µ i /∂ρ j .

The surface term is generated by the gradient term of the energy functional. The space average of H ∇ leads to the surface contribution:

δE ∇ = 1 Ω H ∇ (r)dr = q 2 i,j A i A * j + A * i A j C ∇ ij . (12)
Considering finite-size fluctuations, one shall also take into account the effect of the Coulomb interaction. This was not the case in the thermodynamic framework discussed in Sec. II, since it is well known that the Coulomb energy diverges in infinite charged matter. In star matter, electrons equilibrate the positive charges and remove the divergence. In finite nuclei, the Coulomb energy might be large but finite. In both cases, the physically important quantity is not the total Coulomb energy but its variation δE c induced by density fluctuation. In the present calculation, we have neglected the small contribution of the exchange Coulomb interaction so that the Coulomb contribution to the free-energy variation reads:

δE c = 4πe i e j q 2
i,j

A i A * j + A * i A j 2 , (13)
where e p = q e / √ 4πǫ 0 , e n = 0, with q e the elementary electric charge. Finally, the total free-energy variation is given by the summation of the three terms, δf = δf b + δE ∇ + δE c :

δf = i,j A i A * j + A * i A j 2 ∂µ i (ρ 0,n , ρ 0,p ) ∂ρ j + 2C ∇ ij q 2 + 4πe i e j q 2 , ( 14
)
which can be expressed in a matricial form in the space of the density fluctuations à = (A n , A p ): δf = à * C f Ã, where

C f = C h + 2q 2 C ∇ nn C ∇ np C ∇ pn C ∇ pp + 4πe 2 q 2 0 0 0 1 . ( 15
)
The curvature matrix C f is a function of the temperature T , the average densities ρ 0 i , and the wave number of the fluctuation q. It contains the bulk term of the curvature, C h , the density-gradient term of the nuclear interaction proportional to q 2 , and the Coulomb-interaction term proportional to 1/q 2 . These two additional terms being positive tend to reduce the instability induced by the bulk term. Thus, the instability will be quenched at large values of q by the surface contribution, and at small values of q by the Coulomb interaction. It is then interesting to evaluate the curvature matrix for intermediate values of q, where it may still acquire negative values. In the following, the calculations show that there is still a region where the lower eigen-value of C f is negative for a finite interval of q values, leading to the development of spinodal instabilities.

B. Hydrodynamical evolution

The free-energy curvature analysis presented so far provides a static description of the finite-size instabilities. In order to describe the early dynamics of matter clusterization induced by such instabilities, we use a hydrodynamical framework, which corresponds to the local equilibrium limit and so is consistent with the above Thomas-Fermi approximation of δf. Let us first consider a single fluid described by the Euler and continuity equations [7]:

d(mρv)/dt = -∇P , ( 16
) ∂ρ/∂t + ∇(ρv) = 0 , ( 17
)
where P is the system's pressure, m is the particle mass and ρ the particle density. In the limit of small-amplitude fluctuations ρ = ρ 0 + δρ, Eqs. (16)(17) can be transformed into a diffusion equation:

∂ 2 ρ ∂t 2 = 1 m dP(ρ 0 ) dρ ∇ 2 ρ = ρ 0 m d 2 f(ρ 0 ) dρ 2 ∇ 2 ρ . ( 18
)
which is obtained using the relation between the pressure and the free-energy density f, dP/dρ = ρd 2 f/dρ 2 . Expanding Eq. ( 18) on a plane-wave basis gives the following dispersion relation for the density modes:

ω 2 = ρ 0 m d 2 f(ρ 0 ) dρ 2 q 2 . ( 19
)
If the free-energy curvature d 2 f/dρ 2 is negative, the frequency is purely imaginary, ω = iγ, and induces an exponential growth of the fluctuation with typical time τ = 1/γ. In asymmetric unstable nuclear matter, we generalize Eq. ( 19) introducing fluctuations along the unstable direction of the curvature matrix. This leads to :

ω 2 = -γ 2 = ρ m C f < 2 q 2 , ( 20
)
where m is the nucleon mass, ρ the total average density and C f < the lower eigen-value of the free-energy curvature matrix in the space of density fluctuations. The factor 1/2 is due to the change of variables from (ρ, ρ 3 ) to (ρ n , ρ p ). Equation ( 20) characterizes the properties of the finite-size instabilities in terms of the free-energy curvature matrix, in the framework of a semi-classical and local equilibrium approach. In the following, it will be denoted as the hydrodynamical approach [45].

IV. COLLISIONLESS DYNAMICS OF CLUSTER FORMATION

Let us now consider a fully quantal approach for the description of the growth of density instabilities. This approach is based on the linear response theory to unstable fluctuations. In the hydrodynamical approach, a local equilibrium was assumed: in microscopic words, it means that for each particle species a Fermi sphere centered on zero momentum was defined in each point of space. Now, in contrast, each particle species is associated with a single Fermi sphere affected by a momentum shift (particle-hole excitations) causing density fluctuations in coordinate space. These two limits correspond respectively to first and zero sound.

The linear response function is obtained from the small-amplitude limit of the TDHF theory (Time-Dependent Hartree-Fock): this corresponds to the RPA (Random-Phase Approximation) [START_REF] Ring | The Nuclear Many-Body Problem[END_REF]. In the present work we use the techniques developed in Refs. [38,39] to obtain the RPA response at finite temperature, in the peculiar case of Skyrmetype effective interactions. We first solve the RPA in a fully quantal framework. We then introduce the semi-classical limit of the RPA, namely the Vlasov approach, in order to identify the quantum effects in the comparison between hydrodynamics and RPA results.

A. Quantal effects described in the HF+RPA framework

In the quantal framework, we have to linearize the TDHF equation, which describes the time evolution of the onebody density matrix ρ ij = ĉ † j ĉi , ĉ( †) i being the annihilation (creation) operator associated with the single-particle state i:

i ∂ ρ ∂t = [ ĥ + Vext , ρ] , ( 21
)
where ĥ is the mean-field operator and Vext is an external field. In the absence of perturbation, the mean-field ground state is given by ρ0 and ĥ0 , which is the Hartree-Fock solution satisfying [ ĥ0 , ρ0 ] = 0. In the following, we consider spin-saturated systems. For infinite nuclear matter, the operators ĥ0 and ρ0 are diagonal in momentum space and we have:

ρ0 |k = gf k |k ĥ0 |k = ω k |k ( 22
)
where g is the spin degeneracy, and f k the Fermi-Dirac occupation for a level of momentum k.

In this section, we consider a large system of volume Ω and the plane waves |k are defined such that r|k = e ik•r / √ Ω.

Linear response to an external probe

We calculate the response function of nuclear matter at finite temperature to an infinitesimal external field of the form:

Vext = Ee iq•r-i(ω+iη)t (t < 0) ; Vext = 0 (t ≥ 0) , ( 23
)
where η is a positive infinitesimal part ensuring the adiabatic switching of the external field. The excitation carries a momentum q and an energy ω. For simplicity, we first treat nucleons as a single-component fluid. At late enough times and small enough external field E, the expectation value of the density fluctuation has the same space and time dependence as the external field:

δρ(r, t) = αe iq•r-i(ω+iη)t . ( 24
)
The linearization of Eq. ( 21) leads to the polarization function:

Π(ω, q) = α E = Π 0 (ω, q) 1 -Π 0 (ω, q)v res . ( 25
)
In this expression, we have introduced the Lindhard function Π 0 (corresponding to the response function in the absence of residual interaction):

Π 0 (ω, q) = g (2π) 3 dk f k -f k+q ω -ω kq + iη ; ω kq = ω k+q -ω k , (26)
and the residual interaction v res which accounts for the variation of the mean field δ ĥ due to the density fluctuation. Details for the derivation of the response function (25) are given in Appendix A.

A collective mode (ω, q) of the nuclear fluid corresponds to a large enhancement of the response function, or a small value of the denominator of Eq. ( 25). In the spinodal unstable zone, the response function diverges for a given value of ω and q. The dispersion relation ω(q) for the unstable mode is then defined by (1 -Π 0 (ω(q), q)v res = 0).

Let us now take into account the isospin degree of freedom. In this case, we distinguish neutron and proton fluids, and we get the response function and the residual interaction in terms of the matrices:

Π 0 = Π 0,nn 0 0 Π 0,pp , (27)
v res = v res nn v res np v res pn v res pp , (28)
Π = (1 -Π 0 v res ) -1 Π 0 . (29)
The linear response Π(ω, q) presents a divergence for the zero of matrix (1 -Π 0 v res ), corresponding to the condition:

det (1 -Π 0 v res ) = 0 . ( 30
)
2. Residual interaction and free-energy curvature

We now have to express the residual interaction v res entering the definition of the response function (29). We give in Appendix B the steps demonstrating the link between v res and the free-energy curvature matrix C f discussed in Sect. III. To summarize, considering a density fluctuation

δρ i (r) = q ρ i,q ; ρ i,q = A i,q e iq•r (i = n, p) , (31)
we have the average residual-energy variation per unit volume:

δE res = q i,j A i;q A * j;q 2 v res ij (q) (i, j = n, p) . ( 32
)
The following expression is obtained for the residual-interaction elements as a function of the free-energy curvature:

v res ij (q) = C f ij (q) - δ ij N 0,i , (33)
where N 0,i is the density of states at Fermi level for the particle species i. This corresponds to the matrix expression:

C f = (N 0 ) -1 + v res ; N 0 = N 0,n 0 0 N 0,p . (34)
We remind that according to the semi-classical free-energy criterion discussed in the previous section, an instability occurs when the matrix C f has a negative eigen-value. The corresponding spinodal region is circled by a frontier for which the lower eigen-value is zero, responding to the condition:

det C f = det (N 0 ) -1 + v res = 0 . (35)
Now, the RPA modes are defined by the condition (30), which can be expressed in a closely related form as:

det C RPA = det (N 0 ) -1 + P 0 v res = 0 , (36)
where

P 0 (ω, q) = -(N 0 ) -1 Π 0 (ω, q) . ( 37
)
Note that the matrix C RPA which has just been introduced, in contrast with any curvature matrix, is not symmetric in the general case.

In the RPA framework, the spinodal frontier is defined by the occurrence of a zero transfered energy mode. For T = 0, in the limit q → 0 the uncorrelated response function then reduces to [START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF]:

lim q→0 Π 0,i (ω = 0, q) = -N 0,i . (38)
According to Eqs. (35)(36)(37), this implies that, for T = 0 and q → 0, the condition defining the RPA spinodal frontier reduces to the semi-classical one, det C f = 0. Note however that, taking into account the Coulomb interaction, the instability interval of q has a lower limit, so that the RPA and semi-classical spinodal frontiers always remain separated.

Unstable modes in the RPA framework

In the RPA approach, the instability region is characterized by the occurrence of imaginary modes: ω = ω r + iγ. The imaginary term iγ leads to an exponential increase of the fluctuation amplitude with the typical time τ = 1/γ. Considering the general expression of ω in the complex plane, the uncorrelated response function expressed by Eq. ( 26) can be written (for each particle species i):

Π 0,i (ω, q) = g (2π) 3 dk f i,k -f i,k+q ω -ω kq = 2g (2π) 3 dk f i,k ω kq ω 2 -ω 2 kq . ( 39
)
The condition (30) can be satisfied only if the imaginary part of Π 0 cancels, which according to Eq. (39) implies that ω 2 is real. Consequently, unstable modes correspond to a purely imaginary energy such that ω = iγ. The resulting response function reads:

Π 0,i (γ, q) = - 2g (2π) 3 dkf i,k ω kq γ 2 + ω 2 kq . (40)
At zero temperature, this integral has an analytic solution:

Π 0,i (γ, q) = -N 0,i P 0,i (γ, q) (41)
where, introducing the quantity w 2 i = γ/[ /2m * i ] (w i has the dimension of a momentum), P 0,i is given by: P 0,i (γ, q) = P 0,i (w i (γ), q) = 1 2k Fi q

kF i 0 kdk ln (2kq + q 2 ) 2 + w 4 i (2kq -q 2 ) 2 + w 4 i .
Introducing the reduced variables ν i = w 2 i /2qk Fi = m * i γ/(qk Fi ) and q Fi = q/2k Fi , we obtain the analytical expression: P 0,i (γ, q) = P 0,i (ν i (γ, q), q Fi (q))

= 1 2 1 + 1 + ν 2 i -q 2 Fi 4q Fi ln (1 + q 2 Fi ) 2 + ν 2 i (1 -q 2 Fi ) 2 + ν 2 i -ν i arctan 1 + q Fi ν i + arctan 1 -q Fi ν i (42)
which verifies, consistently with Eq. (38), the condition at the spinodal frontier:

lim q→0 P 0,i (γ = 0, q) = 1 . (43)
The RPA dispersion relation for the unstable mode γ(q) is deduced from Eqs. ( 36) and (42). In this quantal approach, for given temperature and densities, uniform matter is unstable against a fluctuation of finite momentum q if det C RPA (q) (γ) presents a zero for a given value of the typical time τ = 1/γ. The corresponding vector gives the direction of the growing fluctuation in the density plane.

At finite temperature, the same procedure can be followed but the response function defined by Eq. ( 40) needs numerical calculation. Two different methods have been considered: one involves the calculation of Fermi integrals [START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF] and the other consists in performing a pondered sum of zero temperature response functions [38]:

Re Π 0,i (γ, q, T ) = - k Re Π 0,i (γ, q, T = 0, k F = k)df i (k, T ) . ( 44
)
At zero temperature, df i (k) = -δ(k -k Fi )dk, yielding the well-known expression for Re Π 0,i . Both methods have been used in order to check the accuracy of our results at finite temperature.

B. Vlasov approach

We now introduce the semi-classical limit of the RPA. The Vlasov approach describes the matter response to a momentum transfer q and energy transfer ω like the RPA approach, but in a semi-classical approximation. It is then an intermediate description between the hydrodynamics and the full quantal RPA. Again, for simplicity, we give expressions for the case of a single fluid. The semi-classical version of the TDHF equation is obtained in the limit → 0 [START_REF] Ring | The Nuclear Many-Body Problem[END_REF]. Using a Wigner transform, the commutator divided by i is then replaced by a Poisson bracket, where the one-body density and potential operators are expressed by functions of time, position and momentum (t, r, k). This gives the Vlasov equation, describing the semi-classical time evolution in the collisionless limit:

∂f(t, r, k) ∂t = {h(t, r, k) + V ext (t, r, k) , f(t, r, k)} , ( 45
)
where f is the occupation function in phase space, h the mean field and V ext an external potential. In the unperturbed state, the mean field is h 0 (k) and the occupation at temperature T = 1/β is given by the Fermi-Dirac distribution

f 0 (k) = [1 + exp(β(h 0 (k) -µ))].
The density fluctuation corresponds to the occupation variation f(t, r, k) = f 0 (k) + δf(t, r, k), inducing the mean-field variation h(t, r, k) = h 0 (k) + δh(t, r, k).

Let us stress here the difference between Thomas-Fermi and Vlasov approaches. Both approaches give a semiclassical description of the state occupations, using the function f(t, r, k): however, they define the perturbation δf(t, r, k) in different ways. In the Thomas-Fermi approach, the variation of the local density is associated with a variation of the chemical potential, according to the local equilibrium description, and the occupation function is given by a Fermi-Dirac distribution: f(t, r, k) = [1 + exp(β(h 0 (k) -µ(t, r)))]. In terms of the density matrix, this corresponds to a variation of the Fermi surface. In contrast, in the Vlasov formalism f(t, r, k) is given by a canonical transformation of f 0 (k): the variation of the distribution can be expressed in terms of a generating function S(t, r, k) such that δf = {S, f 0 } [44]. In terms of the density matrix, this corresponds to the excitation of particle-hole states with a finite momentum transfer.

The linearized Vlasov equation (semi-classical limit of the RPA) can be solved following a method similar to the one exposed in Ref. [44]. The Vlasov eigen-modes for asymmetric nuclear matter obey the matrix relation:

det C Vlasov (ω, q) = det (N 0 ) -1 + Lv res 2 = 0 , (46)
where (with i = n, p):

L(s) = L(s n ) 0 0 L(s p ) ; s i = ω q dh 0i ({k i = k Fi }) dk -1 ; L(s i ) = 2 -s i ln 1 + s i 1 -s i ; (47)
L(s i ) denotes the semi-classical Lindhard function. The low-q limit of the uncorrelated response function Π 0,i (ω, q) is related to L(s i ) in such way that we get in terms of matrixes:

lim q→0 (-N 0 Π 0 (ω, q)) = lim q→0 P 0 (ω, q) = L/2 . ( 48)

RPA

Hydro. SLy230a

Vlasov FIG. 3: (color online) Dispersion relation at a given density and asymmetry, for zero temperature. Dots represent the boundaries of the unstable interval of q. Top: 1/τ is given as a function of the wave number q. Bottom: associated directions of phase separation in the density plane δρ3/δρ. The horizontal line at 0.4 is the direction of constant proton fraction ρ3/ρ. Left: three approaches applied with SLy230a: RPA (thick line), Vlasov (dashed line), hydrodynamics (thin line). Right: RPA and hydrodynamics with different Skyrme forces.

This means that the RPA eigen-modes obeying Eq. ( 36) reduce to the Vlasov ones in the low-q limit (whatever the value of ω). On the other hand, for ω = 0 we have L(0)/2 = 1, so that C Vlasov (ω = 0, q) reduces to the freeenergy curvature matrix C f (q): consequently, the frontier of the spinodal region is the same for the two semi-classical approaches. On this common frontier, the Vlasov and hydrodynamics eigen-modes are exactly identical. Differences shall appear inside the spinodal region where ω = 0.

V. DISCUSSION OF THE RESULTS

In this section, we present the properties of matter clusterization obtained from the different approaches presented above. The quantal RPA is compared to two semi-classical approaches: Vlasov and hydrodynamics, corresponding respectively to the collisionless (zero sound) and local equilibrium (first sound) limits. We first consider instability against density fluctuations for several momenta q, exploring the dependence on the temperature and on the nuclear interaction. Then, we concentrate on the most unstable mode defined as the mode for which the typical instability time is minimum and explore the clusterization properties.

A. Unstable eigen-modes

Let us determine the different instability properties associated with a given wave-number q. The unstable eigenmodes obey the dispersion relations illustrated for fixed densities on Fig. 3. The upper part of this figure represents the inverse typical time 1/τ (q), which characterizes the time evolution of the modes. The quantum effects can be directly appreciated from the comparison between RPA and Vlasov results. Although these results are quasi-similar at the low-q limit, as expected, they separate for higher values of q. We then observe a quantum quenching of the instability, as discussed in Ref. [36]: the RPA unstable modes are associated with lower values of 1/τ (i.e. a slower amplification of the fluctuation) and they end at a lower value of q. Indeed, in the Vlasov approach the uncorrelated response function is approximated by its limit at low q, which amounts to neglecting a part of the kinetic-energy cost of the fluctuation. In the RPA where the exact value of Π 0 (ω, q) is calculated, the high-q density fluctuations are disfavored. Let us now consider the semi-classical approach based on local equilibrium (hydrodynamics). It presents the same qualitative features as the previous ones: in particular, as discussed earlier, the instability region corresponding to a negative value of the free-energy curvature is formally the same as for the Vlasov approach. However, the hydrodynamical dispersion relation calculated from the curvature C f < leads to a much faster evolution, 1/τ being typically twice as large as for the RPA. It also presents a less symmetric behavior, with a sharper slope toward the high-q border: this is due to the relation 1/τ 2 ∝ |C f < |q 2 , in which the factor q 2 enhances the high-q side of the dispersion relation. The lower part of Fig. 3 deals with the direction of phase separation, also included in the dispersion relation determining the eigen-modes. This feature is connected with the isospin distillation, which is a well-known aspect of the liquid-gas phase transition in asymmetric nuclear matter [42]: in order to minimize the symmetry energy, the dense phase is more symmetric than the dilute one. This is observed in the construction of the thermodynamic phase equilibrium, and reflected by the spinodal instability direction: instead of following a line of constant composition in the density plane, phase separation is rotated towards the isoscalar direction. For a neutron rich point of densities (ρ, ρ 3 ), for which phase separation direction is given by the vector coordinates (δρ, δρ 3 ), the normal isospin distillation is reflected in the relation:

0 < δρ 3 /δρ < ρ 3 /ρ (49)
We have checked that this relation is verified by the thermodynamic spinodal instability for all the Skyrme-like interactions we consider in this article. Moreover, it has been shown in Refs. [43] that for a set of realistic Skyrmelike and Gogny interactions, the instability is still isoscalar-like, which means that δρ 3 /δρ < 0.5 in a domain of reasonable isospin asymmetries. Now when finite-size instabilities are considered and the electric charge comes into play, the Coulomb interaction disfavors proton-density fluctuations: this goes against the normal distillation effect (for neutron-rich matter). As can be seen on Fig. 3, isospin distillation is particularly affected at low values of q, where Coulomb effects dominate: an inverse distillation is even found near the low-q border of the instability interval. However, the normal isospin distillation remains the dominant feature, especially if we consider the most unstable modes [21,44]: it leads to the formation of clusters more symmetric than the homogeneous initial system, immersed in a more neutron-rich gas. Figure 3 shows a very close agreement between Vlasov and RPA results as for the value of δρ 3 /δρ. The direction corresponding to the hydrodynamical approach is nothing but the direction of minimal free-energy curvature, i.e. the eigen-vector u < associated with the lower eigen-value of the curvature matrix C f . Qualitatively, it presents a similar evolution; furthermore, once again hydrodynamics and Vlasov curves are constrained to join at the borders of their common instability interval. However, inside this interval the hydrodynamics curve is clearly distinct from the two collisionless approaches. This shows that the phase separation direction followed in a collisionless approach is quantitatively different from the direction u < of the negative free-energy curvature: it leads to a slightly reduced distillation effect.

The right column of Fig. 3 shows the dependence of the dispersion relation on the Skyrme force. RPA and hydrodynamics results are given for SLy230a, SGII and RATP. All present similar features, ensuring that the results elsewhere shown with SLy230a alone are qualitatively representative. The quantitative differences between SLy230a and RATP are weak. SGII however favors modes of larger q-values, and presents an extended instability interval on the high-q side. This difference comes from the surface properties of the forces: as can be seen in Tab. III, the energy cost for the introduction of a density gradient is lower with SGII.

Let us now address the extension of the instability in the density plane. We first consider instability regions associated with given values of the transfered-momentum q. In a condensed way, we name such regions "q-spinodals". Inside a q-spinodal contour, a density fluctuation of momentum q is unstable and amplified according to the dispersion relation discussed above. Within the RPA approach, the q-spinodal contour at temperature T is determined by the divergence of the linear response Π(ω = 0, q, T ). In the semi-classical description (for both Vlasov and hydrodynamics), the q-spinodal contour corresponds to a zero eigen-value of C f (q, T ). As shown by Eq. ( 38), both curves are equivalent at the limit q → 0. Figure 4 represents different RPA and semi-classical q-spinodals at T = 0. For q = 50 MeV/c almost no difference can be seen in the density plane. On the lower part of the figure, the q-spinodals are drawn in the plane of Fermi momenta, k Fi = (3π 2 ρ i ) 1/3 : we can then check that both approaches become very close for q ≪ k Fi . For higher values of q, the RPA criterion always leads to a reduction of instability compared with the semi-classical approach, consistently with the previous discussion of the quantum quenching affecting the dispersion relation (Fig. 3). Figure 5 represents the global region of instability against clusterization. For each approach, and for a given temperature, it corresponds to the envelope of all the q-spinodals: inside this envelope, nuclear matter can be found unstable against finite-size fluctuations for at least one q-value. The right par of Fig. 5 compares the RPA instability envelopes for different Skyrme-type forces: the slight differences obtained reflect those appearing between the thermodynamic spinodal curves of Fig. 2. Since the widest q-spinodal contour is obtained for q around 100 MeV/c, the corresponding graph of Fig. 4 is representative of the difference between the global instability regions obtained in RPA and semi-classical approaches at T = 0. On the left part of Fig. 5 are shown the instability regions obtained at T = 5 and 10 MeV for both approaches. The difference between the two envelopes is sensitively enhanced at higher temperature, as we approach the limiting temperature of the instability: this means that this limit occurs earlier in the RPA due to the quantum quenching. It is then clear that the temperature range concerned by the liquid-gas instabilities is not sufficient for a suppression of the quantum effects.

The temperature dependence of the instability regions is illustrated in Fig. 6 for symmetric matter, considering the thermodynamic and finite-size (semi-classical and RPA) spinodals. We can observe the reduction of the thermodynamic instability, and the smaller temperature range of the RPA instabilities with respect to the semi-classical ones. The temperature and density coordinates of the highest-temperature points appearing on this figure are reported in Tab. IV for the different Skyrme forces. In the thermodynamic case, the phase-transition region ends in a critical point occurring for symmetric matter; its coordinates are (ρ c , T c ). In contrast, the finite-size instability region does not reach its limiting temperature on the axis of symmetric matter, the charge-exchange symmetry being broken by the Coulomb interaction. Instead, the limiting point is given by the coordinates (ρ lim , (Z/A) lim , T lim ) (associated with the last unstable momentum q lim ). These coordinates are given in Tab. V, in the semi-classical and RPA cases, for the different Skyrme forces. Let us first compare the RPA and semi-classical limiting points. The limiting temperature is reduced by several hundreds of keV due to the quantum quenching of the instability in the RPA approach. The position of the limiting point is found to be at slightly higher density for the RPA: however, this feature is not very pronounced and is not systematic at given Z/A, as can be seen on Fig. 6 for Z/A = 0.5. However, it is noticeable that the RPA limiting point occurs at higher asymmetry: this can be attributed to the larger influence of Coulomb effects, since, as expected, the corresponding transfered momentum is lower. Let us now compare the critical and limiting temperatures. For the forces with similar surface properties (e.g., in the presented results, SLy230a and RATP), the difference between the values of T lim reflects the difference between the thermodynamical critical temperatures T c (typically reduced by 2.4 and 2.8 MeV respectively in the semi-classical and RPA cases). A smaller reduction of the critical temperature is observed for the force with reduced surface energy (SGII), for which the last unstable momentum q lim is higher.

We finally explore on Fig. 7 the behavior of the phase separation direction δρ 3 /δρ associated with the unstable eigen-modes, for discrete values of the momentum q. RPA and hydrodynamics results are compared. The viewgraphs represented on this figure are separated into two groups: on the left-hand side, the total density is fixed at ρ = 0.05 fm -3 and the horizontal axis explores the isospin asymmetry ρ 3 /ρ; on the right-hand side, we fix the proton fraction Z/A = 0.3 (i.e. ρ 3 /ρ = 0.4), and the horizontal axis explores the total density. In both cases, the instability direction is represented for different Skyrme forces (top), different q values (middle), and different temperatures (bottom). First, it could be remarked from the left view-graphs that the different curves are mostly situated inside the cone defined by the horizontal line and the first diagonal, satisfying the inequalities (49). This confirms the predominance of the normal isospin distillation, already noticed about Fig. 3.

(Z/A) SC lim (Z/A) RPA lim Tc T SC lim T RPA lim q SC lim q RPA lim [fm -3 ] [fm -3 ] [fm -3 ] [MeV] [MeV] [MeV] [MeV/c] [MeV/c] Sly230a
Let us now consider the right panels of Fig. 7. As expected, the difference between RPA and hydrodynamics results tends to disappear at the high-density border of the unstable zone, as the Fermi momenta become large enough with respect to the transfered momentum q. Indeed, in this case the RPA reduces to the Vlasov approach, which identifies with the hydrodynamics on the borders of the instability interval. Inside the spinodal zone, the differences are enhanced, especially at T = 0: in the hydrodynamical approach, the instability direction δρ 3 /δρ strongly decreases towards lower densities, while in the RPA approach its density dependence is weak. As a consequence, the RPA values of δρ 3 /δρ are larger (up to a factor 2). This accentuates the observation made on Fig. 3, that the distillation effect is reduced in the RPA framework. However, both approaches again confirm the isoscalar behavior of the unstable mode. We finally notice that the differences are smoothed for the results at finite temperature: this is essentially due to the reduction of the instability, limiting the range of directions leading to a negative free-energy curvature, rather than to the disappearance of the quantum effects. Indeed, we know from the comparison between RPA and Vlasov curves on Fig. 3 that these effects do not affect the separation direction at least for ρ ≥ 0.05 fm -3 . Furthermore, as discussed previously, quantum effects are still present in all the range of temperature concerned by the liquid-gas instabilities.

B. Dynamical properties of cluster formation

In the spinodal decomposition scenario, the early dynamics of cluster formation is induced by the most unstable mode, corresponding to the fastest amplification of the fluctuation. For fixed temperature and densities, this is given by the maximum of the dispersion relation 1/τ (q), that we denote (q 0 , 1/τ 0 ). Clusters should then be formed according to the timescale τ 0 , with a size of the order of λ 0 /2 = π/q 0 , and a composition related to the corresponding direction (δρ 3 /δρ) 0 .

We show on Fig. 8 the projection on the density plane of the iso-contours associated with the typical time constant τ 0 (top) and fluctuation size λ 0 /2 (bottom). Hydrodynamics and RPA results are shown, respectively on the left and right panels. The instability envelope corresponds to τ 0 = ∞. Inside this envelope, the region of fastest time evolution corresponds to the lowest values of τ 0 and is found for low-density symmetric matter (ρ ≃ 0.03 fm -3 for both approaches). The RPA approach has typical times twice larger than the hydrodynamical ones. The half-wavelength FIG. 7: (color online) Directions of phase separation for fixed values of the wave number q studied along two different axis of the density plane, for the hydrodynamical (thin lines) and RPA (thick lines) approaches. Dots give the limit of instability regions. Left: constant density. The horizontal line corresponds to the isoscalar direction, and the first diagonal to the constant-Z/A direction. Right: constant proton fraction. Top figures compare different Skyrme forces; middle figures different values of q; bottom figures different temperatures. λ 0 /2 is larger in the quantal approach than in the semi-classical one. These results generalize the features observed on the dispersion relation 1/τ (q) illustrated on Fig. 3, i.e. the quantum quenching of instabilities and the higher instability in q and 1/τ of the first sound compared to zero sound.

Let us stress here that, although larger time constants and larger cluster sizes are repeatedly associated in the comparison between the different approaches, these two features are actually decoupled. Indeed, τ 0 and λ 0 are not linked by a proportionality relation, since they are determined by the minimum-τ relation (dτ /dq) 0 = 0. In fact, the observed correlation is due to the particular shapes of the respective dispersion curves 1/τ (q). It can even be observed on the low-density part of the hydrodynamics plots that longer times can be associated with shorter wavelengths, stressing the absence of proportionality relation between those two quantities.

It is interesting to give an estimate of the cluster masses associated with the favored wavelengths represented on the bottom part of Fig. 8. Let us first notice that the precise determination of this quantity is beyond the scope of the present approach, which describes the initial dynamics of fragment formation; the properties of the fully separated fragments in the final state coming from the spinodal decomposition should be obtained following the evolution of the density fluctuation beyond the linear regime. However, an order of magnitude can be obtained in the present framework assuming that the half-wavelength λ 0 /2 of the favored mode could give an estimate of the nucleus diameter at saturation density. In Tab. VI are listed some values of the half-wavelength λ 0 /2 of the favored mode and the deduced mass of produced nuclei. Comparing the typical values of the half-wavelength λ 0 /2 shown on Fig. 8 with the masses obtained in Tab. VI, it could be inferred that the spinodal decomposition scenario predicts the formation of light nuclei (up to A≃ 20) through most of the unstable region. Medium-mass nuclei could be formed in the region close to the frontiers of instability (large global density or neutron excess): their formation would thus require much more time than that of lighter elements.

On Fig. 9, we explore the instability direction (δρ 3 /δρ) 0 through the density plane, for RPA and hydrodynamics. The top view-graphs display level-lines of the quantity (δρ 3 /δρ) 0 . A different kind of representation is employed for the bottom view-graphs: they give a direct picture of the instability direction in the plane, by the so-called tangential lines. In each point of a tangential line, the tangent to this curve gives the direction of the most unstable mode. This representation clearly shows the isoscalar character of the modes obtained in both approaches.

Let us now comment the level-line panels of Fig. 9. The (δρ 3 /δρ) 0 level corresponding to the value 0 is the purely isoscalar mode. In the absence of Coulomb interaction, this level strictly corresponds to the axis of symmetric matter [43], as imposed by the isospin-exchange symmetry. In the present case, we have already noticed that the Coulomb interaction disfavors proton-density fluctuations. This results in a global increase of the (δρ 3 /δρ) 0 values, thus the levels are shifted towards the proton-rich side. This effect also appears on the left part of Fig. 7: as expected, the Coulomb shift is increased for larger values of q (middle graph). Coming back to Fig. 9, let us compare the hydrodynamics and RPA upper panels. Two differences can be noticed: i) in RPA the absolute values of (δρ 3 /δρ) 0 are larger, and ii) level-lines on the neutron-rich side approximately follow constant-Z/A lines while hydrodynamics levels are incurved. These panels thus give a global view of the features which have already appeared on Figs. 3 and7, namely a reduced distillation in RPA, with a weaker density-dependence of (δρ 3 /δρ) 0 at fixed Z/A.

Figure 10 displays the clusterization properties obtained along selected axis, for different Skyrme forces, at T = 0. Time, size and direction characterizing the most unstable mode are addressed, comparing RPA and hydrodynamics results. The inverse typical times 1/τ 0 differ typically by a factor of 2 between the two approaches. The cluster sizes are generally increased by about two fermis in the RPA approach. At very low density however, much shorter wavelength are favored in the hydrodynamical approach: note that the validity of the Thomas-Fermi approximation is under question in this domain. Finally, the curves giving the direction of the most-unstable mode are very similar to those obtained on Fig. 7 for q = 100 MeV/c: this is due to the weak dependence of the instability direction δρ 3 /δρ on the momentum q for large enough q values (see Fig. 3). Let us now compare the different Skyrme forces. Similar results are obtained with the three parametrizations. SGII presents faster time evolutions and favors shorter cluster sizes, as was already visible on Fig. 3. As for the direction of phase separation, we can notice that the isospin-distillation effect is slightly weaker with SLy230a.

Finally, we explore on Fig. 11 the temperature dependence of the clusterization properties. Time, size and direction characteristics are shown for T = 0, 5 and 10 MeV, comparing now the three approaches: RPA, Vlasov and hydrodynamics. The differences between RPA and Vlasov results give an estimate of the quantum effects. As can be seen from the size characteristics, quantum effects are still sensible at T = 10 MeV: the Vlasov curve for q 0 is clearly distinct from the RPA one. The smoothing of the differences of clusterization properties between the three approaches with increasing temperature is mainly due to the reduction of the instability. The range of q and δρ 3 /δρ values corresponding to a negative free-energy curvature is narrowed, constraining the different approaches to a closer agreement as for the size and direction associated with the most unstable mode.

The time-evolution part of Fig. 11 deserves a special comment. As T increases, all values of 1/τ 0 go towards zero, which defines the limit of the instability region. At T = 10 MeV, Vlasov and RPA results for this quantity are quite close: however, if we still increase the temperature, the RPA curve will disappear first, then Vlasov and hydrodynamics curves will reach zero together. While the 1/τ 0 values converge to zero, the corresponding divergence of τ 0 on the contrary magnifies the differences between the three approaches. As a result, the range of prediction as for the time-scale of cluster formation is strongly increased as we approach the limit of the instability region.

VI. CONCLUSION

We have presented the clusterization properties of asymmetric nuclear matter in three approaches, comparing the quantal RPA results with those obtains in a semi-classical framework, with Vlasov and hydrodynamical calculations. In all cases, starting from homogeneous infinite matter, we have studied the consequences of introducing a planewave density fluctuation of finite momentum q, working in the small-amplitude limit. Liquid-gas instabilities have thus been identified, corresponding to an amplification of the fluctuation. The unstable eigen-modes (defined for each approach) are characterized by a time constant τ , a wavelength λ = 2π/q, and a direction δρ 3 /δρ. Assuming a scenario of spinodal decomposition, the early dynamics of cluster formation is dominated by the most unstable mode (for given temperature and densities). The clusterization properties are then characterized by the quantities associated with this mode: a time-scale τ 0 , a typical cluster size λ 0 /2 = π/q 0 , and a cluster composition induced by the direction (δρ 3 /δρ) 0 . We have noticed that this study, which assumes fragment formation out of equilibrium, is related to the instability properties of homogeneous matter at sub-saturation density; it is complementary to a description of equilibrated, clusterized matter.

Considering infinite nuclear matter is an approximation which permits to treat analytically in a common framework the liquid-gas instabilities concerning both compact stars and the nuclear multifragmentation observed in heavy-ion collisions. In this framework, the differences between the presented approaches can be attributed to two well-identified origins: comparing Vlasov and RPA gives a direct estimate of the quantum effects, while comparing Vlasov and hydrodynamics shows the difference between a collisionless and local equilibrium description.

The results obtained with the different approaches show similar qualitative features, reproduced with the three Skyrme-like forces that we have used (SLy230a, SGII, RATP). In particular, the regions of instability against finitesize density fluctuations are very close, and the direction of phase separation is always found to be essentially isoscalar, favoring the phenomenon of isospin distillation which leads to more symmetric clusters. An order of magnitude of the cluster masses, estimated as spherical nuclei of diameter λ 0 /2, indicates that light to medium nuclei (A < 40) should be favored through most of the instability region for the various approaches. However, some quantitative differences are noticeable. As for the comparison between the Skyrme forces, the main distinction concerns SGII, associated with shorter size and time characteristics. Let us now summarize the observations concerning the three different approaches. The comparison between RPA and Vlasov shows a quantum quenching leading to larger size and time characteristics in the RPA framework. This effect survives for all the temperature range of the instability. The Vlasov results are framed between the RPA and hydrodynamics ones, which can be considered as giving two opposite limits of predictions. The physical processes should correspond to an intermediate situation, involving different possible relaxation rates. This is an important aspect to take into account in realistic calculations aiming at a quantitative description of nuclear fragmentation. Indeed the two limits we have obtained in the present work indicate that the interval of predictions can be quite large. At T = 0, the inverse time constant 1/τ 0 is typically enhanced by a factor 2 between RPA (slower) and hydrodynamics (faster). The cluster size is about 2 fermis smaller for the hydrodynamics, with a stronger reduction at low density where the validity of the semi-classical approaches gets more questionable (due to the low Fermi momenta). The direction of the unstable modes, for a given value of q, is little sensitive to the quantum effects. The main difference is observed between collisionless and local equilibrium approaches: inside the instability region, the direction corresponding to a divergence of the response function is slightly less isoscalar than the direction of minimal free-energy curvature.

Approaching the frontiers of the instability region (in densities or temperature), the predictions between all approaches get closer as for the size and direction induced by the most unstable mode, due to a reduction of the instability interval of q and δρ 3 /δρ values. The time constant τ 0 characterizing the rapidity of cluster formation has a specific behavior, since it diverges (by definition) on the borders of the instability region: the range of values obtained with the different approaches is then strongly enhanced. This effect is important to notice, since in heavy-ion collisions as well as core-collapse supernovae the time-scale of cluster formation should be compared with that of various equilibration processes. Conversely, a deeper study of the different time scales involved could indicate the respective validity of collisionless and local equilibrium descriptions of matter clusterization.

Let us now consider a finite-size density fluctuation. In order to obtain the expression of the residual interaction in terms of the free-energy curvature matrix, we still work in the local equilibrium approach. Following Eq. (B5), the bulk part of the residual energy variation per unit volume reads (dropping the terms beyond second order in δρ): δρ(r) = q ρ q ; ρ q = A q e iq•r ; A -q = A * q , (B12)

we have:

δE b,res = 1 2 q |A q | 2 v b,res . (B13)
Adding a transfered-momentum dependence to the residual interaction, we generalize the previous expression to:

δE res = 1 2 q |A q | 2 v res (q) . ( B14
)
The q-dependence concerns the gradient and Coulomb terms of the energy variation:

δE ∇ = 1 Ω drH ∇ (r) = 1 Ω drC ∇ (∇ρ(r)) 2 = q |A q | 2 q 2 C ∇ (B15) δE c = 1 2 q |A q | 2 4πe 2 i q 2 (B16)
where e i = 4πq 2 i /ǫ 0 , q i being the particle electric charge for the considered single-component fluid. We can now express v res (q): δE res = δE b,res + δE ∇ + δE c = 1 2 q |A q | 2 v b,res + 2C ∇ q 2 + 4πe 2 i q 2 (B17) v res (q) = v b,res + 2C ∇ q 2 + 4πe 2 q 2 = ∂µ ∂ρ -1 N 0 + 2C ∇ q 2 + 4πe 2 i q 2 . (B18)

The whole procedure can be easily generalized to the two-component nuclear matter, defining the residual interaction operator by the relation (with i, j = n, p):

δE res = 1 2 i,j k1k2k3k4
∂ 2 E ∂ρ i;k1k2 ∂ρ j;k3k4 δρ i;k1k2 δρ j;k3k4 = 1 2 i,j Tr (1) Tr (2) Vres (i1,j2) δ ρ(i1) δ ρ(j2) . (B19) Equations (B9), (B10), (B14) and (B17) then become respectively: A i;q A * j;q 2 v res ij (q) = q i,j

∂µ i ∂ρ j = v b,res ij + δ ij N 0,i (
A i;q A * j;q 2 v b,res ij + 2C ∇ ij q 2 + 4πe i e j q 2 , (B22) and since v res ij (q) = v b,res ij + 2C ∇ ij q 2 + 4πe i e j q 2 = C f ij (q) -

δ ij N 0,i , (B23)
we obtain the following relation between the residual interaction and free-energy density curvature matrixes:

C f = (N 0 ) -1 + v res ; N 0 = N 0,
= k1k ′ 1 k2k ′ 2 k ′ 1 k ′ 2 | Vres |k 1 k 2 k 1 |δ ρ|k ′ 1 k 2 |δ ρ|k ′ 2 (C1)
In the discrete basis, the one-body density matrix elements are:

k i |ρ|k j = n(k i )δ ij , (C2)
whose dependence on the particle number N is contained in the occupation number n(k i ) = n(k i , ρ = N/Ω):

k i |δ ρ|k j = k i | ∂ ρ ∂N |k j δN = 1 Ω ∂n(k i , ρ) ∂ρ δ ij δN . ( C3
)
We have then for the residual interaction term:

2 δE res δN 2 = 1 Ω 2 k1k2 k 1 k 2 | Vres |k 1 k 2 ∂n(k 1 , ρ) ∂ρ ∂n(k 2 , ρ) ∂ρ . (C4)
Let us now take the continuous limit. Denoting k 1 k 2 | Vres |k 1 k 2 = V res k1,k2;k1,k2 , the residual interaction term can be expressed as: For non-relativistic nucleons, with a Skyrme interaction, we have:

X(k, ρ) = X(k 2 , ρ) = µ(ρ) -U(ρ) + k 2 2m * (ρ) = ν(ρ) - k 2 2m * (ρ) . ( C8

FIG. 1 :

 1 FIG.1:(color online) Binding energy in pure neutron matter (top), in symmetric nuclear matter (middle) and symmetry energy (bottom) versus the density for the three Skyrme interactions considered in this paper (SLy230a, SGII and RATP). The saturation point in symmetric nuclear matter is also indicated (dots).

FIG. 2 :

 2 FIG.2:(color online) Contours of the thermodynamic spinodal region for two values of the temperature (T = 0 and T = 10 MeV) and for the Skyrme-like interactions: SLy230a, SGII and RATP.

  ρs

FIG. 4 :

 4 FIG.4: (color online) Instability regions associated with fixed values of the wave-number q, at zero temperature, according to RPA (thick curves) and semi-classical (thin curves) criterions. Upper part: representation in the density plane. Lower part: representation in the plane of Fermi momenta.

FIG. 5 :FIG. 6 :

 56 FIG. 5: (color online) Envelopes of instability against finite-size density fluctuations. Left: comparison between semi-classical and RPA criterions, for temperatures T = 0; 5; 10 MeV. Right: comparison between the RPA instability envelopes obtained with different Skyrme forces, for T = 0; 10 MeV.

FIG. 8 :

 8 FIG. 8: (color online) Time and space characteristics of the most unstable mode, represented in level-lines inside the instability envelope (thick curve). RPA and hydrodynamics results are shown. Top: time constant τ0 giving the fastest amplification of the density fluctuation. Bottom: half-wavelength λ0/2 = π/q0 of the dominant fluctuation mode. Levels are separated by an interval of 0.5 fm.

FIG. 9 :

 9 FIG.9:(color online) Phase-separation directions according to the most unstable mode, reprented by levels and tangential lines (see text for detailed explanation) inside the instability envelope (thick curve). RPA and hydrodynamics results are presented. Levels are separated by intervals of 0.1.

FIG. 10 :

 10 FIG. 10: (color online) Time, size and isospin fractionation characterizing the most unstable mode: dependence on the Skyrme parametrization. RPA (thick line) and hydrodynamics (thin line) results are shown. Left: constant density. Right: constant proton fraction.

3 FIG. 11 :

 311 FIG.11:(color online) Time, size and isospin fractionation characterizing the most unstable mode: temperature evolution. RPA (thick line), Vlasov (dashed line) and hydrodynamics (thin line) results are shown. Constant proton fraction: Z/A = 0.5 for time and size characteristics (four upper lines), Z/A = 0.3 for isospin fractionation.

  fluctuation δρ(r), the average residual energy variation per unit volume is (for the bulk part): Fourier development of δρ(r):

  to explicit n(k, ρ). For a fermion gas at thermodynamic equilibrium, at temperature T = 1/β, with the chemical potential µ(ρ) and the individual energies ǫ(k, ρ), we have:n(k, ρ) = g 1 + e -βX(k,ρ) -1 ; X(k, ρ) = µ(ρ) -ǫ(k, ρ) ,(C6)with ǫ(k, ρ) the individual particle level of momentum k and g the level degeneracy. The residual interaction term now reads

  whatever the form of n(k, ρ), since ρ = dk (2π) 3 n(k, ρ).

TABLE II :

 II Contributions to the quantity Cs coming from the symmetry energy and its derivatives at point ρs (upper border of the spinodal in symmetric matter). The value of ρs is also indicated. The Skyrme-like interactions SLy230a, SGII and RATP have been considered.

TABLE IV :

 IV Limiting density and temperature for symmetric nuclear matter (subscript "s"). The thermodynamic critical point is compared with the values obtained from the semi-classical (SC) and RPA approaches.

				ρc [fm -3 ] [fm -3 ] [fm -3 ] [MeV] [MeV] [MeV] ρ s,SC lim ρ s,RPA lim Tc T s,SC lim T s,RPA lim
		Sly230a 0.0538 0.0461 0.0458 14.55 12.09 11.63
			SGII	0.0518 0.0451 0.0461 14.46 12.33 11.78
		RATP 0.0550 0.0473 0.0474 14.72 12.25 11.81
	ρc	ρ SC lim	ρ RPA lim

TABLE V :

 V Thermodynamic critical point (total density ρc, proton fraction Z/A = 0.5 and temperature Tc) compared with the limiting values obtained from semi-classical (SC) and RPA for finite-size instabilities.

		0.0538 0.0467 0.0474	0.426	0.416	14.55 12.20 11.77	76.7	70.1
	SGII	0.0518 0.0456 0.0460	0.447	0.431	14.46 12.39 11.87	85.0	75.1
	RATP 0.0550 0.0478 0.0478	0.436	0.427	14.72 12.34 11.92	76.6	70.2

  APPENDIX C: LINK BETWEEN RESIDUAL INTERACTION Vres AND ENERGY VARIATION δE resWe first work in a discrete basis |k i , where the residual energy variation is expressed by:2δE res = Tr (1) Tr (2)Vres (1,2) δ ρ(1) δ ρ(2)

	n 0 0 N 0,p	.	(B24)

)

  At zero temperature, the density derivative of the occupation number is a Dirac function:d cos θ 1 d cos θ 2 V res F (cos θ 1 , cos θ 2 ) . (C11)Neglecting the angle dependence of the residual interaction, we have V res F (cos θ 1 , cos θ 2 ) = V res F . Using the relationIn order to generalize this result at finite temperature, let us now assume that the residual interaction elements V res k1,k2;k1,k2 are equal to a constant V res . Starting again from Eq. (C5) we simply have:

	Since X = 0 for k 2 = k 2 F , the integral becomes
		2	δE res δN 2 =	gk F m * (2π) 2	∂X(k 2 F , ρ) ∂ρ
						∂X(ρ, k 2 F ) ∂ρ	=	2π 2 m * gk F	(C12)
	we finally obtain that
							2	δE res δN 2 = V res F .	(C13)
	2	δE res δN 2 = V res 1 (2π) 6	dk 1 dk 2	∂n(k 1 , ρ) ∂ρ	∂n(k 2 , ρ) ∂ρ	= V res ∂ ∂ρ	dk (2π) 3 n(k, ρ)
				n(k 2 , ρ) = gΘ(X(k 2 , ρ))	(C9)
				∂n(k 2 , ρ) ∂ρ	= g	∂X(k 2 , ρ) ∂ρ	δ(X(k 2 , ρ)) .	(C10)

2

ACKNOLEDGMENTS C.D. and J.M. want to thank N. Van Giai for interesting discussions concerning the RPA framework in the unstable phase diagram region. C.D. thanks C. Providência for fruitful discussions on the Vlasov formalism. This work has received the support of the project "modélisation des étoiles à neutrons" supported by CNRS-IN2P3.

APPENDIX A: LINEAR RESPONSE TO AN EXTERNAL PROBE

We detail the calculation of the response function for the single-fluid case. Let us consider the response to an infinitesimal external field of the form: Vext = Ee iq•r-i(ω+iη)t (t < 0) ; Vext = 0 (t ≥ 0) , (A1) carrying a momentum q and an energy ω (the positive infinitesimal part η ensuring the adiabatic switching of the external field). At late enough times and small enough external field E, the expectation value of the density fluctuation has the same space and time dependence as the external field:

where α includes the time dependence for a more compact notation. The linear response Π(ω, q) is defined by the ratio of the density change to the external field strength, for a small external field E:

A density fluctuation δ ρ(t) defined such as ρ(t) = ρ0 + δ ρ(t) induces a variation δ ĥ(t) of the mean-field operator, such that ĥ(t) = ĥ0 + δ ĥ(t). The linearized TDHF equation

reads in momentum space:

The elements of Vext are:

and the elements of δ ρ have the form:

The density factors α and κ(k, q) are related by:

According to Eqs. (A5-A7), the matrix elements of δ ĥ are of the form:

where ν(k, q) is connected to the residual interaction Vres . Denoting E[ρ] the total energy of the system, this residual interaction operator is defined by [START_REF] Ring | The Nuclear Many-Body Problem[END_REF]:

where we have limited the development to the second order in δρ ij . From these relations, we also have the expression of the mean-field variation in terms of the residual interaction (at first order in δρ ij ):

so that from Eq. (A9) we have:

Assuming that the residual interaction element V res k+q,k2 ; k,k2+q does not depend on (k, k 2 ) -Landau monopolar approximation -it can be written:

and we obtain:

with the notation v res = ΩV res . The matrix elements of δ ĥ are then:

Inserting the expressions obtained for the matrix elements of δ ρ, δ ĥ and Vext , the linearized TDHF equation (A5) now reads:

where ω kq = ω k+q -ω k . The second line is obtained by summing over k and k ′ . Taking the continuous limit of Eq. (A18) and multiplying by the time factor e -i(ω-iη)t , we obtain:

where Π 0 is the uncorrelated response function (Lindhard function). A rearrangement of Eq. (A19) gives the following expression for the response function:

This result could also be deduced from the Bethe-Salpeter equation [START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF], keeping the transfered-momentum dependence of the residual approximation and treating the exchange terms at the Landau approximation: this is the LAFET approximation detailed in Ref. [START_REF] Margueron | [END_REF].

APPENDIX B: RESIDUAL INTERACTION AND FREE-ENERGY CURVATURE

In order to establish the link between the residual interaction and the free-energy curvature, we start with the thermodynamic case considering a homogeneous system with large particle number N and volume Ω. For simplicity we first expose the calculation details in the single-fluid case. Fixing Ω, at equilibrium the density matrix ρ is uniquely defined by N and a fluctuation of the particle number δN leads to a fluctuation δ ρ which can be expressed by the Taylor expansion:

and the free-energy variation is:

where S is the entropy and µ the chemical potential. The superscripts "b" (for "bulk") specify that we are in the thermodynamic description. δF b can be separated into an non-interacting term δF b,NI and a residual interaction term δF b,res . The non-interacting part gives the free-energy variation for non-interacting particles of mass m * in a constant mean field. It contains the entropy variation, since the entropy S is determined by the density independently of the particle interaction. The residual interaction term then reduces to δF b,res = δE b,res , which is the term due to the variation of the mean field δ ĥ, so that:

Referring to the expression (A10), δE b,res is expressed in terms of the residual-interaction operator as:

∂ ρ( 1)

Taking constant matrix elements for Vb,res , such that V b,res jl;ik = V b,res , we obtain (cf Appendix C):

so that

Identifying the second-order term of Eqs. (B2) and (B6) gives:

The non-interacting part of the chemical-potential derivative is related to the density of states at the Fermi level, denoted N 0 . Taking the derivative with respect to the density ρ = N/Ω, we have:

where k F is the Fermi momentum. Using the notation ΩV b,res = v b,res , we can re-express Eq. (B7) as:

which establishes the link between the thermodynamic free-energy curvature and the bulk part of the residual interaction.