Camille Ducoin

Jérôme Margueron

Constança Providência

Isaac Vidaña

Core-crust transition in neutron stars: predictivity of density developments

The possibility to draw links between the isospin properties of nuclei and the structure of compact stars is a stimulating perspective. In order to pursue this objective on a sound basis, the correlations from which such links can be deduced have to be carefully checked against model dependence. Using a variety of nuclear effective models and a microscopic approach, we study the relation between the predictions of a given model and those of a Taylor density development of the corresponding equation of state: this establishes to what extent a limited set of phenomenological constraints can determine the core-crust transition properties. From a correlation analysis, we show that (a) the transition density ρt is mainly correlated with the symmetry energy slope L, (b) the proton fraction Yp,t with the symmetry energy and symmetry energy slope (J, L) defined at saturation density, or, even better, with the same quantities defined at ρ = 0.1 fm -3 , and (c) the transition pressure Pt with the symmetry energy slope and curvature (L, Ksym) defined at ρ = 0.1 fm -3 .

I. INTRODUCTION

The study of nuclear systems under extreme conditions is an expanding field of research involving astrophysical processes, laboratory experiments with rare-isotope beams, and the development of more realistic nuclear models. The density dependence of the symmetry energy S(ρ) is one of the central issues in this field [1][2][3]. Although the symmetry energy at saturation density ρ 0 is considered to be well known (J = S(ρ 0 ) = 33 ± 3 MeV), the different nuclear models present a wide range of predictions for the symmetry energy slope L = 3ρ 0 [∂S/∂ρ](ρ 0 ). Several experimental observables have been proposed to obtain a measure of L, for instance neutron-skin thickness [4], isovector dipolar resonances [5,6], isoscaling in multifragmentation [7] and isospin diffusion [8]. The perspective to obtain more stringent constraints on the value of L has motivated several studies concerning the impact of this quantity on compact star structure, especially on the characterization of the core-crust transition in neutron stars [9][10][11]. In a previous work [12], we have addressed the model dependence of the link between L and the core-crust transition properties: density ρ t , proton fraction Y p,t and pressure P t . It appeared that ρ t and Y p,t are unambiguously correlated with L, while the link between L and P t is very sensitive to the model. The transition pressure is the dominant input for the prediction of one of the most important crust properties: namely, the moment of inertia of the crust that may affect pulsar glitches [13,14]. It is then an important challenge to find a better relation between P t and laboratory data.

In this paper, we investigate the role of symmetry-energy properties other than L that could also, hopefully, be related to laboratory data. To perform this study, we use as previously several Skyrme and relativistic effective models, and a microscopic Brueckner-Hartree-Fock (BHF) approach. We also employ a schematic nuclear equation of state (EOS), based on a density development around a reference density, that we call the generalized liquid-drop model (GLDM): our objective is to establish to what extent the core-crust transition obtained with a given nuclear model can be reproduced with a limited set of GLDM coefficients. In other words, we investigate how the specificities of each model can affect the link between phenomenological constraints (i.e., the set of GLDM coefficients, which could be determined from laboratory data) and the core-crust transition.

In this paper, the core-crust transition properties are first calculated in the thermodynamic framework, with the transition defined as the crossing between the β equilibrium line and the thermodynamic spinodal border. The thermodynamic spinodal corresponds to the bulk instability of homogeneous nuclear matter. In fact, in a compact star at the core-crust transition the ground state of stellar matter changes from a clusterized configuration (lattice of nuclei, or more exotic structures named pasta phases [2,15,16]) to homogeneous matter. It has been shown that this transition can be very well approximated by the crossing between β equilibrium and the dynamic spinodal border [17], which corresponds to the finite-size instability region of homogeneous matter: Coulomb and surface contributions make the dynamic spinodal smaller than the thermodynamic one [18][19][20]. Therefore, we also address the dynamic core-crust transition properties, which are strongly correlated with the thermodynamic ones.

This paper is organized as follows. The generalized liquid-drop model is presented in Sec. II. We discuss the correlations existing between various GLDM coefficients and evaluate the accuracy of the core-crust transition predictions obtained with different GLDM expansions. In Sec. III, we present in more details the analysis of the link between the symmetry-energy slope at saturation L and the core-crust transition, which was first presented in Ref. [12]. In Sec. IV, we perform a correlation analysis involving coefficients others than L. This allows us to establish relations between selected GLDM coefficients and realistic transition properties, within a reduced range of model dispersion.

Conclusions are presented in Sec. V.

II. GENERALIZED LIQUID-DROP MODEL

Nuclear density functionals allow us to calculate the nuclear matter equation of state for any total density ρ = ρ n +ρ p and asymmetry y = (ρ nρ p )/(ρ n + ρ p ). They may also present a complex dependence on the density, including, for instance, kinetic densities, spin densities, and density gradients. Such functionals can be derived from different nuclear models. In this paper, we consider Skyrme and relativistic effective models, as well as an equation of state based on BHF calculations. These models, presented in the Appendix, are referred to as 'complete functionals', in contrast with the GLDM that is presented in this section.

The GLDM corresponds to a series expansion of the EOS around a given reference density. It is determined by three choices : (i) the nuclear model giving the full equation of state E(ρ, y), (ii) the reference density ρ ref , and (iii) the order N of the development. The order N and the reference density (ρ ref ) can be chosen at will. In this paper, the GLDM coefficients are derived from the complete EOS E(ρ, y) given by various nuclear models. In a different context, the GLDM coefficients could be nuclear properties determined experimentally at ρ ref , from which we would extrapolate the nuclear EOS at lower or higher densities. We introduce the GLDM in order to investigate how well such an extrapolation allows us to reproduce the core-crust transition predicted by the complete functional; in other words, to what extent the role of higher-order terms can be neglected. Comparing the GLDM predictions with the results from the corresponding complete functional shows how much a more detailed density dependence of the EOS may affect these predictions.

A. GLDM equation of state

For a given reference density ρ ref and order of development N , the GLDM energy per particle reads as:

E GLDM (ρ, y) = N n=0 c IS,n + c IV,n y 2 x n n! + (E kin -E para kin ) , with x = ρ -ρ ref 3ρ ref . (1)
The first term on the right-hand side of Eq.( 1) contains both the kinetic and the potential contributions to the energy in the parabolic approximation with respect to the asymmetry y. The second term gives the contribution of the kinetic term beyond the parabolic approximation, as it will be explained below. We have introduced in this expression the GLDM coefficients c IS,n and c IV,n , respectively associated with the derivatives of the energy E(ρ, y = 0) and of the symmetry energy S(ρ): the index 'IS' ('IV') stands for isoscalar (isovector). They are expressed as:

c IS,n (ρ ref ) = (3ρ ref ) n ∂ n E ∂ρ n (ρ ref , 0) (2) c IV,n (ρ ref ) = (3ρ ref ) n ∂ n S ∂ρ n (ρ ref ) with S(ρ) = 1 2 ∂ 2 E ∂y 2 (ρ, 0) (3)
In the case ρ ref = ρ 0 , the lower-order coefficients are usual nuclear matter properties: c IS,0 = E 0 (saturation energy), c IS,2 = K 0 (incompressibility), c IS,3 = Q 0 , c IV,0 = J (symmetry energy), c IV,1 = L (symmetry-energy slope), c IV,2 = K sym (symmetry incompressibility), and c IV,3 = Q sym . The parabolic approximation, which restricts Eq. ( 1) to the first term, is known to be quite accurate to describe the EOS even at high isospin asymmetry. However, it fails to reproduce the spinodal contour in the neutron-rich region. The reason is that the energy-density curvature in the proton-density direction must diverge at small proton density because of the kinetic term (see the appendix of Ref. [21]); as a result, the spinodal contour can not reach pure neutron matter. Instead, in the parabolic approximation, the curvature in the proton-density direction is constant, and leads to the unphysical prediction of unstable neutron matter. To avoid this discrepancy, we have introduced in Eq. (1) a model-independent correction based on the non-relativistic, free Fermi gas kinetic term:

E kin = 1 ρ 2 2m (τ n + τ p ) , (4)
where τ q = (3π 2 ρ q ) 5/3 /(5π 2 ) is the kinetic density of the nucleon species q = n, p and m is the nucleon mass. As a function of density and asymmetry, we have:

E kin = (3π 2 /2) 5/3 10mπ 2 ρ 2/3 (1 + y) 5/3 + (1 -y) 5/3 (5) E para kin = (3π 2 /2) 5/3 10mπ 2 ρ 2/3 2 + 10 9 y 2 , (6)
where E para kin is the parabolic part of E kin . In the GLDM defined by Eq. ( 1), the extra-parabolic behavior of the functional is sketched by the extra-parabolic behavior of E kin , which brings the model-independent correction E kin -E para kin (second term on the r.h.s. of Eq. ( 1)).

