Symmetry energy and neutron star properties constrained by chiral effective field theory calculations

We investigate the nuclear symmetry energy and neutron star properties using a Bayesian analysis based on constraints from different chiral effective field theory calculations using new energy density functionals that allow for large variations at high densities. Constraints at high densities are included from observations of GW170817 and from NICER. In particular, we show that both NICER analyses lead to very similar posterior results for the symmetry energy and neutron star properties when folded into our equation-of-state framework. Using the posteriors, we provide results for the symmetry energy and the slope parameter, as well as for the proton fraction, the speed of sound, and the central density in neutron stars. Moreover, we explore correlations of neutron star radii with the pressure and the speed of sound in neutron stars. Our 95\% credibility ranges for the symmetry energy $S_v$, the slope parameter $L$, and the radius of a 1.4$\msun$ neutron star, $R_{1.4}$, are $S_v=(30.6\text{--}33.9)$\,MeV, $L=(43.7\text{--}70.0)$\,MeV, and $R_{1.4}=(11.6\text{--}13.2)$\,km. Our analysis for the proton fraction shows that larger and/or heavier neutron stars are more likely to cool rapidly via the direct Urca process. Within our equation-of-state framework a maximum mass of neutron stars $M_{\rm max}>2.1\msun$ indicates that the speed of sound needs to exceed the conformal limit.


I. INTRODUCTION
Understanding dense matter is a central challenge in nuclear physics and astrophysics.In nature, dense matter exists in the cores of neutron stars under extreme neutron-rich conditions.The properties of neutron-rich matter around nuclear densities are described by the nuclear symmetry energy and its density dependence.While there have been impressive constraints from nuclear theory, nuclear experiments, and astrophysics (see, e.g., Refs.[1][2][3][4]), more precise determinations of the symmetry energy and its slope parameter L at saturation density, n 0 = 0.16 fm −3 , are still an open problem.
Recently, new EDFs for the nuclear EOS were introduced by Huth et al. [3], and have the advantage of providing high-density extrapolations that are consistent with causality and with a maximum of the speed of sound.These functionals allow for EOS calculations for the broad ranges of conditions reached in core-collapse supernovae and neutron star mergers.In this work, we use these new EDF EOSs to constrain the symmetry energy and neutron star properties based on a prior informed by chiral EFT calculations of neutron matter.
From the astrophysics side, the strongest constraint on the nuclear EOS comes from the observation of heavy two-solar-mass neutron stars [19][20][21].Moreover, the heaviest well measured neutron star, PSR J0740+6620, was recently also observed by NICER to provide constraints on its radius [22,23].In addition, NICER observed the mass and radius of a typical-mass neutron star, PSR J0030+0451 [24,25].The NICER analyses for both neutron stars by Riley et al. [22,24] and by Miller et al. [23,25] give different mass-radius posteriors, but agree within their uncertainties.The differences in the posteriors are reduced by including realistic assumptions for the EOS, and in this work we explicitly show that in our EDF EOS ensembles the results from both NICER analyses are very similar.In addition to the NICER constraints, we include in our Bayesian inference the tidal deformability information from GW170817 inferred by results from LIGO/Virgo [26].Using the chiral EFT informed priors with the astro posteriors, we provide results for the symmetry energy and neutron star properties.
This paper is organized as follows.In Sec.II we introduce our EOS framework using the new EDFs from Huth et al. [3].These are fit to a range of chiral EFT calculations of neutron matter.Building on this EOS prior, we include constraints at high densities from observations of GW170817 and from NICER using a Bayesian analysis.In Sec.III, we investigate the posterior distributions for the symmetry energy and the slope parameter, as well as for the proton fraction, the speed of sound, and the central density in neutron stars.Moreover, we explore correlations of neutron star radii with the pressure and the speed of sound in neutron stars.Finally, we summarize our results and conclude in Sec.IV.

