Determination of the Equation of State from Nuclear Experiments and Neutron Star Observations

With recent advances in neutron star observations, major progress has been made in determining the pressure of neutron star matter at high density. This pressure is constrained by the neutron star deformability, determined from gravitational waves emitted in a neutron-star merger, and measurements of radii for two neutron stars, using a new X-ray observatory on the International Space Station. Previous studies have relied on nuclear theory calculations to provide the equation of state at low density. Here we use a combination of 15 constraints composed of three astronomical observations and twelve nuclear experimental constraints that extend over a wide range of densities. Bayesian Inference is then used to obtain a comprehensive nuclear equation of state. This datacentric result provides benchmarks for theoretical calculations and modeling of nuclear matter and neutron stars. Furthermore, it provides insights on the composition of neutron stars and their cooling via neutrino radiation.


I. INTRODUCTION
With masses that can be larger than twice the mass of the sun and with radii of approximately 13 km, neutron stars (NS) and their observed properties raise some compelling questions [1].Are pions, strange particles or quark matter important for understanding NS matter and its internal pressure, which supports a NS and prevents it from collapsing into a black hole?If that matter is nucleonic, what roles do repulsive short-range three-neutron or fourneutron forces play in supporting the star?Understanding the connections between the pressure and the composition and the structure of the stellar matter clearly constitutes a key objective for both astrophysics and nuclear physics [2].To better understand the connections between nuclei and neutron stars [3], we combine nuclear measurements and astrophysical observations to obtain an equation of state (EOS) of nuclear matter.
The breakthrough in equation of state research in the past 5 years has been in Astronomy.First the era of Multi-Messenger astronomy, involving simultaneous observations of a merging binary NS system with a wide array of astronomical instruments [4], has enabled more detailed studies of merging neutron stars.Then, these two merging neutron stars (GW170817) observed by the Gravitational-Wave Observatories LIGO/Virgo Collaboration provided information on the deformability of the neutron stars [5].Next, the measurements of the radii of two pulsars with very different masses by the Neutron Star Interior Composition Explorer (NICER) [6][7][8][9] provided constraints on the NS mass-radius correlation.Both the deformability and mass-radius relationship constrain the pressure-density relationships or EOS of NS matter above twice "saturation" density, n > 2n 0 [10,11], where n 0 ≈ 0.16 nucleons/f m 3 (≈ 2.6 × 10 14 g/cm 3 ) is the density at the core of any heavy nucleus.To augment these astrophysical constraints, we add constraints from nuclear physics experiments [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] obtained in the past two decades.When combined, these constraints provide a consistent description of the EOS that describes matter in nuclei and neutron stars at densities of 0.5n 0 < n < 3n 0 .
The methodology used in the current work is composed of three parts which will be described in detail below.Using these constraints, it obtains an EOS with uncertainties that enables consistent descriptions of the observables obtained from nuclear experiments and a neutron star model that uses the same EOS function to calculate neutron star properties for comparison to astronomical observations.A Bayesian analysis framework is used to infer the EOS parameters and to quantify the uncertainties of the results.The resulting EOS provides benchmarks for the microscopic nuclear theory, insights regarding the nature of strongly interacting matter in the outer core of the NS, its composition and the onset of direct Urca cooling processes.In a different context, a comprehensive EOS would also yield insights about the collapse of supernovae and neutrino emission.

II. EQUATION OF STATE FUNCTION
We assume nucleonic matter in nuclei and in neutron stars at zero termperature share a common nuclear EOS that can be chosen to be the energy per nucleon.We parameterize the EOS as a function, ϵ(n n , n p ), of the neutron and proton number densities n n and n p , where n=n n +n p is the total nucleonic density and δ = nn−np n is the asymmetry parameter.To describe the evolution of the EOS with density and asymmetry, δ, ϵ(n n , n p ) =ϵ SN M (n)+S(n)δ 2 consists of two terms where ϵ SN M (n) is the energy per nucleon for symmetric nuclear matter with δ=0 and S(n)δ 2 is the symmetry energy term needed when δ ̸ = 0.This allows the EOS to describe nuclear matter with any neutronproton composition.Such division is useful because the asymmetry δ can be experimentally controlled to enhance the sensitivity to either the energy of symmetric nuclear matter or to the symmetry energy term.Detailed description of the equations of state used in the Bayesian analysis can be found under Methods.

