The nuclear matter density functional under the nucleonic hypothesis

A Bayesian analysis of the possible behaviors of the dense matter equation of state informed through recent LIGO-Virgo as well as NICER measurements reveals that all the present observations are compatible with a fully nucleonic hypothesis for the composition of dense matter, even in the core of the most massive pulsar PSR J0740+6620. Under the hypothesis of a nucleonic composition, we extract the most general behavior of the energy per particle of symmetric matter and density dependence of the symmetry energy, compatible with the astrophysical observations as well as our present knowledge of low energy nuclear physics from effective field theory predictions and experimental nuclear mass data. These results can be used as a null hypothesis to be confronted with future constraints on dense matter to search for possible exotic degrees of freedom.


Introduction
The exceptional progress of multi-messenger astronomy on different astrophysical sources of dense matter has very recently led to quantitative measurements of various properties of neutron stars (NS), such as the correlation between mass and radius (M-R) from X-ray timing with NICER [1][2][3][4] and the tidal polarizability from gravitational wave (GW) LIGO/Virgo data [5][6][7][8][9].These observations, together with the plethora of upcoming data [10], are expected to unveil in the next future exciting open questions such as the structure and degrees of freedom of baryonic matter in extreme conditions, and in particular the presence of phase transitions and the existence of deconfined matter in the core of neutron stars [11].
This direct connection between astrophysical measurements and the microphysics of dense matter is due to the well-known fact that, under the realm of general relativity, there is a one-to-one correspondence between any static observable and the dense matter equation of state (EoS) [12].However, the task is complicated by the fact that there is no ab-initio calculation of ultra-dense matter neither in the hadronic nor in the partonic sectors, and therefore effective models are used.The information about the composition of high density matter is blurred by the uncertainty on the effective energy functional, and similar equations of state can be obtained under different hypotheses on the underlying microphysics [13,14].
A tension was reported between the GW observational data, that tend to favor stiffer EoS, and abinitio nuclear physics calculations, which point towards a slightly softer density dependence [15].Such a tension could in principle suggest the emergence of new degrees of freedom at high density.However, the statistical significance of the dispersion is not sufficient to lead to strong conclusions, and could even be reduced if the new measurement of the neutron-skin thickness of 208 P b by the PREX-II collaboration [16] will confirm a higher value for the skin than previously estimated [17,18].In addition to that, the most recent M-R estimations from the two objects PSR J0740+6620 and PSR J0030+0451 do not report any significant reduction of the NS radius with increasing mass [3,4], in qualitative agreement with the expectations for purely hadronic models for the EoS [19].
For these reasons, the hypothesis of a purely nucleonic composition of the NS cores cannot be ruled out.To identify the observables pointing towards more exotic constituents, it is important to quantitatively evaluate the space of parameters and observables compatible with the nucleonic hypothesis.To this aim, meta-modelling techniques were proposed [19][20][21][22][23][24][25], which allow exploring the complete parameter space of hadronic equations of state, and predicting the astrophysical observables with uncertainties controlled by our present theoretical and experimental knowledge of nuclear physics.This approach can be viewed as a way to transform experimental and observational constraints into nuclear physics empirical parameters to guide the elaboration of phenomenological and microscopic nuclear models, and it can also be used as a null hypothesis to search for exotic degrees of freedom.
In this paper, we address this timely issue by performing a Bayesian statistical analysis of the semiagnostic meta-modelling technique of Refs.[21,24], including both nuclear physics and astrophysical constraints.With respect to previous works by different groups [19][20][21][22][23][24][25], we include the most recent NICER results [4] which give constraints in the density region where many-body perturbation theory (MBPT) cannot be applied, and use a fully unified EoS approach [26] allowing to include in the posterior probabilities the constraints coming from nuclear mass measurements [27].We have not included in the considered constraints the recent skin measurement by PREX II [16] because our model is not presently able to calculate nuclear radii.An extension in this direction would be of foremost interest, and it is left for future work.
The paper is organized as follows.In section 2, we summarize the basic ideas of nucleonic metamodelling developed in Ref. [19,21].We explain the different filters from low energy nuclear physics and astrophysical observations used for the Bayesian analysis in Section 3. The results obtained in the present work are described in Section 4. We make our concluding remarks in Section 5.

