Bayesian analysis of the crust-core transition with a compressible liquid-drop model

TL;DR

This paper uses Bayesian methods and a compressible liquid-drop model to analyze the neutron star crust-core transition, highlighting the dominant role of surface properties and providing constrained estimates of transition density and pressure.

Contribution

It introduces a Bayesian framework combined with a CLD model to quantify uncertainties in the crust-core transition parameters of neutron stars, emphasizing the impact of surface properties.

Findings

Transition density estimated as 0.072±0.011 fm⁻³.

Transition pressure estimated as 0.339±0.115 MeV fm⁻³.

Surface properties strongly influence transition characteristics.

Abstract

The crust-core phase transition of neutron stars is quantitatively studied within a unified meta-modelling of the nuclear Equation of State (EoS). The variational equations in the crust are solved within a Compressible Liquid Drop (CLD) approach, with surface parameters consistently optimized for each EoS set on experimental nuclear mass data. When EoS parameters are taken from known Skyrme or RMF functionals, the transition point of those models is nicely reproduced. A model-independent probability distribution of EoS parameters and of the transition density and pressure is determined with a Bayesian analysis, where the prior is given by an uncorrelated distribution of parameters within the present empirical uncertainties, and constraints are applied both from neutron star physics and ab-initio modelling. We show that the characteristics of the transition point are largely independent of the high density properties of the EoS, while ab-initio EoS calculations of neutron and symmetric matter are far more constraining. The most influential parameter for the determination of the transition point governs the surface properties of extremely neutron rich matter, and it is strongly unconstrained. This explains the large dispersion of existing predictions of the transition point. Only if the surface tension is fixed to a reasonable but somewhat arbitrary value, strong correlations with isovector parameters ($L_{sym},K_{sym}$ and $Q_{sym}$) are recovered. Within the present experimental and theoretical uncertainties on those parameters, we estimate the transition density as $n_t= 0.072\pm 0.011$ fm$^{-3}$ and the transition pressure as $P_t=0.339\pm0.115$ MeV fm$^{-3}$.

Tables

Table 1. Table 1: Transition density n t subscript 𝑛 𝑡 n_{t} (top) and transition pressure P t subscript 𝑃 𝑡 P_{t} (bottom) for several interactions. Meta-modelling unified EoS calculations with p = 2.5 𝑝 2.5 p=2.5 , p = 3 𝑝 3 p=3 , p = 3.5 𝑝 3.5 p=3.5 are given together with the compilation by Ducoin et al. from Ref.  Ducoin2011 using the dynamical ( n t ​ d , P t ​ d subscript 𝑛 𝑡 𝑑 subscript 𝑃 𝑡 𝑑 n_{td},P_{td} ) and thermodynamical ( n t ​ t , P t ​ t subscript 𝑛 𝑡 𝑡 subscript 𝑃 𝑡 𝑡 n_{tt},P_{tt} ) methods.

	$n_t$ (fm-3)
Model	$n_{td}$	$n_{tt}$	$p=2.5$	$p=3.0$	$p=3.5$
BSk14	0.081	0.090	0.073	0.079	0.085
BSk16	0.087	0.096	0.080	0.087	0.092
BSk17	0.086	0.095	0.078	0.085	0.091
NRAPR	0.073	0.083	0.062	0.071	0.076
RATP	0.086	0.097	0.078	0.087	0.092
SkO	0.062	0.073	0.050	0.061	0.064
SLy230a	0.081	0.090	0.076	0.080	0.085
SLy230b	0.080	0.089	0.073	0.079	0.083
SLy4	0.080	0.089	0.073	0.079	0.083
NL3	0.054	0.065	0.046	0.054	0.061
TM1	0.060	0.070	0.049	0.058	0.062
DD-ME1	0.070	0.085	0.064	0.076	0.083
DD-ME2	0.072	0.087	0.071	0.081	0.087

Table 2. Table 2: Minimum value and maximum value of each of the empirical parameters for the prior distribution (prior) and average and standard deviation of each of the empirical parameters of the posterior distribution after application of the HD(LD) filter (see text).

Parameter	Unit	Prior		HD		LD
		Min	Max	Average	$\sigma$	Average	$\sigma$
$n_{sat}$	fm-3	0.15	0.17	0.1600	0.0060	0.1641	0.0049
$E_{sat}$	MeV	-17	-15	-16.01	0.61	-15.29	0.25
$K_{sat}$	MeV	190	270	229	24	234	23
$Q_{sat}$	MeV	-1000	1000	200	535	-31	362
$Z_{sat}$	MeV	-3000	3000	1038	1233	-146	1728
$E_{sym}$	MeV	26	38	33.53	3.48	30.71	0.76
$L_{sym}$	MeV	10	80	45.45	17.97	43.66	3.68
$K_{sym}$	MeV	-400	200	-92	136	-202	42
$Q_{sym}$	MeV	-2000	2000	913	740	-253	673
$Z_{sym}$	MeV	-5000	5000	1463	2216	-114	2868
$m_{sat}^{\ast }/m$		0.6	0.8	0.70	0.06	0.70	0.06
$\Delta m_{sat}^{\ast }/m$		0.0	0.2	0.10	0.06	0.10	0.06
$b$		1	10	5.3	2.7	5.2	2.6

Table 3. Table 3: Average value and standard deviation of the transition density (pressure) with the prior parameter distribution, the posterior using the HD filter, and the posterior using the LD filter for p = 3 𝑝 3 p=3 and p = { 2.5 , 3 , 3.5 } 𝑝 2.5 3 3.5 p=\{2.5,3,3.5\} .

	$n_t$ (fm-3)		$P_t$ (MeV/fm3)
	Average	$\sigma$	Average	$\sigma$
Prior ($p=3$)	0.078	0.040	0.342	0.426
HD ($p=3$)	0.076	0.032	0.394	0.327
LD ($p=3$)	0.074	0.011	0.360	0.122
LD	0.065	0.021	0.307	0.167