III. CONSTRAINTS
Most experimental and astronomical observables are sensitive to the EOS over a limited range of density.For example, radii, masses and deformabilities of neutron stars largely reflect the pressure inside the neutron stars at densities of n > 2n 0 [10,28].On the other hand, nuclear masses and neutron skins of neutron-rich nuclei are mostly sensitive to the nuclear matter EOS at (2/3)n 0 [27].In addition to these two density regions, observables constructed from the distribution patterns of particles emitted in nucleus-nucleus or heavy ion collisions (HIC) have provided constraints on the EOS at densities range 0.22 < n/n 0 < 4.6 [12,13,21,22,29,30].The diversity of the available data is given in Table I and also illustrated in Fig. 1 which shows the effective density region of mostly nuclear physics constraints up to about 3n 0 .Constraints prefaced with "HIC" are obtained from heavy ion collision experiments.Brief description of each of the constraints will be discussed below.For astrophysics, the relevant density range is above 2n 0 .Recall that the nuclear equation of state is parameterized into two terms, the symmetric matter term and the symmetry energy term that accounts for the imbalance of neutrons and protons.We separate these nuclear physics constraints into ones that are sensitive to the symmetry energy and others that probe the symmetric matter EOS.31,32] and a symmetric matter constraint on KSAT obtained from Giant Monopole Resonance (GMR) studies [14].The purple horizontal arrow at the top of the figure indicates the range of density from symmetric matter constraints coming from heavy-ion collision studies [12,13].The remaining terrestrial constraints for nuclear symmetry energy or pressure are labelled in blue and positioned at their sensitive densities [27].The colored band indicates the approximate dependence of pressure on density for pure neutron matter from the findings of current work.

A. Symmetric Matter Constraints
Certain properties of the symmetric nuclear matter such as the saturation density, n 0 , and saturation energy, E SAT are reasonably well determined.We adopt the values of n 0 ≈ 0.16 fm −3 and E SAT ≈ -16 MeV from [31,32].For incompressibility parameter, K SAT , we use the values of 230±30 MeV following Ref.[14].Even though these three parameters are labelled in Fig. 1 at n 0 , only the Taylor expansion coefficient, K SAT , is used as a constraint in the Bayesian analysis.
Above saturation density, measurements of collective flow from energetic Au+Au collisions have constrained the EOS for symmetric matter at densities ranging from n 0 to 4.6n 0 [12,13].The density range is indicated by a horizontal magenta arrow on the upper right corner of Fig. 1.The pressure-density plot for symmetric matter in Fig. 2b shows these two flow constraints [12,13].The black slanted-hatched region (HIC(DLL)) shows the constraint on the symmetric matter EOS published in Ref. [12] from the analysis of both collective transverse and elliptical flow data measured at incident energies of 0.2-10 GeV/u.The extraction and model dependence of the HIC(DLL) constraint are discussed in the text of [12] and its online supplemental material.The magenta slanted-hatched region (HIC(FOPI)) located at lower densities results from an independent analysis of a different set of elliptical flow measurements at 0.4 to 1.5 GeV/u [13].Known sources of theoretical uncertainty are modeled and the symmetric nuclear pressure plotted in panel b reflects this uncertainty estimation.Both contours agree very well in the region around 2n 0 where they overlap.We take the pressure values at 2n 0 as a constraint in our analyses.In view of availability of the new and better measured data [33][34][35], more robust constraints on symmetric matter are expected in the near future.

