From neutron skins and neutron matter to the neutron star crust

We present the first Bayesian inference of neutron star crust properties to incorporate neutron skin data, including the recent PREX measurement of the neutron skin of $^{208}$Pb, combined with recent chiral effective field theory predictions of pure neutron matter with statistical errors. Using a compressible liquid drop model with an extended Skyrme energy-density functional, we obtain the most stringent constraints to date on the transition pressure $P_{\rm cc}=0.33^{+0.07}_{-0.07}$ MeV fm$^{-3}$ and chemical potential $\mu_{\rm cc}=12.6^{+1.8}_{-1.9}$ MeV (which control the mass, moment of inertia and thickness of a neutron star crust), the proton fractions that bracket the pasta phases $y_{\rm p}=0.115^{+0.016}_{-0.017}$ and $y_{\rm cc}=0.041^{+0.007}_{-0.006}$, as well as the relative mass and moment of inertia $\Delta M_{\rm p} / \Delta M_{\rm c}\approx \Delta I_{\rm p} / \Delta I_{\rm c} = 0.54^{+0.05}_{-0.09}$ and thickness $\Delta R_{\rm p} / \Delta R_{\rm c}=0.129^{+0.019}_{-0.030}$ of the layers of non-spherical nuclei (nuclear pasta) in the crust.


Introduction
In the deepest layers of a neutron star crust, nuclear sizes and separations become comparable. Coulomb and surface energy competition predict that nuclei deform and fuse through a sequence of shapes termed nuclear pasta that cushion the base of the crust [1,2,3]. The extent of the crust and pasta phases therein could have observable consequences [4], affecting crustal oscillation modes [5,6,7], crust cooling [8,9,10], crust shattering [11,12], persistent gravitational waves from mountains [6], magnetic field evolution [13] and damping of core modes [14,15]. In order to make progress understanding these phenomena, it is important to constrain the size of the crust and the amount of pasta therein.
The compressible liquid drop model (CLDM) [16,17,18,19] is an efficient method of calculating crustal composition, mass and thickness. A Wigner-Seitz cell containing one nucleus or segment of nuclear pasta surrounded by the neutron gas is modeled. Theoretical uncertainty arises on two fronts: the pure neutron matter (PNM) equation of state (EOS) that governs the hydrostatic pressure and chemical potential in the deep layers of the crust, and the interface energy of the nuclei and pasta shapes, which is specified separately to the bulk EOS [6,20,21,22,23,24,25].
In this letter, we perform, for the first time using neutron skin data, a Bayesian inference of the crust-core transition pressure P cc which controls the mass of the crust, crust-core transition baryon chemical potential µ cc which controls the thickness of the crust, the thickness of the pasta phases relative to the crust ∆R p /∆R c , the relative mass and moment of inertia of the pasta phases ∆M p /∆M c ≈ ∆I p /∆I c , and the proton fraction at the top (y p ) and bottom (y cc ) of the pasta phases. We combine this with theoretical PNM calculations [26] with well quantified errors. Neutron skin measurements and PNM calculations are known to inform the location of the crust-core boundary and the amount of nuclear pasta [27,28,29,21,23]. We account for all uncertainties in the CLDM model used, and detail the remaining sources of uncertainty still to be addressed.