II. EQUATION-OF-STATE FRAMEWORK
The EOS describes the energy density and pressure of matter for given baryon density, composition, and temperature.Since we focus on cold neutron stars, we consider zero temperature.For a given EOS, the mass and radius of neutron stars follow by solving the Tolman-Oppenheimer-Volkoff (TOV) equations [27,28].Our starting point will be the EOS of homogeneous matter, which we constrain by empirical ranges of the properties of symmetric nuclear matter around saturation density and by neutron matter calculations.Based on each EOS, we calculate consistently the structure of the neutron star crust.
Since neutron stars are extremely neutron rich with proton fractions ≈5%, the most important constraints for the EOS come from neutron matter calculations.In this work, we focus on neutron matter calculations based on chiral EFT interactions, which has the advantage that chiral EFT predicts consistent many-body interactions and enables systematic uncertainty estimates based on the EFT expansion [2,29].Neutron matter has been calculated based on chiral two-and threenucleon interactions using many-body perturbation theory (MBPT) [5,6,10,11,13], quantum Monte Carlo (QMC) methods [9,30], self consistent Green's function (SCGF) methods [7], and coupled cluster (CC) theory [8,12].These calculations are able to include all interactions up to next-to-next-to-next-to-leading order (N 3 LO) [6,11,13] and include uncertainty estimates from the EFT truncation [11][12][13]31].

A. Energy density functionals
To extend the EOS to high density we use nonrelativistic EDFs, which depend on the baryon number density n and proton fraction x of uniform matter.The baryonic energy density ε(n, x) is expressed as where τ n /2m N and τ p /2m N are the neutron and proton kinetic densities, with nucleon mass m N .It was shown that the dependence on isospin asymmetry is to a very good approximation quadratic [32,33], with the dominant nonquadratic contributions stemming from the kinetic densities, so that Eq. ( 1) provides a very good approximation for asymmetric nuclear matter.The func- ).For comparison we also show the unitary gas constraint [34].
tionals f n (n) and f s (n) can be chosen to satisfy the constraints from neutron matter calculations and symmetric nuclear matter properties, respectively.For the interaction density functionals, we take the form introduced recently by Huth et al. [3], (2) where a j , b j are fit parameters and d j = d fm −1−j with parameter d = 3 [3].This corresponds to an expansion of the interaction energy density in powers of the Fermi momentum k F ∼ n 1/3 , and the denominator ensures that the interaction part becomes proportional to n 5/3 at higher densities.Note that without the denominator, the interaction part generally causes the speed of sound to exceed the speed of light beyond some baryon density.For a detailed discussion of these new functionals and the parameter choices, see Ref. [3].