Table 4. Table 4: Neutron drip density n d subscript 𝑛 𝑑 n_{d} , crust-core transition density n t subscript 𝑛 𝑡 n_{t} and pressure P t subscript 𝑃 𝑡 P_{t} , and density n μ ∗ superscript subscript 𝑛 𝜇 n_{\mu}^{*} at which muons appear for the meta-model realization of the seven Drischler et al. functionals. Empty cells when no solution is found.

Model	$n_d$ (fm-3)	$n_t$ (fm-3)	$P_t$ (MeV/fm3)	$n_{\mu }^{\ast }$ (fm-3)
1	$2.55\times {10}^{-4}$	0.0785	0.5833	0.1275
2	$2.82\times {10}^{-4}$	0.0769	0.4982	0.1289
3	$2.88\times {10}^{-4}$	0.0758	0.4575	0.1318
4	$2.88\times {10}^{-4}$	0.0740	0.4266	0.1300
5	$1.59\times {10}^{-4}$	0.0320	0.0847
6	$2.82\times {10}^{-4}$
7	$2.68\times {10}^{-4}$	0.0374	0.1045	0.1244

Equations

$$\text{mathpzc}F(A,I,n_0,n_p,n_g)=\frac{E_{nuc}}{V_{WS}}+\left(1-\frac{A}{n_0{V}_{WS}}\right){ϵ}_g+{ϵ}_{el}-\mu n_B,$$

$$E_{nuc}=\frac{ϵ(n_{0p},n_{0n})}{n_0}A+E_c+E_s,$$

$$n_B=n_g+\frac{A}{V_{WS}}\left(1-\frac{n_g}{n_0}\right),$$

$$\displaystyle\frac{\partial(E_{nuc}/A)}{\partial A}\bigg{|}_{I,n_{0},n_{p},n_{g}}=0,$$

$$\displaystyle\frac{2}{A}\frac{\partial E_{nuc}}{\partial I}\bigg{|}_{A,n_{0},n_{p},n_{g}}=\mu_{el}-n_{p}\frac{\partial(E_{c}/A)}{\partial n_{p}}\bigg{|}_{A,I,n_{0}},$$

$$\displaystyle\frac{E_{nuc}}{A}+\frac{1-I}{A}\frac{\partial E^{nuc}}{\partial I}\bigg{|}_{A,n_{0},n_{p},n_{g}}-\frac{\epsilon_{g}}{n_{0}}=\mu\left(1-\frac{n_{g}}{n_{0}}\right),$$

$$\displaystyle{n_{0}}^{2}\frac{\partial(E_{nuc}/A)}{\partial n_{0}}\bigg{|}_{A,I,n_{p},n_{g}}=n_{g}\mu-\epsilon_{g},$$

$$\mu=\frac{2n_{p}}{n_{0}A(1-I)-2n_{p}}\frac{\partial E_{s}}{\partial n_{g}}\bigg{|}_{A,I,n_{0}}+\frac{d\epsilon_{g}}{dn_{g}}.$$

$${\mu }_n^{nuc}$$

$$P_{nuc}$$

$$E_s$$

$$\mu_{n}^{nuc}=\frac{\partial E^{nuc}}{\partial N}\bigg{|}_{Z}+\frac{P_{g}}{n_{0}},$$

$$ϵ(n,\delta )$$

$$e(n,\delta )$$

$$f_1(\delta )$$

$$f_2(\delta )$$

$$E_s(A,Y_p)=4\pi r_{sat}^2{A}^{2/3}\sigma (Y_p),$$

$$\sigma (Y_p)={\sigma }_0\frac{2^{p+1}+b_s}{Y_p^{-p}+b_s+(1-Y_p{)}^{-p}}.$$

$$p_{prior}(\vec{X})=\prod _{k=1}^{2(N+1)+3}f(X_k^{min},X_k^{max};X_k)$$

$$p_{post}(\vec{X})=\text{mathpzc}N{w}_{LD(HD)}(\vec{X})e^{-{\chi }^2(\vec{X})/2}{p}_{prior}(\vec{X}).$$

$$p(Y)=\prod _{k=1}^{2(N+1)+3}{\int }_{X_k^{min}}^{X_k^{max}}d{X}_{k}Y(\vec{X})p_{post}(\vec{X}),$$

Full text

I Introduction

Neutron stars (NS) are a unique observable laboratory of the different phases of hadronic matter Haensel_book . The structure and composition of the inner core is still not completely clear report_Micaela , but it is however well established that a phase transition occurs from a solid crust to a liquid core at some $\approx 1$ km from the surface of the star. The presence of a solid crust plays an important role in a number of phenomena involving NS, from the cooling of proto-neutron stars Page2013 ; Shternin2007 ; Brown2009 to the irregularities (”glitches”) in their rotational motion Haskell2015 ; Link1999 ; Steiner2015 . The precise location of the transition in the star is also important for the theoretical determination of its static properties, notably the star radius Fortin2016 , which accurate measurement will be soon available with upcoming x-ray observations Watts2016 . To understand all these aspects of NS physics, a reliable theoretical estimation of the crustal thickness and its uncertainty is necessary.

For slowly rotating NS, the crustal thickness can be computed with the Tolman-Oppenheimer-Volkov (TOV) equation of hydrostatic equilibrium, if the equation of state (EoS) and the core-crust (CC) transition point is known in density $n_{t}$ and pressure $P_{t}$. A large number of studies have been devoted to the determination of the transition point with different relativistic Horowitz2001 ; Klahn2006 ; Moustakidis2010 ; Fattoyev2010 ; Cai2012 ; Pais2016 and non relativistic Vidana2009 ; Li2016 ; Ducoin2011 ; Routray2016 ; Gonzales2017 models. Most of these studies compute the transition ”from the core”: the transition is defined as the density point where homogeneous nuclear matter becomes unstable with respect to density fluctuations. The simplest version of this technique consists in evaluating the thermodynamical spinodal, that is the instability point of neutral nuclear matter with respect to the nuclear liquid-gas (LG) phase transition. The advantage of this method is that it only requires the knowledge of the energy functional of homogeneous nuclear matter: this functional, at the density of the CC point which is close to nuclear saturation density, is characterized by a small set of parameters which can be strongly constrained through nuclear experiments and/or microscopic ab-initio calculations Tsang2012 ; Lattimer2013 ; Dutra2014 ; Fortin2016 .