Model
Neutron skins, PNM, and bulk nuclear matter in the CLDM are modeled using the extended Skyrme energy-density functional (EDF) [30,31] which allows us to explore a wide range of density dependences of the EOS in the range important for nuclear pasta, the crust-core transition and finite nuclei: 0.25 -1 n 0 where n 0 = 0.16 fm −3 is nuclear saturation density.
The full expression for the extended Skyrme energy density functional (EDF) used here can be found in [32]. The parts of the Skyrme EDF H that we use to control the symmetry energy are the zero-range and density-dependent terms where y p and y n are the proton and neutron fractions and α 3 , α 4 , x 3 , x 4 , t 3 and t 4 are parameters. x 0 , x 3 and x 4 can be adjusted to control the density dependence of the PNM EOS, while the remaining parameters can be readjusted to maintain symmetric nuclear matter properties at the values of the baseline Skyrme Skχ450 [31] within the accepted values from nuclear experiment. The interface tension, consisting of surface σ s (y p ) and curvature σ c (y p ) terms, is given as a function of proton fraction y p [33,34,35]: where the following four model parameters are identified: σ 0 and σ 0,c control the strength of the surface and curvature tension in symmetric matter y p = 0.5 and p and b control the isospin dependence. In particular, the parameter p controls the behavior of the surface energy in very neutron rich environments y p → 0 [36,20], and is particularly important in the modeling of the crust. The curvature tension has been extended in some works to include its own proton fraction dependence [36,20] with additional model parameters, which we do not use here.

Model parameters
The nuclear symmetry energy S (n), a function encoding the energy cost of replacing protons with neutrons in a nuclear system, is a convenient intermediary between nuclear and neutron star properties. It is defined in the expansion of the energy per particle of uniform nuclear matter about a proton fraction of one half. Similarly one can define an expression for the surface symmetry tension σ δ , so we have: where δ = 1 − 2y p . The symmetry energy can then be expanded as a power series in density about nuclear saturation density n 0 = 0.16 fm −3 : using the parameter χ = (n − n 0 )/3n 0 , The strong correlations between the magnitude, slope and curvature of the symmetry energy J, L and K sym and neutron star properties such as their radii and proton fraction, have led to a sustained experimental effort to infer them from nuclear observables [37,38,39,40,41,42].
J, L and K sym can be written in terms of the parameters x 0 , x 3 and x 4 in the extended Skyrme model [32] allowing us to explore the symmetry energy parameter space by translating ranges of J, L and K sym into Skyrme models, keeping all symmetric nuclear matter parameters fixed.
To relate the parameters of our surface energy model Eqn 3 to its isospin expansion Eqn 4 we note that the surface symmetry tension is related to the model parameters via [36] Relating the model parameters to the surface symmetry energy allows us to include a generic correlation between the surface and bulk symmetry energies that emerges from nuclear mass fits and semi-infinite nuclear matter (SINM) calculations [43,44] which can be written [36] Prior is the uniform prior we start with. The regions for which viable crust models exist are shown as the hatched region, and the region consistent with unitary gas constraints in red. The uniform and unitary gas constraint coincide in L-K sym space since the latter does not constrain K sym [29]. The 1 − σ error ellipse from chiral-EFT calculations [45] is shown in grey (a).
c controls the slope of the correlation [44,36] which depends on the exact methods used to extract the surface energy of which there are a variety [43,44]. We take the parameter c to vary widely over a range that encloses all relations extracted from empirical and theoretical fits.
We thus replace the model parameter b in Eqn(2) with the parameter c defined above; given a value of c, b can then be calculated combining Eqns (5) and (6).
It is important to note that we used the full bulk energy density functional (the extended Skyrme) and surface energy functional (Eqn. 3) in our calculations, and not the truncated expansions 4-5 which provide us convenient model parameters.
To summarize, we have seven model parameters: J, L and K sym controlling the bulk EOS and σ 0 , σ 0 , c and p controlling the surface energy.

