Bayesian Survey of the Dense Matter Equation of State Built upon Skyrme Effective Interactions

The nonrelativistic model of nuclear matter (NM) with zero-range Skyrme interactions is employed within a Bayesian approach in order to study the behavior of the neutron star (NS) equation of state (EOS). A minimal number of constraints from nuclear physics and ab initio calculations of pure neutron matter (PNM) are imposed together with causality and a lower limit on the maximum mass of an NS to all our models. Our key result is that accounting for correlations among the values that the energy per neutron in PNM takes at various densities and that are typically disregarded efficiently constrains the behavior of the EOS at high densities. A series of global NS properties, e.g., maximum mass, central density of the maximum mass configuration, minimum NS mass that allows for direct URCA, and radii of intermediate and massive NSs, appear to be correlated with the value of effective neutron mass in PNM at 0.16 fm−3. Together with similar studies in the literature our work contributes to a better understanding of the NS EOS as well as its link with the properties of dense NM.


INTRODUCTION
Hydrostatic equilibrium states of cold mature neutron stars (NSs) are determined by a one parameter equation of state (EOS) that relates pressure (P) to energy density (e).This means that any astrophysical measurement of a NS structurerelated parameter, e.g., mass, radius, tidal deformability, moment of inertia, oscillation modes frequencies, probes in one way or another the EOS.Inference of the NS EOS from astrophysical observations is, nevertheless, not an easy task.First of all, only a few such measurements are available so far.Next, depending on the NS mass and the unknown EOS, different density domains are explored throughout different stars.Then, every global parameter depends in a specific and convoluted way on the behavior the EOS has over a certain density domain.Finally, over limited density domains, different particle compositions may provide similar P(e) dependencies.
Indeed, in spite of considerable effort, the NS EOS is still largely unknown (Oertel et al. 2017;Margueron et al. 2018; mikhail.beznogov@nipne.roaraduta@nipne.roBurgio et al. 2021).In addition to limited astrophysics constraints, this situation is also due to limited constraints from nuclear physics.Information from nuclear physics experiments mainly concerns a narrow density domain around nuclear saturation density n sat ≈ 0.16 fm −3 ≈ 2.7 × 10 14 g/cm 3  and matter with low isospin asymmetry, that is with similar numbers of neutrons and protons, and is customarily translated into values of nuclear empirical parameters (NEPs).Complementary information on the density behavior of matter with extreme isospin asymmetries is provided by ab initio calculations that employ accurate nucleon-nucleon potentials and sophisticated many-body techniques.Nevertheless, so far these calculations are successful in efficiently constraining only the behavior of pure neutron matter (PNM) up to densities of the order of n sat (Drischler et al. 2016(Drischler et al. , 2019)).
From the astrophysical point of view, a tremendous progress was done in the last decade.As such, NS masses in the range 0.8 ≲ M/M ⊙ ≲ 2.1 (Demorest et al. 2010;Antoniadis et al. 2013;Arzoumanian et al. 2018;Cromartie et al. 2020;Fonseca et al. 2021;Doroshenko et al. 2022;Brodie & Haber 2023) have been measured with unprecedented precision.Measurements of massive NSs are informative of the NS EOS stiffness and of utmost relevance for the onset of non-nucleonic particle degrees of freedom (Oertel et al. 2017;Sedrakian et al. 2023).At the other extremity, NS masses lower than 1 M ⊙ challenge our understanding of NS formation in core-collapse supernovae (Suwa et al. 2018).This is, in particular, the case of the NS in the supernova remnant (SNR) HESS J1731-347, whose mass was recently estimated to be M = 0.77 +0.20  −0.17 M ⊙ (or M = 0.83 +0.17 −0.13 M ⊙ ) (Doroshenko et al. 2022;Brodie & Haber 2023).Moreover, its very low radius R = 10.4 +0.86  −0.78 km (or R = 11.25 +0.53 −0.37 km) (Doroshenko et al. 2022;Brodie & Haber 2023) triggers questions about the composition of low mass NS cores or the behavior of NS EOS at low densities.In relation with the first aspect, the occurrence of quark matter (Di Clemente et al. 2022) and ∆-resonances (Li & Sedrakian 2023) has been speculated.According to Li & Sedrakian (2023), the particularly compact configuration could be alternatively explained by low values of the slope of the symmetry energy.Detection of gravitational waves (GW) emitted from the coalescence of two NS in the GW170817 (Abbott et al. 2017;Abbott et al. 2017;Abbott et al. 2019) event or a NS-black hole binary in the GW190425 (Abbott et al. 2020) have allowed the first inferences of NS tidal deformabilities.The relatively low masses of the NSs involved in the GW170817 event, with an estimated total mass M T = 2.73 +0.04  −0.01 M ⊙ and a mass ratio 0.73 ≤ q = M 2 /M 1 ≤ 1, makes that the GW signal mainly probed the intermediate density domain of the EOS.Observations of associated electromagnetic near-infrared event AT2017gfo (Tanvir et al. 2017), gamma ray burst GRB170817A (von Kienlin et al. 2017;Savchenko et al. 2017a,b) and optical (Soares-Santos et al. 2017) counterparts of GW170817 have been, in turn, used to provide an estimate of the maximum mass of non-rotating cold configurations, M * G ≈ 2.3 M ⊙ (Margalit & Metzger 2017;Rezzolla et al. 2018;Ruiz et al. 2018;Shibata et al. 2019).Finally, pulse-profile modeling of X-ray data obtained by NASA's Neutron Star Interior Composition Explorer (NICER) mission allowed simultaneous masses and radii estimates for the canonical mass pulsar PSR J0030+0451 (Miller et al. 2019b;Riley et al. 2019) and the massive millisecond pulsar J0740+6620 (Miller et al. 2021;Riley et al. 2021).In spite of the still large uncertainties on radii, of the order of a few km, these measurements can rule out a number of phenomenological EOSs.
Constraints from multimessenger astrophysics of NSs on the dense matter EOS are best exploited by statistical analyses that allow for systematic exploration of the broad parameter spaces.Many such studies have been performed in the last few years.Most of them employ parametrized EOSs models, e.g., piece-wise polytropes, spectral parameterizations and models of the speed of sound, and focus on the sensitivity of posterior distributions on the prior distributions; narrowing down of the parameter space upon progressive incorporation of constraints; dependence of the posterior distributions on the EOS model; potential benefit from measurements of extra NS observables like the moment of inertia or the gravitational binding energy; the advantage of accounting for the full posterior probability distribution of measurements rather than "cuts"; the role that different kind of measurements play on different density ranges; compatibility between various constraints.Examples in this sense are offered by Fasano et al. (2019); Raaijmakers et al. (2020); Miller et al. (2019a); Dietrich et al. (2020); Essick et al. (2020); Miller et al. (2021); Raaijmakers et al. (2021).
Motivated by the urge for accounting for properties of nuclear matter (NM), an increasing number of works started to use phenomenological energy density functionals (Lim & Holt 2019;Zhang & Li 2019;Ferreira et al. 2020;Güven et al. 2020;Ferreira & Providência 2021;Malik et al. 2022;Ghosh et al. 2022;Patra et al. 2022;Beznogov & Raduta 2023;Malik et al. 2023;Zhou et al. 2023;Huang et al. 2023).The obvious advantage of using phenomenological models consists in having access to matter composition.Models based on the Taylor expansion of the energy per nucleon in terms of deviation from saturation and isospin symmetry (Zhang & Li 2019;Ferreira et al. 2020;Ferreira & Providência 2021;Patra et al. 2022), agnostic metamodels (Güven et al. 2020), mean field models with Skyrme effective interactions (Zhou et al. 2023) as well as models that employ a χEFT inspired expansion in terms of neutron and proton Fermi momenta (Lim & Holt 2019) are non-relativistic.As such, causality has to be enforced up to the central density of the maximum mass configuration.Models by Malik et al. (2022); Ghosh et al. (2022); Beznogov & Raduta (2023); Malik et al. (2023); Huang et al. (2023); Zhou et al. (2023) rely on different formulations of the covariant density functional (CDF) theory of dense matter.Results of Malik et al. (2022Malik et al. ( , 2023) ) -who employ density dependent and nonlinear meson interactions, respectivelyindicate that NS EOS, properties, and composition depend on the EOS model and, in the case of the paper by Malik et al. (2023) also on the assumed strength of the nonlinear scalar vector field contribution.Comparison between posterior distributions of Malik et al. (2022) and of Beznogov & Raduta (2023) shows that even small technicalities in determining the saturation density, which enters the definition of the density dependent coupling constants, leads to small discrepancies between the results.Finally, Beznogov & Raduta (2023) show that constraints on the low density behavior of PNM act differently when conditions are posed on the energy per particle or, alternatively, pressure and that the way in which the EOS is constrained impacts on the correlations between properties of NM on the one hand and properties of NS on the other hand.
The aim of the present paper is to provide better insight into the equation of state of cold catalyzed dense matter.To this end, the non-relativistic mean field approach of nuclear matter with Skyrme effective interactions is employed.This theoretical framework has been intensively used for the study of atomic nuclei but, to our knowledge, not yet utilized in a full Bayesian inference of the NS EOS.EOS models are built by considering constraints from nuclear physics cast in terms of values of NEPs; the density dependence of the energy per neutron in PNM as calculated by the chiral effective field theory (χEFT); nucleon and neutron Landau effective masses in symmetric nuclear matter (SNM) and PNM, respectively, at the density of 0.16 fm −3 , as calculated by χEFT; a lower bound on the maximum mass of a NS as well as compliance with causality at the central density of the most massive NS configuration.
There are several reasons why we opted for this set of constraints.First, it corresponds to the strategy adopted in (Malik et al. 2022(Malik et al. , 2023;;Beznogov & Raduta 2023), meaning that the model dependence of the results can be highlighted by comparing the predictions made by these various studies.Second, by giving more weight to nuclear physics constraints, it offers a perspective complementary to that offered by studies which account for existing observational and/or potential future measurements of radii, tidal deformabilities, moments of inertia, etc.Finally, it keeps the results free from statistical and systematical uncertainties that affect some of the recent and very sophisticated astrophysical measurements on NSs.For instance, the tidal deformability extracted from GW170817 (Abbott et al. 2017;Abbott et al. 2018;Abbott et al. 2019) depends, among other, on whether the same EOS is used for each component star; the assumed EOS model; the spin ranges; the waveform model.Then, according to Vinciguerra et al. (2023), posteriors from NICER are multi-modal and subtle changes in the setup make the inference procedure prefer one mode or another, which translates into different constraints on masses and radii.
Similarly to our previous work (Beznogov & Raduta 2023), different sets of constraints are examined here.Comparison between the thus obtained posterior distributions allows to judge the effectiveness of each of these constraints as well as their compatibility.Compliance with NICER and GW170817 constraints is checked a posteriori.
A purely nucleonic composition is hypothesized for the baryonic component.
The paper is organized as follows.Sec. 2 gives a brief overview of the non-relativistic mean field model with standard Skyrme interactions.Analytic expressions for several quantities, including NEPs up to the fourth order, are also provided.Basic information on the set of constraints posed to our models and the Bayesian setup are presented in Sec. 3. Posterior distributions of NEPs and selected key properties of NS obtained for different sets of constraints together with some correlations between them are presented and analyzed in Sec. 4. The conclusions are drawn in Sec. 6.The sensitivity of the posterior distributions to the parametrization of the model and to the prior distributions is addressed in Appendix B.

THE MODEL
The standard Skyrme forces, that we use in this paper, have the form is the relative momentum operator acting on the right and k ′ is its conjugate acting on the left; is the total local density; n i (r) with i = n, p are the neutron and proton local densities.This standard form of the Skyrme interaction has ten parameters.Eight correspond to x i and t i with i = 0, 1, 2, 3, one corresponds to the power σ of the density dependence and one measures the strength W 0 of the zero-range spin-orbit term.
Within the independent particle approximation the total energy of a system of nucleons corresponds to the average value of the system's Hamiltonian Ĥ: where H (r) represents the local energy density.In the case of NM, that is homogeneous, spin-saturated and with no Coulomb interaction, the energy density H consists of four terms: where k is the kinetic energy term, h 0 is a density-independent two-body term, h 3 is a density-dependent term and h eff is a momentum-dependent term.Their expression are the following where n = n n + n p and n 3 = n n − n p stand for the isoscalar and isovector particle number densities; τ = τ n + τ p and τ 3 = τ n − τ p denote the isoscalar and isovector densities of kinetic energy; 2/m = 1/m n + 1/m p , where m i with i = n, p denotes the bare mass of nucleons.
The kinetic energy in Eq. ( 4) takes the same form as in the Fermi gas model for noninteracting fermions.The density behavior of the potential energy terms, Eqs. ( 5), ( 6) and ( 7), is governed by the coefficients C 0 , D 0 , C 3 , D 3 , C eff and D eff .Each of these coefficients can be expressed in terms of traditional Skyrme parameters via a linear expression (Ducoin et al. 2006): Eqs. ( 5)-( 7) show that in the particular case of NM, the number of independent Skyrme parameters is reduced to seven.
For the limiting cases of SNM and PNM the dependence is reduced further.Specifically, the energy density of SNM only depends on σ, C 0 , C 3 and C eff while that of PNM is governed by σ, (C 0 + D 0 ), (C 3 + D 3 ) and (C eff + D eff ).The isovector and isoscalar sectors are not independent.We note that for σ = 0 and σ = 2/3 the functional dependence on density in Eq. ( 3) gets simplified.Indeed, for the first value Eq. ( 6) has the same dependence as Eq. ( 5) while for the latter value Eq. ( 6) has the same dependence as Eq. ( 7).Throughout this paper, effective interactions are expressed in terms of C 0 , D 0 , C 3 , D 3 , C eff , D eff and σ.This results in a parameter space of lower dimension with respect to the one associated with the traditional coefficients.
In the following we provide expressions for basic thermodynamic quantities and NEPs.
From the thermodynamic definition of pressure, P(n, δ) = n 2 ∂ (H/n) /∂n| δ , where δ = n 3 /n denotes the isospin asymmetry, one obtains From Eq. ( 9) it is straightforward to compute the saturation density corresponding to an arbitrary value of δ by imposing P(n δ sat ) = 0 as well as to compute the value of energy per particle at the saturation, E sat (δ) = H(n, δ)/n| n=n δ sat .We remind that saturation points exist only for |δ| ≤ ∆, where ∆ is an effective interaction dependent value.
EOSs generated by various sets of parameters of the effective interaction are customarily analyzed in terms of the values NEPs have at zero temperature.There are two possible ways to achieve this goal.The first one consists in Taylor expanding the energy per particle E/A = H/n of matter with the particle density n and isospin asymmetry δ with respect to the deviation X = n − n δ sat /3n δ sat from the saturation density n δ sat (Dutra et al. 2012), with The most important coefficients of this series are the low order ones.As any other quantity pertaining to cold matter, these coefficients can be calculated analytically as a function of Skyrme parameters.The corresponding expressions are provided in Appendix A.
The second possibility is to express E/A as a power series in δ (Chen et al. 2009) and then further Taylor expand each coefficient in terms of the deviation x = n − n 0 sat /3n 0 sat from the saturation density of SNM, where The lowest order coefficient in δ 2 present in Eq. ( 12) is called symmetry energy and its value at n 0 sat is denoted by J sym .The expressions of E sym; 2 as well as those corresponding to its slope L sym , curvature K sym , skewness Q sym and kurtosis Z sym are offered in Appendix A.
The asymmetry energy, E asym , is defined as the per-nucleon cost of converting SNM (δ = 0) in PNM (δ = 1), and can further be expressed as Numerical comparison between Eq. ( 16) and Eq.(A5) reveals that, for densities lower than For the sake of completeness we also provide expressions for the Landau effective masses.They write: where "+" ("-") corresponds to neutrons (protons).The conditions to have 0 ≤ m eff; i /m i ≤ 1 lead to C eff ≥ 0 and  1.

THE BAYESIAN SETUP
In this section we discuss the various constraints posed to our EOS models and the procedure adopted to explore the parameter space.

Constraints
A minimum number of constraints is imposed to models built in this work.They correspond to the behavior of NM around (n = n sat , δ = 0); the energy per neutron in PNM [(E/A) PNM ] at densities lower or equal to the saturation density; nucleon (neutron) Landau effective mass in SNM (PNM) at n = 0.16 fm −3 ; lower bound on the maximum gravitational mass of a NS and causality of NS matter.For values, see Table 1.
The NM-related conditions refer to the values the best known NEPs, n sat , E sat , J sym and K sat , have.For each of the considered NEPs, median values and standard deviations (SDs) are fixed to the values obtained by Margueron et al. (2018) based on a collection of 35 standard Skyrme interactions that are frequently employed in the literature (SDs are rounded up).This choice ensures that the domains of values used as a reference correspond to well-behaved energy density functionals identical to those assumed here.We remind at this point that the values obtained for NEPs depend on both the assumed energy density functional and the analyses itself, for a discussion see (Oertel et al. 2017;Margueron et al. 2018).Out for the four NEPs in Table 1, K sat is the least controlled one (Oertel et al. 2017).
The density behavior of PNM is controlled by the values the energy per particle has in a minimum of three and a maximum of four points with densities in the range 0.04 fm −3 ≤ n ≤ 0.16 fm −3 .For all these energies we use the values provided by Somasundaram et al. (2021), where results of χEFT calculations corresponding to six Hamiltonians (Drischler et al. 2016) are reported.The nucleon-nucleon (NN) interactions computed at N 3 LO are supplemented by three-nucleon (3N) interactions computed at N 2 LO.The spread in the energy per particle obtained from these interactions can serve as an uncertainty estimate.More importantly, access to the individual energy per particle values provided by each Hamiltonian, rather than to the values of the means and SDs, allows to calculate correlations between the values corresponding to different densities, see Sec. 3.2.Some of the models built in this paper also impose constraints on the values the nucleon (neutron) effective mass has in SNM (PNM) at n = 0.16 fm −3 .Similarly to the density behavior of the energy per particle in PNM, the values obtained by Somasundaram et al. (2021) are used.Due to the limited flexibility of Landau effective masses as functions of density in models with standard Skyrme interactions (see Eq. ( 17) and Sec.4.1), constraints on m SNM eff; N and m PNM eff; n are not imposed over a wider density domain despite being available from Somasundaram et al. (2021).
The fact that the values used as a reference for NEPs on the one hand and (E/A) PNM , m PNM eff; n and m SNM eff; N on the other hand stem from different sources and have been derived within dissimilar theoretical frameworks might, in principle, raise some tension.In particular, we note that Margueron et al. (2018) provides for m SNM eff; N /m N at saturation a value, 0.77 ± 0.14, that is by 20% higher that the one of Somasundaram et al. (2021).Besides, according to Fig. 3 by Somasundaram et al. (2021), n sat = 0.163 fm −3 and E sat = −15.07MeV.The first value agrees, within 1σ, with the value in Table IV of Margueron et al. (2018), while the latter is by 6% lower than the value put forward in Table IV of Margueron et al. (2018).The compatibility between constraints derived from phenomenological and ab initio models will be addressed elsewhere.
NS matter is minimally constrained.First, we ask that any model complies with the 2 M ⊙ lower limit on the maximum NS mass.Second, we ask that causality is fulfilled at the density corresponding to the central density of the maximum mass configuration.
We note that each of the NM constraints in In a "deterministic approach" four sharp constraints on (E/A) PN M would translate into well defined values for σ, (C 0 + D 0 ), (C 3 + D 3 ), (C eff + D eff ); should an extra sharp constraint on m PNM eff; n be posed at any non-vanishing density, with a high probability, no compatible model will be found.If this is the case, values for the (linear combinations of) effective interaction parameters can be determined through a χ 2 -fitting procedure.This is exactly the procedure used to construct phenomenological effective interactions based on reproduction of nuclear physics observables.As a matter of fact, in most of these cases the number of constraints exceeds the number of the degrees of freedom of the interaction.The fact that, with the exception of run 4, the energy of PNM is constrained in three density values only, confers flexibility to our sets of models; extra flexibility comes from the fact that, instead of providing sharp values, we allow the quantities on which constraints are imposed to span some domains.
Finally, note that constraints imposed on NS properties cannot be connected in a straightforward manner with either parameter of the effective interaction, the reason being that the value that each of these quantity takes depends on the behavior of matter over wide domains of density and proton fraction.

Likelihood
Similarly to the situation we have considered previously (Beznogov & Raduta 2023), only the lower limit on M * G and the upper limit on c s can be regarded as "purely experimental" information.All other constraints are obtained employing sophisticated theoretical modeling.As such, a question of their independence/uncorrelatedness/compatibility arises.We will not go into details here.Note that the topic of independence/uncorrelatedness was discussed at lengths by Beznogov & Raduta (2023) (cf. also Supplemental Materials of Somasundaram et al. (2021)).Some issues regarding compatibility will be addressed in Sec. 4.Here we only focus on the differences compared to the situation of Beznogov & Raduta (2023).
First, since we have employed for the energy per particle of PNM a set of data (Drischler et al. 2016;Somasundaram et al. 2021), where values are provided for each of the six considered Hamiltonians, it is possible to take into account correlations when constraining (E/A) PNM at different values of density.Following the procedure described by Somasundaram et al. (2021) in the Supplemental Materials, we can compute the covariance matrix as follows: Here n i and n j are the values of density at which we compute the correlations between the values of the energy per particle E = E/A and k runs over the six Hamiltonians investigated by Somasundaram et al. (2021).Considering that the mean is also estimated from the same sample, we have introduced Bessel's correction (i.e., dived by 5 instead of 6) to get an unbiased estimate of the covariance matrix.
The log-likelihood function for runs 1-4 has the "typical" form for uncorrelated constraints (Beznogov & Raduta 2023): where q labels the run, N q is the number of constraints for that run except for M * G and c * 2 s constraints, d i and Z i stand for the constraint and its SD as mentioned in Table 1; ξ i (Θ) corresponds to the value the model defined by the parameter set Θ provides for the quantity i.However, for run 0 the log-likelihood reads: where first term corresponds to the constraints that we always treat as uncorrelated (e.g., n sat , E sat , K sat , J sym ) and the second term deals with energy per particle constraints that we consider to be correlated in run 0; δE stands for the difference between the value of the constraint and the value provided by the model.
In Eqs. ( 19) and ( 20) we have discarded the normalization factors as they are not relevant for sampling from the posterior distributions (as long as they do not depend on the input parameters, which they do not in our case).In the same way as in (Beznogov & Raduta 2023), our primary method of Bayesian inference was affine invariant Markov Chain Monte Carlo (MCMC).See the details of implementation, employed software, analysis of convergence and of the effects of the prior in Appendix B.
Second, now we have two "impenetrable wall" constraints.One is the lower limit on M * G and the other is the upper limit on c * 2 s .They are dealt with in a similar fashion to (Beznogov & Raduta 2023), i.e., if one or another (or both) is violated, the log-likelihood function is assigned a large negative value (−10 10 ).This insures that such a point is never accepted into the Markov chain.By employing a similar rejection procedure, we have also discarded all configurations where NS matter was thermodynamically unstable, i.e., where pressure was either negative (at some densities) or decreasing with density.

Prior
The last crucial ingredient of any Bayesian-like analysis is the prior.For technical reason described in Appendix B, we have employed "mixed" input parameters, i.e. instead of using all seven effective interaction parameters (C 0 , D 0 , C 3 , D 3 , C eff , D eff and σ) as the input ones, we have actually used three NEPs (n sat , E sat and J sym ) and four effective interaction parameters (D 3 , C eff , D eff and σ) as the input.The remaining parameters of the effective interaction, C 3 , C 0 and D 0 , were computed according to with n 2/3 sat .The domains of the priors are presented in Table 3.For the three NEPs inputs, the domains cover 10 SDs range symmetrically around the median value (cf.Table 1).For D 3 , the domain was chosen somewhat arbitrarily based on trial and error procedure.The domains of C eff and D eff take into account that the effective masses should lie within 0.02 ≲ m eff /m N ≤ 1, with the lower limit chosen arbitrarily.The domain we choose for σ is slightly larger than the one span by the ensemble of Skyrme interactions, see Table B1 in Dutra (2011), where parameters of the effective interactions are compiled for 240 models.Indeed, according to Table B1, 0.08 ≤ σ ≤ 1.00.Thus, our choice encompasses all behaviors allowed by stan-dard parametrizations.See also Appendix B regarding the importance of σ domain range.
We have chosen uniform (uninformative) priors in the ranges provided in Table 3. Notice, however, that the constraint |D eff | ≤ C eff (see Sec. 2) means that the actual prior distribution of C eff and D eff will not be uniform after this constraint is enforced.Also notice that uniform priors in n sat , E sat , J sym , D 3 , C eff , D eff and σ are not equivalent to uniform priors in the seven effective interaction parameters.See details in Appendix B.

RESULTS
In this section we analyze posterior probability density distributions of various quantities pertaining to SNM, PNM, NS matter and NS properties, as obtained for each of the five runs.For the sake of the comparison, we also show the envelope of 75 standard Skyrme models from Dutra et al. (2012), which satisfy the conditions: i) NS matter is thermodynamically stable, i.e., P NS is positive and non-decreasing as a function of density, ii) MeV.The latter condition corresponds to a dispersion within 5 SDs from the central value reported in Table 1.Hereafter, we refer to these 75 models as Comparison Set models.For each situation we discuss the impact of each type of constraints.Finally, a number of correlations among NM and NS parameters is discussed.
Before analyzing the results it is worth mentioning that different runs in this work have different "fitting" quality.While all the runs of Beznogov & Raduta (2023) have very similar values of χ 2 for their respective maximum a posteriori (MAP) points, in the present work values of χ 2 for MAP points are noticeably different.They are: This is our first indication that standard Skyrme model might not be flexible enough to account for all of the imposed constraints.Further evidence is discussed in the following subsections.