However, the CC transition is very different from the LG one Ducoin2007 ; Chomaz2007 and the thermodynamical spinodal gives only a qualitative (and overestimated) estimation of the transition point. The CC transition occurs at the pressure where the energy density of clusterized matter (the solid cristal embedded in a electron and neutron gas) overcomes the energy density of uniform nuclear matter Baym1971 . The clusterized phase is inhomogeneous and locally charged, and its equilibrium energy is determined by the competition between the Coulomb and the surface energy. These energy components identically vanish in uniform matter, and therefore do not enter in the determination of the thermodynamical spinodal.

A better estimation can be obtained looking at the dynamical response of the homogeneous system with respect to finite size fluctuations Pethick1995 . The resulting dynamical spinodal is available for a limited set of models Ducoin2007 ; Pais2016 ; Xu2009 ; Pearson2012 , and was recently studied by our group within a Bayesian meta-modelling technique Chatterjee ; Antic . One should remark that for the determination of the dynamical spinodal, in addition to the EoS of uniform matter, the isovector gradient terms - which are not well known - play a non negligible role, and its precise location slightly depends on the many body formalism adopted (linear response, Vlasov, or RPA) Ducoin2008 . Even if the dynamical spinodal certainly gives a better estimation of the transition than the thermodynamical one, it is worth mentioning that in the actual dynamical process of supernova collapse that gives birth to neutron stars, matter at subsaturation densities is never uniform but composed of atomic nuclei Lattimer1991 ; Shen1998 , meaning that spinodal decomposition is not the dynamical process leading to the formation of the crust.

For this reason, the most theoretically sound determination of the CC point consists in determining the transition ”from the crust”, by directly comparing the energy density of the two competing phases Baym1971 . This method demands an explicit modelling of clusterized matter in beta equilibrium, which is a complex quantum many body problem. Therefore, after the seminal work by Baym, Bethe and Pethick Baym1971 , only few works using modern energy functionals have been developed along this line Avancini2008 ; Goriely2010 ; Sharma2015 ; Douchin2001 ; Newton2013 ; Gulminelli2015 , and the problem of model dependence clearly arises. In particular, many works have been devoted to the EoS dependence of the transition point and in particular the effect of the $L_{sym}$ parameter Pearson2012 . However, the interplay between the EoS parameters and the the isovector surface tension coefficient, which from the microscopic viewpoint is determined by the isovector gradient terms in the energy functional mentioned above, has been seldom addressed Newton2013 .

In this paper, we will calculate the CC transition point ”from the crust”, by solving the variational equations for non-uniform matter within a compressible liquid drop (CLD) approach, in the same lines of Refs. Baym1971 ; Douchin2001 ; Gulminelli2015 . The distribution probability of the EoS parameters of uniform matter and the extra parameters associated to the cluster surface properties will be determined with a Bayesian analysis, using a fully uncorrelated flat prior and constraints from the low density effective field theory (EFT) modelling by Drischler et al. Drischler2016 . This will allow us presenting model independent estimations of the density and pressure of the transition point, and determining which are the most influential parameters governing the phase transition.

The theoretical uncertainties on the EoS and on the transition point propagate to global observabes of the neutron star such as the crust thickness and moment of inertia, which in turn may be linked to astrophysical observables such as the amplitude of neutron star glitches. Quantative predictions on these global observables were presented in a recent paper Thomas_prl . In the present work, we concentrate on the nuclear physics ingredients, namely the density and pressure of the transition point, and examine in greater details the influence of the different parameters, and the effectiveness of the different constraints.

A very similar Bayesian study was very recently and independently performed in Ref. Holt , for the computation of different quantities than the ones of the present work, namely radii and tidal polarizabilities. The functional expression chosen for the homogeneous matter EoS is not the same as in our work, but our posterior distributions for the EOS parameters are in very good agreement with the results of Ref. Holt , showing the reliability and generality of the meta-modelling technique.

The plan of the paper is as follows. In section II we describe the variational equations which are solved to determine the crust composition and the transition to the core. The meta-modelling technique from Ref. Margueron2018a ; Margueron2018b used for the uniform matter EoS and the expression of the cluster surface tension from Ref. Ravenhall1983 ; Lorentz1993 will also be shortly summarized. Section III demonstrates the ability of our meta-modelling technique to reproduce the published results of specific models. The crust composition and the CC point obtained with specific choices for the model parameters will be compared to the litterature, showing that the parameter space of our meta-modelling is large enough to cover existing functionals, and can thus be used for a Bayesian determination of the EoS parameters. Section IV presents a sensitivity analysis to the different EoS parameters, and we will show that, together with the slope of the symmetry energy $L_{sym}$, the curvature $K_{sym}$ is strongly influential in the determination of the CC pressure. The full Bayesian analysis is reported in section V, where the probability distribution for the transition observables is computed imposing to our uncorrelated prior to reproduce the band predictions in isospin-symmetric and neutron matter of Ref. Drischler2016 obtained from a many-body perturbation theory (MBPT) based on two and three-nucleon chiral EFT interactions at N3LO. We will show that these ab-initio calculations at low density are far more constraining than the astrophysical constraint at high density concerning the maximum mass of NS. A complete correlation study will also be presented, where the importance of the isovector surface energy will be underlined. Conclusions are drawn in section VI.

II Modelling inhomogeneous matter

II.1 Variational equations

The equilibrium configuration of inhomogeneous catalyzed matter is obtained following the standard variational formalism of Refs. Baym1971 ; Douchin2001 ; Gulminelli2015 . Using the Lagrange multipliers technique, the energy density in a Wigner-Seitz cell of volume $V_{WS}$ is minimized with the constraint of a given baryonic density $n_{B}=n_{p}+n_{n}$. The auxiliary function to be minimized reads:

[TABLE]

where $\epsilon_{g}=\epsilon(n_{p}=0,n_{n}=n_{g})$ ($\epsilon_{el}$) is the energy density of a pure uniform neutron (electron) gas at density $n_{g}$ ($n_{e}$), and the bulk interaction between the cluster and the neutron gas is treated in the excluded volume approximation. The cluster energy $E_{nuc}$ depends on the cluster atomic number $A$, isospin asymmetry $I=(N-Z)/A$ and density $n_{0}$, and also on the total electron density $n_{e}=n_{p}$ because of the electrostatic interaction with the electron gas, according to:

[TABLE]

where $n_{0n(p)}=n_{0}(1\pm I/2)$, $E_{s}$ is the surface energy to be discussed in section II.3 below, and we use the standard expression for the Coulomb energy $E_{c}$ from Ref. Baym1971 ; Gulminelli2015 . Minimizing with respect to the five independent variables $A,I,n_{0},n_{p},n_{g}$, and using the baryonic density constraint,

[TABLE]

leads to the following system of coupled differential equations:

[TABLE]

where the baryonic chemical potential $\mu$ results:

[TABLE]

We can see that, in the absence of a possible in-medium modification of the surface energy because of the external gas, the baryonic chemical potential can be identified with the chemical potential of the gas $\mu_{g}\equiv d\epsilon_{g}/dn_{g}$.

It is easy to show that equations (4)-(7) can be equivalently written as chemical and mechanical equilibrium equations between the cluster and the neutron and electron gas, supplemented by the Baym virial theorem Baym1971 :

[TABLE]

where $P_{g}=n_{g}\mu_{g}-\epsilon_{g}$, the cluster chemical potential is modified by the external neutron gas as:

[TABLE]

and the electrostatic interaction between protons in the cluster and the background electrons leads to a modification of the $\beta$-equilibrium condition, $\Delta\mu=-n_{p}\partial(E_{c}/A)/\partial n_{p}$.

Equations (4)-(7) can be numerically solved if the energy functional for homogeneous baryonic matter $\epsilon(n_{p},n_{n})$ and the surface energy $E_{s}$ are specified. These quantities are affected by strong uncertainties, especially in the isovector sector report_Micaela ; Steiner2013 ; Newton2013 ; Newton2014 . To quantify the uncertainty on the crust-core transition induced by our imperfect knowledge of the baryonic bulk and surface energy, we use for both quantities simple and flexible parametrized expressions with parameters whose variation embeds the present uncertainty on nuclear energetics. These functionals are presented in sections II.2 and II.3 below.

II.2 EoS meta-modelling

We parametrize the energy density of homogeneous nuclear matter with baryonic density $n=n_{n}+n_{p}$ and isospin asymmetry $\delta=(n_{n}-n_{p})/n$, using the meta-modelling technique of Refs. Margueron2018a ; Margueron2018b , here briefly summarized.

It is theoretically known that the energy density of homogeneous nuclear matter is analytic at least up to $\approx 2-3n_{sat}$, where $n_{sat}$ is the saturation density of symmetric matter. The integrality of the possible behaviors of the functional can therefore be explored using a Taylor expansion around the saturation point $(n=n_{sat},\delta=0)$, and largely varying the parameters of the expansion, that correspond to the well-known EoS empirical parameters Piekarewicz2009 :

[TABLE]

with $x=(n-n_{sat})/3n_{sat}$. To fasten the series convergence, the $\delta^{5/3}$ term coming from the fermionic zero point energy is explicitly added, as well as an exponential correction insuring the correct limiting behavior at zero density, see Ref. Margueron2018a for more details.

The final form for the energy per particle $e(n,\delta)=\epsilon/n$ at order $N$ of the expansion is given by:

[TABLE]

where the functions $f_{1},f_{2}$ give an effective correction to the parabolic approximation for the symmetry energy:

[TABLE]

The parameters $v_{m}^{is},v_{m}^{iv},a_{N}^{is},a_{N}^{iv}$ of the meta-functional eq.(15) are linear combinations of the successive derivatives $E_{q},L_{q},K_{q},Q_{q},Z_{q},\dots$ in the isoscalar ($q=sat$) and isovector ($q=sym$) sector. The saturation density $n_{sat}$, two parameters linked to the isoscalar effective mass ($\kappa_{sat}$) and effective mass splitting ($\kappa_{sym}$), and the $b$ parameter governing the functional behavior close to the $n\to 0$ limit Margueron2018a , complete the parameter set, which we will note in a compact form as $\vec{X}\equiv\{X_{n},n=1,\dots,2(N+1)+3\}$. It was verified in Ref. Margueron2018a that with an expansion up to $N=4$, the parameter space is sufficiently large to give an excellent reproduction of different popular Skyrme and RMF functionals. By varying the values of $\vec{X}$, eq.(15) thus provides a meta-functional that can continuously interpolate among existing functionals, and possibly explore novel density dependencies which have not yet been proposed in the literature.

A flat distribution of $\vec{X}$ values within physically reasonable intervals determined by empirical evidence Margueron2018a , will give our prior distribution for the EoS parameters. Ab-initio MBPT calculations of uniform nuclear matter Drischler2016 will be used in section V as a constraint, to determine a posterior distribution that will be used to estimate the model dependence of the crust-core transition calculation. A similar strategy was employed to compute different observables of neutron stars in Ref. Margueron2018b ; Thomas_prl ; Chatterjee ; Antic ; Holt , and of finite nuclei in Ref. Chatterjee2017 ; Holt .

II.3 Surface tension

To determine the crust-core phase transition, the energy density of homogeneous matter must be compared with the energy density of clusterized matter. The difference between the two is essentially given by the Coulomb $E_{c}$ and surface $E_{s}$ contribution (see eq.(2)), which identically vanish in homogeneous matter. The Coulomb energy only depends on the cluster charge $Z$ and density $n_{0}$, which are consistently obtained from the variational equations (4)-(7). On the other hand, to calculate the surface energy from the energy functional, strong approximations are needed in the many-body treatment Baym1971 ; Ravenhall1972 . Within these approximations, gradient terms must be specified, involving extra parameters which are added in the fitting protocol of phenomenological functionals Bender2011 . General physical considerations indicate that the surface tension $\sigma(A,I)=E_{s}/S$, with $S$ nuclear surface, should strongly depend on the isospin asymmetry $I$ (with $\sigma\to 0$ as $I\to 1$) Centelles1998 ; Danielewicz2003 and only moderately on the nuclear mass $A$ Durand1993 , but the exact value of $\sigma(I)$ is model dependent, and the uncertainty is particularly important for the extreme $I$ values encountered in the inner crust: indeed no experimental data exist on the surface energy of nuclei beyond the drip-line, which is in-medium modified by the presence of the neutron gas Centelles1998 ; Douchin2000 ; Sil2002 ; Aymard2014 . For this reason, parametrized expressions are generally employed Steiner2005 .