Priors
Priors are given in Fig. 1. We start with uniform priors on both the symmetry energy and surface parameters, over a wide enough range to incorporate all experimental constraints. Not all of the parameter space gives viable crust models: in particular, the large L, small J EOSs lead to PNM instabilities at crust densities. Although these models are effectively filtered out, we still refer to the remaining hatched region of Fig. 1 as our uniform prior and label it Sym:Unif. We can additionally include theoretical constraints from the unitary gas limit of PNM at very low densities [47,29,48,49] leaving the red region in Fig. 1 (Sym: Unif+Unit prior). Note that in L − K sym space the uniform and unitary gas priors coincide.
We refer to the uniform surface priors with the label Surf:Unif. These cover a range within which determinations (r 208 np ) (fm) Figure 2: A histogram (a) of the differences ∆(r 208 np ) between the neutron skins of 208 Pb predicted by the 91 of the most recent Skyrmes from [46], and extended Skyrmes that have been refit to each value of J,L and K sym from that set. We show the values of J versus L (b) and L versus K sym (c) common to two sets of Skyrme parameters: The colors indicate ∆(r 208 np ). There is no obvious correlation between each symmetry energy parameter and the difference between the neutron skin predictions, so we may treat ∆(r 208 np ) as a random variable. The distribution in (a) can be mimicked by a Gaussian of standard deviation 0.0045fm. of surface parameters from mass fits and SINM calculations reside. As a more restrictive set of surface priors we use a fit of the surface parameters to three-dimensional Hartree-Fock (3DHF) calculations of crust nuclei near the onset of nuclear pasta n ∼ 0.25n 0 (Surf:Fit) [21]. Such calculations are computationally intensive and we can only perform a limited number of them, and so it is difficult to obtain robust model uncertainties on the best fit values of most of the parameters: the fits are most sensitive to the parameter p and that is the only parameter we can currently report an error on; we fix the other parameters at their best fit value. The main purpose of these priors is to provide a way to quantify the effect of our lack of knowledge about the surface parameters.

Modeling neutron skins
We use ensembles of extended Skyrme energy-density functionals (EDFs) in 1D Hartree-Fock calculations of neutron skins using the Skyrme RPA code [56] to infer posteriors of J,L and K sym from the neutron skin and PNM data; the Bayesian inference is detailed in full in [57].
Our ensembles of Skyrme parameterizations are created by varying x 0 ,x 3 , and x 4 while holding the other parameters constant. However, this neglects correlations between the symmetry energy parameters and other nuclear matter parameters that would arise if the whole functional was fit to, for example, mass data for each choice of x 0 ,x 3 , and x 4 , something that is beyond the scope of this study but is the subject of ongoing work. These other nuclear matter parameters -for example the isovector gradient coefficient G v -might also influence the neutron skin. There is evidence that, in the case of neutron skins, the error this introduces is small; for example, the isoscalar and isovector gradient coefficients show little correlation with the neutron skins [58,59].
In order to assess this source of modeling uncertainty, we calculate r 208 np for 91 different Skyrme interactions created in the past 20 years taken from [46], fit to various subsets of nuclear data. We then construct 91 extended Skyrme interactions with x 0 ,x 3 , and x 4 adjusted to produce the value of J, L and K sym of each of the 91 original Skyrmes. We calculate the difference in r 208 np for each pair of Skyrmes -the original and its corresponding extended Skyrme counterpart -and plot the results in Fig 2. Figs 2a and 2b show the distribution of symmetry energy parameters for the models, with the differences in r 208 np between the original Skyrmes and the equivalent extended Skyrmes indicated by the color, from lower differences (blue) to higher (yellow). There is no evidence of correlation between the difference in r 208 np and the values of J, L and K sym , so we model the errors as random. In Fig 1c we show the distribution of ∆(r 208 np ). The distribution can be reasonable modeled by a Gaussian of standard deviation 0.0045 fm. We thus account for this uncertainty by adding an extra ±0.0045fm error to our neutron skin data points.