B. Symmetry Energy Constraints
In the past decade, many studies have been conducted to extract the symmetry energy and symmetry pressure [27] mostly at lower densities.Precise symmetry energy constraints shown in Fig. 2d have been obtained at about (2/3)n 0 from nuclear masses using density functional theory in two independent analyses involving different nuclei labeled as Mass(Skyrme) [18] and Mass(DFT) [19].In this density region we also have precise constraints obtained from the energies of isobaric analogue states [20] indicated by a data point labeled as IAS.Measurements of the dipole polarizability of 208 Pb provide a constraint at a lower density of n = 0.05 fm −3 , which is labeled as α D [15].Values for the neutron skin thickness and the slope of the symmetry energy, L(2n 0 /3), at (2/3)n 0 have been published for the PREX-II experiment from which we obtain the pressure P (2n 0 /3).This is labeled as PREX-II [16,17] in Fig. 2a.
Heavy-ion collisions have probed the symmetry energy and pressure at densities far from (2/3)n 0 .At incident energies below 100 MeV per nucleon, lower densities are probed when matter expands after initial impact and compression of the projectile and target.There, experimental observables primarily reflect the symmetry energy at sub-saturation densities with n << n 0 [21,22].Constraints from isospin diffusion in Sn+Sn collisions, labelled as HIC(isodiff) [21] provide a constraint at n/n 0 = 0.22 ± 0.07.Ratios of neutron and proton energy spectra in central collisions, labelled as HIC(n/p) provide a constraint on the symmetry energy at n/n 0 = 0.43 ± 0.05 [22].
Higher energy (>200 MeV/u) central HIC probe the EOS at n > n 0 .The HIC(n/p flow) data point comes from the elliptical collective flows of neutrons and hydrogen nuclei in Au+Au collisions [24][25][26].HIC(π) comes from a recent measurement of the spectral ratios of charged pions in very neutron-rich 132 Sn+ 124 Sn collisions and nearly symmetric 108 Sn+ 112 Sn collisions [23].The use of ( 132 Sn) and ( 108 Sn) radioactive beams enhanced the asymmetry variation and the experimental sensitivity to the symmetry energy.The uncertainty of the pion constraint reflects significant uncertainties in the difference between neutron and proton effective masses and in the contributions from non-resonant pion production [23].Extractions and discussions of the model dependent uncertainties in the constraint marked as HIC(FOPI), HIC(n/p) and HIC(n/p flow) can be found in refs.[13,23,24].Finally, we note that the uncertainties in the sensitive densities listed in Table I are too small to significantly influence the final EOS and therefore not used in the Bayesian analysis.

C. Constraints from astronomical observations
When two neighboring neutron stars begin to merge, the gravitational field of each neutron star induces a tidal deformation in the other.The influence of the EOS of neutron stars on the gravitational-wave signal during the inspiral phase is encoded in the dimensionless, neutron-star tidal deformability, Λ [5].In the inset of Fig. 2c, the value of Λ for 1.4M ⊙ obtained from the binary neutron star merger event of GW170817 [11] is plotted.The observation of a later merger event (GW190425) did not improve the accuracy of Λ due to poorer observation conditions [36].

IV. DESCRIPTION OF THE NEUTRON STAR MODEL CALCULATIONS
We adopt the NS model described in Ref. [39] to calculate properties of the neutron star, such as the deformability, mass and radius for each EOS in the prior distribution, by solving the Tolman-Oppenheimer-Volkoff (TOV) equation [40,41].In particular, we focus on the impact of the EOS on the outer core of the neutron star where neutrons comprise the bulk of the matter and inclusion of small admixtures of protons, electrons and muons is required to attain β-equilibrium.We describe very low density matter in the NS crust with the crustal EOS of Ref. [42].The crust is connected to the outer core with a spline function for smooth connection.Although this connection is simplistic, the uncertainty introduced by this crustal EOS has been estimated to be negligible at densities greater than 0.5n 0 [39].For simplicity, we include no sharp phase transitions to occur within the neutron star core.Both the inner and outer core can be described by the same equation unless the speed of sound in outer core EOS begins to exceed the speed of light.If this occurs, we transition it to the stiffest possible EOS where speed of sound equals to speed of light.When this happens, the EOS can be very repulsive.Such transition occurs mostly above 3n 0 where our study provides insights into the EOS at densities of 0.5-3n 0 and where the current experimental constraints are relevant.More experimental constraints at higher density will be forth coming in the near future.[43].

A. Priors
In this data driven approach, we constrain the Taylor expansion coefficients of the equation of state and the posterior EOS distributions solely by the 15 discrete constraints obtained from astronomical observations and nuclear physics experiments.These constraints are listed in Fig. 1 (n 0 and E SAT are fixed parameters, not constraints) and Table I.The availability of so many constraints precludes the dominance of a single experiment or observation and provides data over a wide range of densities.This approach is different from previous studies that rely on chiral effective field 1XPEHUGHQVLW\n (n0) Density dependence of the pressure for neutron star matter.The purplish-blue area is the EOS obtained in Ref. [10] based on astronomical constraints using only non-parametric EOS.The hatched area is from Ref. [44] with astro and high density heavy ion collision constraints at high energy and χEFT as priors at low density.theory (χEFT) to describe the EOS at low density [2,44,45].To apply the Bayesian analysis, we have to first create priors which is discussed in more details in Method section.Here, we want to emphasize that the general principle is to adopt a wide range of parameters for the EOS so that the prior distributions more than encompass all the experimental constraints.

B. Posteriors
Bayesian inference is used to simultaneously determine all the parameters used in the EOS (see Methods for details).The posterior distributions of the EOS calculated from Bayes theorem are shown in Fig. 2. The 68% and 95% CI of the posteriors are represented by the dark green and light blue shaded contours.In each case, the posterior uncertainties are much narrower than that of the prior (indicated by the blue dotted curves) and nearly all data fall within the 95% CI of the posterior regions.We obtain a value of 530 +115 −138 for the deformability of a 1.4M ⊙ NS which is higher but still within the uncertainties of the value provided by LIGO [11].The analysis also gives a radius value of 12.9 +0.4 −0.5 km for a 1.4M ⊙ NS.
The posterior of the density dependence of the symmetry energy shown in Fig. 2d is similar to that obtained in [27] but the most probable values are somewhat less reflecting the additional astrophysical constraints included in present work.When compared to χEFT, the experimental constrained EOS tends to be stiffer as shown in Fig. 2e.
In Fig. 2f, the purplish blue shaded area shows the EOS for neutron star matter at 90% CI obtained by employing only astronomical constraints [10].The uncertainties of our extracted EOS is narrower especially at the low density region where we employ experimental data.At high density the difference in uncertainties may come from the different priors used in the analysis.Ref. [10] adopts non-parametric EOSs as priors while the current work uses parametric priors based on an expansion widely used in nuclear physics.The recent constraints [44] obtained by incorporating the χEFT at low density and Au+Au collective flow data from [13,24,26] are shown as the hatched regions (black dotted lines) corresponding to 68% (95%) CI.The differences between the current work and Ref. [44] most likely arise from its strict reliance on the χEFT EOS below 1.5n 0 which produces a softer EOS.More comparisons to χEFT results are discussed below and in the methods section.

VI. IMPLICATIONS A. Benchmarking microscopic theories
With recent advances made in the development of χEFT especially in the quantification of the uncertainties [46], it has been extended to 2n 0 in order to help describe neutron star properties.Since our analysis does not include the use of χEFT, it is possible to compare our results with the extended χEFT EOS at these higher densities as shown in Fig. 2e.Such comparisons provides a useful benchmark to assess how well the short-range part of the nuclear interaction is modeled [47].The magenta dotted lines represent 95% CI and the hatched area correspond to 68% CI of this χEFT EOS.(For reference, the experimental data are plotted in Fig. 2d).There are broad agreements between the extracted EOS and the χEFT calculations below 2/3 of the saturation density.The deviations from the experimental symmetry energy increase with increasing density and the χEFT calculations tends to be softer than our extracted EOS and falls below the Urca threshold as discussed in next section.More detailed comparisons to χEFT predictions are described in the Methods section.

B. Composition and the Urca cooling of neutron stars
Data on the pressure of both symmetric matter and neutron rich matter provides crucial insight into the composition of dense matter.The symmetry energy is very important to understand the composition and dynamics of matter within neutron stars because it contributes to the chemical potentials of the particles that compose the stellar matter.A large symmetry energy may increase the fraction of protons, muons, hyperons or other particles that are present in dense matter.In β-equilibrium, the proton fraction (protons per baryon) y p increases with increasing symmetry energy.(Details of calculating y p is given in Methods section).If y p is large enough, the Urca process of rapid neutrino emission may quickly cool a neutron star [48].Many isolated neutron stars are observed to cool relatively slowly where the Urca process is likely not operating [49].However, some neutron stars, perhaps the more massive ones, are observed to cool quickly consistent with Urca [50].
1XPEHUGHQVLW\n (n 0 ) The density difference between central density, ncen of a neutron star of mass M and the lowest density, nUcra where the Urca process is allowed.If this density difference is greater than zero the star may potentially cool quickly.
In the Urca process a neutron beta-decays followed by electron capture on a proton with the net effect of radiating a νν neutrino pair that cools the star.To conserve both momentum and energy in these weak interactions one needs the Fermi momenta of protons k p F , neutrons k n F and electrons k e F to satisfy k p F + k e F > k n F .At a given density n, a symmetry energy above the black solid line in Fig. 2e will have a large enough y p to satisfy this condition so that the Urca process is potentially allowed.Note that the composition, shown in Fig. 3a, is a new result made possible by combining symmetric nuclear matter data and symmetry energy data from heavy ion collisions.
When the central density n cen of a neutron star of mass M (M ⊙ ) exceeds the Urca threshold density, n U rca , the star can potentially cool quickly via the Urca process.Fig. 3b shows the mass dependence of (n cen − n U rca ), we find that Urca cooling is likely for stars with mass larger than 1.8M ⊙ .Note, neutron star cooling depends on possible super fluid pairing gaps in addition to y p , see for example [51].Furthermore, matter in the deep interior of a neutron star may not mainly consist of nucleons.In this simple discussion, we do not include such effects.