In the following, we will use the notation D N (ρ ref ) to identify a development of order N around the density ρ ref .

Developments of this kind are usually considered up to N =3, around the saturation density ρ 0 ≃ 0.16 fm -3 (see, e.g., Refs. [22,23]). In this work, we will also consider an extreme situation labeled D ∞ . Performing an infinite development gives the exact value of E(ρ, 0) and S(ρ) for any density ρ, whatever the choice of ρ ref . In practice, the D ∞ equation of state is simply built using the exact expressions of E(ρ, 0) and S(ρ) to obtain E(ρ, y):

E D∞ (ρ, y) = E(ρ, 0) + S(ρ)y 2 + (E kin -E para kin ) . (7)
Thus, the difference between a complete functional and its associated D ∞ model is just the extra-parabolic content of the nuclear interaction, which has not been taken into account by the correction E kin -E para kin defined above. Let us note that a parabolic approximation is often assumed in BHF calculations, to interpolate between the symmetricmatter and neutron-matter EOS. We will call this approach BHF para . However, by default, our BHF results include the correction E kin -E para kin ; in this case, the complete functional is exactly equivalent to the corresponding D ∞ model.

B. Correlations between GLDM coefficients

The isoscalar and isovector GLDM coefficients obtained with all the nuclear models used in this paper are reported in Table I. These coefficients are intercorrelated, due to common constraints that the various parametrized forces have to verify. In particular, the fitting procedure on finite-nuclei properties, for which the average density is lower than ρ 0 , is likely to provide effective constraints in the density region ρ 0 ≃0.1-0.12 fm -3 (this point is raised, e.g., in Refs. [24,25]); as a result, the coefficients defined at ρ 0 for various effective models should show a tendency to compensate each other in order to focus the various functional predictions at lower density. To understand how a given coefficient (such as L) can affect the EOS properties away from the reference-density region, it is important to evaluate its possible connections with other coefficients. Since the equation of state of asymmetric nuclear matter is our main concern here, we will concentrate on the correlations between isovector properties.

We represent in Fig. 1 different relations between GLDM isovector coefficients. On the left panel, we consider the coefficients taken at saturation density, c IV,n (ρ 0 ); the different plots shows the correlation between L and J, K sym , and Q sym . As we have emphasized in Ref. [12], there is a strong L-K sym correlation. The two eccentric points correspond to the relativistic models DDHδI-25 and DDHδII-30: these models include the δ meson, which is generally associated with an atypical density dependence of the symmetry energy (see, e.g., Ref. [26]). As for the symmetry energy J, it tends to increase with L, although the dispersion between different models is important. The behaviour of the third-derivative coefficient Q sym is less universal: a clear decreasing L-Q sym correlation is obtained among the Skyrme models, but the relativistic models of lower L and the BHF point are completely out of this trend. Linear fits are represented on these plots; in the case of K sym and Q sym , they have been performed using only the models in agreement with the main trend.

The right panel of Fig. 1 shows the relation existing between coefficients defined at two different densities, namely between c IV,n (ρ 0 ) and c IV,n (ρ=0.1 fm -3 ), for n=0, 1, 2. The correlations obtained in the cases of c IV,1 (symmetry energy slope) and c IV,2 (symmetry-energy curvature) reflect the shape similarities between many of the density functionals. The absence of correlation in the case of c IV,0 (symmetry energy) is due to the fact that c IV,0 (0.1) is strongly constrained in nuclear models. This discussion is illustrated in Figs. 2 (Skyrme models) and 3 (relativistic models).

The left panel of each figure displays the density dependence of the symmetry energy, along with its slope and curvature: S(ρ), L(ρ) = (3ρ 0 )[dS/dρ], and K sym (ρ) = (3ρ 0 ) 2 [d 2 S/dρ 2 ]. Note that L(ρ) and K sym (ρ) are defined here with a constant factor involving ρ 0 : they are not equivalent to the coefficients c IV,1 (ρ) and c IV,2 (ρ), except at ρ = ρ 0 . We use here a constant factor in order to represent quantities proportional to the derivatives of S(ρ). The curves are shown for all models under consideration, using the complete functionals. The comparison between Figs. 2 and3 shows that the relativistic models present more variability in the shape of the density functional than Skyrme models,

Model ρ0 E0 K0 Q0 J L Ksym Qsym E01 LIS,01 K01 J01 L01 Ksym,01 [fm -3 ] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] Microscopicc BHF-
J = cIV,0, L = cIV,1, Ksym = cIV,2, Qsym = cIV,3. Isoscalar coefficients at ρ ref = 0.1 fm -3 : E01 = cIS,0, LIS,01 = cIS,1, K01 = cIS,2. Isovector coefficients at ρ ref = 0.1 fm -3 : J01 = cIV,0, L01 = cIV,1, Ksym,01 = cIV,2.
as we have shown in Ref. [26]. This comes from the different ways of describing the nuclear interaction by meson exchange in the effective relativistic framework; in particular, the inclusion of the δ meson has a strong effect on the density dependence of L. We can also notice that relativistic models such as GM1 and GM3, which were built for astrophysical purpose, have not been explicitly fitted on laboratory data. Considering globally the Skyrme and relativistic functionals, we can notice two common convergence regions. One concerns the symmetry-energy: all curves tend to cross around the density ρ ≃ 0.11 fm -3 , taking values S(0.11 fm -3 ) = 24 ± 4 MeV. This behavior was expected, following the previous remark that finite nuclei provide fitting constraints at density slightly below saturation; similar observations have been made e.g. in Ref. [27]. The second convergence region concerns the symmetry-energy slope L(ρ): the different curves show a marked tendency to cross at about ρ 0 /3. In contrast, no convergence effect appears for the second derivative, K sym (ρ). On the right panel of Figs. 2 and3, we represent the density dependence of the variance between the values taken by different models. Denoting X(ρ) = {S(ρ), L(ρ), K sym (ρ)}, the variance is given by:

σ X (ρ) = i [X i (ρ) -X(ρ)] 2 , X(ρ) = 1 N N i=1 X i (ρ) , (8)
where the index i runs over the N models considered. The convergent trends are reflected by the density dependence of the respective variances. In the case of Skyrme functionals, a clear minimum of σ S occurs at ρ ≃ 0.13 fm -3 . With relativistic models, we observe a plateau around an inflexion point at ρ ≃ 0.1 fm -3 . Note that σ S is constrained to cancel at zero density; in such a condition, the inflexion point can be considered as a criterion of convergent trend. The convergence effect is even more clear in the case of the symmetry-energy slope L(ρ); both Skyrme and relativistic models present a sharp minimum of σ L at ρ ≃ 0.06 fm -3 . These results indicate that the correlations existing between J, L, K sym , and Q sym may be associated with effective constraints on the values of S(ρ ≃ 0.12 fm -3 ) and L(ρ ≃ 0.06 fm -3 ). To verify this, we show on the same figures the variances obtained when the curves X(ρ) are calculated using the GLDM at different orders, with reference density ρ 0 . A minimum of σ S in the density range 0.11-0.13 fm -3 is obtained with versions D 1 (ρ 0 ) and D 2 (ρ 0 ) of the GLDM. The version D 3 (ρ 0 ) gives a clear minimum of σ L at ρ ≃ 0.05 fm -3 . We can conclude that the correlations between the GLDM coefficients are related with convergence effects for S(ρ) and L(ρ) in the subsaturation density region.

Let us remark that the convergence of the L values at ρ ≃ 0.06 fm -3 is easily interpreted as a geometrical consequence of a constraint on S(ρ ≃ 0.12 fm -3 ). Indeed, let us imagine a constraint fixing the values (ρ c , E c ) such that all models have to verify S(ρ c ) = E c . Since we also have the condition S(0) = 0, it turns out that for all models, L(ρ) takes the value L c = E c × 3ρ 0 /ρ c at least once on the interval [0, ρ c ]. For instance, if ρ c = 0.12 fm -3 , E c = 25 MeV and ρ 0 = 0.16 fm -3 , we have L c = 100 MeV.

C. Comparison between GLDM and complete functionals

We now explore the accuracy of the GLDM approximation on reproducing the core-crust transition obtained with the complete model on which it has been built. This procedure allows to estimate whether the core-crust transition can be efficiently characterized by a limited set of EOS properties determined at a fixed density.

The discussion of the core-crust transition is performed in the thermodynamic framework, which means that, for simplicity, the transition is defined as the crossing between β equilibrium and the thermodynamic spinodal border. The corresponding transition density, proton fraction, and pressure are denoted, respectively, ρ tt , Y p,tt , and P tt , where the index tt stands for thermodynamic transition. In Sec. IV, we will also discuss the dynamic transition defined using the finite-size spinodal border, in which case the index td will be used. The transition properties in both thermodynamic and dynamic cases are given in Table II.

We compare in Fig. 4 the transition properties (ρ tt , Y p,tt , and P tt ) obtained with the complete functionals and different versions of their associated GLDMs. On the left panel, we consider second-and third-order developments around saturation density, namely, D 2 (ρ 0 ) and D 3 (ρ 0 ); on the right panel, we consider the second-order development around a lower reference density ρ ref = 0.01 fm -3 , namely, D 2 (0.1). On each part, we add the results from the fully developped GLDM, D ∞ , for which the only difference with the complete functional is due to extra-parabolic terms in the nuclear interaction. For the BHF model, the functional is equivalent by construction to its associated D ∞ model.

We can see that D 2 (ρ 0 ) leads generally to an important underestimation of ρ tt , Y p,tt , and P tt . This means that, in a development around saturation density, terms beyond order 2 (i.e., beyond K sym ) have a large impact on the properties of the core-crust transition; the correlations observed between L and these properties can occur only if higher-order corrections are either correlated with L, or similar, for most of the functionals. The underestimation is strongly attenuated, but still present, with D 3 (ρ 0 ), which involves the knowledge of isovector coefficients until Q sym . The situation is much improved if we use a development around ρ ref = 0.1 fm -3 . The accuracy is globally better with D 2 (0.1) than with D 3 (ρ 0 ); furthermore, it is nearly as good with D 2 (0.1) as with D ∞ , which means that, at ρ ref = 0.1 fm -3 , it is enough to perform the development up to second order. It is then important to relate experimental observables with the symmetry-energy density dependence directly in this low-density region. This is also appropriate, since the nucleus properties used to constrain the symmetry energy are often associated with subsaturation densities (neutron-skin thickness, resonances, multifragmentation, isospin diffusion,...).

Let us notice that effective nuclear models have an important role to play for the fine-tuning of the curvature properties of the EOS, crucial for the prediction of the core-crust transition. Indeed, although microscopic methods such as BHF are necessary to obtain reliable EOS values far away from the phenomenological constraints (in particular, for neutron matter and at high density), they can not be used as a reference for the EOS curvature. The numerical BHF results for E(ρ) do not allow the direct determination of second derivatives; the EOS has to be fitted (see, for instance, Ref. [28]), and the resulting curvature properties are sensitive to the fitting conditions. We show in Fig. 5 how this affects the predictions for the core-crust transition. This figure displays the values of ρ tt , Y p,tt , and P tt obtained by using the BHF calculations in different ways. In the first version, BHF-1, the fit is performed on the density interval ρ = [0.1; 0.35] fm -3 : this is the version that is used to establish the GLDM coefficients at ρ 0 . In order to focus on the subsaturation region, we have considered a second version, BHF-2, for which the fit is performed on the density interval ρ = [0.05; 0.18] fm -3 : this is the version we will use afterwards to define the BHF core-crust transition. For now, let us compare the BHF-1 and BHF-2 predictions. The GLDM coefficients at ρ ref = 0.1 fm -3 have been determined for both versions, and there are significant differences in the transition properties predicted by the D 2 (0.1) expansion in each case. A similar contrast appears between the results obtained with the full BHF-1 and BHF-2 EOS. In addition to these two versions of BHF calculations, the figure also shows the complete and GLDM results for the Skyrme force LNS [47], the fitting procedure of which involves the BHF equation of state. From the span of results we obtain, it is clear that a microscopic calculation does not lead to a unique prediction for the corecrust transition properties; phenomenological constraints from finite nuclei will be essential to improve our knowledge of the low-density EOS.

III. CORRELATION BETWEEN L AND THE CORE-CRUST TRANSITION POINT

We consider in this section the specific role of the symmetry energy slope L in the determination of the core-crust transition, defined here as the crossing between the line of β equilibrium and the thermodynamic spinodal contour. We summarize the study performed in Ref. [12], and present more details that support this previous analysis.

A. Position of the transition point

The most direct impact of the symmetry energy on neutron-star structure concerns the proton fraction Y p (ρ) in stellar matter, which is fixed by β equilibrium. For a given density, a lower symmetry energy corresponds to a lower proton fraction. Thus, as far as a high value of L can be correlated with a low value of the symmetry energy at sub-saturation, we expect higher-L models to provide lower values of Y p,tt . This point is illustrated by Fig. 6. At fixed density ρ = 0.08 fm -3 , for increasing L, the β-equilibrium proton fraction decreases. This trend is confirmed and even accentuated if, instead of fixing the density, we consider the proton fraction at the transition point. The dispersion observed in both cases is essentially due to different values of the symmetry energy at saturation density, as will be discussed in the following.

Let us now consider the impact of L on the transition density ρ tt , illustrated by Fig. 7. The correlation between L and ρ tt is a well-known feature [9]; however, its explanation is less intuitive than in the case of the L-Y p,tt correlation. Furthermore, it can not be explained just as a consequence of the behavior of Y p,tt , as we see on the left panel of the figure: even for a fixed proton fraction Y p = 0.025, the spinodal border shows a clear decreasing correlation with L. This feature can be understood as a consequence of the strong link existing between L and the energy-density curvature of neutron matter taken at symmetric spinodal density ρ s :

C NM,s = 2 3ρ 0 L + 1 3ρ 0 n≥2 c IV,n x n-2 s (n -2)! n + 1 n -1 x s + 1 3 + ∂ 2 [ρ(E kin -E para kin )] ∂ρ 2 (9)
with x s = (ρ sρ 0 )/(3ρ 0 ). Since all models yield a symmetric spinodal density close to 0.1 fm -3 , so that x s ≃ -1/9, the term n = 2 is nearly canceled: this reinforces the dominance of L in the determination of C NM,s . Figure 8 illustrates how L affects independently the proton fraction at β equilibrium and the spinodal contour in the neutron-rich region. Due to the typical geometry of these respective lines, the two effects reinforce each-other, leading to a robust correlation between L, ρ tt and Y p,tt .

B. Pressure at the transition point

The link between L and the core-crust transition pressure P tt is more problematic. In order to make this link explicit, let us write the pressure in the GLDM framework: From this expression, we expect that for a given density the pressure of neutron-rich matter should increase with L, which is the leading coefficient. This trend appears on the left panel of Fig. 9, representing the relation between L and the pressure of pure neutron matter, P NM , at ρ=0.08 fm -3 . Thus, a positive correlation between L and P tt should be obtained if we could neglect the density shift due to the L-ρ tt correlation, as well as the effect of higher-order coefficients. However, as it can be seen on the right panel of the figure, the results for P tt (L) present an important dispersion and we cannot extract a clear correlation, although a decreasing trend can be observed among Skyrme models. Four eccentric points close to L=60 MeV weaken this correlation between Skyrme models: they correspond to atypical relations between L and K sym , which also affect the plot P NM (0.08 fm -3 ). The lack of correlation between L and P tt when independent models are considered results from a delicate balance between opposite effects, as we have discussed in Ref. [12]. This is shown by separating the different contributions we can estimate from the GLDM formula. We distinguish two kinds of contributions to the variation dP tt /dL: (i) variations occurring at a fixed density (ρ, y), resulting only from the modifications of the coefficients in Eq. ( 10), which defines P (ρ, y); and (ii) variations due to a shift (δρ, δy), for a fixed expression of P (ρ, y), i.e. frozen values of the coefficients in Eq. (10). The contributions of the first kind come from the explicit L dependence of Eq. (10), and from correlations between L and higher-order coefficients c IV,n . In practice, we will consider the following terms:

P (ρ, y) = ρ 2 3ρ 0   Ly 2 + n≥2 c IS,n + c IV,n y 2 x n-1 (n -1)!   + ρ 2 ∂(E kin -E para kin ) ∂ρ . (10)
L (MeV) 50
δP tt δL e1 = ∂P ∂L (ρ tt , y tt ) ( 11
)
δP tt δL e2 = ∂P ∂K sym (ρ tt , y tt ) δK sym δL ( 12
)
δP tt δL e3 = ∂P ∂Q s (ρ tt , y tt ) δQ s δL (13)
where the index e means that the modification concerns the expression of the pressure, and the number gives the order of the modified coefficient. The contributions of the second kind, resulting from the density position of the transition point, are characterized by the index p. We distinguish the respective effects of total density and asymmetry:

δP tt δL p1 = ∂P ∂ρ (ρ tt , y tt ) δρ tt δL ( 14
)
δP tt δL p2 = ∂P ∂y (ρ tt , y tt ) δy tt δL (15)
The variations of the quantities depending on L are fixed empirically, using as a reference the correlations that are observed between different models. From linear fits, we extract δρ tt /δL=-3.84 × 10 -4 MeV -1 fm -3 , δy tt /δL=6.08 × ) and to the variation of the GLDM coefficients cIV,n in the expression of P (indices e1, e2, and e3 for n=1, 2, and 3 respectively). Right: sum of the contributions, considering the contribution of δcIV,n/δL up to order n=1, 2 and 3 (respective index: tot1, tot2, tot3). The horizontal lines indicate the region where the pressure variation is compatible with zero within the estimated uncertainty. Full symbols: Skyrme models; empty symbols: relativistic models.

10 -4 MeV -1 , δK sym /δL=3.33, δQ s /δL=-6.63. Note that, in the case of Q s , the correlation with L is observed only within the Skyrme models, which are used to perform the linear fit.

These various contributions are represented on Fig. 10. It appears that the contribution of the asymmetry shift is quite marginal, and the overall δP tt /δL essentially results from the balance between three terms: [δP tt /δL] e1 , which is large and positive, is compensated by the conjugated effect of [δP tt /δL] e2 and [δP tt /δL] p1 . These two negative contributions, e2 and p1, are of the same order of magnitude: this means that the correlation L-K sym has the same importance as the correlation L-ρ tt in explaining why we do not observe an increasing correlation P tt (L).

In addition, we can see that the term (δP tt )/δL) e3 due to the L-Q sym relation brings an additional negative contribution, but of much lower magnitude than the term e2. This result means that, although the third-order term of the GLDM has a strong impact on the absolute value of P tt , as it was observed on Fig. 4, the L-Q sym correlation is not crucial in the determination of δP tt /δL; in other words, the third-order correction does not depend strongly on L. On the other hand, let us notice that the strong dispersion of Q sym values in the case of relativistic models is bound to cause a strong dispersion in P tt (L).

To summarize, if we characterize the L-P tt relation using a GLDM development around saturation density, we can identify three effects that are crucial for the determination of δP/δL: (i) the explicit L dependence of the pressure given by Eq. ( 10); (ii) the L-ρ tt correlation and (iii) the L-K sym correlation. These different contributions compensate each other. For some models (those of higher L), the GLDM predicts a decreasing P tt (L); for others (those of lower L), an increase would be obtained. It is interesting to note that δP/δL cancels in the interval of the most realistic L values, namely 50-80 MeV. By estimating an uncertainty of about 20% on the slopes of the L-ρ tt and L-K sym linear fits, we obtain an error bar of ±0.02 fm -3 on δP/δL, which appears compatible with zero throughout this interval. These results are not a quantitative prediction on the evolution of the transition pressure with L; however, they show that the link between L and P tt cannot be deduced from qualitative arguments, and therefore it is not soundly based. The relation between L and P tt is in fact very sensitive to model-dependence, as it will be further discussed in the following.

C. Predictions of a standard GLDM

In order to study how the different GLDM coefficients can affect the core-crust transition, we will make use of a schematic model corresponding to a D 3 (ρ 0 ) expansion with typical values for the different coefficients. The choice of these values is illustrated on Fig. 11, which gives a graphical representation of the saturation properties of the different nuclear models considered in this paper. The lines indicate the intervals of coefficients attributed to the typical GLDM that we are now constructing. We define a reference D 3 (ρ 0 ) model characterized by the following parameters:

ρ 0 = 0.16 fm -3 ; K 0 = 240 MeV ; Q 0 = -350 MeV ; J = 31 MeV .
The coefficient L varies in the interval [40; 100] MeV, and determines K sym and Q sym according to the relations:

K sym (L) = a K × L + b K ; (a K ; b K ) = (3.33; -281) Q s (L) = a Q × L + b Q ; (a Q ; b Q ) = (-6.63; 765)
The predictions of this standard model for ρ tt (L), Y p,tt (L), and P tt (L) are represented in Fig. 12. As expected, we obtain a clear decrease of ρ tt and Y p,tt with L, while the evolution of P tt is quite flat. In the following, we will observe how these curves evolve when some of the standard EOS properties are modified within a realistic interval. We will modify separately the symmetry energy J, the incompressibility K 0 , and the L-K sym relation. On the left panel of Fig. 12, we show the effect of varying J between 29 and 33 MeV. On the right panel, we show the effect of varying K 0 between 220 and 260 MeV. On the central panel, we modify the linear relation between L and K sym by adopting different values of (a K ; b K ): a softer version (a K ; b K ) soft = (2.8; -235), and a stiffer one (a K ; b K ) stiff = (3.8; -320).

The main features to be noted concerning the results of the standard schematic model and how they are affected by typical variations of the GLDM coefficients are the following. (i) The qualitative behavior of ρ tt and Y p,tt is maintained: although the absolute value can be affected by different aspects of the functional, they always unambiguously decrease with increasing L. (ii) The qualitative behavior of P tt is very sensitive to the values of the GLDM coefficients: the application of a very moderate variation, inside a realistic model uncertainty, leads to opposite predictions: P tt (L) either increases or decreases, and most often it is quite flat. These two conclusions confirm the previous analysis.

L (MeV

IV. PARAMETRIZATION OF THE CORE-CRUST TRANSITION

In this last section, we explore the possibility to reduce the model dispersion in the prediction of the core-crust transition, by taking into account the effect of coefficients others than L. First, we check to what extent the dispersion in the L dependence of ρ tt , Y p,tt , and P tt can be attributed to specific GLDM coefficients, such as the symmetry energy at saturation J, the incompressibility K 0 , or the quantity ∆K sym (L) = K sym -(a K × L + b K ). This last quantity characterizes the eccentricity of the model with respect to the typical relation K sym (L) = (a K × L + b K ). On the basis of these results, we propose to fit the core-crust transition properties by a linear dependence on pairs of GLDM coefficients. This idea is applied to the thermodynamical transition, which has been the framework of our analysis, and to the dynamic transition, which is the best approximation to the realistic core-crust transition.

A. Role of the first GLDM coefficients in the dispersion

As we have seen in the previous section, the correlations between L and the core-crust transition properties ρ tt , Y p,tt , and P tt suffer from a certain amount of dispersion when different kinds of models are considered. This effect is particularly harmful for the prediction of the transition pressure; with the link between L and P tt being very sensitive to the details of the functional, the model dispersion destroys the possibility to deduce the value of P tt from a measurement of L. In order to look for model properties that may be responsible for this dispersion, we first use the following procedure:

L (MeV) 50 100 ) -3 (fm
• We call M the GLDM model property whose effects are investigated: namely, J, K 0 or ∆K sym .

• We represent the quantities X(L) calculated with the different models, with X=Y p,tt (Fig. 13), ρ tt (Fig. 14), and P tt (Fig. 15).

For each M property (J, K 0 or ∆K sym ), we check whether the diagram X(L) can be separated in two regions associated with larger/smaller values of M . For this:

• we use a trial frontier, namely a straight line Λ:

X Λ = a × L + b
• we calculate the distance of each data point i to this frontier:

[∆X] Λ,i = X i -(a × L i + b).
• we vary the frontier Λ in order to obtain the best correlation for the diagram [∆X] Λ (M ).

In this way, a different line Λ is defined for each model property M under study. This is clearly seen in Figs. 13,14,15, where the frontier line in the three top graphs, associated, respectively, with M = J, ∆K sym , K 0 , differs from one graph to the other. If the dispersion of the data points can be mainly attributed to differences in their respective M values, the diagram [∆X] Λ (M ) must present a clear correlation. In this case, it would be sufficient to know the values of L and M of a given model to predict accurately the corresponding X value. On the contrary, if the diagram [∆X] Λ (M ) appears completely uncorrelated, we can deduce that the coefficient M under study is not responsible for the dispersion of the data X i (L i ). In several cases, we find an intermediate situation, in which the correlation of [∆X] Λ (M ) is weak but allows to associate eccentric values of M with large values of ∆X.

L (MeV)

The most favorable situation appears in Fig. 13, with the effect of J on ∆Y p,tt . This effect was expected : Y p,tt depends on the value of the symmetry energy at sub-saturation, which is well correlated with L as long as the different models have a similar symmetry energy at ρ 0 . J also appears to affect the values of ρ tt (Fig. 14) and P tt (Fig. 15), although the correlation [∆X] Λ (J) is weaker in these two cases. Let us now consider the eccentricity of the K sym behavior, namely, ∆K sym . It has a clear effect on the relation P tt (L), as we see in Fig. 15; this confirms the analysis of the previous section, where we have underlined the role of the L-K sym correlation in the link between L and P tt . On the other hand, ∆K sym is uncorrelated with the position of the transition, ρ tt and Y p,tt . Finally, the isoscalar incompressibility K 0 has no clear effect on the core-crust transition. It appears completely uncorrelated with the values of Y p,tt and P tt . A weak correlation appears with ρ tt , but this does not affect significantly the quality of the L-ρ tt correlation. In the following, we will concentrate exclusively on the role of isovector coefficients.

B. Prediction of the dynamical core-crust transition

Until now, we have studied the quantities X tt = {ρ tt , Y p,tt , P tt }, which are defined by the crossing between the βequilibrium condition and the thermodynamic spinodal. This framework allowed us to emphasize the analytical role of bulk GLDM coefficients in the transition. However, realistic descriptions of the core-crust transition involve stability comparison between homogeneous matter and clusterized matter. Equilibrium calculations have been performed, e.g., in Refs. [2,[15][16][17]; it has been verified that the resulting transition can be very well approximated by the crossing between the β-equilibrium condition and the dynamic spinodal. In Fig. 16, we illustrate the differences between the thermodynamic and dynamic spinodals and identify with a square (dynamic) and dot (thermodynamic) the crust-core transition, defined by the crossing of the β-equilibrium EOS and the spinodal. We denote X td = {ρ td , Y p,td , P td } as the quantities taken at this dynamic spinodal border. We have obtained these quantities within the effective Skyrme and relativistic approaches. The dynamic spinodal has not been calculated in the BHF framework, since the BHFbased density functional does not include the density-gradient terms needed to modelize the surface effects associated with finite-size density fluctuations.

To check that the thermodynamic framework effectively reflects the correlations between the GLDM coefficients and the core-crust transition, we have to make sure that the transformation from X tt to X td does not destroy these correlations. This is verified in Fig. 17, where we plot the dynamic results as a function of the thermodynamic ones: we observe that these quantities are strongly correlated. In the following, we will extract some empirical relations between the GLDM coefficients and the dynamic core-crust properties, {ρ td , Y p,td , P td }.

A systematic analysis of the effect of the isovector coefficients J, L, and K sym on the transition properties is presented in Table III. Our previous study has shown that the relation X tt (L) is affected by atypical values of J and K sym associated with the various models. To explore this effect, we have performed two-dimensional fits of the transition data:

X ti (M 1 , M 2 ) = a 1 M 1 + a 2 M 2 + b , (16)
where X ti = {ρ tt , Y p,tt , P tt ; ρ td , Y p,td , P td } (dynamic or thermodynamic transition), and M i are two of the isovector GLDM coefficients:

J ref = S(ρ ref ), L ref = 3ρ ref [∂S/∂ρ](ρ ref ), and K sym,ref = (3ρ ref ) 2 [∂ 2 S/∂ρ 2 ](ρ ref ).
We have considered the saturation coefficients J, L, and K sym , as well as coefficients at the reference density ρ ref = 0.1 fm -3 , denoted J 01 , L 01 , and K sym,01 . Table III gives the root mean square (rms) of residuals associated with the different fits, indicating the relevance of the respective combinations of coefficients in the determination of X ti . It appears that ρ ti is well correlated with L; no significant improvement can be obtained by considering pairs of coefficients. The quality of the L-ρ ti correlation can be understood as a consequence of Eq. ( 9), as discussed in Section III A. In the cases of Y p,ti and P ti however, the predictions can be considerably improved by using combinations of coefficients. As expected, Y p,ti is very well correlated with a combination of J and L, and this result is still improved using a combination of J 01 and L 01 . On the other hand, the values of Tab. III indicate that combinations of L and K sym are not relevant to determine Y p,ti . The transition pressure, instead, presents improved correlations with two kinds of parameter combinations: either J and L at saturation density, or L 01 and K sym,01 ; the latter leads to the smallest rms of residuals.

We show in Fig. 18 the most significant correlations obtained between the dynamic transition properties and different combinations of the GLDM coefficients. The quantities ρ td , Y p,td , and P td are shown as a function of L and as a function of selected coefficient combinations. For the transition density, we verify the good L-ρ td correlation; we also consider two coefficient combinations, [J -0.558 × L] and [L 01 + 0.426 × K sym,01 ], which lead to a similar dispersion, lower than ±0.004 fm -3 . The situation is different for the proton fraction Y p,td and the pressure P td , for which the correlation with L is not so good. The selected coefficient combinations lead to clearly improved correlations. Concerning the proton fraction, the combinations [J -0.172 × L] and especially [J 01 -0.160 × L 01 ] present excellent correlations with Y p,td . As for the transition pressure, the combinations [J -0.127 × L] and [L 01 -0.343 × K sym,01 ] allow to considerably reduce the data dispersion, and show unambiguous correlations with P td ; this is especially true for the second combination (coefficients extracted at 0.1 fm -3 ), for which the typical model dispersion for P td becomes ± 0.033 MeV.fm 3 instead of ± 0.085 MeV.fm 3 in the case of P td (L). The linear fits represented in Fig. 18

The relations involving GLDM coefficients at ρ=0.1 fm -3 allow to predict the core-crust transition density, asymmetry and pressure within a reasonable model uncertainty. This indicates that exploring the sub-saturation properties of nuclei in order to constrain directly these low-density coefficients could allow to modelize core-crust transition properties that do not depend strongly on the type of nuclear functional that is used.

V. CONCLUSION

In this paper, we have studied to what extent the core-crust transition properties can be predicted by using a reduced set of phenomenological constraints.

Different kinds of nuclear models have been compared: Skyrme and relativistic effective models, and a microscopic BHF approach. It is seen that the EOS obtained in the relativistic cases present much more variability in their density dependence than the Skyrme ones, due to different ways of describing the interaction; notably, a softer symmetry energy is obtained with density-dependent couplings, and the inclusion of the delta meson leads to an atypic density evolution of the symmetry energy slope, as was noticed in Ref. [26]. A more regular behavior is observed within the Skyrme sets of parameters; this has the drawback to bring possibly spurious correlations between the EOS properties at saturation and subsaturation density. As for the BHF calculations, although microscopic approaches are indispensable to provide realistic predictions for the EOS away from the phenomenological constraints, their predictions for the core-crust transition properties are very sensitive to the numerical fit of the EOS, which is necessary to determine its curvature. However, other ways could be followed to constrain the curvature properties from microscopic approaches, such as the study of Landau parameters [29,30].

To relate the predictions of the different models with their properties at a fixed density, we have introduced a generalized liquid-drop model (GLDM) which consists of a density development of the EOS around a reference density ρ ref , up to a chosen order. When ρ ref is the saturation density, we have seen that a development up to order 3 is necessary to get reasonably close to the thermodynamic core-crust transition properties predicted by the complete functionals. This means that the correlations that can be observed between the symmetry-energy slope at saturation L and the core-crust transition properties are subjected to further correlations existing between the various coefficients of the GLDM. Such correlations reflect two kinds of effects: (i) the possible existence of effective constraints at subsaturation densities, arising from the nuclear data used in the fit of most of the effective models, and (ii) the regularity of the functional shapes, which depend on the construction of the model. We have also considered a development at a lower reference density, ρ ref = 0.1 fm -3 . This approach has the advantage of reducing the model dependence arising from specific functional shapes, and of focusing on a density region closer to most nucleus observations. A development at order 2 around ρ = 0.1 fm -3 allows us to characterize the thermodynamic core-crust transition within nearly the smallest uncertainty allowed by a GLDM approach. This smallest uncertainty is given by the infinite development D ∞ , which gives the best GLDM approximation of the complete functional by neglecting only the extra-parabolic terms in the isospin dependence of the nuclear interaction.

We have also presented a more detailed study of the relation between L and the core-crust transition, the conclusions of which confirm our previous analysis [12]. Namely, the core-crust transition density ρ tt and proton fraction Y p,tt appear clearly correlated with L, despite the variety of models, while the link between L and the transition pressure P tt is much more sensitive to model dependence. Indeed, the impact of L on the transition pressure involves several opposite contributions, which tend to compensate each other; thus, it is not possible to establish a qualitative prediction for the evolution P tt (L), which can change sign depending on the model (see, for instance, the opposite predictions presented in Refs. [10,11]).

To explore the possibility to overcome the model dispersion and predict the core-crust transition properties from a reduced set of nuclear constraints, we have searched those GLDM coefficients others than L that play a major role in the determination of this transition: J and K sym were found to have a significant responsibility in the observed dispersion. Finally, we have addressed the case of the dynamic core-crust transition, given by the crossing between the dynamic spinodal and the β equilibrium. This corresponds to a realistic approximation of the actual core-crust transition, and takes place at lower density than in the thermodynamic approach. We have verified that the dynamic transition is related to the GLDM coefficients by similar correlations. The L-ρ td correlation is quite good, and cannot be significantly improved by considering other coefficients; however, the predictivity of Y p,td and P td is considerably better in terms of selected pairs of coefficients. An excellent correlation appears between Y p,td and a combination of J and L, and it is even better using a combination of J 01 and L 01 (coefficients defined at ρ ref = 0.1 fm -3 ). Furthermore, the model dependence in the prediction of the transition pressure can be considerably reduced if we consider a combination of L 01 and K sym,01 . In this case, an unambiguous correlation is obtained within all the variety of Skyrme and relativistic models considered.

To conclude, it appears that an accurate determination of the first three GLDM coefficients at ρ ref = 0.1 fm -3 would allow a prediction of the core-crust transition properties that do not depend much on the model construction. This gives a strong motivation to focus on the relation between nuclear observables and GLDM coefficients at subsaturation density. It will become possible to use phenomenological nuclear models to restrict the range of the core-crust transition properties in neutron stars; this would have an impact on the interpretation of astrophysical observations, and on the possible scenarios to explain phenomena such as pulsar glitches [13,14]. As discussed in [13], the transition pressure is an essential input to infer the neutron-star mass-radius relation from glitch observations. Since the mass-radius relation predicted by a given EOS is mainly determined by its high-density region, an accurate prediction of the transition pressure would also constrain the high-density EOS.

where, denoting q the third component of the isospin (n for neutrons and p for protons), the kinetic densities are defined by τ q = k2 q . The coefficients C and D, associated respectively with the isoscalar and isovector contributions, are linear combinations of the traditional Skyrme parameters: C 0 = 3t 0 /8 , D 0 = -t 0 (2x 0 + 1)/8 , C 3 = t 3 /16 , D 3 = -t 3 (2x 3 + 1)/48 , C eff = [3t 1 + t 2 (4x 2 + 5)]/16 , D eff = [t 2 (2x 2 + 1)t 1 (2x 1 + 1)]/16 .

(A8)

In this paper, we have considered 21 Skyrme parametrizations, commonly used in the literature, chosen in order to cover a wide range of L values while presenting acceptable saturation properties. Traditionally, Skyrme parameters are fitted in order to reproduce selected nuclear properties measured in a set of nuclei: basically masses and radii, plus several other input with increasing level of sophistication. SV [33] is among the earlier parametrizations. SGII [34], for which spin properties have also been used as constraints, can reproduce isospin effects in giant dipole resonances. R σ and G σ [35] consider spin-orbit splitting in 16 O and surface widths. SkMP [36] was built to improve the fit of 208 Pb charge distribution. The series SkI2, SkI3, SkI4, SkI5 [37] and SkI6 [38] include constraints on the isotope shifts of the charge radius in Pb and Ca. SkO [39] further considers isotopic evolution of two-neutron separation energies in Pb. In addition to nuclear data constraints, many Skyrme forces include in their fitting procedure the neutron matter EOS from microscopic calculations: the objective is to obtain a reliable behavior of the density functional at high isospin asymmetry, especially for astrophysical applications. RATP [40] was the first parametrization using this procedure, including the neutron matter calculation by Friedman and Pandharipande [41]. The Skyrme-Lyon forces SLy230a, SLy230b [32] , SLy4 [42] and SLy10 [43] use the pure neutron matter equation of state UV14+UVII by R.B. Wiringa et al [44]. NRAPR [45] (Non-Relativistic APR) stands for the Skyrme interaction parameters obtained from a fit of the APR equation of state (Akmal-Pandharipande-Ravenhall, Ref. [46]). LNS [47] is based on Brueckner-Hartree-Fock calculations of infinite nuclear matter at different values of isospin asymmetry. The Bruxelles-Skyrme forces BSk14 [48], BSk16 [49] and BSk17 [50] include the Friedman and Pandharipande calculation of neutron matter [41], and a HFB treatment of pairing effects in order to improve mass predictions in the neutron-drip region.

Relativistic effective models

In this paper, we consider two kinds of relativistic effective approaches: RMF models, which have constant coupling parameters described by the Lagrangian density of non-linear Walecka models (NLWM), and DDH models with density-dependent coupling parameters. In each case, we consider models including or not the δ meson, which have been introduced to include in the isovector channel the same symmetry existing already in the isoscalar channel with the meson pair (σ, ω) responsible for saturation in RMF models [51]. The presence of the δ meson softens the symmetry energy at subsaturation densities and hardens it above saturation density. The RMF parametrizations we use are NL3 [52], TM1 [53], GM1, GM3 [54], FSU [55], NLωρ [56], and NLρδ [57]. The DDH parametrizations are TW [58], DD-ME1, DD-ME2 [59], and DDHδ [60]. The models NLρδ(0) and NLρδ(2.5) introduced in [57] have the same isoscalar properties and the same symmetry energy at saturation; however, the last model includes the δ meson with g δ = 2.5, while in the first one, the δ coupling was set to zero. We have also introduced the parametrization NLρδ(1.7) with a weaker δ-meson coupling (1.7 instead of 2.5). The model NLωρ(025) includes a ωρ non-linear term in the Lagrangian as in [56] with strength Λ v = 0.025. The parametrization DDHδI-25 introduced in [60] has a quite low symmetry energy at saturation (25 MeV); therefore, we also consider the parametrization DDHδII-30 where the ρ-meson coupling was adjusted so that, at saturation, the symmetry energy is 30 MeV, and all the isoscalar properties are kept fixed.

The relativistic approach is based on a Lagrangian density given by:

L = i=p,n L i +L σ +L ω +L ρ + L ωρ + L δ . ( A9
)
The nucleon Lagrangians read as

L i = ψi [γ µ iD µ -M * ] ψ i , (A10)
with

iD µ = i∂ µ -Γ v V µ - Γ ρ 2 τ • b µ , (A11)
M * = m -Γ s φ -Γ δ τ • δ , (A12)

FIG. 1 :

 1 FIG.1:(Color online) Correlations between different GLDM coefficients cIV,n. Left : at saturation density, relation between L and other coefficients cIV,n(ρ0), namely, J (top), Ksym (center) and Qsym (bottom). Right : relation between coefficients cIV,n defined at reference density ρ0 or ρ=0.1 fm -3 . Results are shown for different Skyrme models (full symbols), relativistic models (empty symbols), and BHF (star). BHF-1 is used at saturation density, and BHF-2 is used at ρ=0.1 fm -3 .

FIG. 2 :FIG. 3 :

 23 FIG.2:(Color online) Skyrme models: Left: symmetry-energy value (top), slope (center), and curvature (bottom) as a function of the density, calculated with complete functionals. BHF results are indicated for comparison (dots). Right: the corresponding variance for a calculation using the complete functionals (full line) and the GLDM up to different orders: D1(ρ0) (dotted-dashed line), D2(ρ0) (2dotted-dashed line), and D3(ρ0) (3dotted-dashed line).

7 FIG. 4 : 1 FIG. 5 :

 7415 FIG. 4: (Color online) Predictive power of the GLDM development. For different Skyrme and relativistic models, we compare the values of ρtt (top), Yp,tt (center) and Ptt (bottom) versus L, calculated with the complete functionals and with several versions of the corresponding GLDM (see text). Left: D∞, second and third order development around saturation density (D2,3(ρ0)). Right: D∞ and second order development around ρ = 0.1 fm -3 (D2(0.1)).

FIG. 6 :

 6 FIG.6: (Color online) Effect of L on the proton fraction at β equilibrium for a fixed density ρ=0.08 fm -3 (left) and at thermodynamic spinodal crossing (right), for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2para (asterisk).

FIG. 7 :

 7 FIG.7:(Color online) Effect of L on the density at thermodynamic spinodal border for a fixed proton fraction Yp = 0.025 (left) and at β equilibrium (center). Right: correspondence between L and the energy-density curvature of neutron matter at symmetric spinodal density, CNM,s (see text). Results are shown for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2para (asterisk) (both BHF results are identical for neutron matter).

FIG. 8 :

 8 FIG.8:(Color online) Simultaneous effect of L on the proton fraction at β equilibrium and on the thermodynamic spinodal border of neutron-rich matter. Left: spinodal contours and β equilibrium for three Skyrme models: BSk17 (L=36 MeV), SkI6 (L=60 MeV), Rσ (L=86 MeV). Right: correlation between ρtt and Yp,tt for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2para (asterisk).

1 FIG. 9 :

 19 FIG.9:(Color online) Impact of L on the pressure of neutron-rich matter at sub-saturation density, for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2para (asterisk). Both BHF-2 results are identical for neutron matter. Left: pressure of pure neutron matter, for a fixed density typical of the transition: ρ=0.08 fm -3 . Right: thermodynamic transition pressure.

FIG. 10 :

 10 FIG.10:(Color online) Estimation of different contributions to the variation of Ptt with L (see text). Left: separated contributions, due to the transition position shift (indices p1, p2) and to the variation of the GLDM coefficients cIV,n in the expression of P (indices e1, e2, and e3 for n=1, 2, and 3 respectively). Right: sum of the contributions, considering the contribution of δcIV,n/δL up to order n=1, 2 and 3 (respective index: tot1, tot2, tot3). The horizontal lines indicate the region where the pressure variation is compatible with zero within the estimated uncertainty. Full symbols: Skyrme models; empty symbols: relativistic models.

FIG. 11 :

 11 FIG.11:(Color online) Typical variation for the values of J, Ksym(L) and K0. We represent these quantities for different Skyrme models (full symbols), relativistic models (empty symbols), and BHF-1 (star). The lines show the different values which are used for the standard schematic model (see text).

FIG. 12 :

 12 FIG.12:(Color online) Predictions of the standard schematic model for the transition density, proton fraction and pressure, compared with complete effective models: Skyrme (full symbols) and relativistic (empty symbols) models. Left: varying the symmetry energy at saturation, J. Center: varying the relation Ksym(L). Right: varying the incompressibility at saturation K0.

FIG. 13 :

 13 FIG.13: Relation between the dispersion of Yp,tt(L) and the dispersion of GLDM properties (J, K0, ∆Ksym(L)) associated with the different models (see text). Nuclear models: Skyrme (full symbols), relativistic (empty symbols), and BHF-2 (star).

FIG. 14 :

 14 FIG.14: Relation between the dispersion of ρtt(L) and the dispersion of GLDM properties (J, K0, ∆Ksym(L)) associated with the different models (see text). Nuclear models: Skyrme (full symbols), relativistic (empty symbols), and BHF-2 (star).

FIG. 15 :

 15 FIG.15: Relation between the dispersion of Ptt(L) and the dispersion of GLDM properties (J, K0, ∆Ksym(L)) associated with the different models (see text). Nuclear models: Skyrme (full symbols), relativistic (empty symbols), and BHF-2 (star).

FIG. 16 :

 16 FIG.16: Comparison between the thermodynamic (dashed line) and dynamic (full line) spinodals. The dotted line represents the β-equilibrium EOS and the square and dot define the crust-core transition within, respectively, the dynamical and thermodynamical spinodal.

7 FIG. 17 :

 717 FIG.17:(Color online) Comparison between the transition taken at dynamic (td) and thermodynamic (tt) spinodal for different Skyrme (full symbols) and relativistic (empty symbols) models: density (left), proton fraction (center), and pressure (right).

  are: ρ td (L) = (-3.75 × 10 -4 × L + 0.0963) fm-3 (17)ρ td (J, L) = (7.46 × 10 -4 × [J -0.558 × L] + 0.0754) fm -3 (18) ρ td (L 01 , K sym,01 ) = (3.23 × 10 -4 × [L 01 + 0.426 × K sym,01 ] + 0.0802) fm -3(19)Y p,td (L) = -3.15 × 10 -4 × L + 0.0461(20)Y p,td (J, L) = 2.69 × 10 -3 × [J -0.172 × L] -0.0290(21)Y p,td (J 01 , L 01 ) = 3.30 × 10 -3 × [J 01 -0.160 × L 01 ] -0.0253(22)P td (L) = (-2.42 × 10 -3 × L + 0.465) MeV fm -3 (23) P td (J, L) = (3.36 × 10 -2 × [J -0.127 × L] -0.474) MeV fm -3(24) P td (L 01 , K sym,01 ) = (9.59 × 10 -3 × [L 01 -0.343 × K sym,01 ] -0.328) MeV fm-3

FIG. 18 :

 18 FIG.18:(Color online) Correlations between the dynamic core-crust transition properties and GLDM coefficients (see text).

functional fit includes calculation points in the density range

	1	0.187 -15.23 195.50 -280.90 34.30 66.55 -31.30 -112.80 -12.72 -17.58 62.62 22.77 40.14 -28.45
	BHF-2		-12.74 -17.59 67.01 22.77 41.35 -44.30
	Skyrme
	BSk14	0.159 -15.86 239.38 -358.78 30.00 43.91 -152.03 388.30 -13.92 -20.65 118.58 23.29 41.94 -89.79
	BSk16	0.159 -16.06 241.73 -363.69 30.00 34.87 -187.39 461.93 -14.10 -20.87 119.86 24.11 39.44 -109.16
	BSk17	0.159 -16.06 241.74 -363.73 30.00 36.28 -181.86 450.52 -14.10 -20.87 119.87 23.98 39.83 -106.13
	Gσ	0.158 -15.59 237.29 -348.82 31.37 94.02 13.99 -26.77 -13.71 -20.31 118.13 20.03 58.46	6.62
	Rσ	0.158 -15.59 237.41 -348.50 30.58 85.70 -9.13 22.23 -13.71 -20.34 117.99 20.05 55.22	-6.21
	LNS	0.175 -15.32 210.83 -382.67 33.43 61.45 -127.37 302.48 -12.96 -20.07 95.92 23.20 48.07 -66.20
	NRAPR	0.161 -15.86 225.70 -362.65 32.78 59.63 -123.33 311.63 -13.93 -19.89 112.26 24.18 48.86 -72.09
	RATP	0.160 -16.05 239.58 -349.94 29.26 32.39 -191.25 440.74 -14.05 -20.78 116.94 23.55 38.04 -108.29
	SV	0.155 -16.05 305.75 -175.86 32.82 96.10 24.19 47.97 -13.86 -24.15 135.63 21.60 60.45	5.26
	SGII	0.158 -15.60 214.70 -381.02 26.83 37.62 -145.92 330.44 -13.84 -18.91 111.74 20.98 37.16 -83.35
	SkI2	0.158 -15.78 240.99 -339.81 33.38 104.35 70.71 51.60 -13.88 -20.54 119.16 21.17 61.18	22.73
	SkI3	0.158 -15.99 258.25 -303.96 34.83 100.53 73.07 211.53 -13.96 -21.66 122.84 23.03 59.52	11.04
	SkI4	0.160 -15.95 247.98 -331.26 29.50 60.40 -40.52 351.09 -13.88 -21.33 118.11 21.48 43.25 -43.97
	SkI5	0.156 -15.85 255.85 -302.05 36.64 129.34 159.60 11.71 -13.93 -21.22 124.33 22.33 71.01	62.14
	SkI6	0.159 -15.92 248.65 -327.44 30.09 59.70 -47.27 378.96 -13.90 -21.23 119.31 22.19 43.70 -48.82
	SkMP	0.157 -15.57 230.93 -338.15 29.89 70.31 -49.82 159.44 -13.76 -19.72 116.03 20.95 49.62 -32.93
	SkO	0.160 -15.84 223.39 -392.98 31.97 79.14 -43.17 131.12 -13.93 -19.87 113.35 21.63 53.60 -27.51
	Sly230a	0.160 -15.99 229.94 -364.29 31.98 44.31 -98.21 602.92 -14.06 -20.16 114.48 25.43 39.40 -86.38
	Sly230b	0.160 -15.98 229.96 -363.21 32.01 45.96 -119.72 521.54 -14.06 -20.11 114.87 25.15 41.58 -88.09
	SLy4	0.160 -15.98 229.97 -363.22 32.00 45.94 -119.74 521.58 -14.06 -20.11 114.87 25.15 41.56 -88.09
	SLy10	0.156 -15.91 229.74 -358.43 31.98 38.74 -142.19 591.28 -14.16 -19.58 118.84 26.15 39.37 -104.21
	Relativistic
	NL3	0.148 -16.24 270.70 188.80 37.34 118.30 100.50 182.60 -14.64 -20.31 136.42 25.08 73.73	29.99
	TM1	0.145 -16.26 280.40 -295.40 36.84 110.60 33.55 -65.20 -14.68 -21.43 152.80 25.58 73.86	14.72
	GM1	0.153 -16.32 299.70 -222.10 32.48 93.87 17.89 25.77 -14.25 -23.78 141.49 21.74 60.30	2.94
	GM3	0.153 -16.32 239.90 -515.50 32.48 89.66 -6.47 55.86 -14.56 -20.79 135.45 22.07 59.46	-8.12
	FSU	0.148 -16.30 229.20 -537.40 32.54 60.40 -51.41 426.60 -14.81 -19.47 140.97 25.57 47.54 -71.94
	NLωρ(025)	0.148 -16.24 270.70 188.80 32.35 61.05 -34.36 1322.00 -14.64 -20.31 136.42 25.23 50.12 -99.76
	TW	0.153 -16.25 240.20 -541.00 32.76 55.30 -124.70 539.00 -14.48 -21.16 140.67 25.38 48.59 -92.51
	DD-ME1	0.152 -16.23 244.50 307.60 33.06 55.42 -101.00 706.30 -14.66 -18.29 114.43 25.86 48.06 -96.35
	DD-ME2	0.152 -16.14 250.90 478.30 32.30 51.24 -87.19 777.10 -14.58 -17.98 109.23 25.65 44.83 -98.43
	DDHδI-25	0.153 -16.25 240.20 -540.30 25.62 48.56 81.10 928.30 -14.48 -21.16 140.70 20.23 31.73 -48.55
	DDHδII-30	0.153 -16.25 240.20 -540.30 31.89 57.52 80.74 1005.00 -14.48 -21.16 140.70 25.42 38.36 -60.57
	NLρδ(2.5)	0.160 -16.05 240.40 -470.20 30.71 102.70 127.20 282.90 -13.98 -21.65 125.20 18.77 55.76	33.28
	NLρδ(1.7)	0.160 -16.05 240.40 -470.20 30.70 97.14 86.46 202.80 -13.98 -21.65 125.20 19.16 55.12	21.05
	NLρδ(0)	0.160 -16.05 240.40 -470.20 30.34 84.51	3.33 61.40 -13.98 -21.65 125.20 19.77 53.05	-4.90
	TABLE I: GLDM coefficients at saturation density ρ0 and at ρ ref = 0.1 fm -3 for different nuclear models. Details and
	references concerning these models are given in Appendix. BHF-1: functional fit includes calculation points in the density
	range ρ = [0.1; 0.35] fm -3 ; BHF-2:

ρ = [0.05; 0.18] fm -3 . The coefficients are named following traditional notations. Isoscalar coefficients at saturation density: E0 = cIS,0, K0 = cIS,2, Q0 = cIS,3. Isovector coefficients at saturation density:

  Ptt L01 Ksym,01 0.036 MeV.fm -3 P td L01 Ksym,01 0.033 MeV.fm-3

	Xtt M1 M2	rms	X td M1 M2	rms
	ρtt L	0.0038 fm -3	ρ td L	0.0037 fm -3
	ρtt J L	0.0028 fm -3	ρ td J L	0.0035 fm -3
	ρtt L Ksym	0.0038 fm -3	ρ td L Ksym	0.0037 fm -3
	ρtt L01	0.0062 fm -3	ρ td L01	0.0054 fm -3
	ρtt J01 L01	0.0037 fm -3	ρ td J01 L01	0.0038 fm -3
	ρtt L01 Ksym,01	0.0032 fm -3	ρ td L01 Ksym,01	0.0037 fm -3
	Yp,tt L	0.0063	Y p,td L	0.0057
	Yp,tt J L	0.0022	Y p,td J L	0.0024
	Yp,tt L Ksym	0.0063	Y p,td L Ksym	0.0058
	Yp,tt L01	0.0081	Y p,td L01	0.0074
	Yp,tt J01 L01	0.0016	Y p,td J01 L01	0.0014
	Yp,tt L01 Ksym,01	0.0040	Y p,td L01 Ksym,01	0.0040
	Ptt L	0.117 MeV.fm -3 P td L	0.085 MeV.fm -3
	Ptt J L	0.076 MeV.fm -3 P td J L	0.055 MeV.fm -3
	Ptt L Ksym	0.092 MeV.fm -3 P td L Ksym	0.067 MeV.fm -3
	Ptt L01	0.129 MeV.fm -3 P td L01	0.101 MeV.fm -3
	Ptt J01 L01	0.088 MeV.fm -3 P td J01 L01	0.069 MeV.fm -3

TABLE III :

 III Root mean square of residuals associated with two-variable linear fits Xti = a1 × M1 + a2 × M2 + b, where Xti = {ρtt, Yp,tt, Ptt; ρ td , Y p,td , P td } and Mi are two of the isovector GLDM coefficients:J ref = S(ρ ref ), L ref = 3ρ ref [∂S/∂ρ](ρ ref ) and K sym,ref = (3ρ ref ) 2 [∂ 2 S/∂ρ 2 ](ρ ref ). J, Land Ksym are taken at ρ ref = ρ0; J01, L01 and Ksym,01 are taken at ρ ref = 0.1 fm -3 . The results for single-variable linear fits in function of L and L01 are also shown.

Acknowledgments This work was partially supported by the ANR NExEN contract, FCT (Portugal) under grants SFRH/BPD/46802/2008, FCOMP-01-0124-FEDER-008393 with FCT reference CERN/FP/109316/2009, PTDC/FIS/64707/2006, PTDC/FIS/113292/2009 and COMPSTAR, an ESF Research Networking Programme.

Appendix A: Nuclear models

In this appendix, we give an overview of the different types of nuclear models whose properties are compared in the present work. Each of these models gives at every density the energy of symmetric matter and the symmetry energy (in contrast with the GLDM exposed in the Section II A, which is based on a density expansion).

Skyrme-like effective models

The local Skyrme interaction [31] allows us to define an energy density functional H(r) such that the total energy for a system of nucleons in a Slater determinant | ψ > reads as

For homogeneous, spin-saturated matter with no Coulomb interaction, the Skyrme energy density functional is composed of four terms:

In this expression, K is the kinetic-energy term, H 0 is a density-independent two-body term, H 3 is a density-dependent term, and H eff is a momentum-dependent term [32]:

We have introduced the isoscalar and isovector particle densities, ρ and ρ 3 , as well as kinetic densities, τ and τ 3 :

where τ is the isospin operator. We use the vector symbol to designate a vector in isospin space.

The isoscalar part is associated with the scalar sigma (σ) field φ and the vector omega (ω) field V µ , while the isospin dependence comes from the isovector-scalar delta (δ) field δ i and the isovector-vector rho (ρ) field b i µ (where µ is a space-time index and i an isospin-direction index). The associated Lagrangians are:

where

, and Γ j and m j are, respectively, the coupling parameters of the mesons j = s, v, δ, ρ with the nucleons and their masses. The self-interacting terms for the σ meson are included only for the NL3 and NLδ parametrizations, with κ and λ denoting the corresponding coupling constants.

The density-dependent coupling parameters Γ s , Γ v , and Γ ρ , are adjusted in order to reproduce some of the nuclear matter bulk properties, using the following parametrization:

with

where x = ρ/ρ sat and

The values of the parameters m i , Γ i , a i , b i , c i , and d i , i = s, v, ρ, for TW and DD-ME2 are, respectively, given in [58] and [61] and for DDHδ in [60,62]. In this last case, the parametrization for the δ-and ρ-coupling parameters is also given by (A13) with

The Γ i coupling parameters are replaced by the g i coupling constants in the NL3 and NLδ models.

The BHF approach of asymmetric nuclear matter

The BHF approach of asymmetric nuclear matter [63,64] starts with the construction of all the G matrices describing the effective interaction between two nucleons in the presence of a surrounding medium. They are obtained by solving the well-known Bethe-Goldstone equation

where τ = n, p indicates the isospin projection of the two nucleons in the initial, intermediate, and final states, V denotes the bare N N interaction, Q τiτj is the Pauli operator that allows only intermediate states compatible with the Pauli principle, and ω, the so-called starting energy, corresponds to the sum of non-relativistic energies of the interacting nucleons. The single-particle energy ǫ τ of a nucleon with momentum k is given by

where the single-particle potential U τ ( k) represents the mean field "felt" by a nucleon due to its interaction with the other nucleons of the medium. In the BHF approximation, U ( k) is calculated through the "on-shell energy" G-matrix, and is given by

where the sum runs over all neutron and proton occupied states and where the matrix elements are properly antisymmetrized. We note here that the so-called continuous prescription has been adopted for the single-particle potential when solving the Bethe-Goldstone equation. As shown in Refs. [65,66], the contribution to the energy per particle from three-hole line diagrams is minimized in this prescription. Once a self-consistent solution of Eqs. (A16) and (A18) is achieved, the energy per particle can be calculated as

The BHF calculation carried out in this paper uses the realistic Argonne V18 (Av18) [67] nucleon-nucleon interaction supplemented with a three-body force of Urbana type, which (for the use in BHF calculations) was reduced to a twobody density-dependent force by averaging over the third nucleon in the medium [68]. This three-body force contains two parameters that are fixed by requiring that the BHF calculation reproduces the energy and saturation density of symmetric nuclear matter. We note that more microscopically based three-body forces without adjustable parameters have been recently constructed (see Refs. [69][70][71] for a recent analysis of the use of three-body forces in nuclear and neutron matter). We note also that the Av18 interaction contains terms that break explicitly isospin symmetry. Therefore, in principle, we should consider also odd powers of y in the expansion (1) for the Brueckner calculation. However, we have neglected such terms since, as shown by Müther et al. in Ref. [72], the effects of isospin symmetry breaking on the symmetry energy are quite weak (less than 0.5 MeV for a wide range of N N interactions).