A Nuclear models

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 is mainly correlated with the symmetry energy slope , (b) the proton fraction with the symmetry energy and symmetry energy slope defined at saturation density, or, even better, with the same quantities defined at fm, and (c) the transition pressure with the symmetry energy slope and curvature defined at fm.

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 is one of the central issues in this field (1); (2); (3). Although the symmetry energy at saturation density is considered to be well known ( MeV), the different nuclear models present a wide range of predictions for the symmetry energy slope . Several experimental observables have been proposed to obtain a measure of , 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 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 and the core-crust transition properties: density , proton fraction and pressure . It appeared that and are unambiguously correlated with , while the link between and 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 and laboratory data.

In this paper, we investigate the role of symmetry-energy properties other than 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 (15); (2); (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 and the core-crust transition, which was first presented in Ref. (12). In Sec. IV, we perform a correlation analysis involving coefficients others than . 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 and asymmetry . 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 , (ii) the reference density , and (iii) the order of the development. The order and the reference density () can be chosen at will. In this paper, the GLDM coefficients are derived from the complete EOS given by various nuclear models. In a different context, the GLDM coefficients could be nuclear properties determined experimentally at , 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.

ii.1 GLDM equation of state

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


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 . 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 and , respectively associated with the derivatives of the energy and of the symmetry energy : the index ’IS’ (’IV’) stands for isoscalar (isovector). They are expressed as:


In the case , the lower-order coefficients are usual nuclear matter properties: (saturation energy), (incompressibility), , (symmetry energy), (symmetry-energy slope), (symmetry incompressibility), and .

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:


where is the kinetic density of the nucleon species and is the nucleon mass. As a function of density and asymmetry, we have:


where is the parabolic part of . In the GLDM defined by Eq. (1), the extra-parabolic behavior of the functional is sketched by the extra-parabolic behavior of , which brings the model-independent correction (second term on the r.h.s. of Eq. (1)).

In the following, we will use the notation to identify a development of order around the density . Developments of this kind are usually considered up to =3, around the saturation density fm (see, e.g., Refs. (22); (23)). In this work, we will also consider an extreme situation labeled . Performing an infinite development gives the exact value of and for any density , whatever the choice of . In practice, the equation of state is simply built using the exact expressions of and to obtain :


Thus, the difference between a complete functional and its associated model is just the extra-parabolic content of the nuclear interaction, which has not been taken into account by the correction defined above. Let us note that a parabolic approximation is often assumed in BHF calculations, to interpolate between the symmetric-matter and neutron-matter EOS. We will call this approach BHF. However, by default, our BHF results include the correction ; in this case, the complete functional is exactly equivalent to the corresponding model.

ii.2 Correlations between GLDM coefficients

[fm] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV] [MeV]
BHF-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
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
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
DDHI-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
DDHII-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 1: GLDM coefficients at saturation density and at fm 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 fm; BHF-2: functional fit includes calculation points in the density range fm. The coefficients are named following traditional notations. Isoscalar coefficients at saturation density: , , . Isovector coefficients at saturation density: , , , . Isoscalar coefficients at fm: , , . Isovector coefficients at fm: , , .
Figure 1: (Color online) Correlations between different GLDM coefficients . Left : at saturation density, relation between and other coefficients , namely, (top), (center) and (bottom). Right : relation between coefficients defined at reference density or =0.1 fm. 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.