We use the expression originally proposed by Ravenhall et al. Ravenhall1983 on the basis of Thomas-Fermi calculations at extreme isospin ratios, and later employed in different works on neutron star and supernova modelling within the compressible liquid-drop model Lorentz1993 ; Lattimer1991 ; Newton2013 :

[TABLE]

with $r_{sat}=(3/4\pi n_{sat})^{1/3}$, $Y_{p}=Z/A$, and

[TABLE]

In this expression, the parameter $\sigma_{0}$ determines the value of the surface tension of symmetric nuclei, while $b_{s}$ governs the isospin dependence for moderate asymmetries. These parameters are fitted from experimental masses of the spherical magic and semi-magic nuclei 40,48Ca, 48,58Ni, 88Sr, 90Zr, 114,132Sn, and 208Pb. Enlarging the set of mass data does not modify the results presented in this paper. This fit provides optimal values for $\sigma_{0}$ and $b_{s}$ different for each set of uniform matter parameters $\vec{X}$. For illustration, Figure 1 displays the surface tension when $\sigma_{0}$ and $b_{s}$ are fixed for the parameter set $\vec{X}$ corresponding to the SLy4 interaction Chabanat1997 . Different values for the parameter $p$ are also displayed in that figure. We can see that the $p$ parameter determines the behavior of the surface tension for extreme isospin values Newton2013 , and cannot be accessed from empirical nuclear physics data, which are limited to values around $I\leq 0.3$. For this reason, the parameter $p$ will be added to the $\vec{X}$ set as an extra dimension in the parameter space of our meta-modelling, and a reasonable variation interval for its prior distribution will be determined in section IV.

III Results for a representative EoS

In the previous section we have seen that, even in a unified equation of state approach as the one employed in this paper, many different parameters have to be specified to calculate the crust composition and the crust-core phase transition. The importance of the equation of state empirical parameters like $E_{sym},L_{sym}$ has been pointed out by many authors in the context of specific EoS models Ducoin2011 and their influence will be studied in section IV. In addition to that, other parameters, much less constrained in the neutron star and neutron physics literature, enter explicitly in the variational equations eqs.(4)-(7). These parameters comprise the high order derivatives $Q_{q},Z_{q},\dots$ ($q=sat,sym$), the $b$ parameter that determines the limiting low density value where the Taylor expansion around saturation breaks down (see eq.(15)), and the $p$ parameter discussed in section II.3 governing the isovector behavior of the surface energy at extreme isospin values. In this section we will therefore consider a representative EoS model, and investigate whether the transition observables are modified by varying these relatively uncontrolled parameters. The representative EoS model is taken to be the parameter set $\vec{X}$ with empirical parameters fixed from the SLy4 functional, and will be noted as meta-SLy4 in the following.

III.1 Crust composition and CC transition point

The influence of the uncontrolled parameters on the crust composition is analyzed in Fig. 2. The upper part of the figure explores the importance of the high order terms in the density development. We can see that convergence is approximately reached at $N=3$, and a truncation at $N=2$ leads to an important underestimation of the transition density. This means that the $Q_{q}$ (and to a lesser degree $Z_{q}$) parameters cannot be neglected. $N=4$ will be used in the rest of the paper.

The medium panel in Fig. 2 explores the influence of the low density $b$ parameter. A very good reproduction of the SLy4 functional is obtained if $b=10\ln(2)$ is used. Taking different values of $b$ corresponds to considering EoS models which would have exactly the same empirical parameters (including effective masses and high order parameters up to $N=4$) as the SLy4 functional, but would differ in the treatment of the extreme low density domain. This is indeed what happens with the inclusion of deuteron or cluster correlations Roepke2009 ; Roepke2015 . Since again the effect is non-negligible, we keep $b$ as an extra EoS parameter. The low density correction to the Taylor expansion around saturation induced by the $b$ term in eq.(15) gets suppressed of a factor two at a density $n_{min}/n_{sat}=\ln 2/b$. We will consider in the following possible variations of $b/\ln 2$ in the interval $[1,10]$, which corresponds to a breaking down of the Taylor expansion at a density varying between $n_{min}=0.1n_{sat}$ and $n_{min}=n_{sat}$. This latter value corresponding to the minimum $b$ value might look quite extreme, but we will see in section V that the influence of this parameter turns out to be negligible. Finally, the effect of varying the isovector surface $p$ parameter is reported in the lower panel of Fig. 2. Here the variation is done around the value $p=3$ used in the popular Lattimer and Swesty EoS Lattimer1991 . We can see that both the average cluster size and the transition density are strongly influenced by the $p$ parameter, which cannot be directly linked to the uniform matter energy functional. From this observation we can already anticipate that the correlation between the CC transition and the EoS will be considerably blurred by our lack of knowledge of the isovector surface tension.

III.2 Comparison with previous calculations

Our results with the meta-functional optimized to give a good reproduction of the SLy4 model are compared with the other results available in the literature from a direct modelling of the inner crust using the SLy4 functional and different approximations for the surface tension Douchin2001 ; Vinas2017 . In the case of the popular Douchin and Haensel (DH) model Douchin2001 , we also report the average cluster size and charge for comparison with our results. We can see that the value of the cluster charge $Z$ is resonably compatible, but strong differences exist between the different calculations concerning the cluster size as well as the transition density. The difference in the cluster size between our approach and the DH one may be simply due to the different definitions of a cluster in a dense medium, which is affected by a certain degree of arbitrariness Newton2013 ; Panagiota2013 . Conversely, the difference in the transition point is most probably due to the slightly different treatments of the isovector surface tension. We can see that choosing $p=2.61$ ($p=2.34$) allows reproducing the previous results by DH (Vinas et al. Vinas2017 ), and the difference between Ref. Douchin2001 and Ref. Vinas2017 is of the same order as the difference between our results using $p=2.61$ or the canonical value $p=3$ from Ref. Lattimer1991 .