B. Constraints from neutron matter calculations based on chiral effective field theory
For neutron matter constraints we use the MBPT calculations from Ref. [11] based on different chiral N N +3N Hamiltonians, including the Hebeler+ interactions [35], the NNLOsim potentials [36], as well as the N 3 LO 450 MeV and 500 MeV uncertainty bands [11] (using the N N Entem-Machleidt-Nosyk (EMN) interactions [37]).The different neutron matter results and their uncertainties are given by the individual lines shown in Fig. 1.We use the individual lines to fit the a j of the EDF for neutron matter, f n (n) in Eq. ( 2), based on the k F expansion and d = 3.The b j of the corresponding symmetric matter part, f s (n), are determined from empirical properties.We fit to the binding energy E/A(n 0 ) = −16 MeV at saturation density n 0 = 0.16 fm −3 , the incompressibility K = 235 MeV, with K = 9n 2 ∂ 2 (E/A)/∂n 2 (n 0 , x = 1/2), and the skewness Q = −300 MeV, with Q = 27n 3 ∂ 3 (E/A)/∂n 3 (n 0 , x = 1/2).These values are extracted from Skyrme EDFs and constraints for nuclear matter properties [38]; see also Ref. [39].Since neutron star properties are not very sensitive to symmetric nuclear matter, we do not vary all nuclear matter properties, but only explore the most uncertain value of Q in the following, see Sec.III C.
The uncertainties in our EDF EOSs are reflected in the covariance matrix of ⃗ x = (⃗ a, ⃗ b), defined as where x i j is the set of fit parameters (a j , b j ) for the ith individual EOS, ⟨x j ⟩ represents the average of x j , and w i is the weight for each EOS.Since we do not vary the symmetric nuclear matter properties, in this work C jk is a 4 × 4 matrix for the a j from the neutron matter EOSs only.In the initial set given by the 17 neutron matter EOSs, the weights are w i = 1, but when we implement Bayes statistics and inferences, w i < 1.With the average ⟨x j ⟩ and the covariance matrix C jk , a multivariate normal distribution can be used to generate an EOS ensemble based on our EDF EOSs.We note that the statistical uncertainties from this EOS ensemble have of course a prior sensitivity to the initial set of individual EOSs.
The resulting EDF EOS ensemble based on the multivariate normal distribution is shown in Fig. 1 with the 95% credibility region in comparison to the individual EOSs based on MBPT calculations of neutron matter.The ensemble is based on 100,000 EOSs generated using the EDF, Eqs. ( 1) and ( 2), from the average ⟨x j ⟩ and the covariance matrix C jk based on the individual neutron matter MBPT EOSs.The agreement between the band and the individual lines in Fig. 1 indicates that the EDF EOS ensemble employed in this work can generalize chiral EFT results within their uncertainties.Moreover, we compare the EDF EOS ensemble to the unitary gas constraint [34] and observe in Fig. 1 that this is nicely fulfilled by our EOSs.For completeness, we show in Fig. 2 the prior of our EDF ensemble extended up to 8n 0 .

C. Bayesian modelling
We incorporate the astrophysics constraints on the EOS by applying Bayes' theorem, from which the posterior distribution results from the combination of the prior and likelihood, Here, P (⃗ a) represents the EOS prior given by the EDF parameter space obtained from the neutron matter calculations and symmetric nuclear matter properties, and D stands for the astrophysical data so that P (D|⃗ a) is the likelihood or conditional probability of obtaining D for a given EDF with parameter set ⃗ a.In our study, we include the astrophysical observations of GW170817 and results from NICER to constrain the EOS at higher densities.
For the NICER mass-radius constraints for PSR J0030+0451 and PSR J0740+6620 we consider separately either the Amsterdam analysis of Riley et al. [22,24] or the Illinois/Maryland analysis of Miller et al. [23,25].The heaviest neutron star mass of 2.08±0.07,M ⊙ [21] is thus directly implemented through the NICER M -R information of PSR J0740+6620.Folding in the NICER constraints based on our prior leads to the likelihood for the EDF parameters [40,41] where P (M ) is the probability domain of the NICER analysis and P (NICER|M, R(M, ⃗ a)) is the likelihood of each of the two NICER results for a given EDF EOS with parameter set ⃗ a [42].The integral is carried out be discretizing the M − space, summing over all bins which are passed by the M (⃗ a)-R(⃗ a) relation, and weighting those bins with the NICER posterior for each of the sources successively.
In addition to NICER, we use the tidal deformability information from GW170817 inferred by LIGO/Virgo [26], where P (M 1 , M 2 ) is the probability domain from the LIGO/Virgo analysis and P (LIGO|M 1 M 2 Λ 1 Λ 2 ) is the likelihood of the LIGO/Virgo analysis a given EDF EOS with parameter set ⃗ a [42].We assume that the NICER and GW170817 analyses are independent of each other so that, combining both constraints, the likelihood is given by Multiplying the combined likelihood with the prior P (⃗ a) and a normalization constant considering the integral in the denominator, we obtain the posterior distribution P (⃗ a|D) for a given EDF EOS with parameter set ⃗ a.