The isoscalar and isovector GLDM coefficients obtained with all the nuclear models used in this paper are reported in Table  1. 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 , is likely to provide effective constraints in the density region 0.1–0.12 fm (this point is raised, e.g., in Refs. (24); (25)); as a result, the coefficients defined at 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 ) 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, ; the different plots shows the correlation between and , , and . As we have emphasized in Ref. (12), there is a strong - correlation. The two eccentric points correspond to the relativistic models DDHI-25 and DDHII-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 , it tends to increase with , although the dispersion between different models is important. The behaviour of the third-derivative coefficient is less universal: a clear decreasing - correlation is obtained among the Skyrme models, but the relativistic models of lower and the BHF point are completely out of this trend. Linear fits are represented on these plots; in the case of and , 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 and =0.1 fm, for =0, 1, 2. The correlations obtained in the cases of (symmetry energy slope) and (symmetry-energy curvature) reflect the shape similarities between many of the density functionals. The absence of correlation in the case of (symmetry energy) is due to the fact that 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: , , and . Note that and are defined here with a constant factor involving : they are not equivalent to the coefficients and , except at . We use here a constant factor in order to represent quantities proportional to the derivatives of . The curves are shown for all models under consideration, using the complete functionals. The comparison between Figs. 2 and 3 shows that the relativistic models present more variability in the shape of the density functional than Skyrme models, 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 . 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.

Figure 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: (dotted-dashed line), (2dotted-dashed line), and (3dotted-dashed line).
Figure 3: (Color online) Relativistic 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: (dotted-dashed line), (2dotted-dashed line), and (3dotted-dashed line).

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 fm, taking values 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 : the different curves show a marked tendency to cross at about . In contrast, no convergence effect appears for the second derivative, .

On the right panel of Figs. 2 and 3, we represent the density dependence of the variance between the values taken by different models. Denoting , the variance is given by:


where the index runs over the 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 occurs at fm. With relativistic models, we observe a plateau around an inflexion point at fm. Note that 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 ; both Skyrme and relativistic models present a sharp minimum of at fm.

These results indicate that the correlations existing between , , , and may be associated with effective constraints on the values of and . To verify this, we show on the same figures the variances obtained when the curves are calculated using the GLDM at different orders, with reference density . A minimum of in the density range 0.11–0.13 fm is obtained with versions and of the GLDM. The version gives a clear minimum of at fm. We can conclude that the correlations between the GLDM coefficients are related with convergence effects for and in the subsaturation density region.

Let us remark that the convergence of the values at 0.06 fm is easily interpreted as a geometrical consequence of a constraint on . Indeed, let us imagine a constraint fixing the values such that all models have to verify . Since we also have the condition , it turns out that for all models, takes the value at least once on the interval . For instance, if fm, MeV and fm, we have MeV.

ii.3 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, , , and , where the index 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 will be used. The transition properties in both thermodynamic and dynamic cases are given in Table 2.

[fm] [MeVfm] [fm] [MeVfm]
BHF-1 0.061 0.023 0.193
BHF-1 0.083 0.026 0.400
BHF-2 0.078 0.027 0.370
BHF-2 0.094 0.028 0.571
BSk14 0.090 0.033 0.483 0.081 0.030 0.381
BSk16 0.096 0.037 0.502 0.087 0.035 0.402
BSk17 0.095 0.036 0.499 0.086 0.034 0.397
G 0.063 0.013 0.278 0.054 0.010 0.172
R 0.067 0.014 0.312 0.058 0.012 0.202
LNS 0.088 0.031 0.614 0.077 0.028 0.469
NRAPR 0.083 0.034 0.545 0.073 0.030 0.413
RATP 0.097 0.037 0.500 0.086 0.034 0.390
SV 0.071 0.021 0.372 0.061 0.016 0.235
SGII 0.086 0.026 0.401 0.077 0.024 0.311
SkI2 0.064 0.014 0.291 0.054 0.011 0.170
SkI3 0.071 0.022 0.363 0.062 0.018 0.244
SkI4 0.081 0.024 0.332 0.072 0.021 0.234
SkI5 0.061 0.014 0.271 0.051 0.010 0.149
SkI6 0.082 0.026 0.352 0.073 0.024 0.257
SkMP 0.072 0.020 0.357 0.062 0.017 0.241
SkO 0.073 0.020 0.413 0.062 0.017 0.270
Sly230a 0.090 0.039 0.404 0.081 0.037 0.319
Sly230b 0.089 0.038 0.462 0.080 0.036 0.362
SLy4 0.089 0.038 0.461 0.080 0.036 0.361
SLy10 0.091 0.042 0.447 0.083 0.041 0.369
NL3 0.065 0.021 0.422 0.054 0.016 0.236
TM1 0.070 0.025 0.511 0.060 0.020 0.324
GM1 0.074 0.019 0.408 0.067 0.016 0.290
GM3 0.069 0.018 0.356 0.063 0.016 0.267
FSU 0.082 0.037 0.487 0.074 0.035 0.385
NL(025) 0.089 0.038 0.689 0.080 0.034 0.530
TW 0.084 0.037 0.544 0.075 0.033 0.384
DD-ME1 0.085 0.038 0.605 0.070 0.033 0.404
DD-ME2 0.087 0.039 0.594 0.072 0.034 0.409
DDHI-25 0.085 0.022 0.144 0.079 0.021 0.100
DDHII-30 0.086 0.038 0.285 0.080 0.037 0.231
NL(2.5) 0.062 0.009 0.173 0.057 0.008 0.116
NL(1.7) 0.064 0.011 0.197 0.059 0.009 0.143
NL(0) 0.069 0.014 0.276 0.063 0.012 0.197
Table 2: Density, proton fraction, and pressure at the thermodynamic transition () and dynamic transition (), for different nuclear models. BHF-1: functional fit includes calculation points in the density range fm. BHF-2: functional fit includes calculation points in the density range fm.
Figure 4: (Color online) Predictive power of the GLDM development. For different Skyrme and relativistic models, we compare the values of (top), (center) and (bottom) versus L, calculated with the complete functionals and with several versions of the corresponding GLDM (see text). Left: , second and third order development around saturation density (). Right: and second order development around fm ().