As we have already discussed in the introduction, a very limited number of works exists computing the CC transition from a direct modelling of the inner crust, because of the complexity of the simulations. A greater effort has been devoted to the calculations from the core side, using the thermodynamical or dynamical spinodal technique. The compilation of Ref. Ducoin2011 of the transition density and pressure calculated for different models is reported in Table 1 . For each model, the results compiled in Ref. Ducoin2011 are compared to our calculation where the parameter set $\vec{X}$ is fixed such as to reproduce the considered model. A part of these results can also be found in Ref. Thomas_prl .

We can see that the transition density results of the literature for the dynamical spinodal are globally nicely reproduced by our calculation ”from the crust” with the choice $p=3$. Higher values of the isovector surface tension parameter are needed if we want to reproduce the estimations of the thermodynamical spinodal. The transition pressure is more fluctuating, but still a variation of $p$ in the interval $[2.5,3.5]$ allows reproducing all the considered models. A flat probability distribution within this interval will be our prior for the statistical analysis of section V.

IV Sensitivity analysis

One of the advantages of the meta-modelling technique is that, since all the model parameters are a-priori uncorrelated, it is possible to vary each one of them independently of the others, which is not possible using specific functional behaviors such as Skyrme or Gogny or the different versions of RMF. This allows determining the most influential parameters on any given observable. Such a sensitivity analysis is presented in Fig. 4 for the transition density and in Fig. 5 for the transition pressure.

The one-by-one variation of all the EoS parameters is performed around two different reference parameter set $\vec{X}_{ref}$, namely the parameter set corresponding to the SLy4 model (upper part of Figs. 4,5), and the set of average values of the different parameters from the compilation of empirical constraints in Ref. Margueron2018a (lower part of Figs. 4,5).

The minimum and maximum value chosen for each parameter are taken from Ref. Margueron2018b and they are given in Table 2. These values reflect the degree of uncertainty on the different parameters, as measured by their observed variation in the different functionals that have been successfully confronted to low energy nuclear physics data. The vertical lines in Figs. 4,5) give the transition density and pressure domain obtained when the EoS parameters are one by one varied around the reference model, within the interval of Table 2. Since the uncertainty on the different parameters is not the same, the length of the segments is a qualitative measure of the propagation of the uncertainty on the transition point brought by each parameter.

We can see that the sensitivity of each parameter depends on the value of the other parameters, that is on the chosen reference set $\vec{X}_{ref}$. Still, universal trends clearly emerge. We can see that the CC phase transition, at variance with the standard liquid-gas of symmetric matter, is virtually insensitive to isoscalar parameters. Even if extremely large variations of $Q_{sat}$ and $Z_{sat}$ are considered (see Table 2), the prediction of the transition point is almost unaffected. This underlines the importance of the energetics of the neutron gas on the transition point. Concerning the isovector sector, we can see that $L_{sym}$ is the most important parameter. This result is in agreement with previous findings by many authors Ducoin2011 . The symmetry energy at saturation $E_{sym}$ and the effective mass splitting $\kappa_{v}$ do not play any role on the transition, which can be partially explained by the fact that these parameters are already relatively well constrained. Depending on the chosen reference point, the transition pressure shows also a great sensitivity to the isovector compressibility $K_{sym}$. This can explain why the transition pressure exhibits an irregular behavior when plotted as a function of $L_{sym}$ Ducoin2011 : the different functionals considered in the litterature have very different values of $K_{sym}$, which blurs the correlation with $L_{sym}$. This effect is also amplified by the fact that, depending on the reference point, the dependence of $P_{t}$ with $L_{sym}$ is not monotonic. Finally, we can remark that the influence of the fully unknown high order derivatives $Q_{sym}$ and $Z_{sym}$, though less important than the one of $L_{sym}$, is not negligible and comparable to the one of the isovector surface energy parameter $p$ that can be inferred from Table 1. Similar conclusions can be drawn if the sensitivity analysis is performed using the definition of the transition point from the dynamical spinodal Chatterjee ; Antic .

V Statistical analysis

We now come to the quantitative determination of the CC transition point and its uncertainty, on the basis of our imperfect knowledge of the nuclear energy functional. To this aim, we perform a Bayesian determination of the model parameters on the full $2(N+1)+3$ parameter space, up to the fourth order in the Taylor expansion ($N=4$). The prior distribution of $\vec{X}$ is given by an uncorrelated ansatz and a flat distribution of each parameter within the interval specified in Table 2,

[TABLE]

where $f$ is a uniform distribution between $X_{k}^{min}$ and $X_{k}^{max}$, defined in terms of the Heavyside step function $H(x)$ as $f(a,b;x)=(H(x-a)-H(x-b))/(b-a)$. The posterior distribution is obtained by applying different physical filters to the prior distribution,

[TABLE]

In this expression, both strict ($w$ term) and likelyhood (exponential term) filters are applied, and $\mathpzc{N}$ is a normalization. $\chi^{2}(\vec{X})$ represents the $\chi^{2}$ corresponding to the optimal fit of nuclear masses which is done to determine the surface tension parameters $\sigma_{0}$ and $b_{s}$ for each $\vec{X}$ parameter set (see section II.3). $w$ is a sharp $\delta$-function filter that outputs 1 if the constraint is respected, and 0 otherwise.

Two different constraints are considered for the $w$ filter. The first constraint, noted $w_{LD}$, concentrates on the low density (LD) $n\leq n_{sat}$ behavior of the energy functional. We impose to the different functionals generated following eq.(20) to strictly pass through the uncertainty band of the N3LO effective field theory calculation for symmetric and pure neutron matter by Drischler et al. Drischler2016 . This same condition was applied in previous studies Thomas_prl ; Chatterjee ; Antic ; Holt . This filter is applied in the density interval $[0.05,0.2]$ fm*-3*. The very low density region is not considered because of numerical issues due to the very small uncertainty, but we recall that all the generated models by construction converge to zero energy and pressure in the $n\to 0$ limit. The energy per particle and pressure uncertainty bands applied are displayed in Fig. 6. In this figure, the interval delimited by the dotted lines corresponds to our prior distribution. We can see that the uncertainty band of the prior energy per particle and pressure of symmetric matter around saturation is comparable and even narrower than the one corresponding to the ab-initio calculation. This is due to the strong empirical constraints coming from different low energy nuclear physics experiments, that have been considered to determine realistic intervals for the empirical parameters in Table 2. Conversely, as it is well known, the empirical information embedded in the prior distribution is insufficient to effectively constrain the neutron matter EoS, and the ab-initio predictions are much narrower than our prior distribution. In the isovector sector the LD filter is extremely selective: only 2118 models out of the 100 millions generated to numerically sample the prior parameter distribution, fulfill the LD condition.

