Corecrust transition in neutron stars: predictivity of density developments
Abstract
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 corecrust 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 rareisotope 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 neutronskin 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 corecrust transition in neutron stars (9); (10); (11). In a previous work (12), we have addressed the model dependence of the link between and the corecrust 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 symmetryenergy 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 liquiddrop model (GLDM): our objective is to establish to what extent the corecrust 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 corecrust transition.
In this paper, the corecrust 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 corecrust 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 finitesize 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 corecrust transition properties, which are strongly correlated with the thermodynamic ones.
This paper is organized as follows. The generalized liquiddrop model is presented in Sec. II. We discuss the correlations existing between various GLDM coefficients and evaluate the accuracy of the corecrust transition predictions obtained with different GLDM expansions. In Sec. III, we present in more details the analysis of the link between the symmetryenergy slope at saturation and the corecrust 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 Liquiddrop 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 corecrust transition predicted by the complete functional; in other words, to what extent the role of higherorder 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:
(1) 
The first term on the righthand 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:
(2)  
(3) 
In the case , the lowerorder coefficients are usual nuclear matter properties: (saturation energy), (incompressibility), , (symmetry energy), (symmetryenergy 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 neutronrich region. The reason is that the energydensity curvature in the protondensity 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 protondensity direction is constant, and leads to the unphysical prediction of unstable neutron matter. To avoid this discrepancy, we have introduced in Eq. (1) a modelindependent correction based on the nonrelativistic, free Fermi gas kinetic term:
(4) 
where is the kinetic density of the nucleon species and is the nucleon mass. As a function of density and asymmetry, we have:
(5)  
(6) 
where is the parabolic part of . In the GLDM defined by Eq. (1), the extraparabolic behavior of the functional is sketched by the extraparabolic behavior of , which brings the modelindependent 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 :
(7) 
Thus, the difference between a complete functional and its associated model is just the extraparabolic 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 symmetricmatter and neutronmatter 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
Model  

[fm]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  [MeV]  
Microscopicc  
BHF1  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 
BHF2  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 
DDME1  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 
DDME2  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 
DDHI25  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 
DDHII30  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 
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 finitenuclei 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 referencedensity 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 DDHI25 and DDHII30: 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 thirdderivative 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 (symmetryenergy 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.
Considering globally the Skyrme and relativistic functionals, we can notice two common convergence regions. One concerns the symmetryenergy: 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 symmetryenergy 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:
(8) 
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 symmetryenergy 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 corecrust transition obtained with the complete model on which it has been built. This procedure allows to estimate whether the corecrust transition can be efficiently characterized by a limited set of EOS properties determined at a fixed density.
The discussion of the corecrust 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 finitesize spinodal border, in which case the index will be used. The transition properties in both thermodynamic and dynamic cases are given in Table 2.
Model  

[fm]  [MeVfm]  [fm]  [MeVfm]  
Microscopic  
BHF1  0.061  0.023  0.193  
BHF1  0.083  0.026  0.400  
BHF2  0.078  0.027  0.370  
BHF2  0.094  0.028  0.571  
Skyrme  
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 
Relativistic  
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 
DDME1  0.085  0.038  0.605  0.070  0.033  0.404 
DDME2  0.087  0.039  0.594  0.072  0.034  0.409 
DDHI25  0.085  0.022  0.144  0.079  0.021  0.100 
DDHII30  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 
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 thirdorder developments around saturation density, namely, and ; on the right panel, we consider the secondorder 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 extraparabolic 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 corecrust transition; the correlations observed between and these properties can occur only if higherorder 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 symmetryenergy density dependence directly in this lowdensity region. This is also appropriate, since the nucleus properties used to constrain the symmetry energy are often associated with subsaturation densities (neutronskin thickness, resonances, multifragmentation, isospin diffusion,…).
Let us notice that effective nuclear models have an important role to play for the finetuning of the curvature properties of the EOS, crucial for the prediction of the corecrust 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 corecrust transition. This figure displays the values of , , and obtained by using the BHF calculations in different ways. In the first version, BHF1, 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, BHF2, for which the fit is performed on the density interval fm: this is the version we will use afterwards to define the BHF corecrust transition. For now, let us compare the BHF1 and BHF2 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 BHF1 and BHF2 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 lowdensity EOS.
Iii Correlation between and the corecrust transition point
We consider in this section the specific role of the symmetry energy slope in the determination of the corecrust 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
The most direct impact of the symmetry energy on neutronstar 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 subsaturation, 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 wellknown 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 energydensity curvature of neutron matter taken at symmetric spinodal density :
(9) 
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 neutronrich region. Due to the typical geometry of these respective lines, the two effects reinforce eachother, leading to a robust correlation between , and .
iii.2 Pressure at the transition point
The link between and the corecrust transition pressure is more problematic. In order to make this link explicit, let us write the pressure in the GLDM framework:
(10) 
From this expression, we expect that for a given density the pressure of neutronrich 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 higherorder 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 .
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 higherorder coefficients . In practice, we will consider the following terms:
(11)  
(12)  
(13) 
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:
(14)  
(15) 
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 thirdorder 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 thirdorder 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 5080 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 modeldependence, as it will be further discussed in the following.
iii.3 Predictions of a standard GLDM
In order to study how the different GLDM coefficients can affect the corecrust 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 corecrust transition
In this last section, we explore the possibility to reduce the model dispersion in the prediction of the corecrust 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 corecrust 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 corecrust transition.
iv.1 Role of the first GLDM coefficients in the dispersion
As we have seen in the previous section, the correlations between and the corecrust 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 .
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 subsaturation, 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 corecrust 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 corecrust transition
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 corecrust 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 crustcore 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 BHFbased density functional does not include the densitygradient terms needed to modelize the surface effects associated with finitesize density fluctuations.
To check that the thermodynamic framework effectively reflects the correlations between the GLDM coefficients and the corecrust 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 corecrust 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 twodimensional fits of the transition data:
(16) 
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:
(17)  
(18)  
(19)  
(20)  
(21)  