We compare in Fig. 4 the transition properties (, , and ) 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, and ; on the right panel, we consider the second-order development around a lower reference density fm, namely, . On each part, we add the results from the fully developped GLDM, , 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 model.

We can see that leads generally to an important underestimation of , , and . This means that, in a development around saturation density, terms beyond order 2 (i.e., beyond ) have a large impact on the properties of the core-crust transition; the correlations observed between and these properties can occur only if higher-order corrections are either correlated with , or similar, for most of the functionals. The underestimation is strongly attenuated, but still present, with , which involves the knowledge of isovector coefficients until . The situation is much improved if we use a development around fm. The accuracy is globally better with than with ; furthermore, it is nearly as good with as with , which means that, at fm, 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,…).

Figure 5: (Color online) Comparison between exact and GLDM versions of three functionals based on BHF calculations: LNS (47) (Skyrme-type force including constraints from BHF), BHF-1 (fit of the BHF EOS including densities in the range fm), and BHF-2 (fit of the BHF EOS including lower densities, in the range fm).

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 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 , , and obtained by using the BHF calculations in different ways. In the first version, BHF-1, the fit is performed on the density interval fm: this is the version that is used to establish the GLDM coefficients at . 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 fm: 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 fm have been determined for both versions, and there are significant differences in the transition properties predicted by the 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 core-crust transition properties; phenomenological constraints from finite nuclei will be essential to improve our knowledge of the low-density EOS.

Iii Correlation between and the core-crust transition point

We consider in this section the specific role of the symmetry energy slope 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.

iii.1 Position of the transition point

Figure 6: (Color online) Effect of on the proton fraction at equilibrium for a fixed density =0.08 fm (left) and at thermodynamic spinodal crossing (right), for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2 (asterisk).
Figure 7: (Color online) Effect of on the density at thermodynamic spinodal border for a fixed proton fraction (left) and at equilibrium (center). Right: correspondence between and the energy-density curvature of neutron matter at symmetric spinodal density, (see text). Results are shown for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2 (asterisk) (both BHF results are identical for neutron matter).
Figure 8: (Color online) Simultaneous effect of 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 (=36 MeV), SkI6 (=60 MeV), R (=86 MeV). Right: correlation between and for different nuclear models: Skyrme (full symbols), relativistic (empty symbols), BHF-2 (star) and BHF-2 (asterisk).