The second (HD) filter imposes general physical constraints to the global density behavior of the functional, as follows:

• •

positive symmetry energy at all densities,

• •

$P(n)\geq 0$ for $n\geq n_{sat}$ (stability of the EoS),

• •

$0<v_{s}<c$,

• •

$M_{max}>2M_{\odot}$,

where $M_{max}$ is the maximum neutron star mass obtained through a TOV calculation Haensel_book .

Once the posterior parameter distribution is determined, the probability distribution of any parameter or observable $Y$ can be straightforwardly computed as:

[TABLE]

where $Y(\vec{X})$ is the value of the Y variable as obtained with the $\vec{X}$ parameter set.

V.1 Distribution of parameters and observables

The posterior distribution of the most influential isovector EoS parameters is displayed in Figure 7. We can see that the general physical conditions corresponding to the HD filter almost do not constrain the low order empirical parameters, in agreement with the findings of Ref. Margueron2018b . Conversely, the LD filter allows a very tight determination of the empirical parameters $L_{sym}$ and $K_{sym}$.

The effect on the distribution of the transition density and pressure is displayed in Fig. 8. The two-humped posterior distribution is due to the isovector surface energy parameter $p$, for which we consider for simplicity only three values 2.5, 3, and 3.5 with equal probability: the left peak in the distribution of $n_{t}(P_{t})$ is associated to $p=2.5$ while we can not distinguish between $p=3$ and $p=3.5$. This shows that the highest uncertainty in the determination of the transition point, once the EoS is constrained through the most advanced ab-initio calculations of nuclear matter, is due to our poor knowledge of the surface properies of extremely neutron rich matter.

This statement can be quantified by calculating the first moments of the distribution. We give in Table 3 the average value and the standard deviation of the transition density (pressure) $n_{t}(P_{t})$ for the prior and posterior distributions. In each case, the convergence of the results with the number of sampled models fulfilling the most restrictive LD filter, is also checked. We can see that, if the value of the isovector surface tension parameter $p$ is fixed, the LD filter is much more effective in reducing the uncertainty on the transition point, with respect to the HD filter. This underlines the importance of precise constraints on the low density EoS parameters for a reliable prediction of neutron star crust properties, as it has been often stressed in the literature. If the transition point is computed with reasonably well-behaved models for neutron stars fulfilling the HD conditions and with empirical parameters within the accepted bands from nuclear physics experiments (see Table 2), without the more precise constraints from ab-initio calculations, the transition pressure can be badly predicted even in average.

This discussion ignores the uncertainty that we have on the behavior of the surface tension for extreme isospin values. If we incorporate that uncertainty considering $p$ as an extra parameter (last line in Table 3), we can see that, according to the range assumed for the prior distribution of $p$, the uncertainty on the transition point is considerably increased and even the average value is affected.

V.2 Correlations among the empirical parameters and the transition

We now turn to explore the correlations between the empirical EoS parameters and the transition point, as obtained by applying the different filters. Fig. 9 displays the correlation matrix among the different EoS parameters. The isovector surface tension parameter $p$ is fixed to $p=3$ for this study but the results are unmodified if $p$ is allowed to vary. This is because $p$ is decoupled by construction from the homogeneous EoS parameters, and it additionally does not play any role in the mass fit (see section II.3), meaning that it is independent also of the other parameters of the surface tension. The matrix elements above the diagonal in Fig. 9 give the linear correlation coefficient $r_{ij}=|\sigma_{X_{i}X_{j}}|/\sigma_{X_{i}}\sigma_{X_{j}}$ obtained for the posterior distribution eq.(21) with $w=w_{LD}$, that is after application of the low density filter, while the terms below the diagonal refer to the high density filter. In agreement with the findings of Ref. Margueron2018b we can see that the HD filter does not induce any correlation among the empirical parameters, with the exception of a small positive correlation between $L_{sym}$ and $K_{sym}$, essentially due to the EoS stability requirement at high density. On the other side, many different correlations appear due to the LD constraints. The constraint of mass reproduction induces a clear correlation of the surface tension parameters between themselves, as well as with the zero order isoscalar and isovector parameters ($E_{sat},E_{sym}$), which dominate the global energetics of finite nuclei. More interesting correlations among the different isovector parameters are induced by the constraint of reproducing the ab-initio EFT calculations: besides the well known correlation between $E_{sym}$ and $L_{sym}$ which has been observed by many authors in the context of different models  Kortelainen2012 ; Danielewicz2014 ; Trippa2008 ; Colo2014 ; Holt , we also observe a strong correlation among higher order parameters, notably the isovector curvature $K_{sym}$ with the skewness $Q_{sym}$, and the high order isoscalar parameters $Q_{sat},K_{sat}$ with the corresponding isovector ones $Q_{sym},K_{sym}$. These non-trivial correlations can only be observed within the meta-modelling strategy, because in popular functionals like Skyrme the high order parameters are a-priori correlated by the chosen functional form. The fourth order parameters $Z_{sat},Z_{sym}$ do not show any correlation with any other parameter, showing their negligble influence on the density relatively close to saturation implied in the LD filter. Finally, the isoscalar effective mass and effective mass splitting are also essentially uncorrelated with the others: variations of these parameters, which play a crucial role in the structure of finite nuclei, are fully compensated by variations of the density derivatives as long as only the total energetics (kinetic plus potential) is involved Chatterjee2017 .

The correlation matrix between the model parameters and the density and pressure of the transition point is presented in Fig. 10. When the prior parameter distribution is used, which supposes the EoS parameters fully uncorrelated, the transition density is only (slightly) negatively correlated to the $L_{sym}$ parameter, as it was previously reported Ducoin2011 . The transition density directly depends on the energy of $\beta$- equilibrium matter. The transition pressure being linked to the first derivative of the energy density, it is not surprising that it is correlated to higher order parameters of the symmetry energy, namely $K_{sym}$ and $Q_{sym}$. Since these parameters widely vary in existing functionals, this can explain why the present predictions of the transition pressure are so largely scattered (see Table 1). Once physical correlations among the EoS parameters are accounted for in the posterior distribution, new correlations appear for the transition point.

As we have already discussed, for the computation of the transition point it is important to specify, together with the different parameters analyzed in Fig. 9, also the isovector surface tension parameter $p$. Let us first consider the case where we fix this parameter to its canonical value $p=3$ that best reproduces the homogeneous matter dynamical spinodal (see Table 1). In this case, we can see that the correlation of $n_{t}$ (resp. $P_{t}$) with $L_{sym}$ (resp. $K_{sym})$ is preserved. Further interesting correlations with the high order isovector parameters $K_{sym}$ and $Q_{sym}$ emerge if the EoS sample is restricted to respect the HD filter, and even more if the LD filter is applied. All these correlations fade away if the $p$ parameter is allowed to vary (lines ”LD+masses”). In that case, the transition point is solely correlated to the surface properties (parameters $\sigma_{0}$, $b_{s}$ and $p$). This means that the dominant parameter determining the CC transition is the isovector surface tension. Only if we constrain its behavior ar large isospin values imposing $p=3$, the correlation with the symmetry energy is recovered. In that case, we can see that the knowledge of the largely discussed $L_{sym}$ parameters is not enough, and higher order parameters beyond $L_{sym}$ must be constrained to improve the prediction of the transition point, both in density and in pressure.

Our final result for the value and $1\sigma$ confidence interval for the transition point is displayed in Fig. 11. We can observe in a graphical way that the LD and HD filter lead to compatible predictions for the transition point, but the LD filter is by far more constraining. It is clear from this figure that a reliable determination of the CC transition point demands a better control on the isovector properties of the nuclear surface for extreme isospin values, more than more stringent constraints on the EoS parameters. A Bayesian analysis of the microscopic evaluation of the isovector surface energy within the extended Thomas-Fermi method is under progress.

VI Conclusions

In this paper we have presented a detailed study of the core-crust transition point, in the framework of a unified EoS treatment where the inhomogeneous crust is calculated with the same energy functional as employed for the modelling of the homogeneous core. The full parameter space of the EoS, including the successive isoscalar and isovector derivatives of the energy functional at saturation up to $N=4$, the isoscalar and isovector effective mass, and the isospin dependent surface tension, is evenly explored within a completely uncorrelated flat prior, within intervals compatible with the present empirical constraints. The correlated $(2(N+1)+3)$-dimensional parameter distribution is obtained constraining the parameter space such as to reproduce measured mass of magic and semi-magic nuclei, fulfill basic physical conditions as well as modern ab-initio calculations of symmetric and neutron matter. We find that the most influential parameter for the determination of the transition point is linked to the behavior of the surface tension for extreme isospin values, specifically the isospin value at which the surface tension vanishes. Only if this parameter is fixed to an educated, but somewhat arbitrary value, important correlations with the isovector parameters of the EoS are recovered. We confirm the correlation of the transition density with the slope of the symmetry energy at saturation $L_{sym}$ already observed in previous works, and additionally point out the important correlation of the transition pressure with the isovector compressibility $K_{sym}$.

Ab-initio calculations of the symmetric and neutron matter energy and pressure for densities below saturation from Ref. Drischler2016 provide very stringent constraints on the EoS parameter distribution, and are shown to be much more effective for the determination of the transition point than the available constraints at supersaturation density.

Still, if the isovector surface tension is not further constrained, considerable uncertainties affect the transition point ($n_{t}=0.060\pm 0.027$ fm*-3*, $P_{t}=0.25\pm 0.225$ MeV/fm*-3* at the 1$\sigma$ level). If the educated guess $p=3$ can be confirmed by microscopic calculations, the predictions are sensibly changed to $n_{t}=0.072\pm 0.011$ fm*-3*, $P_{t}=0.339\pm 0.115$ MeV/fm*-3*. These values, with the associated uncertainty intervals, determine the crustal width and momentum of inertia of neutron star, which in turn allows quantifying the role of the NS crust in interesting astrophysical phenomena such as pulsar glitches. Our study is performed in the framework of the compressible liquid drop (CLD) model for the inhomogeneous crustal matter and a specific parametrized functional form for the surface energy, but equivalent parameters to our $p$ parameter can be found in alternative modellings, such as for instance isovector gradient couplings in the DFT or Thomas-Fermi approximation Bender2011 , or the surface stiffness parameter in the context of the Droplet model Warda2009 ; Centelles2010 ; Mondal2016 .

For this reason, further studies of the isovector surface tension will be extremely important.

Acknowledgements.

This work was partially supported by the IN2P3 Master Project MAC, ”NewCompStar” COST Action MP1304, PHAROS COST Action MP16214.

Appendix A More on the properties of EFT EoS

The EFT constraints displayed in Fig. 6 are obtained from Ref. Drischler2016 as the limiting surface containing the seven different functionals proposed, which correspond to different hypotheses for the NN interaction, resolution scale, cut-offs, and long range couplings. It can also be interesting to compute the transition point using directly the functionals of Ref. Drischler2016 . This can be done within our meta-modelling technique considering seven $\vec{X}$ parameter sets corresponding to the different functionals, which allows computing the energy density for different isospin ratios. The results for the beta-equilibrium EoS and the transition point are displayed in Table 4 and Fig. 12.

For some functionals, no solution can be found for the variational equations in the whole density range covered by the crust: this is in particular the case for model 6, which predicts a very low value of $E_{sat}$ not compatible with values extracted from empirical information from nuclear mass. Apart from these slight anomalies, the results are well within the range of our predictions imposing the global constraint of Fig. 6.