Results.
We first report on the symmetry energy constraints from the PREX measurement of r 208 np . The 68% credible intervals  L parameter inferred is high ≈ 50 − 114MeV, but about 20MeV smaller and with slightly smaller uncertainty than the value inferred from [61] from a more restricted range of models. The inferred value of J is consistent with nuclear mass fits [60]. We sample ∼100,000 points from our posteriors for each prior/dataset combination and calculate the resulting crust models. In Fig. 3 we show the effect of the surface priors on the marginalized posterior distributions of ∆M p /∆M c (Fig. 3a) and ∆R p /∆R c (Fig. 3b), both calculated using the approximations outlined in Appendix A. We show the priors (blue) and the posteriors obtained from a representative dataset that gives significant constraints: PNM (orange). In each plot, from right to left, we show: our most restrictive surface prior, the fit to 3DHF calculations (the Surf:Fit prior); the same except for allowing c to assume its full range 2.0-7.0 (Fit+c Full), the same but allowing p to vary over its full range 2.0-4.0 (Fit+p Full), and finally all surface parameters allowed to vary over their full ranges (the Surf:Unif prior).
The change in posterior from changing the prior from Surf:Unif to Surf:Fit is almost entirely due to relaxing the allowed range for p: the remaining surface parameters have much smaller effects. This is not surprising, as it is the parameter p which controls the isospin dependence of the interface energy as the proton fraction approaches zero (see e.g.  To account for differences not captured by the median and CIs, we also calculate a difference index [62] defined as d(A, B) = 1 − min[ρ A (x)ρ B (x)]dx where x is the variable of interest, and ρ A,B is the distribution of interest and the reference distribution (in this case the Sym:Unif; Surf:Unif prior) respectively. The difference indices are indicated by the colors of the points in Fig. 5.
The inferred values of P cc , µ cc , ∆M p /∆M c , ∆R p /∆R c , y cc and y p for each dataset are given in Table 1. The values inferred with the fewest initial modeling assumptions -Sym:Unif;Surf:Unif priors -are given, with results from Sym:Unif; Surf:Fit priors in parentheses.
Finally, in Fig. 6 we show the two-dimensional posteriors for the Sym:Unif; Surf:Unif priors. We show results for the prior, the PREX posterior and the Posterior with all our neutron skins data. The 1-σ credible regions for L, K sym , r 208 np , P cc , µ cc and ∆M p /∆M c and ∆R p /∆R c are displayed.
Comparing the prior (circular markers in Fig. 5), PREX (square markers) and Skins (diamond markers) datasets, the effect of the PNM-free data on the medians is to elevate them compare to the prior, except for the case of PREX data applied to ∆M p /∆M c . The effect is relatively small, however. Then, adding the PNM data favors systematically smaller crusts, less pasta, and a smaller proton fraction throughout the pasta phases: medians with the PNM data included are about 25-30% lower than without PNM data for all quantities except for ∆M p /∆M c which drops by ≈ 10%.
As a result of the correlations with the neutron skin of 208 Pb (see, e.g. Fig. 24 and the top panel of Fig. 27 in [21]), tails in the posterior distributions develop at high values of P cc and low values of ∆M p /∆M c and ∆R p /∆R c , which leads to the PREX data increasing the width of the CIs of P cc , ∆M p /∆M c and ∆R p /∆R c . Indeed, there is a smaller peak at zero pasta which can be understood as follows: as the pressures and chemical potentials of the pasta and crust-core transition approach each other, the amount of pasta shrinks. Then all models that predict that the crustcore transition occurs first will predict zero pasta. The bars at the bottom of Fig. 4 show the fraction of models that predict no pasta in the crust for each dataset. The PREX data alone, the PNM data alone, and both combined give a significant strength for zero pasta ∼ 10%. The rest of the neutron skin data strongly favors pasta in the crust. P cc , µ cc ,y cc and y p become significantly more constrained when we include the PNM data (the up, down and left triangles in Fig. 3), as can be seen by inspection in Fig. 4 and quantified in Fig. 5: the 68% CIs shrink by 30-60%. The PNM data do not systematically shrink the 68% CIs for ∆M p /∆M c and ∆R p /∆R c ; although PNM data tightly constrains on J and L, the correlation between J and L and the pasta parameters is relatively weak (see Fig. 22 and top panel of Fig. 27 in [21] or Fig. 7 of [23]).
As we add full neutron skin data (Prior→Skins, PNM→PNM+Skins) the 68% CIs shrink by around 30-50% for P cc -so the neutron skin data is comparably informative to the PNM data for P cc (therefore crust mass). The effect on the 68% CIs for µ cc , y cc and y p is around 10%: PNM data is significantly more constraining for µ cc (therefore crust thickness), y cc and y p . One can see in the top panel of Fig. 27 in [21] that, starting from uniform priors, the quantities most correlated with P cc are indeed the neutron skins and L (the latter constrained by PNM data). We see this correlation in action in the two dimensional posteriors of Fig. 6: the full neutron skin dataset removes only a little strength at low values of µ cc which shows up in Fig. 4, whereas it removes more strength at both the low and high ends of the distributions of P cc , ∆M p /∆M c and ∆R p /∆R c resulting in tighter constraints. In Fig. 4, the shapes of the Surf:Unif and Surf:Fit priors can be seen to be broadly similar for all datasets, and in Fig. 5 markers of the same shape (data) cluster more than those with the same border (prior), and that changes in the posterior's shape is driven by changing data (compare colors between priors and between data). Note that in Fig. 5 we also include the points corresponding to the unitary gas prior Sym:Unif+Unit which makes only a small difference to the median values. As might be expected, the Surf:Fit prior and the Sym:Unif+Unit prior, which incorporate more information into the priors, have the effect of shrinking the 68% CI. For P cc , µ cc , y cc and y p , the effect of the Surf:Fit prior compared to the Surf:Unif prior is no more than that of the data. Additionally, it is clear that the main trend is for greater difference index (color in Fig. 5) to correlate with smaller CI, and that the data gives the dominant effect.
The exception to this is ∆M p /∆M c and ∆R p /∆R c in the case when we include PNM and all skins data; then Surf:Fit priors Table 1: Table of median and 68% (95%) CIs for crust-core transition pressure, chemical potential,mass/moment of inertia,thickness fractions of pasta, crust-core proton fraction and proton fraction at the onset of pasta are shown using uniform surface priors (Surf:Fit in parentheses) and the 5 different data combinations. result in smaller CIs than Surf:Unif priors by over 50%. In this case, the priors on the surface energy parameters are more informative than the data. This highlights the importance of finding more robust constraints on the surface energy (noted also in [23]). Additionally, the difference index does not appear to be correlated particularly with the change in CI for these two quantities: the data and priors mainly change the shape of the distribution characterized by its higher moments -for example, the development of a second peak at zero pasta.

Summary and Conclusions.
The combination of PNM and neutron skin data gives the tightest constraints, and are most robust when uniform priors are used. A significant amount of pasta -50% by mass and moment of inertia, 12% by thickness -is strongly favored, although the possibility that no pasta exists in the crust cannot be ruled out; 12% of models predict pasta is absent for the PREX+PNM dataset. The median proton fraction in the pasta phases is predicted to decrease from ≈16% at the top of the pasta layer to ≈7% at the crust-core transition when only neutron skin data is taken into account, and from ≈12% down to ≈4.5% when PNM data is added. This is of relevance to, for example, possible direct Urca processes in the pasta layers that may enhance the cooling of the neutron star [10]. We note here that extend-ing the chiral-EFT analysis of [26] to infer K sym should lead to even tighter constraints on the crust.
Although we have attempted a thorough characterization of the uncertainties within our CLDM, there are certainly additional uncertainties we should work towards quantifying. Although our extended EDF allows a wide parameter space exploration of the symmetry energy, the third order symmetry energy parameter, Q sym , also correlates with the crust parameters [20,23]. The surface energy function is not unique, and extensions and alternatives have been studied (see e.g. [63,35,36,20]). More microscopic alternatives to the CLDM such as Thomas-Fermi [64] and Hartree-Fock [65] methods give results that appear consistent with the CLDM, but a more thorough comparison should be done.
Nevertheless, there is a concordance between our results and similar recent statistical predictions, as can be seen by the black crosses in Fig 4. Constraints obtained using a similar CLDM [20,23] constrained by low density PNM calculations are shown. For P cc , a range of p =2.5-3.5 ('C PNM') and a single value of p = 3 ('C PNM p3') were used, so these are most comparable with Surf:Unif and Surf:Fit priors respectively. For ∆R p /∆R c and ∆M p /∆M c the same model with surface parameters fit to nuclear masses [23] ('DT PNM') was used. The above studies used a more schematic EDF that allowed for Q sym to be varied too. Also shown it the constraint on P cc from [60] ('LL PNM') obtained by examining the stability of nuclear matter with respect to small density fluctuations, and was also constrained by PNM calculations (but imposes a correlation between J, L and K sym ). In all cases, the best comparison with our calculations is with the 'PNM' dataset (triangles). Let us also note the classic work of [66] predicted a relative mass and moment of inertia of pasta of around 0.5. Previous studies have not included neutron skin data. The thermodynamic spinodal method of [60] is known to give a upper bound on the transition point. Given this, the models are in remarkable agreement. We can conclude that we are entering an era where statistically meaningful constraints on neutron star crusts from experimental and theoretical data are becoming competitive with modeling uncertainties.
Our results point to there being a large amount of pasta in the neutron star crust, which has a number of notable observational consequences. If disordered, as microscopic simulations of pasta suggest [65], this predicts that the cooling timescale of the neutron star crust will be significantly enhanced at late times. It also suggests that magnetic fields in the crust would decay rapidly in pulsars with ages ∼ 10 5 − 10 6 years [13]. We should expect the shear modulus at the crust-core interface to be modified by the existence of pasta, which has implications for a number of scenarios. If the shear modulus is smaller, for example: pasta could less efficiently damp r-modes in rotating neutron stars, therefore allowing greater spin-up of millisecond pulsars [14,15]; pasta would decrease the resonant frequency of crust i-modes possibly observable in resonant shattering of the crust just prior to two neutron stars merging [12]; the maximum size of mountains, and hence their detectability as persistent gravitational wave sources, would be reduced [6]. A large amount of pasta increases the direct Urca emissivity of the star [67] resulting in enhanced cooling early on and subsequently decreasing the rapid cooling rate brought on by the onset of superfluidity in young neutron stars [10]. The effects of pasta on the dynamics of neutron vortices in the inner crust, a key ingredient of many glitch models, have yet to be modeled. Therefore our work suggests renewed efforts to understand the microphysics of nuclear pasta are required.

Acknowledgements
This work is supported in part by the Physics and Astronomy Scholarship for Success (PASS) project funded by the NSF under grant No. 1643567 and the NASA grant 80NSSC18K1019.
Appendix A. Approximations for mass and thickness of pasta.-Our starting point is equations 27 and 28 from [68]. We define x i = 2(µ i − µ 0 )/m b c 2 where i = cc, p labels the location at which to calculate the quantities, the crust-core transition and onset of pasta respectively, and µ 0 ≈ 9 MeV is the baryon chemical potential at the surface of the star. Then we define from which the thickness relative to the stellar radius R of the whole crust i = cc or the crust above the pasta phases i = p is given as where β = GM/Rc 2 is the stellar compactness. Now, x i 1 so G i ≈ 1 + x i and G i − 1 1.
For most neutron stars β < 0.2 and so (1 − 2β) −1 ≈ 1 + 2β, and the denominators read and The thickness of the pasta phases is so the relative thickness of the pasta layers to the crust is Similar expressions were derived in [69], where it was demonstrated that the approximations made led to an accuracy of < 1% in the crust thickness, so similar accuracy can be expected for equation A.6.

7
The mass of a layer in the crust is proportional to the pressure at the bottom of the layer [66,69], so more straightforwardly As shown by Eq. 3 in [66] the moment of inertia of the crust is proportional to the mass of the crust to a first order of approximation, and so P cc controls the moment of inertia of the crust too, and ∆M p /∆M c ≈ ∆I p /∆I c . This is supported by inspecting Table. I of [66] and Table 5 and Fig. 8 of [23].