The most direct impact of the symmetry energy on neutron-star structure concerns the proton fraction 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 can be correlated with a low value of the symmetry energy at sub-saturation, we expect higher- models to provide lower values of . This point is illustrated by Fig. 6. At fixed density fm, for increasing , 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 on the transition density , illustrated by Fig. 7. The correlation between and is a well-known feature (9); however, its explanation is less intuitive than in the case of the - correlation. Furthermore, it can not be explained just as a consequence of the behavior of , as we see on the left panel of the figure: even for a fixed proton fraction , the spinodal border shows a clear decreasing correlation with . This feature can be understood as a consequence of the strong link existing between and the energy-density curvature of neutron matter taken at symmetric spinodal density :


with . Since all models yield a symmetric spinodal density close to 0.1 fm, so that , the term is nearly canceled: this reinforces the dominance of in the determination of .

Figure 8 illustrates how 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 , and .

iii.2 Pressure at the transition point

The link between and the core-crust transition pressure 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 , which is the leading coefficient. This trend appears on the left panel of Fig. 9, representing the relation between and the pressure of pure neutron matter, , at =0.08 fm. Thus, a positive correlation between and should be obtained if we could neglect the density shift due to the - 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 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 =60 MeV weaken this correlation between Skyrme models: they correspond to atypical relations between and , which also affect the plot .

Figure 9: (Color online) Impact of 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-2 (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. Right: thermodynamic transition pressure.
Figure 10: (Color online) Estimation of different contributions to the variation of with (see text). Left: separated contributions, due to the transition position shift (indices p1, p2) and to the variation of the GLDM coefficients in the expression of (indices e1, e2, and e3 for n1, 2, and 3 respectively). Right: sum of the contributions, considering the contribution of 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.

The lack of correlation between and 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 : (i) variations occurring at a fixed density , resulting only from the modifications of the coefficients in Eq. (10), which defines ; and (ii) variations due to a shift , for a fixed expression of , i.e. frozen values of the coefficients in Eq. (10). The contributions of the first kind come from the explicit dependence of Eq. (10), and from correlations between and higher-order coefficients . In practice, we will consider the following terms:


where the index 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 . We distinguish the respective effects of total density and asymmetry:


The variations of the quantities depending on are fixed empirically, using as a reference the correlations that are observed between different models. From linear fits, we extract =-3.84 10 MeV fm, =6.08 10 MeV, =3.33, =-6.63. Note that, in the case of , the correlation with 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 essentially results from the balance between three terms: , which is large and positive, is compensated by the conjugated effect of and . These two negative contributions, e2 and p1, are of the same order of magnitude: this means that the correlation - has the same importance as the correlation - in explaining why we do not observe an increasing correlation .

In addition, we can see that the term due to the - 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 , as it was observed on Fig. 4, the - correlation is not crucial in the determination of ; in other words, the third-order correction does not depend strongly on . On the other hand, let us notice that the strong dispersion of values in the case of relativistic models is bound to cause a strong dispersion in .

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

iii.3 Predictions of a standard GLDM

Figure 11: (Color online) Typical variation for the values of , and . 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).
Figure 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, . Center: varying the relation . Right: varying the incompressibility at saturation .

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 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 model characterized by the following parameters:

The coefficient varies in the interval MeV, and determines and according to the relations:

The predictions of this standard model for , , and are represented in Fig. 12. As expected, we obtain a clear decrease of and with , while the evolution of 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 , the incompressibility , and the - relation. On the left panel of Fig. 12, we show the effect of varying between 29 and 33 MeV. On the right panel, we show the effect of varying between 220 and 260 MeV. On the central panel, we modify the linear relation between and by adopting different values of : a softer version , and a stiffer one .

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 and is maintained: although the absolute value can be affected by different aspects of the functional, they always unambiguously decrease with increasing . (ii) The qualitative behavior of 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: either increases or decreases, and most often it is quite flat. These two conclusions confirm the previous analysis.

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 . First, we check to what extent the dispersion in the dependence of , , and can be attributed to specific GLDM coefficients, such as the symmetry energy at saturation , the incompressibility , or the quantity . This last quantity characterizes the eccentricity of the model with respect to the typical relation . 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.