Meta-modelling of the EoS
Within the assumption that the core of neutron stars is composed of neutrons, protons, electrons, and muons in weak equilibrium, a prior distribution of viable unified EoS model is generated by Monte-Carlo sampling of a large parameter set of 10 independent, uniformly distributed empirical parameters corresponding to the successive density derivatives at saturation up to order 4 of the uniform matter binding energy in the isoscalar and isovector channels.These parameters characterize the density dependence of the energy in symmetric matter and of the symmetry energy, and their prior distribution is consistent with the present empirical knowledge for a large set of nuclear data [21].They are complemented by 5 additional surface and curvature parameters [28] that are optimized, for each set of uniform matter parameters, to the experimental Atomic Mass Evaluation 2016 (AME2016) nuclear mass table [27].The expression of the surface and curvature energy we employ [29] was optimized on Thomas-Fermi calculations at extreme isospin asymmetries, also subsequently employed in different works on neutron star crust and supernova modelling within the compressible liquid drop approximation [26,28,[30][31][32][33].Two additional parameters rule the density dependence of the effective mass and the effective mass splitting, and an extra parameter enforces the correct behavior at zero density, see Ref. [21] for details.The use of the same functional to describe the inhomogeneous crust [26,28] guarantees a consistent estimation of the crust-core transition and is known to be important for a correct estimation of the NS radius [34].

Bayesian analysis
The posterior distributions of the set X of EoS parameters are conditioned by likelihood models of the different observations and constraints c according to the standard definition: where P (X) is the prior, and N is a normalization factor.The different constraints c k used in the present study are the following: (a) nuclear mass measurements in the AME2016 mass table [27]; (b) the bands of allowed region in symmetric and pure neutron matter produced by many-body perturbation theory (MBPT) calculations from Ref. [35] based on two-and three-nucleon chiral effective field-theory (EFT) interactions at next-to-next-to-next-to leading order (N3LO), that are interpreted as a 90% confidence interval; (c) mass measurement from radio-timing observations of pulsar PSR J0348+0432 [36], M J03 = 2.01 ± 0.04M , where M is the solar mass; (d) constraints on the tidal deformability of the binary NS system associated to the gravitational wave event GW170817, detected by the LIGO/Virgo Collaboration (LVC) [7]; (e) X-ray pulse-profile measurements of PSR J0030+0451's mass, M J00 = 1.44 +0.15 −0.14 M , and radius, R J00 = 13.02+1. 24 −1.06 km from Ref. [2]; (f) the radius measurement with NICER and XMM-Newton data [4] of the PSR J0740+6620 pulsar of mass M J07 = 2.08 ± 0.07M [37], R J07 = 13.7 +2.6 −1.5 km [4].Posterior distributions of different observables Y are calculated by marginalizing over the EoS parameters as: where N = 13 is the number of parameters in the metamodel.Y (X) is the value of any observable Y obtained with the X parameter set, with X min(max) k being the minimum (maximum) value in the prior distribution taken as in Ref. [26].In order to see the impact of different constraints on the nuclear physics informed prior, we consider four distributions, each containing around ∼ 18000 models.They are labeled as follows: 1. Prior: models in this set are required to result in meaningful solutions for the crust, that is, the minimization of the canonical thermodynamic potential at a given baryon density leads to positive gas and cluster densities.In addition, the fit of the surface and curvature parameters {σ 0 , b s , σ 0c , β} to the nuclear masses in the AME2016 table must be convergent.These criteria are characterized by the pass-band filter ω 0 .Given that the mentioned conditions are satisfied, i.e., ω 0 = 1, the probability of each model, associated to a parameter set X, is then quantified by the goodness of the optimal fit, in which the original prior P (X) contains uniformly distributed EoS parameters, and the cost function χ 2 (X) has the following form: The sum in Eq. 4 runs over all the nuclei in the AME2016 [27] mass table ; M AM E and M cl (X) being the experimental and theoretical nuclear masses, respectively, in which the latter is calculated within a compressible liquid drop model (CLDM) approximation using the best-fit surface and curvature parameters for each EOS model; σ n represents the systematic theoretical error; and N dof (= n − 4) is the number of degrees of freedom.The distributions obtained with this prior represent the most general predictions, within a purely nucleonic composition hypothesis, that are compatible with low energy nuclear physics experiments.
2. LD: in this sample the models are selected by the strict filter from the chiral EFT calculation, where the energy per nucleon of symmetric nuclear matter (SNM) and pure neutron matter (PNM) of the model are compared with the corresponding energy bands of Ref. [35], enlarged by 5%.This constraint is applied in the low-density region, from 0.02 fm −3 to 0.2 fm −3 .The posterior probability can be written as: in which ω LD (X) = 1 if the model X is consistent with the EFT bands, and ω LD (X) = 0 otherwise.Implementing this low density (LD) filter amounts to including in the nucleonic hypothesis the information from ab-initio nuclear theory.
3. HD + LVC: the posterior probability of this distribution is written as: Here, ω HD is also a pass-band type filter similar to ω LD in Eq. 5.It only accepts models satisfying all the following conditions: causality, thermodynamic stability, and non-negative symmetry energy at all densities.The second term in Eq. 6, P (J03|X), is the likelihood probability from the mass measurement of PSR J0348+0432 [36], which is M J03 = 2.01±0.04M .This likelihood is defined as the cumulative Gaussian distribution function with a mean value of 2.01 and a standard deviation of 0.04: where M max (X) is the maximum NS mass at equilibrium, determined from the solution of the Tolmann-Oppenheimer-Volkoff (TOV) equations [38,39].
We expect these different conditions not to be selective on the low order EOS parameters, but to constitute stringent constraints on the high density (HD) behavior of the EOS that is essentially governed, within the nucleonic hypothesis, by the third and fourth order effective parameters Q sat , Z sat , Q sym and Z sym [21].
The constraint from the GW170817 event, measured by the LVC, evaluates the weight of a model based on its prediction for the tidal deformability Λ.The likelihood is written as : in which q is the ratio of the lighter component mass m 2 to the heavier component mass m 1 , q = m 2 /m 1 ≤ 1, and P LV C ( Λ(q), q) is the joint posterior distribution of Λ and q taken from Refs.[7,40].In Refs.[7,40], the authors performed a Bayesian inference with four different waveform models.The distribution for Λ and q which we are using for this work is the one obtained with the PhenomPNRT waveform, which is mentioned as their "reference model".The tidal deformability Λ is expressed in the form of the component masses, m 1 and m 2 , and the two corresponding dimensionless tidal deformabilities, Λ 1 and Λ 2 , as : The dimensionless tidal deformability Λ is related to the mass M through the expression: where c, G, R(M ), and k 2 are the speed of light, gravitational constant, NS radius at mass M , and Love number, respectively [41][42][43].In our analysis, q is chosen to be in the one-sided 90% confidence interval obtained in Ref. [7], q ∈ [0.73, 1.00].In Ref. [7], it was shown that the chirp mass M c of the binary NS system associated to the GW170817 event was accurately determined, M c = 1.186 ± 0.001M at the median value with 90% confidence limits.The chirp mass M c can be expressed as a function of m 1 and q as: Since the uncertainty in the chirp mass M c is negligible, for each value of the mass ratio q, we calculate m 1 directly from the median value of M through Eq. ( 11).
4. All: Including the three constraints mentioned above together with the likelihood from the joint mass-radius distributions of the two NICER measurements from Refs.[2,4], the posterior probability for the final distribution is written as: The NICER likelihood probability is given by: where p N ICER1 (M, R) is the two-dimensional probability distribution of mass and radius for the pulsar PSR J0030+0451 obtained using the waveform model with three uniform oval spots by Miller et al. in [2]; and p N ICER2 (M, R) is the probability distribution for PSR J0740+6620 using NICER and XMM-Newton data by Miller et al. [4].The intervals of M 1 and M 2 are chosen sufficiently large so that they cover most of the associated joint mass-radius distributions, To insure that the differences in the posterior distributions are induced by the impact of the different constraints, care is taken to have a comparable statistics from the four distributions, for each plot shown in this paper.Moreover, for all shown observables we have checked that an increase in statistics does not affect the results, within the precision chosen for the numerical values given in this paper.

Empirical parameters
In Figures 1 and 2 we show the probability density distributions (PDFs) of isoscalar and isovector empirical parameters of order N < 4, respectively.As we have described previously, the distributions labeled as "Prior" in these figures are not flat, but they carry the information from the experimental nuclear mass measurement.For example, E sat , the energy per particle in SNM at saturation, already has a peaked shape (see Figure 1 (a)) because of this reason.From the HD+LVC distribution in these figures, we can see that the astrophysical constraints on NS mass and tidal deformability have almost no effect on low-order parameters.The impact of the chiral EFT filter on the isoscalar parameters of order N < 3, i.e.E sat , K sat , along with n sat is not prominent.This knowledge can be further reinforced by looking at Table 1, where the LD filter hardly improves the constraints on the aforementioned isoscalar parameters.It can be explained by the fact that the prior intervals of the empirical parameters are chosen based on the current knowledge provided by nuclear physics, in which the deviations of E sat , n sat , and K sat are already relatively small.Table 1: Medians and 68% confidence limits of EoS empirical parameters of order N < 4 in the four distributions.
Prior −16.25  Unlike the lower order parameters in the isoscalar sector, the isovector counterparts are quite poorly determined by nuclear physics experiments.As a result, once the constraint from the chiral EFT calculation is included, E sym , L sym and K sym are strongly affected (see Figure 2 and Table 1).Interestingly, the LD filter also has a non-negligible impact on the high order parameters Q sat and Q sym .This is because  the chiral EFT calculation gives very precise predictions at very low densities, far from nuclear saturation.In this region, the high-order parameters have a non-negligible contribution to the nuclear matter energy.It was shown by Refs.[28,44] that constraining the EoS at very low densities n ∼ 0.02 − 0.1 fm −3 is crucial in studying the crust-core transition.
As one may expect, the constraints from NS observables (HD+LVC) play an important role on highorder parameters, such as Q sat and Q sym , as well as on the poorly constrained isovector compressibility K sym .One can observe that for these parameters, higher values of the chosen intervals are preferred in the nucleonic hypothesis, having low preference on the softer EoSs.However, note that this is the net effect of both the radio mass and GW180817 measurements.We have checked that without the constraint on the tidal deformability, the resulted nuclear matter energies are even higher, which means that the constraint from GW170817 softens the EoS.
As discussed in details in Ref. [21], the density behavior of realistic functionals can be accurately reproduced up to the central density of massive neutron stars by a Taylor expansion truncated at fourth order, but because of the truncation the parameters of order N ≥ 3 have to be considered as effective parameters that govern the high density behavior of the EOS, and do not need to be equal to the corresponding density derivatives at saturation.On the other hand, in the sub-saturation regime, the deviations from the Taylor expansion are accounted for by the low density corrective term that imposes the correct zero density limit [21].These two effects being completely independent, the meaning of the third and fourth order parameters as explored by the EFT calculation and the astrophysical observations is not the same, and we can expect that low and high density constraints might point to different values for those parameters.Comparing the dashed and dashed-dotted lines in Figure 1 we can see that indeed low-density constraints impose lower values of Q sat with respect to high-density ones.This means that low energy experiments aimed at a better measurement of Q sat will not improve our empirical knowledge of the high density EOS.Interestingly, the same is not true for Q sym , for which the dotted and dashdotted distributions closely overlap.Even if the present constraints are quite loose, it appears that the skewness of the symmetry energy at saturation Q sym gives a fair description of the behavior of the EOS at high density, while a deviation is observed at the level of the compressibility K sym .We do not include the results for the fourth-order parameters Z sat , Z sym , because they have very large uncertainties, and very little impact from the different constraints.Furthermore, we will see later on that they have almost no correlations to other parameters as well as observables.
In Figure 3, we plot the bands for SNM energy per nucleon and symmetry energy at 50% and 90% confidence intervals for the four posterior distributions as explained in the previous section.The impact of LD and HD+LVC filters can be observed by looking at panels (b) and (c) of Figure 3, respectively.Their effects are appreciated at different density regime, as is also evident from the analysis done in Figures 1 and 2

Properties of NS crust
Table 2: Estimations of NS crustal properties in four distributions.The results are presented with medians and 68% confidence limits.In our calculation, the transition from the solid crust to the liquid outer core is determined by comparing the corresponding energy density of clusterized matter to that of homogeneous matter at β-equilibrium with the metamodel.For the crust part, metamodel is extended introducing surface parameters within the compressible liquid drop model (CLDM) approach [26].The precision in the prediction of the crust-core transition point is crucial in estimating crustal observables, such as crustal mass, thickness, and moment of inertia.These quantities are in particular thought to have influence on the origin of the pulsar glitches [45].In the literature, there are various works devoting to determining the crust-core transition density n CC with different many-body methods and nuclear functionals, spanning a large range of values, such as n CC = 0.0548 fm −3 in [46] obtained using Thomas-Fermi calculations for the NL3 functional, or n CC = 0.081 fm −3 in [47] within the full fourth-order extended Thomas-Fermi approach for the BSk24 functional.For this reason, an estimation for the uncertainties of the crustal properties with Bayesian tools using both the current nuclear physics and astrophysical data, provided by LVC and NICER, are of great importance.
In Figure 4, we display the joint distributions of the crust-core transition density n CC and pressure P CC .The chiral EFT calculation plays an important role in the determination of the crust-core transition point, which is evident from the LD distribution in Figure 4 (b).One can observe that, the chiral EFT filter puts stringent limits on both the crust core transition density n CC and pressure P CC ; very high and very low values of n CC and P CC get eliminated.In Figure 4 (c) for the HD+LVC distribution, the most noticeable fact is the suppression of models with high transition pressures.However, the probability densities of these models are tiny, and they are outside of the 95% contour in the prior distribution (see Figure 4 (a)).Moreover, they are associated to models violating at least one of the following conditions required in the HD+LVC posterior: causality, thermodynamics stability, or nonnegative symmetry energy.In other words, the astrophysical constraints on NS maximum mass and tidal deformability have very little effects on the crust-core transition.This is consistent with our observations for the nuclear matter energy in Figure 3   The crust-core transition point determines astrophysical observables, such as crust thickness, or moment of inertia [24].In this study, we have chosen the crust thickness to be the demonstrative quantity.Figure 5 presents the PDFs of NS crust thicknesses for 1.4M and 2.0M NSs.In both cases, the uncertainties in the LD distributions are narrowed down compared to the prior, while the effect in the HD+LVC distribution is only marginal.This agrees with our conclusions for the crust-core transition point, that is, the role of the chiral EFT filter is more dominant in the determination of crustal properties.When all constraints are taken into account, crust thicknesses of both 1.4M and 2.0M NSs are known with the relative uncertainties up to 10%.For a quantitative estimation of the effects of different filters, in Table 2 we present crust-core transition density n CC and pressure P CC , along with crustal thickness of 1.4M and 2.0M NSs accompanied by errors on them at 68% confidence interval.Quite conclusively one can comment that the primary effect comes from the LD chiral EFT filter, which also puts stringent constraints when all the filters are combined together, denoted as "All".

NS equation of state
Unlike the crustal properties, HD+LVC filter is expected to put tighter bounds on global NS properties, which is governed chiefly by the high-density part of the EoS.The effects of different filters on the EoS are shown in Figure 6.The light (dark) orange band indicates 90% (50%) confidence interval.For comparison, we also display the result inferred from the gravitational wave data GW170817 by LVC at 90% level in dashed blue lines [6].We have also used the same units for mass-density g cm −3 as in Ref. [6] for the same reason.In this unit the saturation density n sat is denoted as ρ sat ( 2.8 × 10 14 g cm −3 ).In Ref. [6] , Abbott et al. have sampled their EoSs at high density using the spectral parametrization [48].These EoSs are then matched with SLy EoS [49] at around ∼ ρ sat /2.Incidentally , the authors also put some prior criteria similar to our analysis, which are causality, thermodynamic stability, and consistency of NS maximum mass with the observation.For the last condition, they put a sharp limit (M max ≥ 1.97M ) instead of using a likelihood probability like the one used in our analysis (see Eq. 7).However, we have verified that the difference in the maximum NS treatment does not lead to sizable deviation in the final results.In Figure 6 (a), we can see that our prior distribution perfectly covers the whole posterior band given by GW170817 event [6] with good agreement.In our case, the prior distribution carries information from nuclear physics experiments and theoretical calculations via the chosen prior intervals of empirical parameters as well as the mass fit.This is why the EoS in our prior distribution at low densities is relatively narrow compared to other analyses.Note that the uncertainty below ρ sat , appears to be large due to the visual effect of the logarithmic scale in the pressure.Once the chiral EFT filter is applied, this uncertainty is vastly reduced (see Figure 6 (b)), resulting in a very well-constrained band and excellently compatible with the posterior constrained by GW170817 data [6] .Contrarily, the behavior of the EoS at supra-saturation densities is not constrained by the chiral EFT filter.As a result, a larger dispersion is observed at high densities.This dispersion is not as important as in fully agnostic studies [50] because of the nucleonic hypothesis that imposes an analytic behavior of the EoS at all densities.This strong hypothesis can be challenged by the astrophysical measurements, and any inconsistency with the observations will reveal the presence of exotic degrees of freedom.By incorporating the pass-band filter ω HD as well as the condition on the NS maximum mass in Figure 6 (c), the deviation in the lower limit of the pressure at density ρ 10 15 g cm −3 observed in the prior, is eliminated.In particular, the constraint on the NS maximum mass sets a stringent limit on the lower bound of the pressure, and posterior EoS is shifted significantly towards higher values of pressure.Conversely, the constraint from LVC prefers softer EoS, hence setting the limit on the upper bound of the pressure band.In Figure 6 (d), when all constraints are combined together, we obtain as expected a narrower band for the EoS than the one obtained exclusively from GW170817 data [6].In addition, we observe that our EoS is lightly stiffer than the one of Ref. [6] at around 2 − 3ρ sat .The small width of the EoS and its stiffness are assigned to the semi-agnostic hadronic prior, which represents current nuclear physics knowledge.Nevertheless, the overall agreement is excellent.Thus, it indicates the compatibility of the nucleonic EoS with the gravitational wave GW170817 data.
Comparing the "HD+LVC" and "All" distributions in Figure 6 (c) and (d) , it can be observed that the inclusion of the new NICER measurement does not show any significant impact on the EoS.Similar conclusions have been drawn in other studies in the literature.In particular, Pang et al. [51] has carried out a Bayesian analysis using the data from Riley et al. [3] and Miller et al. [4].In both cases , they found that the effect of the constraint from the radius measurement of PSR J0740+6620 only marginally impact the EoS.In Ref. [52] , Raaijmakers et al. performed the Bayesian inference with two EoS parametrizations: a piece-wise polytropic (PP) model and a speed-of-sound (CS) model drawing similar conclusions.For the constraint on PSR J0740+6620, they employed the data from Riley et al. [3], in which the error bar of the radius is smaller than that obtained in Miller et al. [4].They concluded that for the PP models, the impact on the EoS mainly comes from the high mass value of PSR J0740+6620 because their prior distribution in that mass range is within the 68% level of the radius measurement (see Figure 4 in [52]).In Figure 7, we plot the velocity of sound in medium as a function of mass density ρ obtained with four different filters at the 50% and 90% confidence intervals, together with the behavior of some selected models [53][54][55][56].One can observe that for all the filters the most probable equations of state remain causal up to very high densities (∼ 6ρ sat ), even though we do not put this requirement explicitly in our "Prior" and the "LD" filters in Figure 7 (a) and (b), respectively.As expected, the behavior of the sound speed is globally structureless.However, we can surprisingly see a trend for a peaked structure, which is typically presented in the literature as a signature of a transition to exotic matter.This peak may arise from the shoulder observed in Figure 6 above, which is due to the combined constraints of a relatively soft EoS at low density, and the requirement of reaching the maximal mass.These conditions lead to a peak in the global distribution, note however, not for all models individually (see lines in Figure 7 (d)).A very small fraction of non-causal models is present due to the fact that we plot the EoS only up to densities where the nucleon sound velocity is in the interval between 0 and 1. Residual non-causalities (not visible within 90% confidence interval of Fig. 7) originate from the additional lepton contribution in beta-equilibrated matter.In Figure 8 we plot for different filters three shaded regions (from light to dark) sequentially containing 99%, 95%, and 68% confidence intervals for two dimensional distribution for mass and radius of NSs.The two black contour lines at low mass and high mass respectively indicate 68% of the mass-radius distributions for PSR J0030+0451 [2] and PSR J0740+6620 [4].One can observe in Figure 8 (a) that our prior is already quite compatible with both the recent NICER observations [2,4].This explains why the effect of the constraints from NICER is globally small in all our distributions.Moreover, since in Figure 8 (c) we already include the constraint from the radio mass measurement of the high-mass pulsar PSR J0348+0432 beforehand, the impact from the mass of PSR J0740+6620 is obscured in Figure 8  (d).Additionally, the large uncertainty in the new radius measurement does not help constrain further the EoS.The compatibility of NICER measurement and our distributions implies that a nucleonic EoS is flexible enough to reproduce those dense-matter observations.In Ref. [51], Pang et al. computed the Bayes factor to study the possibility of having a strong first-order phase transition from nuclear to quark matter in NS.If the data from Miller et al. [4] is used, the Bayes factor changes from 0.265 to 0.205.Even though the effect from PSR J0740+6620 is not significant, a decrease in the Bayes factor points to the fact that a first-order phase transition to quark matter is disfavored.Similarly, Legred et al. [57] found that the Bayes factor for EoSs having one stable branch against those with at least one disconnected hybrid star branch is 0.156 (0.220) with (without) the PSR J0740+6620 measurement.Both these studies censure the possibility of a strong phase transition and support the suitability of the hadronic EoS with respect to NS observables, which is in line with our present analysis.

NS observables
Figure 9 displays the marginalized distributions of NS radii, R 1.4 and R 2.0 , of the canonical mass 1.4M (panel (a)) and the typical high mass 2.0M (panel (b)) , respectively.The dashed blue lines represent the PDFs obtained when chiral EFT (LD) filter is applied.We can see from the figure that this filter puts a constraint on the upper bound of the distributions.It rejects models with R 1.4 13.6 km and R 2.0 14.0 km.In the HD+LVC distribution for 1.4M NS, the constraint from GW170817 softens the EoS, hence constraining the upper bound of R 1.4 , while the requirement on the NS maximum mass filters out very soft EoSs, which put a limit to the lower bound of R 1.4 .As a result, these two competing effects provide us with a relatively narrow range on the radius.In particular, R 1.4 ∈ [11.8, 14.0] km (see red dashed-dotted line in panel (a)).In the case of R 2.0 , the constraint from radio mass measurement of PSR J0348+0432 becomes redundant because all distributions must support 2.0M NS resulting in no effect shown in the lower value of R 2.0 .Therefore, in the HD+LVC distribution of R 2.0 , the constraint only comes from the LVC measurement.Furthermore, this figure also tells us that the impacts on R 2.0 from the gravitational signal GW170817 and chiral EFT calculation are very similar, even though they affect two different regions of the EoS.Specifically, the former controls the EoS in the NS core, hence the core radius, while the latter dominates the crust EoS, hence the crust thickness.The prediction in the form of median and 68% credible limits for R 1.4 (R 2.0 ) when all constraints are applied together is 12.78 +0.30 −0.29 (12.96 +0.38 −0.37 ) km.In Miller et al. [4], the authors employed three EoS models, namely Gaussian, spectral, and PP.The values of R 1.4 for these three models are respectively 12.63 +0. 48 −0.46 km, 12.30 +0.54 −0.51 km, and 12.56 +0.45 −0.40 km at 68% confidence limit.Despite the difference in EoS sampling methods, these results are in excellent agreement with the results obtained in the present work.Using also the likelihood from PREX-II measurement of R 208 skin [58], Ref. [59] obtains R 1.4 = 12.61 +0.36 −0.41 km, which is also consistent with our prediction.100 500 900 1300  The dimensionless tidal deformability Λ in Eq. 10 suggests a relation between Λ and R for a NS of given mass M .However, this relationship is not straightforward due to the complex radius dependence of the tidal Love number k 2 [41][42][43].The relation between R and Λ , particularly for the mass M = 1.4M has been investigated in several works [60][61][62][63].Interestingly, Figure 10 shows that the distributions of Λ 1.4 and Λ 2.0 behave in accordance with the corresponding radius distributions in Figure 9.This may indicate a strong positive correlation between these two quantities, which will be discussed later on.In addition, we estimate the 90% confidence boundaries of Λ 1.4 (Λ 2.0 ) to be Λ 1.4 ∈ [463, 757] (Λ 2.0 ∈ [43,94]).This prediction of Λ 1.4 agrees excellently with the upper bound extracted from GW170817 signal in Ref. [5] using TaylorF2 model, that is, Λ 1.4 ≤ 800.The limit in Λ 1.4 has been improved in Ref. [6], in which the more realistic waveform PhenomPNRT was employed, and they obtained Λ 1.4 ∈ [70, 580] at 90% confidence level for the EoS-insensitive analysis [5,6].Our distribution is still compatible with this result, but it suggests a slightly too stiff EoS in the nucleonic hypothesis.

Composition
The determination of the proton fraction is crucial for studying NS cooling.The most efficient cooling mechanism of NS is through the direct Urca (dUrca) neutrino emission process.This process is described by the successive following reactions: where l = {e − , µ − }.From the momentum and charge conservations, one can derive the expression for the threshold, below which the dUrca process is forbidden: where x ep (= x e /x p ) is the ratio between electron and proton fraction.Values of x DU can vary in the range from x DU 1/9 in the case of no muons (x ep = 1) to x DU 0.148 at the limit of massless muons (x ep = 0.5) [30,64].DU .We find that this quantity is independent of the constraint used.Furthermore, x mp DU only depends weakly on NS mass, x mp DU 0.134 (0.138) for M = 1.4 (2.0) M .For both masses, the distributions of x p extend to higher values than the corresponding threshold x mp DU .Therefore, it is possible for the dUrca process to operate even in NS of mass 1.4M .Nevertheless, this fast cooling channel is more likely to happen in heavier NSs due to the higher median and deviation of the x p distribution.By integrating the PDF to find the area under the curve for x p ≥ x mp DU , we estimate the possibility for the dUrca process in NS of mass 1.4M (2.0M ) to be approximately 26% (72%).For a more definitive evaluation, the predictions of NS central proton fractions along with the radius and tidal deformability for NSs of mass 1.4M and 2.0M at 68% confidence limit are listed in Table 3 .Studying correlations among parameters and observables reveals a great deal about the many facets of multi-parametric model calculations [65].The most frequently employed tool for this purpose is the linear Pearson correlation, which is defined for two quantities x and y (x, y can be parameters of the model or any observable calculated from it) as, where cov(x, y) is the covariance between x and y , and σ x (σ y ) is the standard deviation on x (y). Figure 12 displays the Pearson correlation coefficients among all bulk, surface, and curvature parameters in the case where all constraints are applied.Since the bulk parameters are initially by construction uncorrelated in the flat prior distribution, we can easily assign the induced correlations to the different filters employed.It is shown in the figure that there is a perfect negative correlation between the surface tension of symmetric matter σ 0 and the saturation energy E sat , with corr(σ 0 , E sat ) = −1.A similar result was found in Ref. [26].The parameters associated to the curvature (σ 0c and β), on the other hand, exhibit strong positive correlations with E sat .These correlations appear due to the fit of the surface and curvature parameters to the experimental nuclear mass table.In addition, if the prior is only constrained  by the experimental masses of nuclei, we also find a strong correlation between b s and E sym , which are the two main parameters governing the energy of asymmetric nuclear matter.However, once the filter from chiral EFT calculation is applied, E sym is tightly constrained, and hence the correlation get blurred.Similar to Refs.[19,26], no significant correlations are found to be induced by the astrophysical constraints.The correlations among the bulk parameters shown in Figure 12 are resulted from the chiral EFT constraint.In particular, the symmetry energy E sym has a moderate (anti)correlation with (E sat ) n sat .Stronger correlations are found among the isovector parameters, which are corr(E sym , L sym ) = 0.67 and corr(L sym , K sym ) = 0.67.The former is found in several works (see Refs. [23,26,66] and references therein for a review), and the latter is also studied in Refs.[67][68][69][70][71][72][73].Slight correlations between high-order parameters, K sat − Q sat and K sym − Q sym , are also induced due to the narrow EFT energy bands at very low densities.

Conclusions
To conclude, we have jointly analyzed different constraints on the nuclear matter EoS coming from nuclear experiments, ab-initio nuclear theory, and several new astrophysical observational data, including the very recent simultaneous observation of mass-radius of PSR J0348+0432 and PSR J0740+6620 from NICER collaboration as well as LIGO-Virgo observation of tidal deformability in GW170817 event.Imposing all these different constraints in a Bayesian framework, we have challenged the hypothesis of a fully analytical (continuous and derivable at all orders) EoS, as obtained in the case where dense baryonic matter is purely constituted of neutrons and protons without any phase transition or exotic degrees of freedom.
Particularly, we have observed that if we have a nuclear physics informed prior including the binding energy data of the whole nuclear chart and chiral EFT constraints on low density SNM and PNM, the posterior for mass-radius of NSs are already in line with NICER observations.Contrarily, bounds on high density matter from radio astronomy observation of NS of 2 solar mass and GW170817 data on tidal deformability is reasonably well appreciated.With the present knowledge on astrophysical observations, we predict that the direct Urca cooling is possible with non-negligible probability (27%) even in a NS with mass as low as 1.4M , which increases much further ∼ 72% for a NS of 2.0M .This might also be very crucial to (in)validate the nucleonic hypothesis of high density matter.As all current data on astrophysical observations comply with the nucleonic hypothesis within our metamodel approach, we need much more stringent constraints from observations to conclusively establish (reject) the presence of exotic degrees of freedom in high density matter.

Figure 1 :
Figure1: Probability density distributions of isoscalar empirical parameters for the prior distribution informed by experimental nuclear masses (black dotted line) and for posteriors of models passing through the low-density (chiral EFT) constraint (blue dashed line), high-density constraints (causality, stability, e sym ≥ 0, maximum NS mass, and tidal deformability) (red dash-dotted line), and all constraints combined (green shaded region).See texts for details.

Figure 2 :
Figure 2: Same as Figure 1 but for isovector empirical parameters.

Figure 3 :
Figure 3: 50% (darker color) and 90% (lighter color) confidence intervals of energy per nucleon of symmetric nuclear matter (e SN M , color green) and symmetry energy (e sym , color orange) as a function of density n.

Figure 4 :
Figure 4: Joint probability density plots of crust-core transition density n CC and pressure P CC .The dashed black contours in each panel indicate the 68%, 95%, and 99% confidence regions.

Figure 6 :
Figure 6: 50% (dark orange) and 90% (light orange) confidence intervals of pressure P as a function of mass density ρ in comparison with the 90% confidence interval of the posterior obtained in Abbott et al. 2018 [6] (blue dashed lines).See text for details.

4. 4
Speed of sound in medium

Figure 7 :
Figure 7: 50% (dark green) and 90% (light green) confidence intervals of sound speed cs c 2 as a function of mass density ρ.Curves in panel (d) show the sound speed of some selected models [53-56] up to the central density corresponding to the maximum mass.See text for details.

Figure 8 :
Figure 8: 2D density plots of NS mass M as a function of radius R in comparison with two NICER measurements at 68% (black contours).The three shaded regions in each panel contain 68%, 95%, and 99% of the distribution.See text for details.

Figure 11 :
Figure 11: Probability density distributions of central proton fractions of NS at M=1.4 M and M=2.0 M .The arrow in each panel indicates the most probable value of x DU .Panel (a): x mp DU 0.134.Panel (b): x mp DU 0.138.At each value of mass, value of x mp DU are very similar in four distributions.See text for details.

Figure 11
Figure 11  shows the PDFs of proton fractions calculated at the center of NS with M = 1.4M and M = 2.0M .The black arrow in each panel indicates the most probable value of x DU , calculated for the central density, denoted as x mp DU .We find that this quantity is independent of the constraint used.Furthermore, x mp DU only depends weakly on NS mass, x mp -0.04 0.66 0.88 0.92 0.06 0.81 0.10 0.03 0.65 0.81 0.88 -0.05 0.83 0.98 -0.02 -0.07 0.14 0.25 0.22 0.89 -0.08 -0.05 -0.12

Figure 13 :
Figure 13: Pearson correlation matrix among some observables in the case all filters are applied.

Figure 14 :
Figure14: Pearson correlation coefficients between some observables with the empirical and surface parameters in the case all filters are applied.

Table 3 :
Medians and 68% confidence intervals of NS radii, dimensionless tidal deformabilities, and central proton fractions at M = 1.4M and M = 2.0M .
Figure 12: Pearson correlation matrix among bulk and surface empirical parameters in the case all filters are applied.