Imposing multi-physics constraints at different densities on the neutron Star Equation of State

Neutron star matter spans a wide range of densities, from that of nuclei at the surface to exceeding several times normal nuclear matter density in the core. While terrestrial experiments, such as nuclear or heavy-ion collision experiments, provide clues about the behaviour of dense nuclear matter, one must resort to theoretical models of neutron star matter to extrapolate to higher density and finite neutron/proton asymmetry relevant for neutron stars. In this work, we explore the parameter space within the framework of the Relativistic Mean Field model allowed by present uncertainties compatible with state-of-the-art experimental data. We apply a cut-off filter scheme to constrain the parameter space using multi-physics constraints at different density regimes: chiral effective field theory, nuclear and heavy-ion collision data as well as multi-messenger astrophysical observations of neutron stars. Using the results of the study, we investigate possible correlations between nuclear and astrophysical observables.


Motivation
Neutron stars (NS) provide us the unique opportunity to study cold and dense matter under extreme conditions far beyond the reach of terrestrial experiments.Having a mass of up to 2 M ⊙ (solar mass) confined within a radius of roughly 10 km, NSs are among the densest objects in the Universe.As the density increases from the surface towards the interior, one encounters different forms of matter [1,2,3,4,5].The NS crust comprises of aggregates of nuclei that become more and more rich in neutrons, until they dissolve in a homogeneous liquid, with baryonic densities comparable to those reached in heavy-ion collisions at relativistic energies.At such high densities, strangeness containing exotic particles, such as hyperons and kaons, which are observed only briefly in particle accelerators, can become stable components of NS matter due to weak equilibrium [6,7,8].
The behaviour of matter can be represented in terms of its pressure-density relationship, also known as the Equation of State (EoS) [4,5,9].NSs probe a complementary regime of the phase diagram of QCD (Quantum Chromo-Send offprint requests to: a debarati@iucaa.indynamics, the theory of strong interactions), than those of terrestrial experiments [4,10].Laboratory nuclear experiments provide information about the cold (T = 0) EoS close to saturation nuclear density (n 0 ∼ 0.16 fm −3 ) [11,12].Heavy Ion Collision (HIC) experiments, such as at GSI in Darmstadt [13,14], probe hot and dense but more or less isospin symmetric (in neutron/proton ratio) matter at intermediate densities ∼ 2 − 3n 0 [15,16].NSs on the other hand describe nuclear matter as highly isospin asymmetric (ANM), at high densities and low temperatures.While the above terrestrial experiments provide clues to the behaviour of NS matter around or below saturation density, NS models require extrapolation to higher densities and asymmetries, giving rise to uncertainties in the model parameters.
Many theoretical schemes are applied to describe the EoS of dense NS matter [17], such as ab-initio (where one has to solve the many-body problem) [18,19,20,21] and phenomenological (where model parameters are fitted to experimental nuclear observables), including both nonrelativistic and relativistic approaches [22,23,24,25,26].Chiral Effective Field Theory (χEFT) is a microscopic ab-initio technique which provides an effective and reliable description of pure neutron matter (PNM) at densities ∼ 0.5 − 1.4n 0 [27].In phenomenological NS models such as the Relativistic Mean Field (RMF) model [24,25], the parameters are fitted to nuclear empirical observables measured at saturation density (e.g.binding energy at saturation, compressibility, symmetry energy and its slope, effective nucleon mass), which are determined from the measurements of neutron skin thickness of 208 Pb or 48 Ca [28,29,30], from electric dipole polarizability α D , from isoscalar giant monopole and dipole resonances [31,32], nuclear masses [33,34,35,36] and from isobaric analog states [37].
On the other hand, constraints for NS models at high densities come from astrophysical observations [9,11,38].NSs are observable throughout the electromagnetic spectrum, from gamma-rays to radio frequencies.With rapid advance in technology, the current generation of spacebased and ground-based observatories (e.g.Fermi, INTE-GRAL, Parkes pulsar timing array, LOFAR, NICER, XMM-Newton, Chandra) are providing us with a vast pool of Multi-Wavelength (MW) astrophysical data, and a wealth of information about these objects.Many more such MW astronomy projects are under contruction or in conception (Thirty Meter Telescope, E-ELT, SKA, THESEUS, ATHENA, CTA).
There are many different astrophysical observables concerning NSs that may be inferred using MW astronomy.Spin frequency of NSs can be measured very accurately using pulsar timing.General Relativistic (GR) effects in pulsars in binaries can lead to a very accurate determination of NS masses.It is now well known that NSs can have masses as high as 2 M ⊙ [39,40,41].Recent observations of the heaviest known pulsar PSR J0740+6620 indicate that the maximum mass of the neutron stars should exceed 2.14 +0.10  −0.09 M ⊙ [42].This value has recently been revised to 2.08 +0.07 −0.07 M ⊙ [43].NS radii may be determined using hot spots on NS surface, offset from the rotating pole modulated by stellar rotation [44,45,46,47,48].One of the primary goals of the next generation of hard X ray timing instruments as well as the recently launched NICER (Neutron star Interior Composition ExploreR) mission [49] is to measure NS radii with different masses to high accuracy.Recent NICER data provide the mass and radius measurements for pulsars PSR J0030+0451 [50,51] and PSR J0740+6620 [52,53].Using equations of hydrostatic equilibrium, the Tolman-Oppenheimer-Volkoff equation (TOV), one may map the nuclear EoS to observable quantities such as its mass and radius.Then the measurement of the mass and radius effectively allows us to constrain the NS EoS.
With the advent of the LIGO-Virgo-Kagra (LVK) global network of gravitational wave detectors [54,55,56,57], a new window to the Universe has opened up, and the pool of data is expected to be further enriched with the upcoming generation of detectors such as the Cosmic Explorer, Einstein telescope or LIGO-India.The most exciting and revolutionary astrophysical discovery in recent times has been the direct detection of Gravitational Waves (GW) from NSs [58,59,60].The detection of the NS merger event GW170817 is particularly unique because this event seen in gravitational waves and its subsequent followup in electromagnetic counterparts by telescopes such as Fermi, Integral and Chandra opened up a new era of Multi-Messenger (MM) astronomy.The role of dense matter EoS is crucial in the tidal response of NSs during late stage of binary inspiral and post-merger remnant oscillations [61].Recent analyses of the GW170817 event [62,63,64] and updated values in the recent catalog paper [65] put a tight constraint on the effective tidal deformability Λ. Abbott et al. (2019) [64] provide confidence intervals of the effective tidal deformability depending on the spin prior and type of interval.For the low-spin case, the highest posterior density interval yields a value of tidal deformability to be 300 +420 −230 .As the tidal deformability depends on the stellar mass and radius, and therefore this result also leads to a constraint on the mass-radius relation [62,66,67].
Since the pioneering works by Steiner et al. [68], there have been many attempts to impose constraints on the EoS by using multi-messenger observations of neutron stars using a statistical Bayesian scheme [69,70,71,72,73,74,75].The idea of this scheme is to match the low density EOS constrained by theoretical and experimental nuclear physics with parametrized high density EOSs satisfying gravitational wave and electromagnetic data [76,77,78,79,80].Usually EoSs based on different parametrization schemes such as piecewise polytropes [67,81,82,83], spectral representation [84,85] or speed-of-sound parametrization [77,86,87] have been used for such studies, convenient for numerical relativity simulations and parameter estimation from gravitational waveforms.Non-parametric inference schemes have also been developed recently [88,89] and combined with chiral EFT results [90].However in such studies the degrees of freedom at intermediate or high densities are overlooked and the EoS curves are modelled by chosen functional forms or considering general physical conditions such as causality.The role of the underlying nuclear physics, such as the key nuclear saturation parameters that control the behaviour of NS global observables, are not apparent from such studies.Few recent studies have used Bayesian scheme within the RMF model [91] or hybrid (nuclear + piecewise polytope) parametrizations [71] to obtain posterior distributions of empirical parameters but no clear correlation between the nuclear empirical and astrophysical observables was established in such studies.A few recent investigations employed Bayesian techniques to explore correlations among empirical nuclear parameters and a few NS observables [92,93,94,95].Some of these works employed nuclear meta-modelling technique [96], but such a parametrization is valid around saturation density, and their validity at high densities is questionable.
The aim of this investigation is to apply state-of-theart data from multiple experimental and astrophysical channels at different densities to constrain the nuclear EoS and to systematically explore possible correlations between astrophysical and nuclear empirical parameters.We perform this investigation within the framework of a realistic phenomenological model, the Relativistic Mean Field (RMF) model, to describe NS matter.The uncertainty in the behaviour of nuclear empirical quantities is reflected in the uncertainty in the determination of the RMF models parameters.Motivated by the Bayesian approach, we vary the parameters of the realistic nuclear model within their allowed uncertainties, compatible with state-of-the-art nuclear experimental data.We then apply a simple cut-off filter scheme to constrain the parameter space using a combination of current best-known physical constraints at different density regimes: theoretical (chiral effective field theory) at low densities, experimental (nuclear and heavy-ion collision) at intermediate densities and multi-messenger (multi-wavelength electromagnetic as well as GW) astrophysical data at high densities to restrict the parameter space of the nuclear model.A similar scheme was applied by Hornick et al. [97], where limits were imposed on only isovector saturation nuclear parameters within the RMF parameter space using chiral effective field theory at low densities and multi-messenger (multi-wavelength electromagnetic as well as GW) astrophysical data at high densities.We extend this study by considering variation of all the parameters within their uncertainties, apply updated multi-messenger constraints and also impose additional constraints from heavy-ion data at intermediate densities.
The structure of the article is as follows: in Sec. 2, we describe the formalism, particularly the microscopic description (details of the nuclear model) and the macroscopic global modelling.The details of the filter scheme are discussed in Sec. 3. We first test the scheme by reproducing the results of [97] in Sec 4. We then extend the scheme to a full analysis including a wider range of parameters and additional constraints in Sec. 5. We demonstrate the results of the investigation, including comparison with previous results and the study of correlations.Finally, in Sec.6 we discuss the main implications of this work.