III. RESULTS
Next we present our results for the properties of neutron stars and the symmetry energy based on the EOS framework developed in the previous section.This combines the information from neutron matter based on chiral EFT interactions, with empirical properties of symmetric nuclear matter, as well as astrophysical constraints from GW170817 and NICER using a family of EDFs for nucleonic matter.Since matter in neutron stars is very neutron-rich, we have focused more on the propagation of the theoretical uncertainties in our knowledge of neutron matter.An advantage of our EOS framework is that we use the same EDF to construct the crust and core EOSs for neutron stars.In the following, we present our results for the neutron star mass and radius, the proton fraction, the speed of sound, and the central density in neutron stars.We also provide results for the symmetry energy and the slope parameter and explore correlations of neutron star radii with the pressure and the speed of sound in neutron stars.

A. Mass-radius relation
The mass and radius of neutron stars are obtained by solving the TOV equations for nonrotating stars.Figure 3 shows the 95% credibility regions for the mass M and radius R generated from the multivariate normal distribution for the EDF EOSs based on an ensemble of ≈ 10 5 EOSs.The top panel shows the prior distribution for the k F expansion using different values of d = 1, 3, 5, 7, and d = ∞.The middle and lower panels show the posterior distribution including astrophysics information from GW170817 and the NICER analysis of Riley et al. [22,24] or the NICER analysis of Miller et al. [23,25], respectively.Our results show that the posterior distributions obtained from the two different NICER analyses are very Posterior distribution based on the top panel prior and including astrophysics information from GW170817 and NICER (dashed 95% and 68% contours from the Amsterdam analysis, Riley et al. [22,24]).Bottom panel: Same as the middle panel but using the NICER results from Miller et al. [23,25].
similar once the nuclear physics information is encoded in the EOS framework.
Regarding the different EDF choices, we find that the d = 3 distribution is similar to the case of d = 5, 7, and d = ∞.However, large d, and in particular d = ∞ allows    for the speed of sound to become acausal, c 2 s > 1 (in units with the speed of light c = 1), as the density increases, which is not the case in either neutron or symmetric matter for d = 3 by construction.In addition, as d = 1 makes the interaction energy density rapidly behave like n 5/3 , the EOS becomes soft at rather low densities compared to the larger d values [3].As a result the 95% credibility regions for mass and radius only extend slightly above 2M ⊙ .Therefore, in the following, we will show results only for the EDF EOSs with d = 3.Before doing so, we also list the radius ranges of typical 1.4M ⊙ and 2M ⊙ neutron stars to show the rather minor sensitivity to the choice of d (see Table I).
In Table I we give the prior and posterior ranges for the radius R 1.4 of a 1.4M ⊙ neutron star at 95% (±2σ) and 68% (±1σ) credibilities as well as the most likely ra- of EOSs have a maximum mass of neutron stars greater than 2.0M ⊙ , while for the posterior distribution, 97% (98%) of EOSs have a maximum mass above 2.0M ⊙ using the NICER analysis of Riley et al. [22] (Miller et al. [23]).In Fig. 4,we show the color-coded prior and posterior distributions for the case of d = 3.In both posterior distributions, the most probable radii for neutron stars between 1.0M ⊙ and 1.8M ⊙ vary only within 0.3 km.Moreover, the mass and radius distribution for M > 2.0M ⊙ is very similar for the prior and the two posteriors, because the astrophysics information mainly removes EOSs that give low maximum mass and small radii.
Table II gives the prior and posterior ranges for the radius R 2.0 of a 2.0M ⊙ neutron star for the EDF EOSs with d = 3.The prior distribution shows a wider radius range because it does not include information of a massive neutron star.Again the two posterior ranges for R 2.0 are very similar and merely shifted by less than 100 m.In the case of d = 3, the maximum mass of neutron stars among the ∼ 10 5 EOS ensemble reaches up to 2.23M ⊙ , while it can go up to 2.32M ⊙ for d = ∞.