VII. SUMMARY
In summary, we perform a Bayesian analyses of 12 nuclear experimental constraints with 3 complementary constraints from astronomical observation to obtain a consistent understanding of the equation of state of neutron rich matter at number densities from about half to three times saturation density.Our method eliminates over-reliance on the data of one particular experiment or observation or theory.Within the density range we study, the data do not appear to require the existence of a phase transition inside the core of neutron stars nor do they rule them out.We adopt an equation of state constructed from metamodeling in the form of Taylor expansion around saturation density up to fourth order.The extracted values of the expansion coefficients allow construction of the equation of state for symmetric and asymmetric matter including neutron star matter.With the neutron star equation of state we obtain a radius of 12.9 +0.4 −0.5 km and a deformability of 530 +115 −138 for the canonical neutron star with 1.4 solar mass.The tension between the posteriors and the data nominally reveal areas for experimental or observational improvements as in the need for additional measurement of the tidal deformability of the neutron stars and more precise heavy ion data on symmetry energy and pressure above saturation density [43].Such new data are expected in the near future.

SUPPLEMENTARY MATERIAL A. Parameterization of the equation of state and the neutron star model
Neutron star radii and deformabilities combine to provide constraints on matter at densities of 2 < n/n 0 < 4 in the outer core of a cold neutron star [8,9,11] where we assume the hadronic component of matter is largely nucleons (neutrons and protons).Experimental observables that probe the EOS via nucleus-nucleus collisions provide constraints on nucleonic EOS functionals that can be used to describe such densities [12,13,23,24].To perform these calculations, we express the dominant hadronic component of the EOS in terms of the energy per nucleon ϵ(n, δ) or pressure, P = n 2 ∂ϵ/∂n, where, n = n n + n p is the total density, n n and n p are the neutron and proton number densities and δ ≡ (n n − n p )/n is the asymmetry.
To second order in δ, the energy ϵ(n, δ) can be written as [52], where ϵ(n, 0), is the energy per nucleon of symmetric matter with equal neutron and proton densities (i.e.δ = 0) and S(n), is the symmetry energy of pure neutron matter where δ = 1.
We calculate the symmetric matter EOS (ϵ(n, 0)) and the density dependence of the symmetry energy (S(n)) using the Metamodeling ELFc formalism of ref. [53].In metamodeling, both ϵ(n, 0) and S(n) contain kinetic energy, (K.E.) and mean-field potential terms, (U), as well as effective mass (m * ) corrections that reflect the non-localities of the nucleonic mean field potentials.These effective mass corrections appear in the kinetic energy term where the mass term is replaced by effective masses ( Away from n = 0 where the kinetic energy has a branch cut singularity, the metamodeling EOS is dominated by values and derivatives of ϵ(n, 0) and S(n) at saturation density n 0 ≈ 0.16 nucleons/fm 3 which are labeled as follows, where x = (n − n 0 )/(3n 0 ).The exact formulae for these coefficients and how they are combined with additional correction terms to obtain the EOS down to n = 0 can be found in Ref. [53].By using an EOS of the Meta-modelling form [53], instead of more customary Skyrme [14] or Relativistic Mean Field (RMF) [54] density functionals, which are commonly used to describe nuclei and nuclear matter, we can reduce the model dependent correlations in asymptotic expansion coefficients listed in Eqs. ( 3) and ( 4) that are characteristic features of the Skyrme and RMF functionals [53].This allows the influence of the experimental and astrophysical constraints to be modeled more flexibly and clearly.Even though the nucleon effective mass m * /m 0 is an adjustable parameter in the Meta-modeling approach [53], we assign nominal values of m * /m 0 = 0.7 and m * s = m * v where m * s and m * v are the isoscalar and isovector effective masses, respectively.We adopt the metamodeling EOS for nuclear matter in the outer core of a neutron star which is composed of neutrons, protons, electrons and muons, and is β-equilibrated at a temperature close to zero MeV.Consistent with our experiment and observation driven approach, we do not introduce a theoretical model for a different form for matter in an inner core of a neutron star, but look to see if strong signatures for its presence can be observed by changes in the EOS at high densities which might emerge in the analyses.We find none.We do ensure when we extend the EOS to higher densities that the speed of sound, c s (=c ∂P/∂ϵ) remains less than the speed of light, c.
To calculate NS radius, we need to extend the EOS from the core to the crust.We adopt a commonly used EOS for the solid outer crust from Ref. [42] and extrapolate this through the inner crust region of 0.3n T < n < n T via spline interpolations.Here, n T = −3.75× 10 −4 L + 0.0963 fm −3 is the crust to core transition density.In [39] the uncertainty of this crustal EOS description is estimated to be small compared to the present uncertainty of the EOS in the core of the neutron star.Using these EOS expressions, NS properties such as deformability, masses and radii are obtained by solving the Tolman-Oppenheimer-Volkoff (TOV) equations [40,41].More details of the neutron star model and its uncertainties are provided in [39].

B. Constraints
The constraints that are used in the Bayesian analysis have been discussed in detail in the main text.

D. Priors
Bayesian analysis is used to constrain the relevant Taylor expansion coefficients defined in Eqs.(3,4) that govern the EOS at high densities.Table II lists the range of parameters used to generate the priors that will cover most of the experimental and astrophysical constraints listed in Table I.For comparison, the corresponding posterior values are also listed in Table II.Here, we give the rationale in choosing the ranges of the parameters listed in Table II using experimental information as guidance when possible.The general principle is to allow each parameters with a very wide range to ensure that the priors encompass as much of the experimental constraints listed in Table I as possible.
The ranges of the symmetry energy parameters of Eq. ( 4) are guided by the comprehensive study in Ref. [27].However, we do not have a corresponding comprehensive analysis for the symmetric matter EoS parameters to provide guidance regarding the Taylor expansion parameters appropriate for Eq. ( 3).Together with the well-known values of n 0 , E SAT and K SAT , the existence of two consistent constraints at 2n 0 on the symmetric matter EoS provides ample guidance on the EoS up to 2n 0 .Preliminary analyses of recent experimental results at high collision energy and therefore higher density from Beam Energy Scan [34,35] and E895 [55,56], suggest that the HIC(DLL) constraint on the pressure at high density may have to be re-evaluated.We assess both the existing HIC(DLL) constraint and new but still preliminary analysis [57] and concluded that both analysis agree that the symmetric matter pressure at 4n 0 does not exceed 300 MeV/fm 3 .We include this knowledge in the prior distribution by selecting values for Z SAT so that the symmetry pressure P SNM (4n 0 ) is uniformly distributed from 0 − 300 MeV/fm 3 , nearly twice the published HIC(DLL) standard deviation upper limit of 156 MeV/fm 3 [12].Effectively, this means that Z SAT prior is not independent but is correlated with the priors for other Taylor coefficients, especially Q SAT , in the expression for the symmetric matter EoS.Thus, P SNM (4n 0 ) appears in Table II   To illustrate the influence of neutron star on the priors, 5000 randomly chosen priors based on the parameter ranges listed in Table II are shown as grey curves in Fig. 4 where the left panel shows the density dependence of energy per particle of the pure neutron matter (commonly used in nuclear physics community) while the right panels shows the corresponding pressure as a function of number density (more commonly used in astrophysics community).The TOV equations are solved and EOS which do not lead to a maximum neutron star mass of at least 2.17M ⊙ are rejected before applying the Bayesian analysis.(Lowering the least maximum neutron star mass to 1.8M ⊙ does not change the results/conclusions significantly.)The neutron star requirement removes very stiff and very soft EOS especially those with very small or negative pressure at high density, and greatly reduces the number of priors to about 5%.They are shown as blue curves in Fig. 4. Since these are the priors used in the Bayesian analysis, they will be referred to simply as "priors" from here on.

E. Method
Bayesian inference is performed to generate the posterior distribution, which is governed by various expansion coefficients defined in Eqs.(3,4) and the experimental and astronomical measurements.The posterior distribution is given by Bayes theorem [58], which can be written as, In this equation, P (EOS) is the prior and L(i th Constraint) is the likelihood, defined as the probability of observing the experimental and astronomical results assuming that a given set of EOS expansion coeficients is the perfect description of nuclear matter.In this analysis, for an observable O, which the measurements/observations constrain to have a mean of x 0 and a standard deviation σ equal to the error, the likelihood of it having a true value of x is given by a Gaussian exp[−(x − x 0 ) 2 /(2σ 2 )].For gravitational wave constraints with asymmetric uncertainty, an asymmetric Gaussian is used.The latter is constructed from two Gaussian functions, having the same mean value but different widths and scaled to match the height of each other at the mean.These are used to describe the two halves of the distribution that lie on opposite sides of the mean value.
The marginalized posterior distributions are calculated by distributing each parameter on a histogram, weighted by a factor from Eq. (5).A total of 4,500,000 EOS were sampled in two stages.In the first stage, 1,500,000 EOS were sampled uniformly within ranges in Table II.This small set of expansion coefficients informed us of the range of preferred values.The remaining 3,000,000 expansion coefficient sets were sampled uniformly within 99% confidence intervals (CIs) from the marginalized distributions from the initial EOS set to create the final posterior distribution.This two-step process allowed us to more efficiently sample the parameter space that is relevant.The CIs for the energy and pressure values for the EOS, as well as the values of various observables directly related to the nuclear EOS, were calculated by weighing each prediction with the corresponding factor.MeV R(1.4M ) 13.1 +1.5 1.5 12.9 +0.4 0.5 km (1.4M ) 518 +522 261 530 +115 138 FIG. 5. Marginalized probability distributions for EOS parameters.The off-diagonal plots show the pair-wise probability distributions for S0, L(0.67n0), L(n0), L(1.5n0), R(1.4M⊙) and Λ(1.4M⊙).The two red contours, from inside to outside, correspond to 68% and 95% confidence intervals, respectively.The diagonal plots show the marginalized probability distribution of each parameter correspond to posterior after all experimental constraints are applied.The two outer red vertical dashed line shows the 68% confidence interval and the single solid red vertical line shows the median (50 th percentile).The table on the upper right corner of the figure shows the median ± 68% confidence interval values.

F. Posterior distributions
The Bayesian analysis is performed with parameters allowed to vary within ranges listed in Table II.The posterior of EOS parameters and predictions for the radius and deformability of the 1.4M ⊙ neutron star are also listed in the same Table .In the same table, we listed the predictions for the Taylor expansion coefficients for the symmetry energy at the density of 0.1 fm −3 where the experimental constraints are most robust.Parameter is defined as variables that EOS takes as input and prediction is defined as everything else that requires calculation to be performed to get the values.We made this distinction for the tables for clarity, but they will all be called parameters in other sections for brevity.The posterior distribution, which describes the probability of the i th parameter having value m i given the constraints listed in Table I, is calculated from Bayes theorem.Fig. 5 shows the corner correlation plots for six parameters: symmetry energy at saturation density, S 0 , slopes of the symmetry energy at 0.67n 0 , n 0 and 1.5n 0 , and the radii and deformability of a neutron star with the nominal mass of 1.4 M ⊙ .In red, we show the 1σ and 2σ contours of these correlations.Because we construct the prior EOS using meta-modeling, we avoid any apriori correlations between the slope of the symmetry energy at different densities and other neutron star properties.The correlations in these plots reflect strong connections between the EOS at specific densities and neutron star radii and deformabilities.Nonetheless, both the radii and deformability of the neutron star are correlated with the slope of the symmetry energy at 0.67n 0 , n 0 and 1.5n 0 as well as with each other, consistent with previous studies.Along the diagonal are the posterior 1D plots for the parameters with the solid vertical lines corresponding to the median (50 th precentile) and the dashed vertical lines corresponding to the 68% (1σ) CI.The posterior values with 1σ uncertainties are also listed in Table II.In each case, the posterior uncertainties are much narrower than that of the prior.In Fig. 6, the light blue and dark green shaded regions represent the 95% and 68% CI of the posteriors, respectively.

G. Comparison to chiral effective field theory
The EOS extracted in this work provides the energy per nucleon of symmetric matter, ϵ(n, δ = 0) and the symmetry energy term, S(n).The sum of these two terms yields ϵ(n, δ = 1), the EOS for pure neutron matter.Unlike previous studies [44,59,60], our methodology does not require theoretical EoS.As the Chiral Effective Field Theory (χEFT) has become a popular choice to describe neutron star properties, it is interesting to compare our result with the theoretical calculations.In Fig. 6, we show comparisons of the constrained EOS to the Chiral Effective Field Theory (χEFT) calculations with the 500 MeV momentum cutoff [31,46,61] for the symmetry energy term of the EOS (left panels), symmetric nuclear matter term (middle panels) and for pure neutron matter (right panels) both in terms of pressure vs. density (upper panels) and energy per nucleon vs. density (lower panels).Note that there are different implementations of χEFT.The one used here has a truncation that is fully quantified as shown by the solid magenta contours corresponding to 95% CI and the magenta hatched region corresponding to 68% CI.
In general, the extracted EOS is broadly consistent with χEFT calculations especially at densities below n 0 .For the symmetric nuclear matter, the agreement is rather good up to 1.5n 0 even though it is slightly offset from the saturation energy and saturation density as noted in [61].To compare the symmetry term of the EOS, the symmetry energy (pressure) of the χEFT EOS is obtained by subtracting the energy (pressure) of the symmetric matter from that of the pure neutron matter [31,46,61].This leads to large uncertainty bands in (Fig. 6a) and (Fig. 6d).The EOS derived from χEFT is in general softer, consistent with smaller predicted NS radii [44].It is interesting to note that due to the softness of the χEFT, Urca process will be mostly disallowed in NS matter utilizing χEFT as EOS as shown in Fig. 6d where the black line indicates the Urca threshold.(Bottom) Same as top panels except the y-axis is energy per nucleon.The solid green (light blue) region is the 68% (95%) CI of the posteriors.The magenta (hatched) contours represent the calculations from the Chiral Effective Field Theory (χEFT) [31,46,61].

FIG. 1 .
FIG. 1. Constraints from nuclear experiments and astronomy observations with their corresponding sensitive densities.The red horizontal arrow indicates the high density region (> 2n0) probed by neutron star observations.The magenta text indicates chosen values at ESAT ≈ −16M eV for the symmetric matter EOS located the saturation density n0 ≈ 0.16 fm −3 following[31,32] and a symmetric matter constraint on KSAT obtained from Giant Monopole Resonance (GMR) studies[14].The purple horizontal arrow at the top of the figure indicates the range of density from symmetric matter constraints coming from heavy-ion collision studies[12,13].The remaining terrestrial constraints for nuclear symmetry energy or pressure are labelled in blue and positioned at their sensitive densities[27].The colored band indicates the approximate dependence of pressure on density for pure neutron matter from the findings of current work.

FIG. 2 .
FIG. 2. Constraints on nuclear equation of state and neutron star properties.All symbols in these panels are data discussed in the text and listed in Table I.The blue dotted lines are 95% CI boundaries of the prior distributions.The dark green and light blue shaded regions represent 68% and 95% CI boundaries of the posterior distributions.Panels: (a) Density dependence of the symmetry pressure.(b) Density dependence of symmetric nuclear matter pressure.(c) Mass-Radius correlation plot with the NICER results.The insert shows the mass-deformability correlation plot.The symbol represents the experimental deformability value for 1.4M⊙.(d) Density dependence of the symmetry energy.(e) Density dependence of the symmetry energy compared to χEFT calculations (magenta).The onset of Urca cooling occurs above the solid black line.(f)Density dependence of the pressure for neutron star matter.The purplish-blue area is the EOS obtained in Ref.[10] based on astronomical constraints using only non-parametric EOS.The hatched area is from Ref.[44] with astro and high density heavy ion collision constraints at high energy and χEFT as priors at low density.

FIG. 3 .
FIG. 3. Predictions on the proton fraction and the direct Urca cooling from our analysis.(a) The density dependence of the proton fraction.The onset of Urca cooling occurs above the solid black line.(b)The density difference between central density, ncen of a neutron star of mass M and the lowest density, nUcra where the Urca process is allowed.If this density difference is greater than zero the star may potentially cool quickly.

FIG. 4 .
FIG.4.Prior EOS (grey curves) for pure neutron matter and priors that form stable neutron stars with maximum neutron star mass greater than 2.17M⊙ (blue curves).The y-axis is energy per nucleon in the left plot and pressure in the right plot.
TABLE II.(Top) Ranges of Taylor expansion coefficients used to generate priors and posteriors.(Bottom) Predictions forTaylor expansion parameters at n = 0.1 fm −3 .The parameters of the first group, KSAT, QSAT, S0, L, Ksym, Qsym, Zsym, and PSNM(4n0) are varied uniformly within the ranges listed under priors.The symmetric matter constraint at supra-saturation density dictates that ZSAT values are strongly correlated with QSAT.

FIG. 6 .
FIG.6.(Top) Left, middle and right plots show the symmetry pressure, pressure of symmetric nuclear matter, pressure of pure neutron matter, respectively, as a function of number density n. (Bottom) Same as top panels except the y-axis is energy per nucleon.The solid green (light blue) region is the 68% (95%) CI of the posteriors.The magenta (hatched) contours represent the calculations from the Chiral Effective Field Theory (χEFT)[31,46,61].

TABLE I .
List of Constraints used in Bayesian Analysis.The NICER constraints marked with * are given half weight in the Bayesian analysis as the 4 constraints listed here come from two measurements analyzed by two different groups.
instead of Z SAT .