Nuclear matter
The density dependence of the energy per nucleon in SNM and PNM is investigated in Fig. 1 (left and middle columns).Results corresponding to different runs in Table 2 are reported on subsequent rows.The considered density domain is 0.02 ≤ n ≤ 0.5 fm −3 .This and the following (Figs. 2 and  3) conditional probability density plots of the posterior distributions are normalized such that for any given fixed value of n the total probability (of a considered function of n to have any value) is 1.0.The colormap is normalized accordingly, i.e., black color corresponds to the maximum probability density at each fixed n (unlike the joint probability density plots where the maximum of the colormap would correspond to the maximum of probability density in the whole 2D domain, not at each cross-section).In this and the following figures, CR stands for credible region.
Regarding the behavior of SNM we note that: (i) different runs feature rather similar dispersions between (E/A) SNM curves; this is understandable given that the behavior of the isoscalar channel is constrained by conditions imposed on n sat , E sat and K sat and those are common for all runs; curves in each family start to diverge at a density slightly larger than n sat ; this is a remainder that constraints at n sat are not sufficient to control the high density behavior of the EOS, (ii) one more constraint posed on (E/A) PNM over n ≤ n sat does not seem to impact the scattering of curves over the complementary domain n > n sat ; indeed, the difference between the dispersion of curves in runs 1 and 4, which contains one extra constraint with respect to run 1, is small, (iii) accounting for correlations among (E/A) PNM over n ≤ n sat slightly reduces the dispersion of curves over the complementary domain of density (run 0 vs. run 1), (iv) implementation of the constraint on m SNM eff; N slightly broadens the uncertainty band of (E/A) SNM over n > n sat (run 2 vs. run 1), while no such effect is obtained when the constraint is imposed on m PNM eff; n (run 3 vs.run 1), (v) the envelopes corresponding to each of our runs and the one that corresponds to the Comparison Set models are of similar width; overall, models built in this work prefer slightly lower values than those of Comparison Set models, (vi) the median curves of runs 0 and 4 are close to the median curve in run 1; the median curve in run 3 (run 2) explores larger (lower) values than those of the median curve in run 1.Now, let us consider the consequences that various constraints have on (E/A) PNM .The following comments are in order: (i) run 0, which accounts for correlations among (E/A) PNM at 0.08, 0.12 and 0.16 fm −3 , provides by far the smallest dispersion of curves; this means that accounting for correlations at n ≤ n sat is very efficient in controlling the behavior of PNM even at densities much larger than n sat , (ii) any extra constraint imposed on the behavior of PNM at low densities results in reducing the curve dispersion over n > n sat , see the uncertainty bands in runs 4 and 1, (iii) constraints on the effective mass at saturation are efficient in reducing the dispersion (runs 2 and 3 vs.run 1); constraints on m SNM eff; N are more efficient than those on m PNM eff; n ; the reason is that conditions on m SNM eff; N control C eff , while those on m PNM eff; n control (C eff + D eff ); together with (C 0 + D 0 ) and (C 3 + D 3 ), (C eff + D eff ) is also controlled by the conditions on (E/A) i , though in a loose way (see Fig. 6), (iv) the widths of the envelopes that correspond to all runs except run 1 are much lower than the width of the envelope associated with the Comparison Set models; this means that the models in each of our sets manifest less variety in the density dependence of the symmetry energy; the situation can be attributed to the constraints we have imposed on (E/A) i and that are absent in many standard Skyrme interactions, (v) the largest (lowest) value of the median curve at n = 0.5 fm −3 corresponds to run 1 (run 0), that is, to the run with the largest (lowest) dispersion, (vi) the constraints we have imposed play a much larger role on (E/A)(n) in PNM than on (E/A)(n) in SNM; this is easy to understand given that most of these constraints are imposed on the behavior of PNM, (vii) over n ≲ 0.16 fm −3 , the posterior distributions of all our runs are in good accord with the results from χEFT (Somasundaram et al. 2021), illustrated by a light blue band.The right column in Fig. 1 addresses the density dependence of E asym , Eq. ( 15).Largely different behaviors of PNM allowed by models belonging to each run are translated into largely different density dependencies of the asymmetry energy.In particular, runs 0 and 1, which provide the most and the least constrained behavior of (E/A) PNM (n) (see the middle column), respectively, manifest the smallest and largest dispersion of E asym (n) at n ≥ n sat .It is worth noting that all our runs accommodate models where E asym increases with density along with models that show the opposite behavior.Models in the second category allow for negative values of E asym at den- sities exceeding a certain value.The Comparison Set models manifest the same variety of behaviors.The only cases where the median curve increases with density correspond to run 1 and run 2.
Recently, constraints on the symmetry energy and symmetry pressure (P sym = n 2 ∂E sym /∂n) at around (2/3)n sat have become available.They have been obtained either from homogeneous sets of measurements/calculations or compilations of data that span a range of densities.Examples are provided by P sym;0.1 = 2.38 ± 0.75 MeV/fm 3 , extracted from the neutron skin thickness of 208 Pb (Adhikari et al. 2021;Reed et al. 2021); P sym;0.1 = 1.8 ± 0.2 MeV/fm 3 , inferred from a collection of nuclear structure and heavy ion collision data and astrophysical observations of NSs (Lynch & Tsang 2022); E sym;0.1 = 25.6 +1.4  −1.3 MeV, obtained from a correlation between a linear combination of chemical potentials in nuclei neighboring shell closures and the symmetry energy at (2/3)n sat (Qiu et al. 2024).A χEFT E asym;0.1 value can be computed from the energies per nucleon in SNM and PNM provided in (Somasundaram et al. 2021) and amounts to 23.65 ± 0.37 MeV.The values of E sym;0.1 and P sym;0.1 corresponding to each of our runs are reported in Table 5.In all cases, the medians and SDs of these quantities are slightly smaller than those in (Lynch & Tsang 2022;Qiu et al. 2024).The relative similitude with P sym;0.1 in (Lynch & Tsang 2022) can be attributed to the equivalence of the yet different constraints.The relative similitude with E sym;0.1 in (Qiu et al. 2024) is most probably attributable to our usage of realistic constraints on J sym and E/A in PNM.Our E sym;0.1 -values are in excellent agreement with E asym;0.1 extracted from (Somasundaram et al. 2021); this can be seen as an a posteriori proof of effectiveness of constraints on E/A in PNM.Our P sym;0.1 -values are much lower than those in (Reed et al. 2021).This, however, should not come as a surprise given that the values of J sym = 38.1 ± 4.7 MeV and L sym = 106 ± 37 MeV in (Reed et al. 2021) exceed considerably those in this work, see Sec.4.3.2, as well as other values extracted from nuclear physics (Oertel et al. 2017).
Further insight into the properties of NM is offered in Fig. 2, where the effective nucleon mass in SNM (left column), effective neutron mass in PNM (middle column) and neutronproton effective mass splitting in neutron rich matter with Y p = 0.2 (right column) are plotted as a function of density.The results corresponding to the Comparison Set models and the predictions of χEFT calculations of Somasundaram et al. (2021) are illustrated as well.We note that all models generated within all of the runs are characterized by effective masses smaller or equal to the bare mass and that they decrease or stay constant with the density.This is a direct consequence of conditions imposed on C eff and C eff ± D eff (see Sec. 2) in order to produce this behavior.Note that this is nevertheless not the case of all standard Skyrme interactions is the literature.For instance, Ska45s20 (Dutra et al. 2012), Sefm09 (Dutra et al. 2012), Sefm1 (Dutra et al. 2012), SK255 (Agrawal et al. 2003), which are present in the Comparison Set, lead to m PNM eff; n increasing with density.Ska45s20 also leads to m SNM eff; N increasing with density.Other examples are offered by BSk1 (Samyn et al. 2002) and BSk2 (Samyn et al. 2002) for which both m PNM eff; n and m SNM eff; N augment with the density.The density behavior of the upper boundary of these models' envelope is due to Ska45s20 (Sefm1) for SNM (PNM).
The light blue shadowed areas, illustrating the predictions of microscopic calculations of Somasundaram et al. (2021), show that the decrease in m SNM eff; N (n) gets saturated around 0.2 fm −3 and m PNM eff; n (n) is non-monotonic.As obvious from the analytic expressions in Eq. ( 17), the standard Skyrme effective interactions are too simple to account for such a behavior.In order for it to be described a more complex density dependence of the effective mass is needed.Bruxelles extended Skyrme models (Chamel et al. 2009) provide a possible solution and the performances of more sophisticated Skyrme interactions will be considered elsewhere.
The left column in Fig. 2 shows that, with the exception of run 2 where a constraint is posed on the value of the effective mass of nucleons in SNM at the density 0.16 fm −3 , the EOS models in each run manifest a large dispersion of m SNM eff; N (n); the largest dispersion by 90% CR corresponds to EOS models in run 3, where a constraint is posed on the value of neutron effective mass in PNM.The scattering of models in run 2 is small, which means that the constraint is efficient.Moreover, one can see that the target value of m SNM eff; N at 0.16 fm −3 is met rather accurately.The dispersion that characterizes runs 0, 1, 3 and 4 makes that a significant percentage of our models predict for m SNM eff; N values smaller than the those of the Comparison Set.This situation is attributable to the much larger variety of models explored here.We also note that absence of constraints sensitive to the value of m SNM eff; N in the fitting procedure of the Comparison Set models makes that some of them have m SNM eff; N increasing with density and exceeding the bare nucleon mass.
The middle column in Fig. 2 shows that m PNM eff; n (n) is very sensitive to the way in which the constraints are incorporated.This situation is easy to anticipate.Accounting for correlations among the values of (E/A) i in PNM at different densities makes models gather in a relatively narrow valley though the envelope remains as large as when correlations are disregarded (run 0 vs. run 1).Constraints on effective masses at the saturation density result in the accumulation of models in a narrow band (runs 2 and 3 vs.run 1).Posterior distributions from run 2 and run 3 differ much one from the other.As expected, the effect is stronger when the constraint is imposed on the neutron effective mass in PNM than on the nucleon effective mass in SNM.The first of these constraints also leads to values of m PNM eff; n (n) larger than those obtained when the second constraint is imposed (run 3 vs.run 2).The fact that at 0.16 fm −3 the values of m PNM eff; n in run 3 are close to the target values means that the constraint is effectual.The effect obtained when an extra constraint is posed on (E/A) in PNM at low densities (run 4) is weaker than the one obtained when correlations among values of (E/A) i are accounted for (run 0).We remind that, for obvious reasons, this was also the case of (E/A)(n) in PNM, see Fig. 1 (middle column).We note that the domain of maximum curve density in run 4 has a width intermediate to widths of domains of maximum curve density in run 0 and run 1.Finally, from all runs, the upper value border of all models is characterized by m PNM eff; n (n)/m n = 1, which corresponds to The right column in Fig. 2 illustrates the density dependence of δm eff; np = m eff; n /m n − m eff; p /m p in neutron-rich matter with the proton fraction Y p = n p /n = (1 − δ)/2 equal to 0.2 (arbitrarily chosen).For all runs a wide spread of curves is observed; the values range between −0.35 and 0.35.Runs 0 and 3 mostly generate positive values, which means that most models agree with the predictions of microscopic calculations that include three body forces (Sammarruca et al. 2005;Dalen et al. 2005;Baldo et al. 2014;Shang et al. 2020).We remind that all of the latter calculations provide values of neutron Landau effective mass that exceed, in neutron rich matter, the values of the proton Landau effective mass.This agreement does not come as a surprise given that run 0 accounts for correlations among values of (E/A) in PNM at different densities and run 3 explicitly implements a condition on the neutron effective mass in PNM.These two runs also produce the strongest accumulation of models in a narrow domain; they correspond to the situations where the smallest dispersion of m PNM eff; n (n) is obtained.The median curves in these runs explore values round 0.1 (0.2) for run 0 (run 3).Run 2 mostly generates negative values; this shows that conditions posed on m SNM eff; N at 0.16 fm −3 play differently than those posed on m PNM eff; n at 0.16 fm −3 .The median curve in run 2 is characterized by values of the order −0.1.The models in the runs 1 and 4 are more uniformly distributed around δm eff; np = 0; models in the run 4, which implements constraints on a wider density domain, gather in a moderately narrow band centered around 0.05 while no such behavior is shown by models in run 1.Positive (negative) values of δm eff; np are obtained by models with D eff < 0 (D eff > 0).Positive and negative values of δm eff; np are obtained also by models in the Comparison Set; the median curve of this bunch of models is characterized by δm eff; np < 0.

Neutron stars matter
In this subsection we consider the role that each of the constraints we impose play on the NS EOS.NS EOSs are obtained by smoothly matching the core and the crust EOSs.For the outer and inner crusts we adopt the models by Haensel et al. (1989) and Negele & Vautherin (1973), respectively.
The smooth matching between crust and core EOSs is done at a density of about n sat /2.In the leptonic sector both electrons and muons are accounted for.Modifications brought to the NS stiffness are investigated by analyzing the density dependence of pressure (P NS ) and speed of sound squared (c 2 s ).We remind that c 2 s /c 2 = (dP/de) fr , where the subscript "fr" indicates that the derivatives have to be evaluated with the composition frozen.Modifications brought to the chemical composition are judged from the density dependence of the proton fraction Y p .The predictions of the five sets of models considered in this paper are depicted in Fig. 3.
The left and middle columns show that the patterns in (E/A)(n) in PNM are somehow replicated in P NS (n) and c 2 s (n).More precisely, the smallest and largest dispersions among the P NS (n) and c 2 s (n) curves generated within one run correspond to the run 0 and run 1, respectively, which show the smallest and largest dispersions among (E/A)(n) in PNM curves, see Fig. 1 (middle column).For the huge majority of models c 2 s is a monotonic function of n.In this situation, models with values of the speed of sound exceeding, at high n-values, the speed of light should not be considered as violating causality.The situation should be rather understood as the consequence of prolonging the curves beyond the largest density value for which stable NSs configurations are obtained.Indeed, compliance with causality is imposed in our models by requiring that, for a density equal to the central density of the maximum mass NS configuration (n * c ), c 2 s ≤ 1.For a tiny fraction of the models the maximum value of c 2 s is nevertheless reached at a density lower than n * c .Some of those may violate causality.We have checked that, due to their low number, such models have no effect on the statistics.
The sudden change of slope in the lower boundary of the pressure envelope of the runs 1, 3 and 4 is attributable to the softening that the softest models in these runs suffer at high densities.This phenomenon is responsible also for the behavior of the lower boundary of the c 2 s envelope that, starting from 0.6 − 0.7 fm −3 , decreases with the density, see the middle column.With the exception of runs 0 (the most constrained) and 1 (the least constrained), the models generated here manifest roughly the same scattering in P NS (n) as the Comparison Set models, though models in the latter set do not exhibit the extra softening at high densities.As a consequence, the agreement between the values that c 2 s spans in the runs 2, 3, and 4 and the values that c 2 s covers in the Comparison Set holds only up to ≈ 0.7 fm −3 .
The widely different behaviors of E asym (n) in our runs lead to a large dispersion of the proton fraction Y p (n) in βequilibrated matter.Indeed, all the five runs contain models where Y p increases steeply with the density as well as models where Y p vanishes at densities a few times larger than n sat .Runs 1 and 2 show that the models that yield large (small) values of Y p also yield large (small) values of P NS , the strength of this correlation depending on the value of density.Note, nevertheless, that this correlation is not universal, as P NS (n) has an explicit dependence on Y p (n), E sym (n) and dE sym (n)/dn (see Eq. ( 101) of Steiner et al. (2005)).Models with Y p steeply increasing with n will allow the direct URCA process to operate even in low mass stars.It is know that in order for direct URCA to operate, Y p should exceed a threshold value.The latter depends on the electron (n e ) and muon (n µ ) densities and writes Y p; DU = 1/ 1 + 1 + x 1/3 e 3 , with x e = n e / n e + n µ .Models with vanishing Y p not only do not accommodate for direct URCA, but lead to NSs whose inner core is exclusively made of neutrons.The shape of the lower envelope of the Comparison Set models shows that such a behavior is not an artifact of our models.Indeed, vanishing Y p values are caused by negative values of E asym .We note that only the models in the runs 1 and 2 have median curves that increase with the density, which corresponds to E sym (or E asym ) increasing with the density, see Fig. 1

(right column).
M −R diagrams corresponding to the models in each run are represented in Fig. 4 as conditional probability density (also known as curve density) plots P(R|M).Unlike the previous conditional probability density figures, this plot requires an additional normalization.The reason is that for M > 2 M ⊙ the number of curves is no longer constant, it decreases until there is only one curve left corresponding to the maximum of the maximum masses.Thus, we have first to normalize the 2D bins' counts corresponding to a given value of M to the number of EOSs that exist for this value of M. Then we normalize the colormap to the maximum probability density at each fixed value of M similarly to the previous conditional probability density plots (for the latter, see the explanation in the beginning of Sec.4.1).Results corresponding to each run are depicted in a separate panel.It comes out that the run 0 and run 1, which provide the smallest and the largest dispersion among the models, populate the narrowest and the broadest domains of M − R space, respectively.For instance, the difference among the largest and smallest radii of 1.4 M ⊙ stars is 1.2 km for run 0 and 3.0 km for run 1.The corresponding figures for the 2.0 M ⊙ stars are 2.0 km and 4.3 km, respectively.The lowest (largest) values of maximum gravitational mass, M * G , are obtained for the run 0 (run 1), which shows the softest (stiffest) P NS (n) dependence (see Fig. 3 left column).The runs 0 and 1 distinguish one from the other also by the sensitivity NS radii have to NS masses.Indeed, correlations between (E/A) in PNM at different densities for models in the run 0 are translated in NS radii which decrease with the increase of NS masses.The lowest number of constraints implemented in the run 1 is responsible for a more ambiguous evolution of R(M), especially in low mass stars.The largest radii predicted by the Comparison Set models for stars with M/M ⊙ ≤ 1.5 are significantly larger than the largest radii predicted by our models.These large values are due to the models with particularly steep dependence of E asym as a function of n, see the right column of Fig. 1.The fact that this situation does not hold also in more massive stars is a nice illustration that the density dependence of symmetry energy impacts only low and intermediate mass NSs.
Joint mass and radius posterior probability density contours corresponding to PSR J0030+0451 (Miller et al. 2019b), PSR J0740+6620 (Miller et al. 2021) and to the NS in the SNR HESS J1731-34 (Doroshenko et al. 2022) are also plotted in each of the six panels in Fig. 4 at 50% (solid lines) and 90% (dashed lines) CR.We present two variants of constraints for HESS J1731-34: orange contours correspond to the "X-ray only" data of Doroshenko et al. (2022), while red contours correspond to the "full priors" data.Note that the latter data include specific assumptions on the NS EOS (see Doroshenko et al. 2022).Those assumptions can be incompatible with our assumptions, thus rendering the comparison between our results and the "full prior" dataset meaningless.The obvious conclusions are that, i) in each of the five runs there are models that comply at 50% CR with "constraints" from PSR J0030+0451 and PSR J0740+6620, ii) none of those also comply with any of the 50% CR "constraints" from HESS J1731-34.The same holds for the Comparison Set models.These features signal either a tension between the χEFT-inspired conditions that we have imposed on (E/A) PNM and HESS J1731-34, or insufficient flexibility of the density functional we have assumed.This aspect will be addressed elsewhere.
Predictions of the least and most constrained runs in what regards the combined tidal deformability Λ as a function of the mass ratio q are confronted in Fig. 5, where the constraints from GW170817 (Abbott et al. 2019) are also plotted.In addition to the expected shrinking of the uncertainty band upon accounting for correlations among (E/A) i in PNM, we note that median values in run 0 are somewhat lower than those in run 1.The values spanned in run 0 are slightly higher than the median value from (Abbott et al. 2019); the 90% CR in run 1 is very close to the 90% upper bound from (Abbott et al. 2019), but this is a coincidence.One can see that our calculations result in larger values of Λ compared to the GW data.They are, nevertheless, well within 90% confidence interval of the observational constraints.

Model dependence of the results
The role that the different constraints play on models of NM and NSs is discussed in this subsection by analyzing 1D marginalized posterior distributions of different quantities.Considered are: parameters of the effective interactions as well as linear combinations of them; parameters of NM; energy per particle in PNM at various densities lower or equal to 0.16 fm −3 ; nucleon (neutron) effective mass in SNM (PNM); global properties of NS.  2 are considered.X-and Y-axis ranges were chosen to increase readability.As such, some of the very wide or very narrow distributions are cut at X-axis or Y-axis limits, respectively.certain degeneracy is introduced in both isoscalar and isovector channels.This feature makes the conditions imposed on the various quantities easier to meet, which translates into accumulation of points.The posterior distribution of σ in run 3 peaks around 2/3 and towards the maximum value allowed in the simulation, 1.1; also, the values σ ≲ 0.3 are forbidden.The preference for large values of σ as well as the suppression of low values of σ stem from the necessity to compensate for the low values of (C eff + D eff ), which decrease the h eff term of pressure, see Eq. ( 9).Yet, the pressure should be high enough to comply with the 2 M ⊙ constraint.The low population of the narrow domain around σ ≈ 3/4 can be regarded as the outcome of strong competition among σ = 2/3 and σ = 1.1.
Run 1 corresponds to a looser version of run 4.This suggests that posterior distributions of (C 0 + D 0 ), (C 3 + D 3 ), (C eff + D eff ) should be wider than those for run 4. Fig. 6 confirms also this expectation; in particular, the strongest widening is obtained for (C eff + D eff ).This figure also shows that the posterior distribution of σ in run 1 is peaked at the same values as the posterior distribution for run 4, though this time preference is given to low σ values.This corroborates the interdependence between (C eff + D eff ) and σ noticed in the previous paragraph.
Correlations between values taken by (E/A) in PNM at different densities, accounted for in run 0, suggest that the peaks in (C 0 + D 0 ), (C 3 + D 3 ), (C eff + D eff ) should be narrower than those corresponding to run 1.Fig. 6 shows that this is, indeed, the case.It also shows that the distribution of σ resembles qualitatively the distribution in run 4 and run 1.A clear preference for the lowest allowed value 0.01 is to be noticed.
In run 2 the extra constraint is posed on m SNM eff; N , which translates into a constraint on C eff .The posterior distribution of (C eff + D eff ) becomes rather wide and the large value tail of σ distribution is suppressed, which makes that the peak around σ = 2/3 disappears.
Information on the posterior distributions for the C 0 , D 0 , C 3 , D 3 , C eff , D eff parameters of the effective interaction is provided in the other panels.It turns out that the distributions of the C-parameters somewhat mirror the distributions of the D-parameters.The only exception from this rule corresponds to C eff and D eff in run 2, where the "symmetry" is "broken" by the constraint on m SNM eff; N that acts only on C eff .In fact, run 2 is the only case where the distributions of C eff and (to a lesser extent) D eff are narrow.This coherence between C-and Dparameters obviously means that the behaviors in isovector and isoscalar channels are not independent.
Finally, we note that, in some cases, the posterior distributions of some parameters are cut.Specifically, these are D 3 in runs 1 and 2 and σ in run 3.This situation is the direct consequence of the chosen priors, see Sec. 3, and suggests that slightly different distributions would be obtained if parameters were allowed to span wider domains.Note that other parameters that appear to be cut on Fig. 6 are not really cut, they just go outside the axis limits, which were chosen to increase readability.

Posterior distributions of parameters of NM, effective masses and (E/A) in PNM
In order to increase the efficiency, the MCMC exploration is done by controlling n sat , E sat and J sym , see Sec. 3.2 and Appendix B.
The first two quantities are allowed to span narrow domains of values, see Table 1.Fig. 7 shows that marginalized posterior distributions meet very well the target distributions.The fact that distributions corresponding to different runs sit one on the top of the other is explained by the fact that the constraint on K sat is common to all runs and other constraints in Table 1, which are specific to one run or another, have limited effect on the zero order parameters in the isoscalar channel.Indeed, out of the parameters that govern the density dependence in the isoscalar channel, the constraints on (E/A) in PNM only act on σ while that on m SNM eff; N acts only on C eff .The posterior distributions of K sat differ much from one run to another.This shows that the high order parameters in the density dependence of symmetric matter do depend on constraints on PNM as they do depend on extra constraints on SNM.This result is another manifestation of the entanglement of isoscalar and isovector sectors in Skyrme models, already visible in Eq. ( 6).Distributions of K sat in run 0, run 1 and run 4 are double peaked and qualitatively similar; each of the peaks is correlated with a peak in the posterior distribution of σ.Distributions of K sat in run 2 and run 3 differ from distributions of other runs and among themselves; the same was the case of posterior distributions of σ.The fact that all the posterior distributions of K sat miss the target distribution can be regarded as an incompatibility between different conditions that we have imposed.Patterns in the distributions of K sat are replicated in the distributions of higher order parameters Q sat and Z sat .Distributions of these latter quantities also extend over wider domains, which shows that uncertainties in the density behavior of the EOS increase with the density.It appears that low values in K sat are correlated with low values in Q sat and high values in Z sat .Notice that the same holds for EOSs derived within CDF theory, irrespective whether the meson couplings depend on density (Malik et al. 2022;Beznogov & Raduta 2023) or include non-linear terms (Malik et al. 2023), which suggests that those correlations are universal.
Posterior distributions of J sym in runs 0, 1, 2 and 4 are bellshaped, while the corresponding distribution in run 3, where an extra constraint was posed on m PNM eff; n , is asymmetric; posterior distributions in runs 0, 1 and 4 are quite similar and cover the same domain of values while those for runs 2 and 3 are shifted towards larger and lower values, respectively.Runs 2 and 3 also provide the narrowest and largest distributions.Given that runs 0, 1 and 4 contain only constraints on (E/A) PNM , while runs 2 (3) also constrains m SNM eff; N (m PNM eff; n ), the result looks natural.While none of the runs is successful in reproducing the shape of the target distribution, the peak value is met with relatively high accuracy.This was not the case of the K sat distributions.Posterior distributions of L sym are bell-shaped for all runs and differ somewhat from one run to another.We note in particular that run 1, that allows for the largest values of E asym (see Fig. 1), also allows for the largest values of L sym .Discrepancies in L sym mean that the various sets of constraints are effectual in controlling the low density behavior of the isovector channel and that each set of constraints act differently.The situation is very different in what regards high order parameters.Indeed, the distributions of K sym , Q sym and Z sym are wide and increasingly different from one run to the other.We note in particular that the only double peaked distribution of K sym corresponds to run 1, while the only single-peaked distributions of Q sym and Z sym correspond to run 2. Peaks in the posterior distributions of Q sym and Z sym appear to be correlated with the peaks in σ.Patterns in the Q sym distributions are to some extent replicated in the Z sym distributions and low values of one quantity are correlated with high values of the second.This negative correlation between Q sym and Z sym reminds the negative correlation between Q sat and Z sat discussed in the previous paragraph.The similitude between the behaviors of the third and fourth order parameters in the isoscalar and isovector sectors is striking the more if we note that negative correlations between Q sym and Z sym are spotted also in CDF models (Malik et al. 2022;Beznogov & Raduta 2023;Malik et al. 2023) though their strength is model-dependent.
The target distributions of m SNM eff; N and m PNM eff; n at 0.16 fm −3 are met only in runs 2 and 3, where conditions on these quantities are posed.In other runs, the posterior distributions of the two effective masses show complex structures and extend over wide domains.Distributions of m SNM eff; N at 0.16 fm −3 in runs 0, 1 and, to a lesser extent, 4 are similar; there are no two runs that provide, for m PNM eff; n , similar distributions.The double peaked structure of m SNM eff; N (m PNM eff; n ) in run 3 (run 1) obviously originates from the double peaked structure of C eff (C eff + D eff ), see Fig. 6.Small SDs of (E/A) in PNM at low densities make that the corresponding constraints are very efficient.As a consequence, run 4 (runs 1, 2 and, to a lesser extent, 3) meet the target distribution at 0.04 fm −3 (0.08 fm −3 ) while in all other cases target distributions on (E/A) in PNM are missed.The situation of (E/A) 1 in run 0 is interesting as even if its distribution does not agree with the target, it is still fairly close to it and of comparable width.Notice that for run 0 the value of (E/A) in PNM at 0.04 fm −3 is not explicitly constrained.This result is the out-turn of accounting for correlations at different (but higher) densities.We also note that the distribution of (E/A) 1 for run 1 is wide.This means that the uncorrelated constraints posed at higher densities are ineffective at low densities.The fact that also the distribution of (E/A) 1 for run 2 is wide shows that the extra constraint in this run has limited effect, which is understandable given that it mainly acts on SNM.
Medians and 68% CI of marginalized posterior distributions of key NM quantities, including those plotted in Fig. 7, are provided in Table 5 in Appendix C.

Posterior distributions of global parameters of NS
Equilibrium properties of non-rotating and sphericallysymmetric NSs are inferred by solving, for each model of NS EOS, the Tolman-Oppenheimer-Volkoff equations.Tidal deformabilities are computed according to Hinderer (2008); Hinderer et al. (2010).The lowest mass of a NS that accommodates for the direct URCA process (M DU ) is the mass of the NS with ρ c = ρ DU , where the latter quantity is defined by Posterior distributions of selected global properties of NSs are plotted on Fig. 8.The results of all five runs are illustrated.
Runs 0 and 3 provide the two narrowest distributions of In what regards run 0, the situation is expected as this is the case where by far the narrowest distributions of P NS (n) is obtained, see Fig. 3.The fact that for these runs the most probable value of M * G is 2 M ⊙ , that is, the low limit imposed for the maximum mass, means that these runs gather the largest number of soft NS EOSs.As a consequence, NS matter can be strongly compressed, which explains on the one hand the large values obtained in the distributions of n * c and the related ρ * c and, on the other hand, the small values of R 1.4 , Λ 1.4 , R 2.0 and Λ 2.0 .The narrowest (broadest) posteriors of R 1.4 , Λ 1.4 , R 2.0 and Λ 2.0 correspond to run 0 (1).Runs 0 and 2 -4 provide for Λ 1.4 relatively narrow posteriors, localized between the median value and the upper bound of the constraint 190 +390 −120 (at 90% confidence, Abbott et al. 2018) extracted from GW170817.The posterior of run 1 has a long tail whose extremity exceeds by 40% the above mentioned upper limit from GW170817.

The widest distribution of M *
G corresponds to run 1, in agreement with expectations based on the wide distribution of P NS (n) in Fig. 3.As a results, wide distributions characterize all other parameters plotted in Fig. 8, including M DU .The large number of models that allow for massive NSs explains why this is the only run where the c * 2 s distribution is peaked at 1, the largest accepted value.
We have seen in Fig. 3 that an extra condition on m SNM eff; N at 0.16 fm −3 or on (E/A) PNM at a low density leads to NS EOSs stiffer that those of run 0.This explains why, for runs 2 and 4, the probability to have NS with M * G close to 2 M ⊙ is strongly reduced and the M * G -distribution is peaked at larger values.The EOS stiffening gets translated into lower values of n * c and ρ * c and higher values of R 2.0 and Λ 2.0 .
Posterior distributions of M DU are particularly interesting in that in all runs there are models which allow for direct URCA as well as models where direct URCA is not allowed to operate in stable NSs.By far the smallest (largest) number of models which allow for direct URCA corresponds to run 0 (runs 1 and 2).This situation could be anticipated based on the values of E asym (n) in Fig. 1 even if in that figure the largest considered density is only a fraction of the largest value of n * c in Fig. 8. Out of all the models in run 1 that allow for direct URCA, 3.3% accommodate it in NSs with masses smaller than 1.2 M ⊙ .For run 2 the corresponding figure is 1.2%.The variety of distributions of M DU is an extra proof that each and every constraint acts on the density dependence of the symmetry energy and, thus, on NS particle composition.
Medians and 68% CI of marginalized posterior distributions of key NS parameters, including some of those plotted in Fig. 8, are provided in Table 5 in Appendix C.
For comparison, also shown in Fig. 8 are marginalized posteriors corresponding to run 4 of (Beznogov & Raduta 2023), which is roughly equivalent with run 1 here.No posterior distribution is shown for c * 2 s as, with a maximum value of 0.77, the models in run 4 span a domain of values situated outside the X-axis range of our figure.No posterior distribution is shown for M DU , the reason being that direct URCA is not allowed in any of the models in run 4.All CDF distributions are shifted and of different shape with respect to those derived in this paper.The lowest (largest) values of c * 2 s , n * c , ρ * c and P * c (R 1.4 , Λ 1.4 , R 2.0 and Λ 2.0 ) in CDF indicate that these EOS models are considerably stiffer than those built within the non-relativistic theory of nuclear matter.Although not really visible in the figure, the M * G -distribution is noticeably wider for the CDF model than for any of the models in this work.Along with discrepancies in what regards the threshold for direct URCA, this corroborates the idea that the structure of the density functional prevails over constraints posed around n sat .−120 constraint from (Abbott et al. 2018); note that both the median value and the lower boundary sit outside the figure X-range.For M DU the value "−1" corresponds to the models that do not allow this process to operate in stable NSs.For comparison, marginalized posteriors corresponding to run 4 of (Beznogov & Raduta 2023) are illustrated as well (orange shaded areas).Y-axis ranges have been chosen to increase readability.

Correlations between parameters of NM
Correlations among parameters that describe the isoscalar and isovector behavior of the NM, m SNM eff; N and m PNM eff; n are investigated in Figs. 9 and 10, which correspond to run 0 and run 1 in Table 2, respectively.Correlations with n sat and E sat are not included as the tight constraints imposed on these quantities do not allow for any correlation to manifest.
It comes out that some of the correlations manifest themselves irrespective of whether correlations among (E/A) i with i = 2, 3, 4 are accounted for, while other are sensitive to how the constraints from PNM are implemented.
Examples for the first situation are provided by the positive (negative) correlation between K sat and Q sat (K sat and Z sat ); the negative correlation between Q sat and Z sat ; the negative correlation between Q sym and Z sym .As mentioned in Sec.4.3.2,these correlations manifest also in other models, which suggests that they might be universal.A positive correlation exists also between J sym and L sym , and it becomes stronger when correlations among (E/A) i are included.We note that a J sym − L sym correlation shows up also in (Beznogov & Raduta 2023), while it does not in (Malik et al. 2022(Malik et al. , 2023)); this situation is surprising the more considering that the model in (Beznogov & Raduta 2023) is a version of the model previously proposed in (Malik et al. 2022).We also note that, modulo the long tail at low values of m SNM eff; N , both runs 0 and 1 show negative correlations between m SNM eff; N and K sat .This re-sult is easy to understand considering that m SNM eff; N is negatively correlated with C eff and C eff is positively correlated with K sat , see Eq. (A2).Runs 0 and 1 also show negative (positive) correlations between m SNM eff; N on the one hand and Q sat (Z sat ) on the other hand.The latter correlations are straightforward to explain based on the other correlations discussed above.
A spectacular outcome corresponds to correlations between the isoscalar and isovector channels X sat − Y sym , with X, Y = K, Q, Z that exist for run 0 and do not exist for run 1.Notice that no strong cross-channel correlations have been identified by Malik et al. (2022); Beznogov & Raduta (2023), while correlations between Q sat and Q sym appear to be favored in CDF models with strong nonlinear vector ω field contribution (Malik et al. 2023).A possible explanation for the absence of cross-channel correlations in the literature is that correlations among the values a function, which is used as a constraint, takes at different densities, are disregarded.If true, the situation might change if these correlations among the constraints are accounted for.Following the same reasoning, one might speculate that, by accounting for correlations among (E/A) i , a J sym − L sym correlation may appear in the CDF models of Malik et al. (2022Malik et al. ( , 2023) ) and the correlation of Beznogov & Raduta (2023) may become stronger.Alteration of correlation patterns upon inclusion of correlations among (E/A) i will be addressed in detail elsewhere.Until  then, it is certain that a strong model dependence comes from the coupling between isovector and isoscalar channels.Experimental signatures of this "connection" would probably boost our understanding of the density dependence of the nuclear symmetry energy.Other examples of correlations that depend on how the constraints from (E/A) i are implemented are provided by the positive (negative) correlations between m SNM eff; N and J sym , L sym , K sym , Q sym (Z sym ), which exist for run 0 and do not exist for run 1; the correlations m SNM eff; N − J sym and m SNM eff; N − L sym are weak and related one to another by the strong correlation L sym − J sym ; the correlations m SNM eff; N − X sym with X = K, Q, Z can be viewed as the consequence of the correlations m SNM eff; N − K sat and correlations K sat − X sym .Results of this section can be cross-checked with those of Beznogov & Raduta (2023), where a study similar to the one performed here was done using EOSs built within the CDF theory with simplified density-dependent couplings.Note nevertheless that, for all quantities that have been constrained, the values of the constraints differ slightly from those used in this work.Moreover, there is only one case that, modulo these slight differences, is common for our paper (Beznogov & Raduta 2023) and the present work.That is run 4 of Beznogov & Raduta (2023), which is equivalent to run 1 in this paper.The fiducial run of Beznogov & Raduta (2023), to which Figs. 2 and 3 of that paper correspond, additionally accounts for constraints on the pressure in PNM at three values of the density.As such it can be considered somehow intermediate between runs 1 and 0 here.For a complementary discussion about correlation dependence on theoretical framework; effective interaction; size of domains; progressive incorporation of constraints, the reader is referred to Sec.III of Beznogov & Raduta (2023) as well as to Fattoyev et al. (2012); Zhou et al. (2023); Malik et al. (2023).

Correlations between parameters of NM and properties of NSs
Correlations between parameters of NM and properties of NSs are considered in Fig. 11 for run 1.
It comes out that R 1.4 and, to a lesser extent, R 2.0 are positively correlated with L sym .Correlations occur also between R 1.4 and R 2.0 and K sym , but are not shown on this figure.For K sym ≳ −100 MeV, K sym is negatively correlated with P * c .These results agree with the conclusions of Fortin et al. (2016); Zhou et al. (2023), where both phenomenological non-relativistic and relativistic mean field models have been considered, and are expected to manifest also in other models.The rationale is that large values of L sym and K sym engender stiff NS EOSs that prevent matter from being compressed too much.Limited compression translates into low values of n * c and P * c and large values of R 1.4 and R 2.0 , in agreement with panels (a), (b), (c) and (g), respectively.
M DU appears to be negatively correlated with L sym [panel (d)] and K sym [panel (h)].The first of these correlations was previously commented by Fortin et al. (2016).Both correlations are easy to understand and, as before, are expected to be present also in other models.The explanation is that large values of L sym and K sym lead to large values of E sym (n) that prevent β-equilibrated matter from getting depleted from its protons.The fact that none of these correlations are strong is attributable to higher order parameters that are negatively correlated with L sym and K sym , as it is the case of Z sym (see Fig. 10).The strengths of these correlations nevertheless depend on the details of the model.We underscore that, when correlations among E/A in PNM are accounted for, only the K sym − M DU correlation survives.A plausible explanation for this situation is that the domains where the various quantities take values shrink too much too allow correlations to manifest, see Figs. 7 and 8.

SUMMARY
The non-relativistic mean field theory of nuclear matter and the standard Skyrme parametrization of the nucleonic effective interaction have been employed to generate, within a Bayesian framework, models for NM and NS EOS.Several types and different combinations of constraints have been used.In the first place constraints were posed on the four best known NM properties, i.e., saturation density of symmetric matter; energy per nucleon and compressibility of symmetric saturated matter; symmetry energy at saturation.The considered values were those obtained by Margueron et al. (2018), based on a compilation of 35 standard Skyrme interactions frequently used in the literature.Then, the behavior of PNM was controlled by the values that (E/A) PNM takes in a minimum of three and a maximum of four points with densities in the range 0.04 fm −3 ≤ n ≤ 0.16 fm −3 .For (E/A) PNM values obtained within χEFT with NN at N 3 LO and 3N at N 2 LO were used (Somasundaram et al. 2021).Finally, models of NSs were required to comply with the 2 M ⊙ condition on the maximum mass and be causal at the density corresponding to the central density of the maximum mass configuration.The effectiveness of constraints from PNM was tested by accounting for or disregarding the correlations among the values that (E/A) PNM takes at different densities; implementing constraints at three or four densities.Two sets of models also implemented constraints on the values of the nucleon (neutron) Landau effective masses in SNM (PNM) at 0.16 fm −3 ; similarly to (E/A) PNM , the values obtained by Somasundaram et al. (2021) were used.
The wide domain of values associated with the constraint on K sat together with the limited effect the constraints on (E/A) PNM have on SNM get reflected in a significant uncertainty in what regards the high density behavior of E/A (n) in SNM and an expected similitude among the domains of values that the models in each run span for a fixed density.At variance with this, the behavior of PNM differs much from one run to the other.The smallest uncertainties in the behavior of PNM as a function of density correspond to the situations where correlations among (E/A) PNM at various densities have been accounted for.The largest uncertainties correspond to the opposite case of minimal number of constraints on (E/A) PNM .Largely different behaviors of PNM get translated into largely different behaviors of the density dependence of E asym and NS EOS.We note, in particular, that the softest NS EOS, though still compatible, by construction, with the 2 M ⊙ limit, corresponds to the case where correlations among (E/A) i in  PNM were implemented; the largest values of the maximum mass, of the order of 2.5 M ⊙ , were obtained for models in the least constrained run.Large uncertainties in E asym , due to large uncertainties in both SNM and PNM, make that all runs contain models where the (a)symmetry energy increases steeply with the density as well as models where this quantity vanishes at densities a few times larger than the saturation density.Models in the first category allow the direct URCA to operate in stars as light as ≈ 1 M ⊙ , while models in the second category do not allow at all for this fast neutrino emission process.
Conditional probability density distributions P(R|M) show that, for each of our runs, a significant number of models comply with joint mass and radius posterior probability density contours corresponding to PSR J0030+0451 (Miller et al. 2019b) and PSR J0740+6620 (Miller et al. 2021).In contrast, a very small number of models appear to comply with joint mass and radius measurements of the NS in the SNR HESS J1731-34 (Doroshenko et al. 2022); moreover, none of In what regards the combined tidal deformability as a function of mass ratio it turns out that the models in the most constrained run provide values by ≈ 30% higher than the median value extracted from GW170817 (Abbott et al. 2019).At variance with this, models in the least constrained run provide values that span a wide domain; their median value is by 60% − 70% larger than the median value in (Abbott et al. 2019); the 90% upper bound of the CR of models in this run is close to the 90% symmetric upper bound in (Abbott et al. 2019).Overall, our results are fully compatible with GW170817 data.
Marginalized posterior distributions of all NEPs except n sat and E sat ; (E/A) PNM over 0.04 fm −3 ≤ n ≤ 0.16 fm −3 ; m PNM eff; n ; m SNM eff; N differ considerably from one run to the others, which means that every constraint is effectual.Target distributions of K sat and (E/A) PNM at 0.16 fm −3 are never met, which, again, might signal either an incompatibility among constraints or the limited flexibility of the energy density functional.Target distributions of m PNM eff; n and m SNM eff; N are met only when constraints on these quantities are imposed; in all other cases their marginalized posteriors are broad and shifted towards lower values.
Marginalized posterior distributions of various NS properties also differ from one run to the others; the largest (narrowest) distributions are obtained for the least (most) constrained runs (runs 1 and 0, respectively).When compared with equivalent distributions obtained within Bayesian inferences where EOS models are built based on the CDF (Malik et al. 2022;Beznogov & Raduta 2023), considerable discrepancies arise.First, Skyrme-like interactions lead to NS EOS that are softer than their CDF counterparts, which gets reflected in NS with  ones handled via Eqs.( 19) and ( 20) and the "hard wall" ones.The latter constraints can be viewed as post-filtering criteria for the posterior obtained by imposing only the former constraints.Moreover, one call of the log-likelihood requires ≈ 1 ms, while the time needed to compute M * G once is ≈ 1 s (on the same machine).Thus, we have split the calculation into two parts.First, we compute the posterior via MCMC without M * G and c s constraints.While this is as inefficient as before, it is now ≈ 1000 times faster as we do not compute M * G .This makes the first part reasonable in terms of computation time.Then, in principle, we need to post-filter the results by accepting only the models that comply with constraints on M * G and c s .Yet, if we did just that, it would have taken the same (very long) amount of time as running MCMC directly without splitting the constraints.Instead, we have "thinned" the chain by taking only approximately one sample per the integrated autocorrelation length before post-filtering.We aimed for the integrated autocorrelation length of the "thinned" chain to be τ th = 1.0 − 1.2.By doing so we have simultaneously reduced the computational time for post-filtering by a factor of ≈ τ orig and ensured that the "thinned" chain contained only (almost) independent samples."Thinning" also alleviates the issue of low acceptance fraction of the MCMC sampling.The acceptance fraction of the post-filtering itself varied from ≈ 0.8% to ≈ 50% depending on the run.As a result of such two stage "thinning", we can be certain that the final posterior contains only independent samples.
The length of the original chain was chosen such that (i) the final posterior has at least 100 000 samples for all runs except run 0, where we were only able to compute ≈ 32 000 samples.This run had the highest integrated autocorrelation length (τ orig ≈ 10 000) and the lowest acceptance fraction of the post-filtering (≈ 0.8%), (ii) the length is ≥ 200τ orig .This is needed to ensure that the estimation of τ orig is reliable (as mentioned in the emcee documentation, the length of the chain should be at the very least 50τ).Apart from the autocorrelation length analysis, we have performed "bootstrap" tests of the stability of the posterior distributions of the input parameters to ensure proper convergence of the chains.The procedure is the same as described in our paper (Beznogov & Raduta 2023).The results showed no systematic drifts of the posterior.For run 0 (our worst case run in terms of convergence), we have computed the values of 10%, 50% and 90% quantiles for each bootstrap slice.The maximum difference between corresponding values was ≈ 2%.Typical differences were less than 1%.
To additionally verify the correctness of the results, the posterior for the most numerically problematic run (run 0) was calculated independently by means of a nested sampling (Skilling 2004(Skilling , 2006) ) algorithm MLFriends (Buchner 2016(Buchner , 2019)), implemented in the Python package ultranest (v.3.5.7) (Buchner 2021) 2 .After ≈ 3.1 × 10 9 likelihood evaluations we ended up with ≈ 90 000 effective samples of the posterior (we employed region sampling to avoid any possible biases related to step sampling, see ultranest documentation; this approach might be very inefficient, but robust).Again, we have computed the values of 10%, 50% and 90% quantiles and compared ultranest results against emcee results.And again, the maximum difference between corresponding values was ≈ 2.6% and typical differences were less than 1%.Plotted on a figure, the histograms were basically identical.Note that in our "pre-production" tests we have also utilized extensively Nested Sampling Python package dynesty (v.2.1.0)(Speagle 2020;Koposov et al. 2022)  3 .
Thus, we have confirmed that despite above-mentioned difficulties and bi-modality of some parameters (e.g., Figs 6, 9 and 10), our results are reliable.
Another topic worth discussing in more details is the change of the prior due to re-parametrization.To investigate this question we have computed two simplified analogs of run 1.Run 1a is the full equivalent of run 1 except for prior range of σ, which was reduced to [0.1, 1.1] from [0.01, 1.1].Since σ = 0 is a singular point, this change was sufficient to make computation in the effective interactions parameters parametrization feasible.Thus, run 1b is analog of run 1a, but in the effective interaction parameters parametrization.The priors for run 1b were distributed uniformly in the ranges indicated in Table 4.Those ranges were chosen by analyzing with the results of run 1a.
The results are presented on Fig. 12, where runs 1, 1a and 1b are confronted against each other.Let us first compare runs 1 and 1a.They use the same "mixed" parametrization.We can see that reducing the range of σ does not strongly affect the results.C eff and D eff are almost unchanged, while C 0 , D 0 , C 3 and D 3 are much more localized in run 1a.This means that the values of σ close to zero are responsible for long tails of C 0 , D 0 , C 3 and D 3 in run 1 (and other runs as well).This is not surprising, as σ = 0 is not only a singular point of transformations ( 21), but also a degenerate point of energy functionals as discussed in Sec. 2. Comparing runs 1a and 1b, which differ in parametrization and, thus, in priors, one can conclude that re-parametrization has little impact on the posteriors.Indeed, only the σ posterior distribution shows slight difference between the runs, all other parameters are virtually identical.This means that the posterior distributions are strongly dominated by the constraints and are almost insensitive to the priors.This, in turn, justifies our re-parametrization and proves the robustness of the posteriors to the variations of the rather arbitrary chosen priors.
C. POSTERIOR DATA TABLES Information complementary to Figs. 7 and 8 is provided here.
Table 5. Parameters of marginalized posterior distributions corresponding to each run.Considered quantities are NEPs; effective masses of nucleons in SMN and PNM at n = 0.16 fm −3 , respectively; energy per particle in PNM at n = 0.04, 0.08, 0.12 and 0.16 fm −3 ; key parameters of NSs, as discussed in the text; symmetry energy (E sym;0.1 ) and symmetry pressure (P sym;0.1 ) at n = 0.1 fm −3 .Medians (Med.) and 68% CI are shown.The median values and the values of the CI are rounded to 3 and 2 significant digits, respectively; all trailing zeros after the decimal point are removed.

Figure 1 .
Figure 1.Conditional probability density (also known as curve density) plots corresponding to the density dependence of the energy per nucleon in SNM (left column) and PNM (middle column) and asymmetry energy, Eq. (15), (right column).Results of various runs are illustrated on subsequent rows, as indicated on the left panels.Curve density is indicated by colors: dark red (light yellow) corresponds to large (low) densities.Dotted and dashed dodger blue curves demonstrate the median and 90% CR, respectively.Solid black curves mark the envelope of the bunch of models associated with each run.Solid cyan curves mark the envelope of Comparison Set models, cyan dot-dashed curves show their median, see text for details.Light blue shadowed regions mark the uncertainty domain for (E/A) PNM , as computed by Somasundaram et al. (2021).

Figure 2 .
Figure 2. Conditional probability density (also known as curve density) plots corresponding to the density dependence of the nucleon effective mass in saturated SNM (left column), neutron effective mass (middle column) and neutron-proton effective mass splitting in neutron-rich matter with Y p =0.2 (right column).Results of various runs are illustrated on subsequent rows, as indicated on the left panels.The color map shows the curve density.Dotted and dashed dodger blue curves demonstrate the median and 90% CR, respectively.Solid black curves mark the envelope of the bunch of models associated with each run.Solid cyan curves mark the envelope of Comparison Set models, cyan dot-dashed curves show their median, see text for details.Light blue shadowed regions mark uncertainty domains as computed by Somasundaram et al. (2021).

Figure 3 .
Figure 3. Conditional probability density (also known as curve density) plots corresponding to the density dependence of pressure (left column), speed of sound squared (middle column) and proton fraction (right column) in NS matter.Results of various runs are illustrated on subsequent rows, as indicated on the left panels.The color map shows the curve density.Dotted and dashed dodger blue curves demonstrate the median and 90% CR, respectively.Solid black curves mark the envelope of the bunch of models associated with each run.Solid cyan curves mark the envelope of the Comparison Set models, cyan dot-dashed curves show their median.The light blue domains on the right panels correspond to the direct URCA threshold, as predicted by our models; see text for details.

Figure 4 .Figure 5 .−Figure 6 .
Figure 4. Conditional probability density (also known as curves density) plots P(R|M).The color map shows the curve density.Dotted and dashed light cyan curves demonstrate the median and 90% CR, respectively.Solid black curves mark the envelope of the bunch of models associated with each run.The envelope and the median the Comparison Set models are illustrated with cyan solid and dot-dashed curves in the panel "5".Only stable configurations are plotted.Also depicted are joint mass and radius posterior probability density contours corresponding to PSR J0030+0451 (Miller et al. 2019b) (blue contours), PSR J0740+6620 (Miller et al. 2021) (green contours) and to the NS in the SNR HESS J1731-34 (Doroshenko et al. 2022) (orange and red contours); solid and dashed lines stand for 50% and 90% CR, respectively.

Figure 7 .
Figure7.Marginalized posteriors of NEPs, effective nucleon mass in SNM, effective neutron mass in PNM and energy per particle in PNM at densities between 0.04 and 0.16 fm −3 .All runs are considered.When available, target distributions for constrained parameters are also plotted.Y-axis ranges have been chosen to increase readability.As such, some of the very narrow distributions are cut.

Figure 8 .
Figure 8. Marginalized posteriors of selected properties of NSs.Considered are: maximum gravitational mass (M * G ); the central density corresponding to the most massive configuration (n * c ); speed of sound squared (c * 2 s ), energy density (ρ * c ) and pressure (P * c ) at n * c ; radii (R 1.4 , R 2.0 ) and tidal deformabilities (Λ 1.4 , Λ 2.0 ) of NSs with masses equal to 1.4 M ⊙ and 2.0 M ⊙ ; the lowest mass of a NS that accommodates for direct URCA (M DU ).The vertical solid line on the Λ 1.4 plot illustrates the upper bound (at 90% confidence) of the 190 +390−120 constraint from(Abbott et al. 2018); note that both the median value and the lower boundary sit outside the figure X-range.For M DU the value "−1" corresponds to the models that do not allow this process to operate in stable NSs.For comparison, marginalized posteriors corresponding to run 4 of(Beznogov & Raduta 2023) are illustrated as well (orange shaded areas).Y-axis ranges have been chosen to increase readability.

Figure 9 .
Figure 9. Two-dimensional (2D) marginalized posteriors of some of the NM parameters.X sat and X sym parameters are expressed in MeV; m SNM eff; N and m PNM eff; n are expressed as fractions of the nucleon and neutron masses, respectively.The color map indicates the probability density.The light solid and the black dashed contours show 50% and 90% CR, respectively.Results correspond to run 0 in Table2.

Figure 10 .
Figure10.The same as in Fig.9but for run 1 in Table2.

Figure 11 .
Figure 11.2D marginalized posterior distributions for selected NM and NS parameters for which correlations are obtained.The color map indicates the probability density.The light cyan solid and the black dashed contours demonstrate 50% and 90% CR, respectively.Results correspond to run 1 in Table 2. them agrees with constraints from the other two observations above.At the very least, this tension can be attributed to a tension between HESS J1731-34 data and data on (E/A) PNM from χEFT (Somasundaram et al. 2021) or, rather, a limited flexibility of the energy density functional.In what regards the combined tidal deformability as a function of mass ratio it turns out that the models in the most constrained run provide values by ≈ 30% higher than the median value extracted from GW170817(Abbott et al. 2019).At variance with this, models in the least constrained run provide values that span a wide domain; their median value is by 60% − 70% larger than the median value in(Abbott et al. 2019); the 90% upper bound of the CR of models in this run is close to the 90% symmetric upper bound in(Abbott et al. 2019).Overall, our results are fully compatible with GW170817 data.Marginalized posterior distributions of all NEPs except n sat and E sat ; (E/A) PNM over 0.04 fm −3 ≤ n ≤ 0.16 fm −3 ; m PNM eff; n ;

Figure 12 .
Figure 12.Marginalized posterior distributions for the parameters of standard Skyrme effective interactions.Also shown are marginalized posteriors for (C 0 + D 0 ), (C 3 + D 3 ), (C eff + D eff ) on which constraints on (E/A) in PNM and m PNM eff; n act.Runs 1, 1a and 1b.

Table 1 .
Constraints posed on EOS models.E sat and K sat represent the energy per particle and compression modulus of symmetric saturated matter with the density n sat ; J sym stays for the symmetry

Table 2 .
Sets of constraints, other than those on n sat , E sat , K sat , J sym , Table 1 acts on a different number of input parameters of the model.As such, constraints on n sat , E sat , K sat get translated into constraints on σ, C 0 , C 3 and C eff ; constraints on (E/A) PNM lead to constraints on σ, (C 0 + D 0 ), (C 3 + D 3 ), (C eff + D eff ); constraints on m SNM eff; N lead to constraints on C eff ; constraints on m PNM eff; n lead to constraints on (C eff + D eff ); constraints on J sym act on D 0 , D 3 , C eff , D eff , σ; notice that J sym depends on C 0 and C 3 through n sat .

Table 3 .
Domains of the priors.

Table 4 .
The ranges of the priors for run 1b.