Microscopic description
As mentioned in Sec. 1, we perform this study within the framework of the RMF model to calculate the NS EoS.We start with the interaction Lagrangian density [97,98,99]: where Ψ N is the Dirac spinor for nucleons N , m is the vacuum nucleon mass, while γ µ and ⃗ τ are the Dirac and Pauli matrices respectively.The interaction among the nucleons is mediated by the exchange of the scalar (σ), vector (ω) and the isovector (ρ) mesons.The isoscalar nucleon-nucleon couplings g σ and g ω are determined by fixing them to nuclear saturation properties.The σ meson self-interaction terms b and c ensure the correct description of nuclear matter at saturation.The effective nucleon mass is then defined as m * = m − g σ σ.The isovector and mixed ω − ρ couplings g ρ and Λ ω can be related to empirical quantities such as symmetry energy (E sym ) and its slope (L sym ) [98].The quartic ω self-coupling ζ is set to zero in this study, because the corresponding term is know to soften the EoS which is in tension with the 2M ⊙ mass constraint from pulsar data.The choice of model parameters is discussed in Sec.2.3.
From the Lagrangian density Eq. ( 1), one can solve the equations of motion of the constituent particles as well as those of the mesons [97].In the mean-field approach, the meson fields are replaced by their mean-field expectation values.One can then calculate the EoS starting with this RMF model.The energy density is given by [97] where k F N and E F N are the Fermi momentum and energy of the corresponding nucleon N respectively.The pressure P can be derived from the energy density using the Euler relation where n is the number density and the nucleon chemical potentials (µ) are given by

Macroscopic description
Given an EoS, the equilibrium configurations of non-rotating relativistic NSs can be obtained by solving the Tolman-Oppenheimer-Volkoff (TOV) equations of hydrostatic equilibrium [100,101], Integrating the TOV equations from the centre of the star to the surface, one can obtain global NS observables, such as mass, M , radius, R, and compactness, C = M/R that can be deduced from astrophysical data.The boundary conditions that must be satisfied are a vanishing mass, m| r=0 = 0, at the centre of the star, and a vanishing pressure, p| r=R = 0), at the surface.The tidal deformability, Λ, can be obtained by solving a set of differential equations coupled with the TOV equations [61,102,103,104].It was shown in these works that the tidal deformability parameter Λ depends on the mass and radius of the NS (and hence its EoS) [105] where k 2 is the l = 2 dimensionless tidal Love number.

Model parameters
To obtain the EoS, one must know the coupling constants, which are referred to as the model parameters.The nucleon isoscalar coupling constants (g σ , g ω , b, c) defined in Sec.2.1 are calculated by fixing the empirical properties obtained from nuclear experiments, such as nuclear saturation density (n 0 ), binding energy per nucleon (E/A), incompressibility (K sat ) and the effective nucleon mass (m * ) at saturation.On the other hand, isovector coupling constants (g ρ N , Λ ω ) are obtained by fixing the symmetry energy (E sym ) and the slope (L sym ) of symmetry energy at saturation (see e.g., [98,97]).
In order to test our calculations, we first calibrated the model parameters by calculating the EoS, corresponding mass-radius relations and tidal deformability for the wellknown parametrizations GM1 and GM3 [106].We then reproduced the EoS for neutron matter and the mass-radius relations given in Hornick et al. [97] (see e.g.Figs. 1, 2, 5, 6 of that paper).The saturation properties of these parameter sets are recalled in Table 1.

Cut-off filter scheme
In a typical Bayesian analysis, an initial range of parameters is varied to produce a prior probability distribution.The posterior distribution is then obtained by multiplying the prior and likelihood functions.The likelihood functions are appropriately chosen physical conditions.
In this study, we apply a "cut-off filter" scheme, where we impose strict limits on the nuclear parameter space using multi-density constraints.The nuclear empirical parameters are varied randomly within their uncertainty range in Table 1 to generate the EoSs.Among the EoSs, we then allow only the ones which satisfy certain physical requirements or "filters", described in detail in Sec.3.1.
Although Miller et al. (2019) [107] point out the statistical uncertainties in constraining EoS by putting strict limits from multi-messenger observations of neutron stars like we have used here, using the cut-off scheme gives a correct estimate of the nuclear parameter ranges consistent with the observations.Many recent works [108,109,110,111,112] also used similar cut-off schemes for constraints for ultra-dense matter.Further, one must note that chiral effective theory and heavy-ion collision experiments give a band of nuclear observables as constraints and a statistical weighting of the model-dependent heavyion data might lead to a false confidence in the results.So, using a cut-off filter scheme for the constraints seems appropriate in these cases, in order to extract the underlying physics.We have also generated a very large number of prior EoSs to constrain the parameter space.One may verify whether Bayesian analysis with statistical weighting affects the posterior probability distribution, or significantly alters the physical correlations between nuclear empirical parameters and astrophysical observables which is the main aim of this study.In Sec.5.2 and Sec.5.5 we check this explicitly by redoing the analysis with a proper calculation of the likelihood of the data.
As mentioned in the introduction, such limits were also imposed on the RMF parameter space in the recent study by Hornick et al. [97].A similar scheme using the metamodelling approach [113] was also employed in the analysis of the effect of uncertainties in nuclear empirical parameters and the study of their correlations with NS observables [114] and experimental observables of finite nuclei by one of the authors (D.C.) of this work [115].This work [115] found no correlation between the symmetry energy E sym and its first order derivative L sym , which is at variance with several studies in the literature which find a good correlation between E sym and L sym [116,117,118].

Filter functions
The following physical constraints from multi-disciplinary physics are applied in this work at different densities, as described below: -At low densities: χEFT Important constraints on the EoS of NS matter at low baryon densities n b in the range of ∼ 0.5 − 1.4n 0 come from microscopic calculations of pure neutron matter (PNM) using chiral effective field theory [27,81,119].Chiral EFT applies low momentum expansion of nuclear forces related to the symmetries of QCD.It accounts for many-body interactions among nucleons using order-by-order expansions in terms of contact interactions and long-range pion exchange interactions.
Thus the EoS of PNM can be compared with the latest χEFT calculation results [77].
-At high densities: astrophysical observations The constraints at high density come from multi-messenger astrophysical observations, such as high mass NS observations and GW constraints on tidal deformability, as follows : Table 1.Empirical saturation parameters for the RMF models GM1, GM3 [106] and the parameter sets from HTZCS [97]. 1.From the recent observations of the heaviest known pulsar PSR J0740+6620, the maximum mass of the neutron stars should exceed 2.08 +0.07 −0.07 M ⊙ [43].This sets an upper bound on the maximum NS masses corresponding to the EoSs considered in this study.2. Using the low-spin highest posterior density interval for tidal deformability from the recent analyses of the GW170817 event [64], we apply a constraint on the upper bound of the effective tidal deformability Λ < 720 [120].We do not consider the lower limit on tidal deformability in this study.
As explained in Sec. 1, the tidal deformability depends on the mass and radius (see Eq. 6), and therefore this result also leads to a constraint on the mass-radius relation, imposing that NSs should have R 1.4M⊙ < 13.5 km [66,67].
-At intermediate densities: heavy-ion collisions At intermediate densities in the range 1 − 3n 0 , we have additional information from heavy-ion collision experiments.There are several heavy-ion collision experiments, which impose important constraints on the nuclear parameters in the density region of n b /n 0 ∼ 1 − 3 which we discuss in detail below: 1. KaoS experiment The measurements on K + meson production in nuclear collisions at subthreshold energies were performed with the Kaon Spectrometer (KaoS) at GSI, Darmstadt [121].The Kaon multiplicity in Au+Au & C+C collisions is an indicator of the compressibility of dense matter at n b ∼ 2 − 3n 0 .The level of compression is in turn controlled by the stiffness of nuclear matter through the nucleon potential U N .The more attractive U N is, the higher is the produced K + meson abundance, which therefore can serve as a probe for the stiffness of nuclear matter.The analysis of the experimental results using IQMD transport models indicated that the EoS in this density range must be soft, as described by a simple Skyrme ansatz with an incompressibility K of ∼ 200 MeV or less [121,122].[14].From the comparison of the elliptic flow ratio of neutrons with respect to charged particles with UrQMD predictions, a value γ = 0.72 ± 0.19 is obtained for the power-law coefficient describing the density dependence of the potential part in the parametrization of the symmetry energy, as defined in the equation below [14] : We would like to emphasize here that although the constraints from heavy-ion collisions are important and interesting, they are model-dependent.The heavy-ion data are confronted by well established transport codes such as IQMD [123,124] or Ultra relativistic Quantum Molecular Dynamics (UrQMD) using various phenomenological EoSs to assess which EoS is most compatible with the data.There are two important points to be aware of in this regard.First of all, the data that are claimed to constrain the EoS of nuclear matter require confirmation by independent experimental ef-forts.Secondly, it must be proved that the conclusions using a particular transport code (IQMD or UrQMD) are not limited to that particular code.Some efforts have already been made in this direction [125].In view of the points raised earlier we do not consider the constraints from collective flow of nucleons [126].

Correlations
Using the filtered EoSs obtained from the analysis, one can then perform a search for physical correlations of the different empirical parameters among themselves, as well as with the NS astrophysical observables such as radii and tidal deformability of canonical 1.4 M ⊙ and massive 2 M ⊙ NSs.For this study, we use Pearson's linear correlation coefficient defined as : between two random variables X and Y where Cov(X, Y ) denotes the covariance and S(X), S(Y ) denote the standard deviation for the variables X and Y respectively.The correlations studied here are calculated from the confidence bounds used in this study.

Statistical weighting
In Sec. 3, we discussed the need to take into account the statistical uncertainties in the measurement of the nuclear and astrophysical observables while constraining the EoS.
In this section, we summarise how to include the effects of statistical weighting on the correlations.
In this work, we use a χ 2 distribution, as is done for hypotheses testing or goodness of fit [127,115,128].We obtain our priors from the parameter space {P} by uniformly varying the nuclear parameters from Table 2.The posterior is obtained by multiplying the prior and the likelihood functions.For Gaussian variables, the likelihood is proportional to which we assign as weight W to each parameter set {P}.In (9), the χ 2 is a weighted sum of squared deviations defined as below : As discussed in Sec.3.1, for low density χEFT filters, we get bounds for the PNM EoS in the density range of ∼ 0.5 − 1.4n 0 [27].We use these bounds as the uncertainty range to calculate the standard deviation(σ i ).
For astrophysical observations, mass measurement of PSR J0740+6620 [43] (2.08 +0.07 −0.07 M ⊙ ) and tidal deformability measurement from GW170817 [64](300 +420 −230 ) give the observed mean(O i ) and σ i values for the observables.Using the equations ( 10) and ( 9) we then calculate the weights associated with each parameter set {P} for these filters.For χEF T and multi-messenger observations of the NSs, the uncertainties can be determined to a good accuracy, but due to the model dependence of the heavy-ion collision data, statistical weighting may lead to false conclusions.We include weights for HIC filters also in Sec.5.5 with a word of caution.As discussed in Sec.5.3, for KaOS constraint, we calculate the nucleon potential and associated errors(σ i ) from the non-relativistic Skyrme model for K sat =200 MeV and n 0 = 0.15f m −3 and use it as a upper limit [129] to calculate the weight.From FOPI [13] and ASY-EOS [14] experiments, we obtain uncertainty bands for binding energy and Symmetry energy respectively which are used to calculate the associated mean and standard deviations in the corresponding density range.The final weight for a parameter set {P} is the product of the weights associated with each constraint: The definition of the Pearson's linear correlation coefficient in Eqn. ( 8) also gets modified when we include weights.The definition of covariance between two variables X and Y with a weight W is where M (X, W ) is the mean of the variable X defined as M (X, W ) = i XiWi i Wi .The weighted correlation coefficient is defined as 4 Preliminary results: Testing the scheme In order to test whether the applied scheme can be used to impose constraints on the EoS parameter space, we first reproduce the results from the parameter study in Hornick et al. [97], hereafter referred to as HTZCS.In this work, the isoscalar couplings (n 0 , E sat , K sat ), which have comparatively lower uncertainties, were kept fixed while the isovector couplings (E sym , L sym ) along with effective mass m * /m were varied individually within the chosen range motivated by theoretical and experimental uncertainties (see Table 1), i.e. to obtain the EoSs.

Input EoSs
We initially obtain input EoSs by randomly picking up values within the parameter space from HTZCS given in Table 1.For these sets, we calculate the PNM and ANM EoS as described previously.In Fig. 1, we plot the binding energies vs density of low density PNM EoS superposed with the χEFT band [27].

Output EoSs
After generating the random EoSs, we use the χEFT and astrophysical filters described in Sec.3.1 to obtain filtered sets for the nuclear parameters and NS observables.Around 50% of the EoSs pass through both the filters (around 60% for only χEFT and 80% for only NS observables.)

Low density: χEFT
We first apply the cut-off filter applying the constraints from χEFT on the EoS of pure neutron matter [27].In order to do so, we evaluate the binding energies in the density range of n b /n 0 ∼ 0.5 − 1.4 corresponding to the χEFT data and allow only those parameter sets that lie within the band allowed by χEFT calculations (Fig. 2).

High density: NS astrophysical observations
Comparing with the recent observations of mass, radius, moment of inertia and tidal deformability using multimessenger astrophysical and GW data, we can rule out further combinations of parameter sets and allow only those combinations which simultaneously satisfy all constraints on NS observables.
In Fig. 3, we plot the mass-radius relations corresponding to the filtered ANM EoSs, along with the uncertainty in the measurement of maximum mass of the most massive pulsar yet observed PSR J0740+6620 [43].It is evident from this figure that the NS radii span a wide range from 11-14 km.We plot the dimensionless tidal deformability as a function of NS mass in the posterior ANM EoSs.Posterior Fig. 3. Mass-radius relations for posterior ANM EoSs after passing through χEFT and NS observations filters(maximum mass of PSR J0740+6620 [43] and tidal deformability Λ from GW170817 [64]).The light green band indicates the uncertainty uncertainties in the measurement of maximum mass of PSR J0740+6620 [43]

Correlations among nuclear empirical parameters and with NS observables
To verify some of the observations in HTZCS, we probe whether there exist correlations among nuclear empirical Posterior Fig. 4. Dimensionless tidal deformability for posterior ANM EoSs after passing through χEFT and NS observations filters (maximum mass of PSR J0740+6620 [43] and tidal deformability Λ from GW170817 [64]).
parameters after applying the χEFT band as a filter in our analysis.We plot the correlations among the parameters L sym and m * /m in Fig. 5. From the plot, we observe that there are no points in the region below L sym = 55 MeV and m * /m = 0.65.This was also concluded from Fig. 3 and 4 in HTZCS: that it is more difficult to obtain physical solutions for the neutron matter EoS simultaneously for small L sym and m * /m compatible with the χEF T results.This can be explained by the Hugenholtz-van Hove theorem, which states that the binding energy per particle has to be equal to the Fermi energy at saturation.The decrease in m * /m (increase in the scalar potential) would lead to larger values for the vector field and a stiffer EoS.The softening of the EoS for small L sym competes with the stiffening of the EoS as m * /m is decreased, leading to either a solution outside the band of χEFT or the appearance of unstable solutions at sub-saturation density.
In Fig. 6, we display the correlation matrix of the following quantities: nuclear empirical parameters (E sym , L sym , m * /m) and the NS observables (R 1.4M⊙ , Λ 1.4M⊙ , R 2M⊙ , Λ 2M⊙ ).Some of the crucial observations we make from the correlation matrix are listed below: -As given in Table .1, the values of the isoscalar parameters n 0 , E sat and K sat are not varied as in HZTCS (this will be done in the full analysis in Sec. 5) and therefore the correlations among them do not appear in the correlation matrix.-We see a very high correlation between the tidal deformability Λ and radius R which is expected from the formula of tidal deformability parameter (Λ ∝ R 5 ) (see Eq. 6).-The effective nucleon mass is highly correlated with the radius and tidal deformability.As explained in HTZCS, due to the Hugenholtz-van-Hove theorem, the smaller the effective nucleon mass is, the stiffer the EoS is.Thus the effective mass controls both the variation in the maximum mass and the radius, and consequently also the tidal deformability.-Among the nuclear empirical parameters, the symmetry energy E sym and its first order derivative L sym appear to be have a moderate correlation (0.59).-There is a weak correlation (0.36) between the slope of symmetry energy parameter L sym and radius at 1.4M ⊙ .A correlation between L sym and R 1.4M⊙ has been reported in several articles in the literature [130,131,132,133], although in HTZCS R 1.4M⊙ was found to be nearly independent of L sym .It was also concluded in HTZCS, that the variation of m * /m is more important than L sym in the determination of R 1.4M⊙ .This is clearly in agreement with our results.-We also note that the correlation of L sym with R 2M⊙ is very poor (0.04).As explained in HTZCS, at high densities the isoscalar vector field ω grows linearly with density, while the equation of motion for the isovector vector field ρ has a trivial solution for a ρ field growing inversely proportional to the density.For a 2M ⊙ neutron star, these two effects lead to a cancellation of the density dependence of the nuclear symmetry energy, i.e. the slope L sym around saturation density.So, we expect very little correlation of L sym with observables for a 2M ⊙ NS.

Results: Full Analysis
In HTZCS, the isoscalar saturation parameters were kept fixed and only the isovector parameters varied within their uncertainties.It is however unclear whether simultaneous variation of all the parameters would affect the above conclusions.This is crucial, as a study of sensitivity to each individual empirical parameter does not take into account any underlying correlations among the parameters themselves, or with the observables.
In this section, we will vary both isoscalar and isovector nuclear empirical parameters simultaneously within their uncertainty ranges, motivated by state-of-the-art nuclear experimental data.We will then impose additional constraint filters from heavy-ion collision experiments to further restrict the parameter space.Finally we will investigate possible correlations among the empirical nuclear parameters and the astrophysical observables.

Variation of all nuclear parameters
Having established our scheme and verified the results of HZTCS, now we extend our analysis to all the nuclear parameters and also for larger number of parameter sets.The ranges of variation for the nuclear parameters are taken to be compatible with estimations coming from the analysis of a variety of nuclear data from terrestrial experiments, astrophysical observations and theoretical calculations.The compressibility modulus K sat is taken to be in the range 200-300 MeV, as obtained from current experimental data on isoscalar giant monopole and dipole resonances (compression modes) in nuclei [31,134].For the energy per particle sat saturation E sat , we adapt values -16±0.2MeV from binding energy of atomic nuclei [33,34].The symmetry energy E sym is varied between 28-34 MeV, the slope of symmetry energy L sym in the range 40-70 MeV, while the nucleon effective mass m*/m is considered to vary between 0.55-0.75 as in HTZCS.The ranges of all the nuclear parameters considered in the full analysis are listed in Table 2 Fig. 7. Posterior correlation matrix for variation of all nuclear empirical parameters and NS observables, after application of the χEFT and NS observations filter (maximum mass of PSR J0740+6620 [43] and tidal deformability Λ from GW170817 [64]) Following the same cut-off filter scheme described in Sec. 3, we now obtain posteriors of 1000 points for the nuclear parameters and NS observables, for the filters from χEFT and NS observations.In Fig. 7, we display the correlation matrix of the following quantities: nuclear empirical parameters (n 0 , E sat , K sat , E sym , L sym ), the effective mass m * /m and the NS observables (R 1.4M⊙ , Λ 1.4M⊙ , R 2M⊙ , Λ 2M⊙ ).It was found that the correlation values between the variables do not significantly change with increasing the number of posteriors points.Some of the main observations from the correlation matrix are listed below: n 0 and m * /m now show a moderately good correlation (0.42).-L sym and E sym display a strong correlation (0.72), higher than our observation in the previous section n 0 has a weak correlation with the NS observables.
The correlation is noticeable (0.51) for 1.4M ⊙ NS but is negligible for 2 M ⊙ .-The correlation of m * /m with 1.4 M ⊙ NS observables is moderate (≈ 0.5), but is higher for 2 M ⊙ NS. -The correlation between L sym and radius of 1.4M ⊙ NS is also lower (0.32) than our previous findings.-All the NS observables (radius and dimensionless tidal deformability for 1.4 M ⊙ and 2 M ⊙ NS), as expected, show a strong correlation with each other (according to Eq. 6).

Correlations with statistical weighting
As discussed in Sec. 3, one must check whether the inclusion of statistical weighting in the applied scheme with the hard cut-off limit considered significantly affects the conclusions of this study.In order to do that, we recalculate the posterior correlation matrix for variation of all nuclear empirical parameters and NS observables, after application of the χEFT and NS observations filters as in Fig. 7, by including a proper calculation of the likelihood of the data for the given model.The results are shown in Fig. 8.In summary, other than the expected correlations among NS observables and between symmetry energy and its slope, we found strong correlations of the effective nucleon mass with the NS observables and between saturation density and the symmetry energy.

Additional constraints from heavy-ion collision experiments
In the previous section, we imposed constraints on the nuclear parameters at very low-density from χEFT data and at very high density from multi-messenger observations of neutron stars.Now, in this section we will apply the results from several other heavy-ion collision experiments, which impose important constraints on the nuclear parameters in the density region of n b /n 0 ∼ 1 − 3.However, it is well known that the above constraints are modeldependent, and therefore the implications of the results must be taken with caution.

KaoS Experiment
As discussed in Sec.3.1, results from the KaoS collaboration directly test the EoS for densities ∼ 2 − 3n 0 .For Skyrme type models this restriction implies compression moduli of K sat ≤ 200 MeV [135,136].The constraint given by the KaoS data is accomplished by only allowing nucleon potentials which are more attractive than the Skyrme parametrization within the considered density regime.Therefore for RMF models, one should derive the nucleon potential and compare it with the Skyrme parametrization with K sat ≤ 200 MeV.A more attractive potential at high densities will allow matter to be compressed more for the same incident energy, enhancing multiple scattering processes in subthreshold kaon production [129].
In this analysis the non-relativistic Skyrme model is used to produce the EoS for K sat =200 MeV and n 0 = 0.15f m −3 .To impose the KaoS constraint, we calculate the nucleon potentials for the priors (saturation parameters randomly picked from the parameter space in Table 1 as described in Sec.4.1) and allow only those for which the nucleon potential is less than that for the Skyrme EoS.The results are shown in Fig. 9. Around 40% of our prior sets pass through this filter. 2.0

FOPI Constraint
In Sec.3.1, we discussed the fact that using the FOPI data, one can obtain a constraint for the binding energy of SNM in the density region of 1.4-2 n 0 .To impose this constraint, we calculate the binding energy for SNM for the input parameters randomly picked from the parameter space as in Table 2 and allow only those for which the binding energy lies inside the band allowed by the FOPI data in this density range.This is depicted in Fig. 10.Nearly 20% of our prior sets pass through the FOPI filter.

ASY-EOS constraint
It was discussed in Sec.3.1 that information about the symmetry energy for ANM at supra-saturation density can be obtained from the recent ASY-EOS data.We computed the symmetry energy for ANM EoS using Eq.14 below [100]: where k F is the Fermi momentum of the nucleons.Since the value of the symmetry energy at saturation in the interval between 30 MeV and 32 MeV seem to be favoured by a majority of terrestrial experiments and astrophysical observations [137,138], we used the constraints on E sym for the choice of parametrization at E sym (n 0 ) = 31 MeV as in [14].To impose the constraint filter in the Bayesian scheme, we calculate the symmetry energy of ANM EoS using Eq.14 for the input parameters randomly picked from the parameter space as in Table 2 and allow only those for which the symmetry energy lies inside the band allowed by the ASY-EOS data in the range of ∼ 1.1−2.0n0 (Fig. 11).We find that around 40% of our prior sets pass through this filter.In Fig. 7, we displayed the correlation matrix for all the nuclear parameters and the NS observables for the posteriors obtained after passing through the filters from χEFT and NS observations.Now we will add three additional HIC constraints from KaoS, FOPI and ASY-EoS experiments discussed in Sec.3.1.On applying the constraint from KaoS experiment in addition to χEFT and NS observables, the correlation of n 0 with the NS observables increases while that with m * /m decreases.The correlation between K sat and m * /m also increases (not displayed here).Fig. 12. Posterior correlation matrix for all the nuclear empirical parameters and NS observables, after application of the χEFT,HIC experiments (KaOS, FOPI, ASY-EOS) and NS observations filters (maximum mass of PSR J0740+6620 [43] and tidal deformability Λ from GW170817 [64]) After adding all the HIC filters, we reproduce the correlation plot given in Fig. 12.There are a few differences between the plots before (Fig. 7) and after applying the constraints from KaoS, FOPI and ASY-EOS experiments (Fig. 12

Correlations with all (χEFT+astro+HIC) constraints including weighting
In Sec.5.2, we studied the correlations applying filters from χEFT and NS astrophysical observables, including statistical weighting in the applied scheme.In the previous section, we included the constraints from HIC experiments, and studied their effects on the correlations (without statistical weighting).As discussed in Sec. 3, inclusion of weights for HIC may introduce a false confidence in the results, as HIC experiments are known to be model dependent and have large uncertainties.In this section, we proceed to calculate the correlations including all filters (χEFT+astro+HIC) and with statistical weighting, with a word of caution.In Fig. 13, we display the resulting correlations.Upon comparison with Fig. 8, we can draw the following conclusions: 1.The correlation of n 0 -E sym (0.54) still persists, while the n 0 -L sym correlation weakens (0.15) 2. The correlation of L sym and E sym changes significantly on application of HIC filters, from 0.71 to 0.25 3.As in Fig. 8, there is a weak correlation of n 0 with NS observables 4. The correlation of K sat with m * /m has increased from 0.1 to 0.51 5.The (L sym , m * /m) correlation is not significantly affected 6 Discussions

Summary of the results
In this work we applied the current best-available information from multiple channels, namely nuclear theory and experiments, heavy-ion collision experiments as well as multi-messenger (multi-wavelength electromagnetic as well as GW) astrophysical observations to constrain the nuclear EoS at different densities.The phenomenological RMF model applied here satisfies state-of-the-art nuclear matter saturation properties, low-density constraints from neutron matter from χEFT ab-initio approach, heavy-ion collision experimental results at intermediate densities and multi-messenger astrophysical observational data at high densities.We applied the constrained EoS parameter space to study correlations among nuclear empirical parameters and NS observables by applying a cut-off filter scheme.
In order to test the approach, we reproduced the results of the recent investigation by Hornick et al. [97] (HTZCS).We simultaneously varied the isovector nuclear empirical parameters within the allowed parameter space of uncertainties, to generate random EoSs.Then we imposed filters from physical constraints such as χEFT data and recent NS astrophysical observations of mass, radius and tidal deformability to restrict the choice of parameters.We used the filtered parameter sets to study correlations among the parameters themselves and also with the NS observables.We found strong correlations between E sym and L sym , in agreement with previous literature.We confirmed the findings of HTZCS that the correlation of R 1.4M⊙ with m * /m is stronger than that of L sym .Our results also agreed with HTZCS that the effective nucleon mass is highly correlated with the radius and tidal deformability.The expected strong correlation between tidal deformability and radius was also seen for both 1.4 and 2 M ⊙ stars.
The scheme was then extended to a full analysis, considering variation of both isoscalar and isovector saturation parameters and also for a larger sample of data in order to improve the statistics and to test the robustness of the results.We observed the usual strong correlation of L sym with E sym , and a weaker correlation of L sym with R 1.4M⊙ .We also recovered the expected correlations among the NS radii and tidal deformability.In addition, we also found weak correlations of n 0 with m * /m, and with the NS observables.
NICER data provide constraints on the mass and radius of a pulsar.The Bayesian parameter estimation of the mass and equatorial radius of the millisecond pulsar PSR J0030+0451 [50] yield M = 1.44 +0.15 −0.14 M ⊙ and R e = 13.02+1. 24 −1.06 km, consistent with the independent analysis by [51].More recently for PSR J0740+6620, the equatorial circumferential radius obtained was 13.7 +2.6 −1.5 km [52] (68%).In Fig 14, we superimpose the filtered M-R curves with the mass and radius measurements from the recent NICER data, for pulsar PSR J0030+0451 [50,51] and PSR J0740+6620 [52].We verified that our results are consistent with the latest analyses from NICER.All our filtered points pass through the 1-σ limit of the M-R region of the NICER observations for pulsar PSR J0030+0451 [50] and PSR J0740+6620 [52].However, another independent radius measurement for J0740+6620 was performed by Riley et al.(2021) [53], which differs from the Miller et al.( 2021) results [52] by ≈ 0.8 km at the 1σ credibility contour.Some of our M-R posterior curves do not agree with the analysis by Riley et al. [53] within 1-σ but all posteriors agree within 2-σ (not shown here).
In Fig 15, we superimpose the M-R curves with the posterior mass-radius relations of the two neutron stars from GW170817 observations [62,139] obtained by using a low-spin prior and parametrized-EoS with constraint of maximum mass of 1.97M ⊙ for comparison.We find that that all our posteriors points are within the 90% confidence interval.This study does not use posteriors from previous LIGO-Virgo Collaboration publications, and any correlations that may be implicit in such posteriors have not been considered.The correlations studied here are calculated from the confidence bounds used in this study.We then investigated the effect of imposing additional constraints from heavy-ion data and investigated their impact on the EoS at intermediate densities (∼ 1−3n 0 ).Applying only the KaoS constraint increases the correlation of n 0 with the NS observables while decreasing the same for m * /m.It also increases the correlation of incompressibility K sat with m * /m.We further studied the combined effect of imposing all the HIC filters (KaoS, FOPI and ASY-EOS experiments) in addition to the χEFT and NS observational constraints.The two considerable changes in correlation with application of HIC filters are decreasing of E sym and L sym correlation and increase in the correlation of m * /m and K sat .The first one can be attributed to either or both of FOPI and ASY-EOS.
Further, we verified whether the inclusion of statistical weighting for the filters affect the conclusions significantly.Other than the expected correlations among NS observables and between symmetry energy and its slope, we found strong correlations of the effective nucleon mass with the NS observables and between saturation density and the symmetry energy.The latter was found to weaken on inclusion of heavy-ion filters with corresponding statistical weighting, but this data is known to be modeldependent.In summary, the most important nuclear parameters to consider for astrophysical data are the effective nucleon mass m * /m and the nuclear saturation density n 0 .

Comparison with previous studies
NS radii are known to be sensitive to the symmetry energy at around two times the saturation density [140,141].It was first discussed by Fattoyev et al. [142] that the tidal deformability of low mass neutron stars is sensitive to L sym , while the tidal deformability of massive neutron stars is sensitive to the high density behavior of symmetry energy.A possible relation between L sym − R 1.4M⊙ − Λ 1.4M⊙ was searched by Zhang et al. [143,144].Hornick et al. [97] found R 1.4M⊙ is almost independent of L sym on applying constraints from χEFT.Using several non-relativistic and relativistic RMF and Skyrme Hartree-Fock (SHF) models, Malik et al. [145] found that the tidal deformability is weakly or only moderately correlated with the individual nuclear matter parameters of the EoS.They also found stronger correlations of Λ for specific choices of the linear combinations of the isoscalar and isovector EoS parameters.
The strong correlation between L sym and E sym as well as L sym and R 1.4M⊙ were reported in many works [96] combining different constraints such as astrophysical data from LIGO/Virgo or NICER, χEFT theory and neutron skin thickness from PREX II results [146].A mild tension was reported between nuclear empirical parametrization with those from astrophysical and χEFT results [147].Other works which have employed a hybrid parametrization (nuclear+piecewise polytrope) have not found such apparent tension [72,148] and found a strong correlation between R 1.4M⊙ and L sym based on PREX II results.
Recently Huth et al. [149] combined data from astrophysical multi-messenger observations of neutron stars and from heavy-ion collisions (FOPI [13] and ASY-EOS [14] experiments as well as EoS constraint for symmetric nuclear matter [126]) with microscopic χEF T calculations and used a Bayesian inference technique to analyse the nuclear EoS and NS properties.While χEF T was used to constrain the EoSs at densities below 1.5n 0 , an extrapolation using the speed-of-sound in NS matter was used to extend the EoS at higher densities, with the additional criterion to support at least 1.9M ⊙ NS mass.The mass measurements of the PSR J0348+404232 and PSR J1614-22303 were used to obtain lower bound and the NS in GW170817 to obtain an upper bound for the NS maximum mass.Other multi-messenger information such as from NICER and XMM Newton missions, as well as GW data and kilonova AT2017gfo9 associated with the GW signal were also used in the Bayesian framework.The final EOS constraints were obtained through the combination of both the HIC information and astrophysical multi-messenger observations.The study concluded that the constraints from HIC experiments is in excellent agreement with NICER observations.In comparison with our work, Huth et al. [149] start with a soft EOS from chiral EFT up to 1.5 n 0 and the EOS becomes stiffer due to the large L sym values of the ASY-EOS predictions, and the final results only slightly depend on the EOS of SNM.In contrast for the RMF model used in our work, RMF models are more isospin symmetric at high density so that they are less sensitive to parameters of PNM.

Implications and outlook
We must recall here the difference between the astrophysical versus the heavy-ion constraints applied here and their impact on the parameters.The constraints from heavy-ion data are model-dependent and therefore we treat them separately from astrophysical constraints and discuss the implications of the results with a word of caution.We also note here that the EoS beyond 2-3 n 0 can be associated with the emergence of new degrees of freedom, such as the appearance of hyperons or quarks which we have not considered in this study.With the upcoming CBM experiment at FAIR [150] the EoS in this density regime can be probed further in the near future.
In conclusion, nuclear empirical parameters control important properties of nuclear matter from which accurate predictions for dense matter EoS can be made based on our current knowledge, improved by up-to-date constraints from experimental and astrophysical data as well as ab-initio approaches.Numerical relativity simulations of astrophysical scenarios such as supernovae or binary NS mergers require cold EoSs in beta-equilibrium for the initial data.In order to probe the effect of the microphysics (e.g.new degrees of freedom) on NS astrophysical observables, one requires a large variability of EoSs covering the entire parameter space.The knowledge of the constrained EoS parameter space resulting from simultaneous application of constraints, coming from different physics at different density regimes, can help to make an informed choice of the parameters in the simulation codes.
In the present work, we have pointed out the most important empirical parameters which are mainly responsible for the uncertainty in the NS observables based on the nucleonic EoS.According to the results of our investigation, the essential parameters for RMF models to vary are the saturation nuclear density n 0 and effective nucleon mass m * /m.This was also concluded in a recent investigation of the effect of dominant RMF parameters on f -mode oscillations in neutron stars [151].This conclusion results from the fact that n 0 showed a moderate correlation with R 1.4M⊙ , while m * /m showed a strong correlation with NS observables like radius and tidal deformability.However on imposing additional constraints from heavy-ion data, the correlation of n 0 with R 1.4M⊙ increased while that of m * /m decreased, although it is noted that HIC constraints are less robust and model-dependent.These results should therefore be of interest to the community performing astrophysical and numerical relativity simulations.
) It represents a new and more stringent constraint on the symmetry energy for ANM in the regime of supra-saturation density.The densities probed are shown to reach beyond 2 n 0 .

Fig. 1 .
Fig. 1.Binding energy E/A of pure neutron matter EoS as a function of normalized baryon density n b /n0 for the prior distribution, superimposed with the chiral effective field theory (χEFT) band [27].

FigFig. 2 .
Fig. 2. Binding energy E/A of posterior pure neutron matter EoSs as a function of normalized density n b /n0 allowed by chiral effective field theory (χEFT) data [27].

Fig. 5 .
Fig. 5. Correlation between Lsym and m * /m from the posterior obtained after the χEFT filter

): 1 .
The correlation of n 0 with other nuclear parameters except E sat has increased.2. The correlation of L sym with E sym has decreased (almost by half to 0.36) and now L sym has a significant correlation (0.52) with m * /m. 3. The correlation among NS observables is not significantly altered due to the HIC filters.4. m * /m now has negligible correlation with the observables for 1.4 M ⊙ stars but correlation with n 0 (0.78) has increased from previous case.

Fig. 13 .
Fig. 13.Same as in Fig. 12 after inclusion of statistical weighting

Fig. 14 .
Fig.14.Mass-radius relations for ANM EoSs.The elliptical red and deep blue regions denoted by lines correspond to 1-σ limits of the mass and radius measurement of the pulsar PSR J0740+6620[52] and PSR J0030+0451[50] respectively.The light green band indicates the uncertainty in the measurement of maximum mass of PSR J0740+6620[43].

Fig. 15 .
Fig.15.Mass-radius relations for posterior ANM EoSs.The light blue and orange regions above and below the mass value of M=1.365 M⊙ correspond to the mass-radius estimates for the two compact stars of the merger in GW170817[62,139].The light green band indicates the uncertainty in the measurement of maximum mass of PSR J0740+6620[43].

Table 2 .
Range of empirical parameters considered for the full analysis in Sec. 5.