B. Symmetry energy and L parameter
We can also extract the symmetry energy S v and the slope parameter L from our calculations.To this end, we calculate the symmetry energy parameters from our EDF EOS ensemble of (a i , b i ) generated by the covariance matrix C jk [see Eq. ( 3)] from where both symmetry parameters are evaluated at n = n 0 and x = 1 2 .The resulting distributions for S v and L follow Gaussian distributions, with the mean and covariance matrix given in the following equation.This is shown in Fig. 5 for the individual MBPT calculations for the different chiral NN+3N Hamiltonians from Ref. [11] as points, where the dashed (solid) line connects the 500 (450) MeV cutoff N 3 LO results.As discussed, our EDF EOS ensembles are built from all the different chiral NN+3N results.The resulting 95% prior and posterior distributions are shown for the EDF EOS ensemble with the k F expansion and d = 3.Note that our 95% distribution contours are obtained by integrating the 2D (twodimensional) domain with 2D probability for S v and L. We find that the prior range for S v and L is narrowed to larger values with the astrophysics constraints included.For both NICER analyses the posteriors are again very similar.
For the prior distribution these are given by (mean values in MeV and convariance matrix in MeV while the posterior distributions for the astrophysical inferences are given for the Riley et al. [22] and Miller et al. [23] analyses, respectively: (10) We observe that the astrophysics constraints move the posterior distributions to larger S v and L values within the prior range.Moreover, all MBPT calculations for the different chiral N N + 3N Hamiltonians are still largely within the posterior range, but some of them only borderline.This points to that astrophysics prefers EOSs on the stiffer part of the neutron matter EOS band based on chiral EFT.This is consistent with the EOS findings in Ref. [43].
In Fig. 5 we also show the GP-B (Gussian process, BUQEYE Collaboration) results at N 3 LO from Ref. [31].Since the GP-B contours are based on the same N 3 LO 500 (450) MeV results [11] included in our analysis, we can trace the difference between the GP-B contours and the N 3 LO points to the evaluation of S v and L for the correlated range of 95% of the calculated saturation density, while our distributions are at a fixed reference saturation density n 0 = 0.16 fm −3 .Since the L parameter scales linearly with the density, this mainly affects the L value, while the range of symmetry energies is broadened due to the additional uncertainty in the calculated saturation density.
Finally, we compare our 95% posterior distributions in Fig. 5 with the recent results from Essick et al. [4], which are however 90% contours.These are based on a different set of chiral N N + 3N calculations and astrophysics constraints through a more general Gaussian process extension to high densities.Nevertheless both contours (at the same reference saturation density n 0 ) are remarkably consistent.

C. Proton fraction
The ground state of neutron star matter is obtained by solving the condition for beta equilibrium, where the neutron, proton, and electron chemical potentials µ n , µ p , and µ e are given by with total energy density ε.Since the core is composed of uniform nuclear matter, Eq. ( 11) is straightforward for a given EDF.For the crust EOS, where matter exists in inhomogeneous form, we employ the liquid drop model (LDM) [44] using the same EDF to construct the EOSs of the inner and outer crust.
In the inner crust, the total energy density including the electron contribution is given by [44] where u is the volume fraction of the nucleus to the Wigner-Seitz cell, n i is the baryon number density of the heavy nucleus, n no is the density of unbound neutrons, x i is the proton fraction in the heavy nucleus, and f i = f (n i , x i ) and f no = f (n no , x no = 0) are the energies per baryon for the heavy nucleus and unbound neutrons, respectively.σ(x i ) is the surface tension at zero temperature as a function of the proton fraction in heavy nuclei, r N the radius of the heavy nucleus, e the electric charge, d the dimension of the nuclear pasta phase, f d (u) the Coulomb shape function corresponding to the nuclear pasta phase, and ε e is the electron energy density.We use the surface tension from [45] σ where σ 0 , α, and q are parameters fit to the calculation of the surface tension.In this work, we use σ 0 = 1.14 MeV fm −2 , α = 3.4, and q = 30, but note that the crust properties depend only weakly on the surface tension parameters, and also the impact of the crust on the investigated neutron star properties is minor.
Based on the viral theorem, the Coulomb energy is approximately twice the nuclear surface energy.Thus, we can combine the surface and Coulomb energy to a single form of energy contribution, which leads to a simpler equation for the energy density [46], where D(u) is a continuous dimension function introduced in Ref. [46].For total baryon density and proton fraction Y p , and thus electron density n e = Y p n, the conditions u, n i , x i , and n no are found by minimizing the total energy density, Eq. ( 15), using the Lagrange multiplier method for the constraints of baryon density and charge neutrality, For an outer crust EOS, which is defined as the region without unbound neutrons, the outside neutron density n no is neglected.Using the LDM construction, the transitions from the outer to inner crust and to the outer core are thus smooth, since the same EDF is employed to construct the entire neutron star EOS. Figure 6 shows the average proton fraction at the central density ⟨Y c p ⟩ based on the EDF EOS ensemble for the k F expansion and d = 3, as well as the variance over the average σ Y c p /⟨Y c p ⟩.The average proton fraction is dominated by the core, but include the details of the crust calculation discussed above.We note that in Fig. 6 (and in Figs. 9 and 10), the mass and radius domain is restricted to the region where the relative probability to the maximum probability ratio P (M, R)/P max ≥ 10 −2 (as in Fig. 4).As expected, the proton fraction increases as the mass increases, and for a given mass, it increases with radius as the EOS becomes stiffer.Our EOS ensemble assumes for the proton fraction that matter is nucleonic, which may not be valid for massive stars.However, for typical 1.4M ⊙ neutron stars, this may not be such a large extrapolation.
In addition, we plot in Fig. 6 the threshold Y p = 1/9 for the direct Urca process, which leads to fast cooling  neutron stars [47].We find that typical neutron stars around 1.4M ⊙ do not exceed this threshold for radii around 12 km, but only in our largest radius configurations.However, based on our results, we expect that massive neutron stars with M > 2.1M ⊙ would cool via the direct Urca process.
Figure 7 shows the total proton fraction Y tot p of the maximum mass star versus the maximum mass.The total proton fraction increases along a band as the maximum mass increases, due to the stiffer EOS. Figure 7 shows results for four different Q values of symmetric nuclear matter, keeping in mind that negative Q values are favored by nuclear masses, ab initio calculations, and astrophysics [4,33,38].With increasing Q, the total proton fraction for a given mass decreases and also the maximum mass increases, as larger Q stiffens the EOS.Naturally, the sensitivity to Q is much less pronounced for typical neutron stars.Figure 8 shows the proton fraction at the central density Y c p versus the radius of a 1.4M ⊙ star, which exhibits a tight correlation and is only very weakly dependent on Q. Larger radii thus have a larger proton fraction.Again we see that radii around 12 km, as expected based on most recent EOS astrophysical inferences [23,39,43,[48][49][50][51][52][53][54], do not cool via the direct Urca threshold.However for larger radii R 1.4 > 12.6 km (for Q = −300 MeV) even typical neutron stars would be fast coolers.

D. Central density and speed of sound
Next, we study the posterior distribution for the central density and the speed of sound in neutron stars.Figure 9 shows the average central density in units of saturation density, ⟨n c /n 0 ⟩, and its variance over the average, σ nc /⟨n c ⟩.The average central density increases with increasing mass, while it decreases as the radius increases for a given mass of neutron star.This results from stiffer EOSs leading to larger radii.In our EDF EOSs, the maximal central density reaches up to ≈ 7n 0 , which is reached for softer EOSs in the most massive neutron stars with smaller radii.Figure 10 shows the speed of sound squared c 2 s = ∂P/∂ε at the central densities in neutron stars.In our EDFs, the speed of sound increases but remains causal and decreases at high density [3].As we see from Fig. 10, the speed of sound is increasing as the mass increases, so in neutron stars most matter is on the part of the EOS that has an increasing c 2 s in our ensemble of EOSs.In Fig. 10, the red dashed line represents c 2 s = 1/3, which shows that even typical 1.4M ⊙ stars exceed the conformal limit, except when they have radii larger than 13 km (see also the middle panel of Fig. 11).Moreover, information on the radii of massive stars with M ≳ 2.0M ⊙ would inform us about c 2 s at the central density (see also Fig. 12).This could be realized with an improved NICER radius measurement [22,23] of the 2.08 ± 0.07M ⊙ pulsar PSR J0740+6620 [21].

E. Correlations
Finally, we study the correlation of neutron star radii with the pressure and the speed of sound.In Ref. [55] it was suggested that the radius of a 1.4M ⊙ neutron star would follow the emprical relation R 1.4 ∼ p 1/4 2n0 , where p 2n0 is the pressure at twice saturation density.In the top panel of Fig. 11 we show that this correlation is indeed fulfilled in our EDF EOS ensemble within a band.For the radius in km and the pressure in MeV fm −3 , we find R 1.4 = 1.279 + 5.063 P 1/4 2n0 for the mean line of the correlation shown in Fig. 11, with a correlation coefficient r xy = 0.985.While the details of this correlation depend on the EOS model, this indicates that astrophysical observations of neutron star radii provide constraints for the pressure at twice saturation density.
The middle panel of Fig. 11 shows the distribution of R 1.4 versus the speed of sound at the central density of neutron stars.Most of the distribution follows a linear trend, but the correlation coefficient r xy = −0.907 is weaker in this case.We also observe that c 2 s at the central density exceeds the conformal limit c 2 s = 1/3 in our EDF EOS ensemble for R 1.4 smaller than 12.8 km The correlation is even weaker at lower densities when comparing R 1.4 with the L parameter in the bottom panel of Fig. 11, which is proportional to the pressure of pure neutron matter at saturation density.This is as expected because the central density of a 1.4M ⊙ neutron star is ∼ 3n 0 .Nevertheless, there is a general trend that R 1.4 increases as L increases.
Figure 12 shows the correlation of the radius of a 2.0M ⊙ neutron star with the speed of sound at the central density.The strong correlation indicates that the radius measurement of massive neutron stars provides constraints for the speed of sound in dense nuclear matter.For the radius in km, we find R 2.0 = 16.493−7.846c 2 s (n c ), with a correlation coefficient r xy = −0.996.Moreover, we find within our EDF EOS ensemble that the speed of sound at the central density of 2.0M ⊙ stars is always greater than the conformal limit.Table III summarizes the fitting functions for neutron star radii for different nuclear matter properties and the corresponding correlation coefficients discussed above.
Figure 13 shows the mass and radius prior when we impose the conformal limit for the speed of sound.The top panel shows the case when the speed of sound continues to increase up to 1/3 and maintains the conformal   As can be seen from our results for S v and L, present astrophysics constraints prefer larger pressures within the prior ranges.To this end, we have also explored correlations of neutron star radii with the pressure and the speed of sound.The radius of 1.4M ⊙ stars was found to correlate well with the pressure at twice saturation density, and R 2.0 was shown to correlate tightly with the speed of sound at the central density.Therefore, precise measurements of R 1.4 provide key information for density regimes at the limits of chiral EFT calculations, and radii of massive neutron stars will help to constrain the behavior of the speed of sound in dense matter.Finally, by constructing EOS ensembles with imposed conformal limit on the speed of sound, we found that a maximum mass of neutron stars M max > 2.1M ⊙ indicates that the speed of sound needs to exceed the conformal limit.

FIG. 2 :
FIG.2: Same as Fig.1but for the 95% credibility region of the EDF ensemble used in this work (based on the kF expansion and d = 3), but extended up to 8n0.

FIG. 3 :
FIG.3: Top panel: 95% credibility regions for the neutron star mass M and radius R generated from the multivariate normal distribution for the EDF EOSs with the kF expansion using different d values.Middle panel: Posterior distribution based on the top panel prior and including astrophysics information from GW170817 and NICER (dashed 95% and 68% contours from the Amsterdam analysis, Riley et al.[22,24]).Bottom panel: Same as the middle panel but using the NICER results from Miller et al.[23,25].

2 FIG. 4 :
FIG.4: Same as Fig.3but showing the color-coded prior and posterior distributions for the case of d = 3.

FIG. 5 :
FIG.5: Plot of the L parameter versus symmertry energy Sv based on the MBPT calculations for different chiral NN+3N Hamiltonians from Ref.[11] [points, see upper left legend; the dashed (solid) line connects the 500 (450) MeV cutoff N3 LO results] and the resulting prior and posterior distributions based on the EDF EOS ensemble of this work for the kF expansion and d = 3 (see lower right legend).Note that the posteriors for both NICER analyses are very similar.For comparison, we also plot the GP-B results at N 3 LO from Ref.[31] and the recent results from Essick et al.[4].All contours are 95%, except for 90% from Essick et al.Note that all results are for a fixed reference saturation density n0 = 0.16 fm −3 , while the GP-B results are for the correlated range of 95% of the calculated saturation density.

1 FIG. 6 :
FIG. 6: Color-coded posterior distribution of the proton fraction at the central density based on the EDF EOS ensemble for the kF expansion and d = 3 (using the NICER Riley et al. analysis).The left and right panels show the average proton fraction at the central density ⟨Y c p ⟩ and the variance over the average σY c p /⟨Y c p ⟩, respectively.The black dashed line denotes the direct Urca threshold of Yp = 1/9.

2 FIG. 7 :
FIG. 7: Color-coded posterior distribution of the total proton fraction Y tot p of the maximum mass star versus the maximum mass for four different Q values (different panels) based on the EDF EOS ensemble for the kF expansion and d = 3 (using the NICER Riley et al. analysis).The black dashed line denotes the direct Urca threshold of Yp = 1/9.

2 FIG. 8 :
FIG.8: Same as Fig.7for the proton fraction at the central density Y c p versus the radius of a 1.4M⊙ neutron star.

FIG. 9 :
FIG.9: Same as Fig.6but for the average central density in units of saturation density ⟨nc/n0⟩ (left panel) and its variance over the average σn c /⟨nc⟩ (righ panel).

11 FIG. 10 :
FIG.10: Same as Fig.9but for the square of the speed of sound c 2 s at the central densities in neutron stars.The red dashed line represents the conformal limit c 2 s = 1/3.

2 FIG. 13 :
FIG. 13: Neutron star mass and radius prior based on the EDF EOS ensemble for the kF expansion and d = 3 when we impose the conformal limit of the speed of sound squared c 2 s ≤ 1/3 (top panel) or a jump in the speed of sound squared to constant c 2 s = 1/3 at n = 2n0 (bottom panel).

of a 1 .
4M ⊙ neutron star R 1.4 are S v = (30.6-33.9)MeV, L = (43.7-70.0)MeV, and R 1.4 = (11.6-13.2) km.Moreover, we have shown that larger and/or heavier neutron stars have a larger proton fraction and are thus more likely to cool rapidly via the direct Urca process.

TABLE II :
Same as Table I but for the radius of a 2.0M⊙ neutron star and for the case of d = 3 only..

TABLE III :
Fitting functions for neutron star radii for different nuclear matter properties and the corresponding correlation coefficients (as discussed in the text).