© H. Dinh Thi et al. 2021

1. Introduction

The interior of a cold neutron star (NS) is predicted to be made of various phases of matter, from a solid crust to a liquid core. In particular, the outer crust is thought to be made of ions arranged in a lattice, embedded in an electron gas, while in the inner crust the neutron-proton clusters, neutralised by the electron gas, also coexist with a (free) neutron gas, until, at about half the saturation density (n_sat ≈ 0.15 fm⁻³), nuclei dissolve into homogeneous matter, thus marking the transition to the liquid core (Haensel et al. 2007). In a liquid-drop picture of the nucleus, the equilibrium structure of the clusters in the crust results from the competition between the Coulomb and surface energy of nuclei. At lower densities, in the outer crust and in most regions of the inner crust, nuclear clusters are sufficiently far apart from each other that their structure is not influenced by neighbouring clusters. The effect of the surface energy prevails and clusters are thus expected to be spherical. However, at a higher density, nuclei become very close to each other, and, eventually, at the bottom of the crust, matter may arrange itself into various exotic configurations known as ‘pasta’ phases. Although so far there is no direct observational evidence of pasta phases, their existence may have a sizeable impact on different NS phenomena, such as NS cooling (Newton et al. 2013b; Horowitz et al. 2015; Lin et al. 2020), magnetic and rotational evolution of pulsars (Pons et al. 2013), crust oscillations (Gearheart et al. 2011; Sotani et al. 2012), and transport properties (see e.g., Schmitt & Shternin 2018 for a review). In addition, Gearheart et al. (2011) also estimated the effect of non-spherical configurations on the NS ‘mountains’, showing that the presence of pasta phases can decrease by up to an order of magnitude the NS maximum quadrupole ellipticity sustainable by the crust.

Since the pioneering works by Ravenhall et al. (1983), Hashimoto et al. (1984), and Oyamatsu et al. (1984), several studies have been conducted on the nuclear pasta, using different approaches, including compressible liquid-drop (CLD) models, the (extended) Thomas-Fermi method, nuclear energy-density functional theory, and molecular dynamics calculations (see Pethick & Ravenhall 1995; Chamel & Haensel 2008; Watanabe & Maruyama 2012; Blaschke & Chamel 2018 for a review and references therein). Very recently, Balliet et al. (2021) studied the pasta properties in NSs within a CLD model with surface parameters fitted on quantum 3D Hartree-Fock calculations of nuclei in a neutron gas, and calculated the distributions of pasta observables, namely the mass and the thickness of the pasta layer, using a Bayesian analysis. For all models they explored, they found that more than 50% of the crust by mass and 15% by thickness are made up of pasta, in agreement with the results of previous investigations (Lorenz et al. 1993; Newton et al. 2013a). Nevertheless, the presence of these non-spherical structures, such as cylindrical clusters (rods), slabs, cylindrical holes (tubes) and spherical holes (bubbles) in the NS, and the transition density among the different geometries, remain model dependent.

In this work, we studied the properties of the pasta phases in cold, non-accreting NSs under the so-called ‘catalysed matter hypothesis’, that is matter in its absolute ground state at zero temperature. To this aim, we extended the CLD model of Carreau et al. (2019a, 2020b) to account for non-spherical pasta structures in the inner crust. This formalism is described in Sect. 2. The use of the CLD approach allows us to address the question of the model dependence of the results, which is discussed in Sect. 3 considering several representative nuclear functionals. In addition, to obtain quantitative predictions of the uncertainties in the pasta observables, we performed a Bayesian analysis using flat uniform priors to generate posteriors with filters accounting for both our current knowledge of low-density nuclear physics and constraints at high density from general and NS physics. This investigation is assessed in Sect. 4. Finally, we present our conclusions in Sect. 5.

2. Theoretical framework

2.1. Model of the pasta phases

Following the seminal work by Baym et al. (1971), the equilibrium configuration of inhomogeneous catalysed matter (we assume zero temperature throughout the paper) present in the inner crust of NSs is variationally obtained within a CLD approach, along the same lines as in Douchin & Haensel (2001), and Carreau et al. (2019a, 2020b). We consider a periodic lattice configuration consisting of Wigner-Seitz cells of volume V_WS containing a clustered structure (‘pasta’) composed of Z protons of mass m_p and A − Z neutrons of mass m_n (A being the cluster total mass number), immersed in a uniform gas of neutrons and electrons with respective densities n_g and n_e. If inhomogeneities appear as clusters (holes) of volume V, the density distribution in the Wigner-Seitz cell is n_i (n_g) if l < r_N, and n_g (n_i) otherwise, r_N being the linear dimension of the pasta structure and l the linear coordinate of the Wigner-Seitz cell. The denser clustered phase is characterised by a density n_i = A/V (n_i = A/(V_WS − V) in the case of holes) and a proton fraction y_p = Z/A = n_p/n_i. The volume fraction occupied by the cluster (of density n_i) or hole (of density n_g), u = V/V_WS, thus reads

$$\begin{array}{c}\hfill u=\{\begin{array}{cc}(n_{\mathrm{B}}-n_{\mathrm{g}})/(n_{\mathrm{i}}-n_{\mathrm{g}})& \mathrm{for}\mathrm{clusters},\mathrm{and}\hfill \\ (n_{\mathrm{i}}-n_{\mathrm{B}})/(n_{\mathrm{i}}-n_{\mathrm{g}})& \mathrm{for}\mathrm{holes}.\hfill \end{array}\end{array}$$(1)

The energy density of the cell is minimised with the constraint of a given total baryonic density n_B = n_n + n_p, n_p (n_n) being the proton (neutron) density, respectively. Moreover, charge neutrality holds, thus n_e = n_p. The corresponding thermodynamic potential per unit volume is written as

$$\begin{array}{cc}& \mathrm{\Omega }=n_{\mathrm{n}}{m}_{\mathrm{n}}{c}^2+n_{\mathrm{p}}{m}_{\mathrm{p}}{c}^2+{ϵ}_{\mathrm{B}}(n_{\mathrm{i}},1-2y_{\mathrm{p}})f(u)\hfill \\ \hfill & +{ϵ}_{\mathrm{B}}(n_{\mathrm{g}},1)(1-f(u))+{ϵ}_{\mathrm{Coul}}+{ϵ}_{\mathrm{surf}+\mathrm{curv}}\hfill \\ \hfill & +{ϵ}_{\mathrm{e}}-{\mu }_{\mathrm{B}}^{\mathrm{tot}}{n}_{\mathrm{B}},\hfill \end{array}$$(2)

where ϵ_B(n, δ) is the energy density of uniform nuclear matter at baryonic density n and isospin asymmetry δ = (n_n − n_p)/n, ϵ_e is the energy density of a pure uniform electron gas at density n_e = n_p, ${\mu }_{\mathrm{B}}^{\mathrm{tot}}$ is the baryonic chemical potential (including the rest mass), and ϵ_{surf + curv} and ϵ_Coul are the finite-size corrections accounting for the interface tension between the cluster and the neutron gas and the electrostatic energy density, respectively. Finally, the function f(u) is given by

$$\begin{array}{c}\hfill f(u)=\{\begin{array}{cc}u& \mathrm{for}\mathrm{clusters},\mathrm{and}\hfill \\ 1-u& \mathrm{for}\mathrm{holes}.\hfill \end{array}\end{array}$$(3)

As was recognised in early works (see e.g., Ravenhall et al. 1983; Hashimoto et al. 1984), one of the advantages of the decomposition of Eq. (2) in terms of bulk and interface is that the geometry of the pasta structures only affects the finite-size corrections, and the latter can be expressed analytically as a function of the dimensionality of the structure (d = 1 for slabs, d = 2 for cylinders, d = 3 for spheres). We write the interface energy density as in Maruyama et al. (2005) and Newton et al. (2013a):

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{surf}+\mathrm{curv}}=\frac{\mathit{ud}}{r_{\mathrm{N}}}({\sigma }_{\mathrm{s}}+\frac{(d-1){\sigma }_{\mathrm{c}}}{r_{\mathrm{N}}}),\end{array}$$(4)

where the surface tension σ_s and curvature tension σ_c are independent of the dimensionality. For the latter quantities, we used the expressions originally proposed by Ravenhall et al. (1983) on the basis of Thomas-Fermi calculations at extreme isospin asymmetries, also subsequently employed in different works on NS crust and supernova modelling within the CLD approximation (Lattimer & Swesty 1991; Lorenz et al. 1993; Newton et al. 2013a; Carreau et al. 2019a, 2020b; Balliet et al. 2021), namely,

$$\begin{array}{cc}& {\sigma }_{\mathrm{s}}={\sigma }_0\frac{2^{p+1}+b_{\mathrm{s}}}{y_{\mathrm{p}}^{-p}+b_{\mathrm{s}}+{(1-y_{\mathrm{p}})}^{-p}},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill & {\sigma }_{\mathrm{c}}=5.5{\sigma }_{\mathrm{s}}\frac{{\sigma }_{0,\mathrm{c}}}{{\sigma }_0}(\beta -y_{\mathrm{p}}),\hfill \end{array}$$(6)

where the parameters (σ₀, σ_0,c, b_s, β, p) must be optimised on theoretical calculations or experimental data. The Coulomb energy density reads as follows:

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{Coul}}=2\pi {\left(e{y}_{\mathrm{p}}{n}_{\mathrm{i}}{r}_{\mathrm{N}}\right)}^{2}u{\eta }_d,\end{array}$$(7)

with e being the elementary charge, and

$$\begin{array}{cc}& {\eta }_1=\frac{1}{3}[u-2(1-\frac{1}{2u})],\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill & {\eta }_2=\frac{1}{4}[u-lnu-1],\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill & {\eta }_3=\frac{1}{5}[u+2(1-\frac{3}{2}{u}^{1/3})].\hfill \end{array}$$(10)

A nuclear model consists in a choice for the energy functional ϵ_B(n, δ), complemented by a set of values for the surface parameters (σ₀, σ_0,c, b_s, β, p). This point is detailed in the next sub-section. Once the nuclear model is specified, the pasta structure and composition at a given baryonic density n_B is determined by a two-step process. First, a geometry d and a shape s (clusters or holes) is considered, and the thermodynamical potential of Eq. (2) is minimised with respect to the variational parameters (n_i, I = 1 − 2y_p, A, n_p, n_g). This allows us to identify the baryonic chemical potential $\mu ={\mu }_{\mathrm{B}}^{\mathrm{tot}}-m_{\mathrm{n}}{c}^2$ with the chemical potential of the neutron gas,

$$\begin{array}{c}\hfill \mu =\frac{d{ϵ}_{\mathrm{B}}(n_{\mathrm{g}},1)}{\mathrm{d}{n}_{\mathrm{g}}},\end{array}$$(11)

and gives the optimal value of Ω = Ω_opt(d, s) for each geometry. Then, the equilibrium configuration is defined as the one corresponding to the values of d and s that produce the minimum value of Ω_opt.

2.2. Energy functional

The energy density ϵ_B of homogeneous nuclear matter is poorly known out of the saturation point n_sat ≈ 0.15 fm⁻³ of symmetric n_n = n_p matter, and a large number of models have been proposed in the literature based on effective Hamiltonians or Lagrangians. Margueron et al. (2018) showed that those different models can be accurately reproduced using a Taylor expansion in x = (n − n_sat)/3n_sat up to order 4 around the saturation point (n = n_sat, δ = 0), varying the parameters of the expansion, which correspond to the so-called equation-of-state empirical parameters:

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{B}}(n,\delta )\approx n{\displaystyle \sum _{m=0}^4\frac{1}{m!}({\frac{d^m{e}_{\mathrm{sat}}}{d{x}^m}|}_{x=0}+{\frac{d^m{e}_{\mathrm{sym}}}{d{x}^m}|}_{x=0}{\delta }^2)x^m.}\end{array}$$(12)

In Eq. (12), e_sat = ϵ_B(n, 0)/n is the energy per baryon of symmetric matter, and e_sym = (ϵ_B(n, 1)−ϵ_B(n, 0))/n is the symmetry energy per baryon, defined here as the difference between the energy of pure neutron matter and that of symmetric matter. To speed up the series convergence, a δ^5/3 term from the fermionic zero-point energy is added, as well as an exponential correction ensuring the correct limiting behaviour at zero density (see Eq. (17) in Carreau et al. 2019a). Following the common notation in the literature, we denote E_sat(sym) = e_sat(sym)(n = n_sat), while the successive derivatives of e_q, with q= sat,sym, are called L_q, K_q, Q_q, Z_q. The bulk parameters {E_q, L_q, K_q, Q_q, Z_q, q= sat,sym} are complemented with the saturation density parameter, n_sat, two parameters related to the isoscalar effective mass and effective mass splitting, which are denoted as 𝒦_q, and the b parameter governing the functional behaviour close to the zero-density limit (see Sect. 2.2 in Carreau et al. 2019a for details). The complete parameter set thus has 13 parameters and will be noted in a compact form as X_bulk ≡ {n_sat, b, (E_q, L_q, K_q, Q_q, Z_q, 𝒦_q, q= sat,sym)}. Different nuclear models will then correspond to different sets of X_bulk parameters.

Concerning the parameters (σ₀, σ_0,c, b_s, β, p) entering the surface functional, Eq. (4), Carreau et al. (2019a) showed that, in order to have a realistic and consistent treatment of the composition of the crust, it is important that these parameters are fixed consistently with the functional employed for the bulk part of the nuclear energy ϵ_B. Indeed, the bulk energy is model dependent, but in the case of a spherical geometry the total (bulk plus surface plus Coulomb) energy is constrained by the requirement of reproducing the experimentally measured nuclear masses, which creates an obvious correlation between the bulk and surface parameters that should be accounted for.

In the vacuum, the nuclear mass corresponding to a spherical fully ionised atom of charge Z and mass number A can be deduced from Eqs. (2), (4), and (7) as:

$$\begin{array}{cc}& M(A,Z)c^2=m_{\mathrm{p}}{c}^{2}Z+m_{\mathrm{n}}{c}^2(A-Z)\hfill \\ \hfill & +\frac{A}{n_0}{ϵ}_{\mathrm{B}}(n_0,I)+4\pi r_{\mathrm{N}}^2({\sigma }_{\mathrm{s}}+\frac{2{\sigma }_{\mathrm{c}}}{r_{\mathrm{N}}})\hfill \\ \hfill & +\frac{3}{5}\frac{e^2{Z}^2}{r_{\mathrm{N}}},\hfill \end{array}$$(13)

where I = 1 − 2Z/A, the nuclear radius is r_N = (4πn₀/3)^−1/3A^1/3, and the bulk density n₀ is given by the equilibrium density of nuclear matter at isospin asymmetry I, defined by ∂ϵ_B/∂n|_I, n0 = 0.

For each choice of the parameter set X_bulk, we determine the associated surface parameters by a χ²-fit of Eq. (13) to the experimental Atomic Mass Evaluation (AME) 2016 (Wang et al. 2017).

As an example, in Table 1 we give the bulk parameters that reproduce the density behaviour of symmetric nuclear matter and pure neutron matter of a chosen set of popular nuclear models that are analysed in the next section. The associated optimal surface parameters are given in Table 2. The (σ₀, b_s, σ_0,c, β) parameters were obtained from the fit of the experimentally measured nuclear masses (Wang et al. 2017), while the p parameter was optimised to provide a good reproduction of the crust–core transition density of the different functionals, whenever available, or fixed to p = 3 otherwise (see Table 3 and Carreau et al. 2019a, 2020b for a discussion).

3. Comparison of different nuclear functionals

Using the formalism described in Sect. 2, we computed the properties of the pasta layer predicted at the bottom of the inner crust. To explore the model dependence of the results, we examined the equilibrium configurations obtained with different popular nuclear models, corresponding to the different parameter sets listed in Tables 1 and 2, as illustrative examples. Results are shown in Fig. 1, where the different colours correspond to the density regions where different geometries (spheres, rods, slabs, tubes, and possibly bubbles) dominate. The upper edge of each column gives the transition point from the inhomogeneous crust to the homogeneous core, defined as the point where the thermodynamical potential, Eq. (2), corresponding to the optimal geometry, equals the thermodynamical potential of homogeneous nuclear matter in beta equilibrium. The sequence of the different geometries appears to be model independent, although not all the considered models predict bubble configurations, and it is consistent with previous results (see Pethick & Ravenhall 1995 for a review).

Fig. 1. Sequence of equilibrium pasta phases and corresponding transition densities among the geometries for different nuclear functionals (see text for details).

For a given energy functional, it is well known that the precise value of the transition density from the spherical to the non-spherical phases depends on the details of the many-body model used to treat the clustered structure. A comparison with results in the literature is reported in Table 3. The most sophisticated calculation presently available is the one by Pearson et al. (2020) with full fourth-order extended Thomas-Fermi calculations within the BSk24 functional. In that paper, the authors obtained a core-crust transition density very close to the one displayed in Fig. 1 for BSk24 and a slightly higher value for the transition density to the rod shape (see Table 3). We consider this comparison as a very good agreement, especially considering that beyond-mean-field contributions such as pairing, which were fully neglected here, have recently been shown to affect the composition of the crust (Pearson et al. 2020; Shelley & Pastore 2020).

In the case of the SLy4 functional, Martin & Urban (2015) employed a slightly less sophisticated second-order extended Thomas-Fermi approach, and they observed a transition to the cylindrical shape at a higher density than the value obtained with our approach for SLy4. However, the thermodynamical potentials corresponding to d = 3 and d = 2 in Martin & Urban (2015) are almost indistinguishable starting from n = 0.05 fm⁻³, in good agreement with our results. The crust-core transition point is also in excellent agreement with our findings. At variance with this result, Viñas et al. (2017) and Douchin & Haensel (2001), who also employed the SLy4 functional, reported no deviation from the spherical shape within a zeroth-order Thomas-Fermi calculation and a CLD approach, respectively. Moreover, they also obtain a lower density for the crust-core transition. Finally, we can also compare our results concerning the relativistic functionals NL3, DD-ME2, and DD-MEδ with the extensive Thomas-Fermi calculations of Grill et al. (2012); in our approach, we obtain a lower transition density between spheres and cylinders, and we do find non-spherical configurations for the NL3 functional, unlike Grill et al. (2012).

Even if the comparison of our work with previous results can be globally considered as satisfactory, Table 3 shows that some differences exist between the inner-crust composition of different many-body methods used to compute the energy of clusterised matter, for a fixed equation of state (or, equivalently, a fixed set of bulk parameters). This highlights the importance of finite-size contributions to the nuclear energy, which are essential in determining the optimal composition and are not uniquely linked to the bulk matter properties. Only the bulk parameters of our meta-model are adjusted to reproduce a given functional, while we employ a different fitting protocol for the surface (plus curvature) energy. The effects of the surface contribution are studied in greater detail in Sect. 4.2.

A complementary piece of information on the inner-crust composition in the presence of non-spherical geometries is given by Fig. 2, which reports the total proton fraction in the cell, $Y_{\mathrm{p}}^{\mathrm{WS}}$, as a function of the baryonic density for the different geometries, for a selected set of functionals, as illustrative examples. We can see that some of the functionals exhibit the characteristic parabolic shape already reported by Pearson et al. (2020) for the BSk24 model, and the absolute value of the proton fraction is also in reasonable agreement with the findings of Pearson et al. (2020). However, we can also observe that the trend of the proton fraction, the effect of the geometry, and the numerical value throughout the inner crust are very strongly model dependent.

Fig. 2. Proton fraction in the Wigner-Seitz cell as a function of the baryonic density in the neutron-star crust for different geometries of the clustered structures using different popular nuclear models (see text for details).

To summarise the previous discussion, from Figs. 1 and 2 we can conclude that the uncertainty in the predictions due to the choice of the nuclear functional is much more important than the uncertainty due to the many-body method adopted to describe the inhomogeneities. We limited our analysis to a small set of popular models, that span a quite wide range of the empirical parameters, to provide illustrative examples. However, these nine models might not fully encompass the current uncertainties in the nuclear equation of state. For this reason, in the next section we assess the model dependence of the crust composition in a more quantitative way by largely exploring the parameter space of the nuclear functionals compatible with the present theoretical and experimental constraints.

4. Bayesian analysis

We now come to the quantitative determination of the uncertainty on the pasta observables, due to the incomplete knowledge of the nuclear energy functional. To this aim, we performed a Bayesian analysis by largely varying the model parameters X = {X_bulk, p} using flat non-informative priors, and generating posterior distribution with filters that impose both our present low-density (LD) nuclear physics knowledge and high-density (HD) constraints coming from general and NS physics:

$$\begin{array}{c}\hfill p_{\mathrm{post}}(\mathit{X})=\mathcal{N}{w}_{\mathrm{LD}}(\mathit{X})w_{\mathrm{HD}}(\mathit{X})e^{-{\chi }^2(\mathit{X})/2}{p}_{\mathrm{prior}}(\mathit{X}),\end{array}$$(14)

where 𝒩 is the normalisation. A first filter w_LD on the bulk parameters is given by the uncertainty band of the chiral N³LO effective field theory (EFT) calculation for symmetric matter and pure neutron matter by Drischler et al. (2016), which is interpreted as a 90% confidence interval. Moreover, the quality of reproduction of the experimental nuclear masses M_exp is defined from the error estimator ${\chi }^2(\mathit{X},{\sigma }_0,{\sigma }_{0,\mathrm{c}},b_{\mathrm{s}},\beta )$:

$$\begin{array}{c}\hfill {\chi }^2=\frac{1}{N_{\mathrm{dof}}}{\displaystyle \sum _{\mathrm{n}}\frac{{(M(A_{\mathrm{n}},Z_{\mathrm{n}})-M_{\exp }(A_{\mathrm{n}},Z_{\mathrm{n}}))}^2}{{\sigma }_n^2},}\end{array}$$(15)

where N_dof is the number of degrees of freedom, the sum runs over the AME2016 nuclear mass table from Wang et al. (2017), M is calculated from Eq. (13) for each model described by a combined set $\{\mathit{X},{\sigma }_0,{\sigma }_{0,\mathrm{c}},b_{\mathrm{s}},\beta \}$, and σ_n corresponds to the systematic theoretical error. As for the p parameter, unless specified otherwise (see Sect. 4.2), its value was allowed to vary in the range from 2 to 4.

A further constraint is given by the condition that the minimisation of the thermodynamic potential Eq. (2) leads to physically meaningful results for the crust, namely positive values for the optimal gas and cluster densities. Finally, the high-density behaviour of the equation of state is also controlled through the w_HD filter in Eq. (14) by imposing (i) stability; (ii) causality; (iii) a positive symmetry energy at all densities; and (iv) the resulting equation of state to support M_max > 1.97 M_⊙, where M_max is the maximum NS mass at equilibrium determined from the solution of the Tolmann-Oppenheimer-Volkoff (TOV) equations (Haensel et al. 2007), M_⊙ being the solar mass. The resulting posterior equations of state were shown by Carreau et al. (2019b) and Carreau (2020) to be fully compatible with the recent measurement of the tidal polarisability parameter Λ in the gravitational-wave event GW170817 (Abbott et al. 2018) (see Carreau et al. 2019b for details and their Fig. 1).

Out of the 10⁸ models generated to numerically sample the prior parameter distribution, 7008 models are retained by the applied filters¹ to compute marginalised posteriors for the different observables Y whose average value is given by

$$\begin{array}{c}\hfill ⟨Y⟩={\displaystyle \prod _{k=1}^N{\int }_{X_k^{\min }}^{X_k^{\max }}d{X}_{k}Y(\mathit{X})p_{\mathrm{post}}(\mathit{X}),}\end{array}$$(16)

where p_post(X) is the posterior distribution, Y(X) is the value of the Y variable as obtained with the X parameter set, $X_k^{\text{min}(\text{max})}$ is the minimum (maximum) value in the prior distribution of parameter X_k taken from Table 2.2 in Carreau (2020) (see also Margueron et al. 2018 for details)² and of p from 2 to 4, and N = 14 is the number of parameters defining the nuclear model. The marginalised posteriors of given observables can then also be confronted with additional constraints, not originally included in the filters, to check the validity of the assumptions and the flexibility of the prior. An example of such comparison with both recent astrophysical and nuclear-physics data is given in Fig. 3. In the top panel, we show the marginalised posterior for the gravitational mass of the NS as a function of the radius. The shaded green areas correspond to the 68% (1σ), 95% (2σ), and 99% (3σ) (from dark to light green) confidence intervals, while the blue areas represent the constraint inferred from the gravitational-wave event GW170817 by the LIGO/Virgo collaboration (Abbott et al. 2018), and the pink areas show the recent results from NICER (Miller et al. 2019, 2021), all at the 2σ level. For comparison, in the same figure we also show the mass-radius relation obtained with the nine models analysed in Sect. 3. We can see that our posterior distribution is in agreement with the recent observations coming from NICER (Miller et al. 2019), NICER combined with XMM-Newton (Miller et al. 2021), and those inferred from GW170817 (Abbott et al. 2018). On the other hand, some of the popular models displayed lie outside the posterior distribution, meaning that they would be filtered out by our low- and/or high-density constraints. In the bottom panels of Fig. 3, we plot the marginalised posterior distributions for the symmetry energy parameters E_sym and L_sym. Very recently, Reed et al. (2021) inferred the values of E_sym = (38.1 ± 4.7) MeV and L_sym = (106 ± 37) MeV from measurements of the neutron-skin thickness in ²⁰⁸Pb performed by the PREX collaboration (shown as shaded rectangles in the figure). These findings point towards rather high values of these coefficients, overestimating the present limits deduced from both theoretical and experimental measurements (see e.g., Baldo & Burgio 2016; Burgio & Fantina 2018 for a discussion, and references therein). Indeed, our posterior distributions suggest E_sym = (30.8 ± 1.3) MeV and L_sym = (47.3 ± 9.2) MeV, meaning that the average values we obtain are considerably lower than those reported in Reed et al. (2021). However, the estimation of these parameters made in the latter work is not model independent, and it is based on a correlation extracted by employing a limited number of models that might not span all the possible functional dependences. In our study, the low-density constraint imposed by EFT calculations would disfavour such high values of the symmetry energy parameters, and models that predict high E_sym and L_sym are filtered out in our posterior distribution. This is in full agreement with Reed et al. (2021), where the tension between the extracted values of L_sym and the ab initio calculations was already highlighted (see Fig. 1 of Reed et al. 2021). We can also observe that the NL3 and TM1 models that have, among the models considered in Sect. 3, the closest E_sym and L_sym values to those reported in Reed et al. (2021), lie outside our posterior distribution and appear in disagreement with the recent astrophysical constraints from GW170817 (see top panel of Fig. 3).

Fig. 3. Top panel: marginalised posterior distribution for the neutron-star gravitational mass as a function of radius for the models including both the low- and high-density filters. The shaded areas correspond to the 1σ, 2σ, and 3σ (from dark to light green) confidence intervals, while lines show the mass-radius relation for different selected models. The pink (blue) shaded areas represent the constraint from Miller et al. (2019, 2021) (Abbott et al. 2018) at 2σ. Bottom panels: marginalised posterior distributions for the symmetry energy at saturation, Esym (left panel), and its slope, Lsym (right panel). The shaded rectangles show the values of Esym and Lsym from Reed et al. (2021) (see text for details).

The theoretical ab initio chiral EFT calculation and the experimental AME2016 mass measurements constitute two powerful and complementary constraints allowing us to predict pasta observables with controlled uncertainties. The theoretical filter only acts on the bulk parameters, and it is particularly effective in constraining the parameters in the symmetry sector X_sym. Conversely, the mass filter is not very effective in selecting the bulk parameters, because a stronger binding from the bulk term can be compensated by an increased importance of the surface terms via the correlation imposed by Eq. (13). However, the mass filter is crucial in fixing the surface properties of inhomogeneous matter and, consequently, the crustal properties of NSs. These two aspects are analysed separately in the next sub-sections.

4.1. The influence of the bulk functional

The chiral EFT constraint was applied in numerous previous studies of static properties of NSs (Carreau et al. 2019b; Lim et al. 2019; Tews et al. 2019; Essick et al. 2020; Brown et al. 2020). In these studies, the accent is typically put on the high-density equation of state, which is the dominant ingredient of the different astrophysical observables that are integrated over the whole NS volume, such as the tidal polarisability or the NS mass and radius. The compatibility of the functionals with the ab initio predictions in the very low-density region below n ≈ 0.05 fm⁻³ has a negligible influence for these global NS observables, and as such it was often overlooked.

However, as long as crustal observables are concerned, this very low-density region is crucial. This can be seen in Fig. 4, where we represent the posterior distributions of the density and pressure at the crust-core transition, (n_CC, P_CC) (top panels), and the density and pressure at the transition from spherical nuclei to non-spherical pasta structures, (n_p, P_p) (lower panels). For comparison, we also indicate the predictions corresponding to some selected models with arrows: BSk24 (black-filled arrow), RATP (red-filled arrow), DD-ME2 (green-filled arrow), and TM1 (yellow-filled arrow). Two different density intervals are considered for the application of the chiral EFT filter, namely [0.1, 0.2] fm⁻³ and [0.02, 0.2] fm⁻³. We can see that the consideration of the very low-density part of the equation of state has a notable effect on the determination of the transition points. In particular, a significant number of models that correspond to reasonable properties of nuclear matter close to saturation, and therefore fulfill the chiral EFT condition for n ≥ 0.1 fm⁻³, are seen to produce very low values for the crust-core and the sphere-pasta transition point. This implies a very thin crust for the NS and a small or negligible contribution of the pasta phases, as shown in Fig. 5, which displays the prediction for the relative crustal thickness associated with the non-spherical pasta structures, R_pasta/R_crust (top panel), the associated mass, m_pasta/m_crust (middle panel), and the fraction of the moment of inertia, I_pasta/I_crust (lower panel). The two choices for the low-density constraints and the arrows indicating the predictions of some selected models (BSk24, RATP, DD-ME2, and TM1) are the same as those adopted in Fig. 4. These observables are computed from the numerical solution of the TOV equations, the crust radius (mass) being calculated as the difference between the total NS radius (mass) and that of the core (we considered NSs with a total mass equal to the maximum mass predicted by each computed model). As can be seen in Figs. 4 and 5, when the chiral EFT constraint is applied from n ≥ 0.02 fm⁻³, almost all models corresponding to a very small (or even null) pasta contribution are filtered out of the posterior distribution. Knowing that the precision and reliability of the ab initio calculations improve with decreasing density, we can conclude that our present theoretical understanding of nuclear matter implies that pasta phases should exist in the inner crust of NSs.

Fig. 4. Posterior distributions of crust-core (top panels) and sphere-pasta (lower panels) transition density (left panels) and pressure (right panels), obtained considering different density intervals for the application of the chiral EFT constraint. Arrows correspond to the predictions of some selected models: BSk24 (black), RATP (red), DD-ME2 (green), and TM1 (yellow) (see text for details).

Fig. 5. Posterior distributions of the thickness (top panel), mass (middle panel), and fractional moment of inertia (lower panel) of the pasta layer with respect to the whole crust for a neutron star with a mass M = Mmax. The low-density EFT constraints and the arrows corresponding to the predictions of selected models are the same as in Fig. 4 (see text for details).

To investigate this point further, in the top panel of Fig. 6 we plot the energy per baryon distribution of the different models for several baryon densities, for symmetric (δ = 0) and pure neutron matter (δ = 1). The corresponding uncertainty bands of the chiral EFT constraint from Drischler et al. (2016) are also given (delimited by black dash-dotted lines). Models for which EFT constraints are applied from n ≥ 0.02 fm⁻³ (n ≥ 0.1 fm⁻³ but predicting a crust-core transition lower than 0.05 fm⁻³, see Fig. 4), are displayed in coral (light blue) on the left (right) part of the density axis. We can clearly see that models that are not compatible with the EFT constraints in the range 0.02 ≤ n < 0.1 fm⁻³ and yielding a crust-core transition below 0.05 fm⁻³ also predict lower energy per baryon in the very low-density region. This results in a stiffer equation of state in the NS crust in the baryon density range 0.03 ≲ n_B ≲ 0.1 fm⁻³ (or, equivalently, in the mass-energy range 0.5 × 10¹⁴ ≲ ρ_B ≲ 1.8 × 10¹⁴ g cm⁻³), as can been seen by comparing, in the lower panel of Fig. 6, the red and blue bands that represent the total pressure in the Wigner-Seitz cell as a function of the baryon (lower axis) and mass-energy (top axis) density in the inner crust. It is also clear that applying the EFT constraint from n ≥ 0.02 fm⁻³ induces a stricter filter and thus a narrower band for the resulting equation of state. Figures 4–6 thus show that it is essential to verify the compatibility with the ab initio calculations down to very low densities to obtain realistic predictions for the properties of the NS crust. Incidentally, Shelley & Pastore (2021), who performed a systematic study of the NS inner-crust composition using several Skyrme functionals, also pointed out the importance of constraining the low-density part of the equation of state.

Fig. 6. Top panel: bands of the energy per baryon of symmetric (δ = 0) and pure neutron matter (δ = 1) as a function of density representing the chiral EFT constraint from Drischler et al. (2016). The probability distributions of models for which EFT constraints are applied from n ≥ 0.02 fm−3 (n ≥ 0.1 fm−3 and predicting a crust-core transition nCC < 0.05 fm−3) are represented as a violin-like shape. Lower panel: 1σ band of the total pressure versus baryon density nB and versus mass-energy density ρB (in units of 1014 g cm−3) predicted by these models (see text for details).

The relative influence on the pasta observables of the different parameters governing the bulk properties of nuclear matter can be quantified from the Pearson linear correlation coefficients displayed in Fig. 7 for the pasta thickness with respect to the crustal one (top panel), and the relative fractional moment of inertia (lower panel). We do not display the Pearson coefficients for the mass of the pasta layers here since they are equivalent to those of the fractional moment of inertia. In the absence of physical constraints (line labelled ‘prior’ in Fig. 7), no very strong correlation can be observed between the bulk parameters and the pasta observables, as expected from the even exploration of the parameter space of our prior distribution. Some loose correlation is observed with E_sat and the second (K_sym) and third (Q_sym) derivatives of the symmetry energy at saturation n_sat. The prior distributions of the fourth-order parameter Z_sym and of the high-order parameters in the symmetric energy sector, Q_sat and Z_sat, are also very large. However, they are less influential because the inner crust is characterised by a very low proton fraction (see Fig. 2), and its density is close to n_sat, which implies a small contribution of the high-order derivatives in the Taylor expansion Eq. (12). When only the models compatible with our present knowledge of nuclear matter are retained (lines labelled ‘LD+HD’ in Fig. 7), the correlations with the physically significant parameters start to appear. It is interesting to observe that the energy functional must be controlled in the whole sub-saturation region (line labelled ‘n ≥ 0.02 fm⁻³’ in Fig. 7) for the correlations to appear. In the absence of this constraint, large compensations are possible among the different terms of the functional, thus washing up the physical correlations. The most relevant parameters are seen to be the energy of symmetric matter E_sat and symmetry energy E_sym at saturation, as well as the slope of the symmetry energy L_sym. These parameters are already relatively well constrained by nuclear theory and experiments, but we can expect that their uncertainty will be further reduced in upcoming studies, which will lead to an increased precision in the determination of the pasta contribution to the physics of the crust.

Fig. 7. Pearson correlations between pasta properties and the bulk, surface, and curvature parameters. Two different density intervals for the application of the chiral EFT constraint are considered (see text for details).

4.2. The influence of the surface tension

Figure 7 also shows that a better modelling of the pasta properties would be obtained if we could better constrain the energy functional close to saturation. However, the most striking information of Fig. 7 is the fact that the surface parameters are more influential than the bulk properties, as shown by the high correlation coefficients. A part of the uncertainty on the surface tension can certainly be ascribed to the very simplified model of nuclear masses that we employ, Eq. (13), based on a spherical CLD approximation without accounting for deformation, shell, and pairing effects. A more conceptual limitation on the determination of those coefficients is also the fact that the information on the nuclear mass is not enough to fully pin down the surface properties, due to the partial compensation between surface and bulk.

We can thus conclude that the strong influence of the surface parameters underlines the importance of a unified treatment of the equation of state, which should be employed not only in developing specific equation-of-state models, but also in the statistical studies of NS properties.

To assess the importance of a self-consistent calculation of surface properties from the assumed bulk energy functional in the modelling of the NS crust, we performed two calculations where the surface parameters in the Bayesian analysis are fixed to the values obtained from two different accurate fits of extended Thomas-Fermi calculations of a large pool of nuclei. In the work by Carreau et al. (2020b), the full BSk24 extended Thomas-Fermi mass table was fitted from dripline to dripline, while Furtado & Gulminelli (2020) used the SLy4 functional and included in the fit extended Thomas-Fermi calculations beyond drip with global proton fractions as low as $Y_{\mathrm{p}}^{\mathrm{WS}}=0.02$. The resulting surface parameters are reported in Table 4. These fits lead to surface energy functionals slightly different from the ones obtained with the parameters of Table 2 because of the larger pool of nuclei included in the optimisation. Moreover, in both cases, the extension of the mass table to extremely neutron-rich nuclei has allowed also an optimal determination of the p parameter (see Eq. (5)), which governs the extreme isospin behaviour of the surface tension and cannot be fixed from the properties of terrestrial nuclei.

The surface tensions in Table 4 are more accurately determined than those deduced from the global reproduction of measured nuclear masses (see Table 2) and which depend on the empirical X_bulk parameters. However, they are not consistent with the bulk properties of functionals different from BSk24 and SLy4, respectively. The resulting posterior distribution³ for the pasta obervables is displayed in Fig. 8. We can see that the use of a surface tension that is not consistent with the bulk functional (green dashed lines and red dash-dotted lines) leads to a small shift in the most probable values of the represented pasta-layer properties, and an important deformation of the distribution, with a clear underestimation of the uncertainties in the pasta observables.

Fig. 8. Posterior distributions of thickness (top panel), mass (middle panel), and fractional momentum of inertia (bottom panel) of the pasta layer with respect to the whole crust. Blue solid curves: surface parameters optimised on the energy functional; red dash-dotted (green dashed) curves: fixed surface parameters from Furtado & Gulminelli (2020) (Carreau et al. 2020b). The chiral EFT filter is applied from n = 0.02 fm−3.

5. Conclusions

We presented a study of the properties of pasta phases in cold catalysed NSs, extending the work of Carreau et al. (2019a, 2020b). In particular, we applied a CLD model with surface parameters adjusted either on experimental nuclear masses or theoretical calculations. We addressed the model dependence of the results by employing different nuclear models, which are reproduced by fixing the bulk parameters to reproduce the low-density behaviour of symmetric matter and pure neutron matter of the corresponding functional. We find that the value of the transition densities between different geometries is model dependent. The sequence of the predicted geometries, namely spheres, rods, slab, and tubes, appears to be the same for all considered models, while the presence of bubbles before matter becomes homogeneous in the core is only predicted by some of the considered models. Moreover, our findings show that the uncertainty in the predictions due to the choice of the nuclear functional is more significant than that due to the many-body method adopted to describe the inhomogeneities.

To assess the model uncertainties in the pasta observables, we performed a Bayesian analysis by largely varying the model parameters using uniform priors and generating posterior distributions with filters accounting for both our present low-density nuclear physics knowledge and high-density general and NS physics constraints. Our results show that the low-density nuclear physics constraints are crucial in determining the crustal and pasta observables. Here, the low-density filters, given by the uncertainty band of the chiral N³LO EFT calculation for symmetric matter and pure neutron matter by Drischler et al. (2016) are applied from n ≥ 0.02 fm⁻³. Recent works have been devoted to deriving EFT-inspired energy-density functionals benchmarked by ab initio calculations at even lower density (0 ≲ n ≲ 0.1 fm⁻³), and these could constitute additional constraints in the future (Yang et al. 2016). Our final predictions for the different crustal properties, namely the fractional thickness R_pasta/R_crust, moment of inertia I_pasta/I_crust, and mass m_pasta/m_crust of the pasta layer for a NS with a mass equal to its maximum mass, are summarised in Table 5. The left (right) column lists the results obtained when the low-density EFT filters on the energy of homogeneous matter are applied from n ≥ 0.02 fm⁻³ (n ≥ 0.1 fm⁻³).

The most relevant nuclear parameters for predicting the pasta properties are the lower-order empirical parameters, namely the energy of symmetric matter, the symmetry energy, and the slope of the symmetry energy at saturation, as well as the surface and curvature parameters. In particular, we find that the latter parameters are more influential than the bulk parameters for the description of the pasta observables. In addition, we show that the use of a surface tension that is inconsistent with the bulk functional leads to an underestimation of the uncertainties in the pasta properties, thus highlighting the importance of a consistent calculation of the nuclear functional.

This formalism can be extended at finite temperature, similarly to the approach adopted in Fantina et al. (2020) and Carreau et al. (2020a) for the outer and inner crusts, respectively, to account for the so-called impurities. This question is of particular interest since the presence of a pasta layer at the bottom of the NS crust may have sizeable impact on different NS phenomena including NS cooling, and this will be addressed in a future work. Moreover, very recently, upper limits on the NS equatorial ellipticity were provided through gravitational-wave data by the LIGO/Virgo collaboration (Abbott et al. 2020). The possible impact of pasta phases on this observable, which was already pointed out by Gearheart et al. (2011), thus deserves further investigation.

Acknowledgments

This work has been partially supported by the IN2P3 Master Project NewMAC, the CNRS PICS07889, and the CNRS International Research Project (IRP) “Origine des éléments lourds dans l’univers: Astres Compacts et Nucléosynthèse (ACNu)”. We thank Brynmor Haskell for interesting suggestions.

References

All Tables

All Figures

	Fig. 1. Sequence of equilibrium pasta phases and corresponding transition densities among the geometries for different nuclear functionals (see text for details).
In the text

	Fig. 2. Proton fraction in the Wigner-Seitz cell as a function of the baryonic density in the neutron-star crust for different geometries of the clustered structures using different popular nuclear models (see text for details).
In the text

Fig. 3. Top panel: marginalised posterior distribution for the neutron-star gravitational mass as a function of radius for the models including both the low- and high-density filters. The shaded areas correspond to the 1σ, 2σ, and 3σ (from dark to light green) confidence intervals, while lines show the mass-radius relation for different selected models. The pink (blue) shaded areas represent the constraint from Miller et al. (2019, 2021) (Abbott et al. 2018) at 2σ. Bottom panels: marginalised posterior distributions for the symmetry energy at saturation, Esym (left panel), and its slope, Lsym (right panel). The shaded rectangles show the values of Esym and Lsym from Reed et al. (2021) (see text for details).

In the text

	Fig. 4. Posterior distributions of crust-core (top panels) and sphere-pasta (lower panels) transition density (left panels) and pressure (right panels), obtained considering different density intervals for the application of the chiral EFT constraint. Arrows correspond to the predictions of some selected models: BSk24 (black), RATP (red), DD-ME2 (green), and TM1 (yellow) (see text for details).
In the text

	Fig. 5. Posterior distributions of the thickness (top panel), mass (middle panel), and fractional moment of inertia (lower panel) of the pasta layer with respect to the whole crust for a neutron star with a mass M = Mmax. The low-density EFT constraints and the arrows corresponding to the predictions of selected models are the same as in Fig. 4 (see text for details).
In the text

Fig. 6. Top panel: bands of the energy per baryon of symmetric (δ = 0) and pure neutron matter (δ = 1) as a function of density representing the chiral EFT constraint from Drischler et al. (2016). The probability distributions of models for which EFT constraints are applied from n ≥ 0.02 fm−3 (n ≥ 0.1 fm−3 and predicting a crust-core transition nCC < 0.05 fm−3) are represented as a violin-like shape. Lower panel: 1σ band of the total pressure versus baryon density nB and versus mass-energy density ρB (in units of 1014 g cm−3) predicted by these models (see text for details).

In the text

	Fig. 7. Pearson correlations between pasta properties and the bulk, surface, and curvature parameters. Two different density intervals for the application of the chiral EFT constraint are considered (see text for details).
In the text

Fig. 8. Posterior distributions of thickness (top panel), mass (middle panel), and fractional momentum of inertia (bottom panel) of the pasta layer with respect to the whole crust. Blue solid curves: surface parameters optimised on the energy functional; red dash-dotted (green dashed) curves: fixed surface parameters from Furtado & Gulminelli (2020) (Carreau et al. 2020b). The chiral EFT filter is applied from n = 0.02 fm−3.

In the text