iv.1 Role of the first GLDM coefficients in the dispersion

Figure 13: Relation between the dispersion of and the dispersion of GLDM properties (, , ) associated with the different models (see text). Nuclear models: Skyrme (full symbols), relativistic (empty symbols), and BHF-2 (star).
Figure 14: Relation between the dispersion of and the dispersion of GLDM properties (, , ) associated with the different models (see text). Nuclear models: Skyrme (full symbols), relativistic (empty symbols), and BHF-2 (star).
Figure 15: Relation between the dispersion of and the dispersion of GLDM properties (, , ) associated with the different models (see text). Nuclear models: Skyrme (full symbols), relativistic (empty symbols), and BHF-2 (star).

As we have seen in the previous section, the correlations between and the core-crust transition properties , , and 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 and being very sensitive to the details of the functional, the model dispersion destroys the possibility to deduce the value of from a measurement of . In order to look for model properties that may be responsible for this dispersion, we first use the following procedure:

  • We call the GLDM model property whose effects are investigated: namely, , or .

  • We represent the quantities calculated with the different models, with = (Fig. 13), (Fig. 14), and (Fig. 15).

For each property (, or ), we check whether the diagram can be separated in two regions associated with larger/smaller values of . For this:

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

  • we calculate the distance of each data point to this frontier: .

  • we vary the frontier in order to obtain the best correlation for the diagram .

In this way, a different line is defined for each model property under study. This is clearly seen in Figs. 13, 14, 15, where the frontier line in the three top graphs, associated, respectively, with , differs from one graph to the other.

If the dispersion of the data points can be mainly attributed to differences in their respective values, the diagram must present a clear correlation. In this case, it would be sufficient to know the values of and of a given model to predict accurately the corresponding value. On the contrary, if the diagram appears completely uncorrelated, we can deduce that the coefficient under study is not responsible for the dispersion of the data . In several cases, we find an intermediate situation, in which the correlation of is weak but allows to associate eccentric values of with large values of .

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

iv.2 Prediction of the dynamical core-crust transition

Figure 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.
Figure 17: (Color online) Comparison between the transition taken at dynamic () and thermodynamic () spinodal for different Skyrme (full symbols) and relativistic (empty symbols) models: density (left), proton fraction (center), and pressure (right).

Until now, we have studied the quantities , 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. (15); (2); (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 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 BHF-based 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 to 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, .

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


where (dynamic or thermodynamic transition), and are two of the isovector GLDM coefficients: , , and . We have considered the saturation coefficients , , and , as well as coefficients at the reference density fm, denoted , , and . Table 3 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 .

It appears that is well correlated with ; no significant improvement can be obtained by considering pairs of coefficients. The quality of the - correlation can be understood as a consequence of Eq. (9), as discussed in Section III.1. In the cases of and however, the predictions can be considerably improved by using combinations of coefficients. As expected, is very well correlated with a combination of and , and this result is still improved using a combination of and . On the other hand, the values of Tab. 3 indicate that combinations of and are not relevant to determine . The transition pressure, instead, presents improved correlations with two kinds of parameter combinations: either and at saturation density, or and ; 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 , , and are shown as a function of and as a function of selected coefficient combinations. For the transition density, we verify the good - correlation; we also consider two coefficient combinations, and , which lead to a similar dispersion, lower than fm. The situation is different for the proton fraction and the pressure , for which the correlation with is not so good. The selected coefficient combinations lead to clearly improved correlations. Concerning the proton fraction, the combinations and especially present excellent correlations with . As for the transition pressure, the combinations and allow to considerably reduce the data dispersion, and show unambiguous correlations with ; this is especially true for the second combination (coefficients extracted at 0.1 fm), for which the typical model dispersion for becomes 0.033 MeV.fm instead of 0.085 MeV.fm in the case of . The linear fits represented in Fig. 18 are: