Implication of multimessenger observations on the relativistic mean-field equation of state of dense nuclear matter and skin thickness of nuclei

Abstract

The composition and properties of infinite nuclear matter under extreme conditions of temperature and pressure remain incompletely understood. In this work, we constrain the equation of state (EoS) of nuclear matter—constructed within the framework of the Relativistic Mean Field (RMF) model—by combining results from chiral effective field theory and multimessenger observations of neutron stars. Using the saturation properties of nuclear matter, we generate a wide ensemble of EoS, which are subsequently constrained within a Bayesian framework. The resulting posterior distributions provide tight bounds on both the saturation parameters and the coupling constants of the RMF model. Our results indicate that the GW170817 event and the latest NICER observation favor a relatively soft EoS, leading to lower crust–core transition densities and thinner neutron star crusts. The radius of a $1.4\,M_{\odot}$ neutron star is tightly constrained to $12.508_{-0.241}^{+0.257}\,\mathrm{km}$, while the maximum mass reaches $2.174_{-0.123}^{+0.174}\,M_{\odot}$. Furthermore, our analysis reveals that the $\omega$–$\rho$ coupling, which governs the density dependence of the symmetry energy, becomes increasingly significant under successive astrophysical constraints. Finally, the predicted neutron skin thickness of ${}^{48}\mathrm{Ca}$ agrees well with the CREX measurement, whereas that of ${}^{208}\mathrm{Pb}$ remains in tension with PREX-II. In contrast to earlier studies, we do not observe a clear correlation between the neutron skin thickness of ²⁰⁸Pb and the symmetry energy slope parameter $L$.

I Introduction

Neutron stars (NS) are among the most extreme and fascinating objects in the universe, providing a natural laboratory for studying matter under conditions of immense density and pressure. The mass of a canonical NS is around $1.4\,M_{\odot}$, contained within a radius of about 10 km. This results in matter being compressed to densities several times higher than the normal nuclear density ($\rho_{0}$) Glendenning (1997). The composition and properties of matter at such extreme densities are governed by the equation of state (EoS), which remains uncertain because of the lack of precise theoretical and observational constraints. In principle, the behavior of nuclear matter can be derived directly from quantum chromodynamics (QCD), the fundamental theory of strong interactions. However, QCD becomes highly non-perturbative at low energies and finite baryon densities, making direct, first-principles calculations of nuclear systems too difficult and complex with current techniques.

To overcome this challenge, a wide range of theoretical approaches have been developed to model the EoS of nuclear matter across densities relevant to finite nuclei, neutron stars, and core-collapse supernovae Oertel et al. (2017). Ab initio methods aim to solve the nuclear many-body problem starting from realistic two- and three-nucleon forces derived from chiral effective field theory ($\chi$-EFT) Tews et al. (2013); Hebeler et al. (2013); Lynn et al. (2016); Drischler et al. (2019); Huth et al. (2021). These methods have been used to compute the EoS of neutron-rich matter at densities up to $\sim 2\rho_{0}$. However, they are less reliable at higher densities because of the breakdown of the EFT expansion and the increasing importance of unresolved short-range physics.

To go beyond the limitations of ab initio methods, various phenomenological models have been developed, designed to reproduce experimental observables in finite nuclei while also describing the matter at extreme densities as those found in NS cores. These include both non-relativistic energy density functionals, such as Skyrme and Gogny interactions, and relativistic density functionals based on the meson exchange Lagrangian Dutra et al. (2012, 2014); Sun et al. (2024). Phenomenological models also require significantly less computational resources compared to ab initio methods.

Among the relativistic frameworks, the Relativistic Mean-Field (RMF) model, originating from the Walecka model and extended to include nonlinear and cross-coupling terms, has been especially successfulGlendenning (1997); Walecka (1974); Serot and Walecka (1992). RMF models are Lorentz covariant and naturally extendable to finite temperatures, suitable for supernova and NS merger simulations.

In some modern variants, density-dependent couplings are introduced as an alternative to adding higher-order meson terms Fuchs et al. (1995); Typel and Wolter (1999); Lalazissis et al. (2005); Typel (2018); Gogelein et al. (2008). Then there are meta-modeling approaches that employ a model-agnostic expansion of the EoS around nuclear saturation density Margueron et al. (2018a, b). Both non-relativistic and relativistic variants of meta-modeling exist and are often used within Bayesian frameworks to study correlations between nuclear matter parameters and neutron star observables Dinh Thi et al. (2021); Mondal and Gulminelli (2023); Traversi et al. (2020); Zhu et al. (2023); Malik et al. (2022); Beznogov and Raduta (2023, 2024); Char et al. (2023); Scurto et al. (2024); Passarella et al. (2025); Biswas et al. (2019, 2021).

Moreover, there are several approaches that utilize flexible, semi-empirical parameterizations of the EoS, designed to explore a wide range of behaviors while satisfying physical constraints. These include piecewise polytropes Read et al. (2009), spectral representations Lindblom (2010); Lindblom and Indik (2012, 2014) and speed of sound parameterizations Greif et al. (2019). Machine learning techniques have also gained significant attention in recent years for EoS inference, offering flexible, data-driven models that can capture complex behaviors without relying on fixed functional forms Landry and Essick (2019); Essick et al. (2020); Soma et al. (2022); Han et al. (2023); Zhou et al. (2023); Cuceu and Robles (2025). Although still in development, these methods show promise in capturing complex EoS behaviors that may not be easily represented in traditional parametric forms.

A variety of observational and theoretical inputs are essential for constraining the nuclear equation of state. At low densities, the chiral effective field theory ($\chi$-EFT) provides reliable predictions, anchoring the EoS near saturation density. Recent advancements in astrophysical observations have significantly improved our ability to probe the internal structure of NS and, consequently, the EoS governing dense matter. At higher densities, gravitational wave observations from binary neutron star mergers (e.g., GW170817) constrain the tidal deformability Abbott et al. (2017, 2019, 2018a), while NICER measurements of pulsars such as PSR J0030+0451Riley et al. (2019); Miller et al. (2019) and PSR J0740+6620Miller et al. (2021); Riley et al. (2021) yield simultaneous mass–radius estimates. Enhanced observational techniques and advanced data analysis methods have accelerated NICER’s mass-radius measurements of pulsars. In the last year alone, NICER has delivered three new radius measurements (PSR J0437-4715Choudhury et al. (2024), PSR J1231-1411Salmi et al. (2024), and PSR J0614-3329Mauviard et al. (2025)).

Because mass-radius and tidal deformability map directly onto the EoS, these measurements tightly constrain its behavior over the relevant density range. Additionally, kilonova emissions from neutron star mergers provide emerging constraints by linking the EoS to the properties of ejected matterLund et al. (2025), and radio pulsar timing has confirmed neutron stars with masses above 2 $M_{\odot}$, setting firm lower limits on the stiffness of the EoS Antoniadis et al. (2013); Fonseca et al. (2016, 2021).

Beyond astrophysical data, theory also provides an anchor at asymptotically high densities, where quark degrees of freedom become relevant. In this regime, perturbative QCD (pQCD) calculations can be used to estimate the EoS of cold, dense matterGorda et al. (2023). Although pQCD is not applicable at neutron star densities, it provides an essential boundary condition for models attempting to extrapolate the EoS toward high-density limits.

An important laboratory observable that connects the nuclear structure to the dense matter EoS is the neutron skin thickness, defined as the difference between the neutron and proton root mean squared radii in a nucleus. It provides important insight into the iso-spin dependence of the nuclear force. This quantity is found to be sensitive to the density dependence of the symmetry energy, especially its slope parameter $L$ at nuclear saturation density Viñas et al. (2014); Reinhard and Nazarewicz (2016). A larger $L$ corresponds to a stiffer symmetry energy, which tends to push the neutrons further out, leading to a thicker neutron skin.

Recent measurements from the PREX-II and CREX experiments have provided model-independent determinations of the neutron skin in ${}^{208}\mathrm{Pb}$ and ${}^{48}\mathrm{Ca}$, respectively, using parity-violating electron scatteringAdhikari et al. (2021, 2022). While PREX-II suggested a relatively thick neutron skin in ${}^{208}\mathrm{Pb}$ pointing towards a larger value of $L$Adhikari et al. (2021), the CREX result for ${}^{48}\mathrm{Ca}$ reported a rather smaller skin thickness, favoring a softer symmetry energy around saturation densityAdhikari et al. (2022). The tension between these two measurements has generated significant theoretical interestReed et al. (2024); Reinhard et al. (2022), as it may point to a more complex density dependence of the symmetry energy or to limitations in existing nuclear models Kumar et al. (2024). Together, PREX-II and CREX provide complementary constraints on the EoS probing the isovector sector of nuclear interactions from both heavy and medium-mass nuclei. Since the symmetry energy and its slope influence the neutron star radii and the pressure of neutron-rich matter, accurate neutron skin measurements from PREX-II and CREX provide constraints on the EoS that are complementary to astrophysical observations.

Given the growing number of theoretical, observational, and laboratory constraints, Bayesian inference has emerged as a powerful and widely adopted tool to extract meaningful constraints on the EoS. In this work, we use Bayesian parameter estimation to explore the parameter space and constrain the EoS in light of the latest astrophysical observations. In particular, we use the uncertainties of the nuclear empirical parameters to generate a large sample of EoSs using a nucleonic RMF model. We impose further constraints coming from $\chi$-EFT calculations, and combined astrophysical observations from GW, radio, and X-ray observations. Finally, we use our posterior samples for the RMF coupling parameters to calculate the neutron skin thickness for ²⁰⁸Pb and ⁴⁸Ca. Then, we compare the skin thickness distributions coming from our calculations with the experimental results from PREX-II and CREX. This is the first time the Bayesian posteriors for the RMF parameters, informed by ab initio calculations and the multimessenger observations, have been used to calculate the distribution of the neutron skin thickness.

This paper is organized as follows. In the next section, we discuss the EoS model and its parameters. Section III describes the details of the Bayesian framework that we have adopted and various theoretical, experimental, and observational data that we have utilized to constrain the EoS parameters. In section IV we discuss our findings, and finally we conclude with section V. In this article, we use natural units, i.e., $\hbar=c=1$.

II Formalism

II.1 Equation of State

We employ a nonlinear finite-range RMF model, which comes under the class of Walecka-type models. In these models, the baryons interact among themselves via the exchange of mesons, which in our case are the scalar-isoscalar ($\sigma$), the vector-isoscalar ($\omega$), and the vector-isovector ($\rho$) mesons. The Lagrangian density of the model is given by Glendenning (1997); Nandi and Pal (2021); Dutra et al. (2014):

$${\cal L}= \sum_{N}\bar{\psi}_{N}\left[\gamma^{\mu}\left(i\partial_{\mu}-g_{\omega}\omega_{\mu}-\frac{1}{2}g_{\rho}\bm{\tau\cdot\rho_{\mu}}\right)-\left(M_{N}-g_{\sigma}\sigma\right)\right]\psi_{N}+\frac{1}{2}\left(\partial_{\mu}\sigma\partial^{\mu}\sigma-m_{\sigma}^{2}\sigma^{2}\right)-\frac{1}{3}bm_{N}(g_{\sigma}\sigma)^{3}-\frac{1}{4}c(g_{\sigma}\sigma)^{4}$$ (1)

$$-\frac{1}{4}\omega_{\mu\nu}\omega^{\mu\nu}+\frac{1}{2}m_{\omega}^{2}\omega_{\mu}\omega^{\mu}-\frac{1}{4}\bm{\rho_{\mu\nu}\cdot\rho^{\mu\nu}}+\frac{1}{2}m_{\rho}^{2}\bm{\rho_{\mu}\cdot\rho^{\mu}}+\Lambda_{\omega\rho}g_{\omega}^{2}g_{\rho}^{2}\omega_{\mu}\omega^{\mu}\bm{\rho_{\mu}\cdot\rho^{\mu}}$$

where, $\psi_{N}$ is the isospin doublet of nucleons. Starting from the above Lagrangian density, we calculate the EoS of the nucleonic matter, relevant for the NS core. However, an NS core also contains leptons ($e^{-}$ and $\mu^{-}$) that maintain beta equilibrium with nucleons. Therefore, we need to add their contributions to the EoS.

For the crust, we use the standard SLy4 EoS Douchin and Haensel (2001). The crust-core transition density is determined by calculating the quantity $K_{\mu}$, the compressibility of npe$\mu$ matter under constant chemical potential:

$$\begin{split}K_{\mu}&=\rho^{2}\frac{d^{2}E_{0}}{d\rho^{2}}+2\rho\frac{dE_{0}}{d\rho}\\ &+\delta^{2}\left[\rho^{2}\frac{d^{2}E_{\text{sym}}}{d\rho^{2}}+2\rho\frac{dE_{\text{sym}}}{d\rho}-2E_{\text{sym}}^{-1}\left(\rho\frac{dE_{\text{sym}}}{d\rho}\right)^{2}\right]\end{split}$$ (2)

The meaning of different symbols in the above equation is given in the next section. As the density decreases, the nuclear matter becomes unstable against the formation of clusters at a certain density, making $K_{\mu}$ negative. This density corresponds to the crust-core transition point. During joining the crust, we also ensure that there is no pressure jump and causality violation at the crust-core boundary.

Table 1: Prior Distributions for the Six EoS Parameters

Parameter	Unit	Prior Distribution
$\rho_{0}$	fm-3	$\mathcal{N}(0.152,\,0.004^{2})$
$E_{0}$	MeV	$\mathcal{N}(-16.17,\,0.36^{2})$
$M^{*}/M$	—	$\mathcal{U}[0.60,\,0.80]$
$K$	MeV	$\mathcal{U}[170,\,320]$
$E_{\rm sym}$	MeV	$\mathcal{U}[25,\,40]$
$L$	MeV	$\mathcal{U}[30,\,130]$

II.2 Parameters

The energy per nucleon $E(\rho,\delta)$ at nucleon density $\rho$, isospin asymmetry $\delta\equiv(\rho_{n}-\rho_{p})/\rho$ is given by,

$$E(\rho,\delta)=E_{0}(\rho)+E_{\rm{sym}}(\rho)\cdot\delta^{2}+\mathcal{O}(\delta^{4}),$$ (3)

Here, $E_{0}(\rho)$ is the energy per nucleon in symmetric nuclear matter (SNM), and $E_{\rm{sym}}(\rho)$ is the symmetry energy which can be expanded in Taylor series about $\rho_{0}$ as:

$$E_{0}(\rho) =E_{0}(\rho_{0})+\frac{K_{0}}{2}\chi^{2}+\dots$$ (4)

$$E_{\rm{sym}}(\rho) =E_{\rm{sym}}(\rho_{0})+L\,\raisebox{1.72218pt}{$\chi$}$$

$$\quad+\frac{K_{\rm{sym}}}{2}\chi^{2}+\dots$$ (5)

where the dimensionless parameter $\chi$ is defined as

$$\chi=\left(\frac{\rho-\rho_{0}}{3\rho_{0}}\right)$$

In the expression for $E_{0}(\rho)$, the quantity $K_{0}$ is the incompressibility of SNM at saturation density which is given by

$$K_{0}=9\rho_{0}^{2}\left(\frac{\partial^{2}E_{0}(\rho)}{\partial\rho^{2}}\right)_{\!0}\,,$$ (6)

where $()_{0}$ indicates that the quantity in bracket is evaluated at saturation density. The other three parameters in the expression for $E_{\rm{sym}}(\rho)$ are the magnitude $E_{\rm{sym}}(\rho_{0})$, slope $L$ and curvature $K_{\rm{sym}}$ of nuclear symmetry energy at saturation density, respectively. The expressions for $L$ and $K_{\rm{sym}}$ are given by

$$L=3\rho_{0}\left(\frac{\partial E_{\rm{sym}}(\rho)}{\partial\rho}\right)_{\!0}$$ (7)

and

$$K_{\rm{sym}}=9\rho_{0}^{2}\left(\frac{\partial^{2}E_{\rm{sym}}(\rho)}{\partial\rho^{2}}\right)_{\!0}$$ (8)

We consider six saturation properties of nuclear matter: the saturation density $\rho_{0}$, energy per nucleon $E_{0}(\rho_{0})$, the effective nucleon mass $M^{*}(=M_{N}-g_{\sigma}\sigma$), the incompressibility $K_{0}$, the symmetry energy $E_{\rm{sym}}(\rho_{0})$ and the slope of the symmetry energy $L$. The six coupling constants of the theory ($g_{\sigma},g_{\omega},g_{\rho},b,c,\Lambda_{\omega\rho}$) can be related algebraically to these six saturation properties. Given a set of $\rho_{0}$, $E_{0}(\rho_{0})$, $M^{*}$, and $K_{0}$, we calculate the couplings related to the isoscalar sector ($g_{\sigma},g_{\omega},b,$ and $c$) by following the prescription given in Ref. Glendenning (1997). Once the couplings in the isoscalar sector have been fixed, the isovector couplings, $g_{\rho}$ and $\Lambda_{\omega\rho}$, can be determined from $E_{\text{sym}}(\rho_{0})$ and $L$ as shown in Ref.Chen and Piekarewicz (2014) (see also Hornick et al. (2018)). To derive analytical expressions for these couplings, one begins with the density-dependent expression for the symmetry energy, which, for the Lagrangian density given in Eq. (1), can be written as:

$$E_{\rm{sym}}(\rho) = \frac{k_{\rm F}^{2}}{6E_{\rm F}}+\frac{g_{\rho}^{2}\,\rho}{8m_{\rho}^{\ast 2}}$$ (9)

$$= E_{{\rm sym},0}(\rho)+E_{{\rm sym},1}(\rho)$$ (10)

where

$$E_{{\rm sym},0}(\rho)=\frac{k_{\rm F}^{2}}{6E_{\rm F}},$$ (11)

and

$$E_{{\rm sym},1}(\rho)=\frac{g_{\rho}^{2}\,\rho}{8m_{\rho}^{\ast 2}}$$ (12)

with

$$\frac{m_{\rho}^{\ast 2}}{g_{\rho}^{2}}\equiv\frac{m_{\rho}^{2}}{g_{\rho}^{2}}+2\Lambda_{\omega\rho}g_{\omega}^{2}\omega_{0}^{2}\,,$$ (13)

represent the isoscalar and isovector contributions to the symmetry energy, respectively. Here $k_{F}=(1.5\pi^{2}\rho)^{1/3}$ is the Fermi momentum of the nucleon and $E_{F}=\sqrt{k_{F}^{2}+M^{*2}}$ is the corresponding Fermi energy.

Now writing $J=E_{\rm sym}(\rho_{0})$ and using the definition of $L$ (eq. (7)) we can split them into isoscalar and isovector parts as:

$$J=J_{0}+J_{1}\;\;{\rm and}\;\;L=L_{0}+L_{1}\,,$$ (14)

where

$$J_{0} =E_{\rm{sym},0}(\rho_{0})=\left(\frac{k_{\rm F}^{2}}{6E_{\rm F}}\right)_{\!0}\,,$$ (15)

$$J_{1} =E_{\rm{sym},1}(\rho_{0})=\left(\frac{g_{\rho}^{2}\rho}{8m_{\rho}^{\ast 2}}\right)_{\!0}\,,$$ (16)

$$L_{0} =3\rho_{0}\!\left(\frac{dE_{\rm{sym},0}}{d\rho}\right)_{\!0}\,,$$ (17)

$$L_{1} =3\rho_{0}\!\left(\frac{dE_{\rm{sym},1}}{d\rho}\right)_{\!0}\,,$$ (18)

Since $J_{0}$ is already known from the saturation properties related to the isoscalar sector, we can readily calculate $J_{1}$ with a given value of $J$ as:

$$J_{1}=J-J_{0}.$$ (19)

A similar procedure can be followed to calculate $L_{1}$ with a given $L$:

$$L_{1}=L-L_{0}$$ (20)

where

$$L_{0} =3\rho_{0}\!\left(\frac{dE_{\rm{sym},0}}{d\rho}\right)_{\!\!0}$$

$$=J_{0}\Bigg(1+\frac{M^{\ast 2}}{E_{\rm F}^{2}}\left[1-\!\frac{3\rho}{\,M^{\ast}}\!\left(\frac{\partial M^{\ast}}{\partial\rho}\right)\right]\Bigg)_{\!0}\,,$$ (21)

with

$$\left(\frac{{\partial M^{\ast}}}{\partial\rho}\right)_{\!\!0} = -\Bigg[\frac{M^{\ast}}{E_{\rm F}}\left(\frac{m_{\sigma}^{\ast 2}}{g_{\sigma}^{2}}+\rho_{\rm s}^{\,\prime}(M^{\ast})\right)^{-1}\Bigg]_{\!0}\,,$$

$$\frac{m_{\sigma}^{\ast 2}}{g_{\sigma}^{2}} \equiv \frac{m_{\sigma}^{2}}{g_{\sigma}^{2}}+2bm_{N}g_{\sigma}\sigma+3cg_{\sigma}^{2}\sigma^{2}\,,$$

$$\rho_{\rm s}^{\,\prime}(M^{\ast}) = \left(\frac{\partial\rho_{\rm s}}{\partial M^{\ast}}\right)$$

$$= \frac{1}{\pi^{2}}\left[\frac{k_{\rm F}}{E_{\rm F}}(E_{\rm F}^{2}+2M^{\ast 2})-3M^{\ast 2}\ln\left(\frac{k_{\rm F}+E_{\rm F}}{M^{\ast}}\right)\right].$$

can again be determined from the isoscalar sector.

Once $L_{1}$ is known, we can calculate $\Lambda_{\omega\rho}$ from eq.(18) as:

$$L_{1} =3\rho_{0}\!\left(\frac{dE_{\rm{sym},1}}{d\rho}\right)_{\!\!0}$$

$$=3J_{1}\left[1-32\left(\frac{g_{\omega}^{2}}{m_{\omega}^{2}}\right)\!g_{\omega}\omega_{0}\Lambda_{\omega\rho}J_{1}\right]_{0}$$

Or,

$$\Lambda_{\omega\rho}=\frac{3J_{1}-L_{1}}{96J_{1}^{2}(g_{\omega}/m_{\omega})^{2}g_{\omega}\omega_{0}}$$ (22)

Finally, we can determine $g_{\rho}$ by putting the value of $\Lambda_{\omega\rho}$ in eq. (13) as:

$$\frac{m_{\rho}^{2}}{g_{\rho}^{2}}=\frac{m_{\rho}^{\ast 2}}{g_{\rho}^{2}}-2\Lambda_{\omega\rho}g_{\omega}^{2}\omega_{0}^{2}=\frac{\rho_{0}}{8J_{1}}-2\Lambda_{\omega\rho}g_{\omega}^{2}\omega_{0}^{2}\,.$$

Or,

$$g_{\rho}^{2}=\frac{m_{\rho}^{2}}{\rho_{0}/8J_{1}-2\Lambda_{\omega\rho}g_{\omega}^{2}\omega_{0}^{2}}\,.$$ (23)

In this study, we adopted the following values for the masses of nucleons and mesons: $M_{N}=939$ MeV $m_{\sigma}=550$ MeV, $m_{\omega}=782.5$ MeV, and $m_{\rho}=763$ MeV.

III Constraints and Bayesian Inference

III.1 Bayesian Inference

In the Bayesian framework, prior knowledge about model parameters is formally incorporated through a prior distribution, which is then updated using the likelihood of the observed data. The result is the posterior distribution, which is obtained via Bayes’ theorem:

$$P(\theta|D)=\frac{P(D|\theta)P(\theta)}{\int P(D|\theta)P(\theta)d\theta}\,.$$

Here, $P(\theta)$ is the prior probability, $P(D|\theta)$ is the likelihood function, $P(D)=\int P(D|\theta)P(\theta)d\theta$ is the normalization constant also called the evidence, and $P(\theta|D)$ is the posterior probability of the model parameters given the dataset $D$. The evidence can be computed by integrating the numerator over the entire parameter space. But as the dimension of the parameter space increases, it becomes computationally inefficient. However, we can overcome this difficulty by adopting statistical computation techniques such as Markov Chain Monte Carlo (MCMC) or Nested Sampling.

In this work, we’re using Bayesian statistics for parameter estimation via PymultinestBuchner et al. (2014), which is based on a nested sampling algorithm. While MCMC can struggle with complex or multi-modal distributions, nested sampling, by design, can more efficiently map out different regions of high likelihood, providing a better global view of the parameter space. This method not only provides the posterior distribution but also accurately computes the Bayesian evidence.

III.2 Prior

As discussed in subsection II.2, we have the six saturation properties of nuclear matter: the saturation density $\rho_{0}$, energy per baryon $E_{0}$, the effective nucleon mass $M^{*}$, the incompressibility $K_{0}$, the symmetry energy $E_{\text{sym}}(\rho_{0})$ and the slope of the symmetry energy $L$ as parameters for our work. The ranges and the type of prior we’re using is given in Table 1. Since the first two parameters control mainly the low-density behavior of EoS, these two are well constrained from nuclear theory and experiments, and therefore, we use gaussian priors for them. For the remaining less-constrained parameters, we use uninformative priors between the specified ranges, i.e., we sample each parameter uniformly between its minimum and maximum value.

III.3 Constraints

III.3.1 $\chi$-EFT

$\chi$-EFT is rooted in Quantum Chromodynamics (QCD) and describes the interactions of nucleons and pions. It is effective at low-energy scales. For our work, we use constraints provided by Drischler et al. (seeDrischler et al. (2016)Drischler et al. (2019)Drischler et al. (2021)). They provide uncertainty bands for the EoS at each density upto 1.12 times the saturation density. This can be used to constrain the low-density part of the EoS.

III.3.2 Maximum Mass

The precise measurement of the neutron star PSR J0740+6620, with a mass of $2.08\pm 0.07\,M_{\odot}$, currently represents the highest reliably determined mass of any known neutron starRiley et al. (2021). This serves as a critical constraint on the EoS, as it sets a firm lower bound on the pressure required to support such a massive object against gravitational collapse. An EoS that becomes too soft at high densities would fail to support a star of this mass, and is therefore ruled out by this observation. As such, the existence of PSR J0740+6620 imposes one of the most stringent constraints on the stiffness of the EoS at supranuclear densities.

III.3.3 GW170817

On August 17, 2017, at 12:41:04 UTC, the Advanced LIGO and Advanced Virgo gravitational-wave detectors made their first observation of a binary neutron star inspiral Abbott et al. (2017, 2018b). The masses and tidal deformabilities of both stars were estimated from this detection. For each EoS generated in the Bayesian analysis, we simulate a hundred binaries by sampling the primary mass $m_{1}$ from a uniform distribution and fixing the companion mass $m_{2}$ satisfying the observed chirp mass of GW170817. The corresponding tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ are obtained from the EoS. These values are then supplied to joint probability distribution for $(m_{1},m_{2},\Lambda_{1},\Lambda_{2})$, which was constructed from the GW170817 data file. Repeating this process for all EoS gives the likelihood corresponding to this event.

III.3.4 NICER

NICER provides mass and radius measurements for pulsars that provide strong constraints on the EoS. The headline 68% credible intervals for the five pulsar measurements are summarized in Table 2.

Although mass-radius data are available for five pulsars to date, we only incorporate the first three and the latest one in our likelihood calculations. We omit NICER 4 due to ambiguity in its mass-radius data. Salmi et al. (2024).

Table 2: NICER mass and radius measurements

NICER	Pulsar	Mass($M_{\odot}$)	Radius(km)	Ref.
1	PSR J0030+0451	$1.44^{+0.15}_{-0.14}$	$13.02^{+1.24}_{-1.06}$	Riley et al. (2019)
2	PSR J0740+6620	$2.08\pm 0.07$	$13.7^{+2.6}_{-1.5}$	Riley et al. (2021)
3	PSR J0437–4715	$1.418\pm 0.037$	$11.36^{+0.95}_{-0.63}$	Choudhury et al. (2024)
4†	PSR J1231–1411	$1.04^{+0.05}_{-0.03}$	$12.6\pm 0.3$	Salmi et al. (2024)
5	PSR J0614–3329	$1.44^{+0.06}_{-0.07}$	$10.29^{+1.01}_{-0.86}$	Mauviard et al. (2025)

^∗NICER numbers are shorthand labels used in this work for convenience; the Pulsar column lists the original name of the pulsars
^† NICER 4 is not used in this work.

IV Results and Discussion

We have used nuclear saturation properties as our parameters instead of the couplings in the Lagrangian. This choice offers several advantages. Saturation properties are physically meaningful and directly constrained by experimental data from finite nuclei, unlike the more abstract meson–nucleon and meson-meson couplings. Moreover, this parametrization allows for a clearer connection between nuclear observables and the EoS, and is better suited for Bayesian inference, where physically motivated priors can be imposed and interpreted more directly.

We impose our observational and theoretical constraints in four sequential steps, building a progressively tighter posterior ensemble of equations of state(EoSs). First, we select all EoSs, that satisfy chiral effective field theory at low densities and reproduce at least the maximum observed neutron‑star mass (PSR J0740+6620). This defines our $\chi$-EFT+M_max baseline. Second, we further restrict this ensemble by applying the tidal deformability bounds from GW170817. Third, we jointly enforce the mass–radius posteriors from the first two NICER targets NICER 1 and 2, as well as the more recent measurement from NICER 3. Finally, we apply the latest NICER measurement NICER 5 to arrive at our most tightly constrained set. At each step the same ordering appears in the figure legends, allowing us to see precisely how each new piece of data refines our knowledge of the neutron‑star EoS.

We do not include NICER 4 (PSR J1231-1411) in our analysis due to the strong dependence of its inferred radius on the assumed priors. The results vary significantly based on whether informative or uninformative priors are used, with the radius ranging from $12.6\pm 0.3$ km to $13.5^{+0.3}_{-0.5}$ km. Moreover, when broader radius priors are applied, the models fail to fit the data well. This sensitivity to prior assumptions, along with the model-dependent nature of the inferred geometry of the hot spots, makes the constraints from PSR J1231 less robust. It is worth noting, however, that some recent studies do incorporate PSR J1231 in their analysis Li et al. (2025a, b), but we choose not to do it in this study.

IV.1 Saturation Properties

Figure 1: Corner plot depicting the posterior distributions of the nuclear saturation parameters. The diagonal panels display the one-dimensional marginalized distributions for each parameter, while the off-diagonal panels show the 90% contours for the two-dimensional joint probability density functions. The colored contours and distributions represent successive constraints applied in the order indicated by the legend. On top of each diagonal panels we also show the 68% credible interval for each parameter in the order indicated by the legend.

We construct the joint probability density functions (PDFs) of the parameters $\rho_{0}$, $E_{0}$, $M^{*}$, $K$, $E_{\text{sym}}$, $L$ using posterior samples obtained from the chosen prior distributions and the likelihood function. Figure 1 shows the joint and marginalized posterior distributions of these parameters, with successive application of constraints from nuclear theory and astrophysical observations. The diagonal panels show the marginalized one-dimensional probability distributions for each parameter, while the off-diagonal panels represent the joint PDFs between pairs of parameters, illustrating their correlations.

Table 3: Posterior 90% credible intervals for the six EoS parameters

Parameters	Unit	$\chi$-EFT+ $\text{M}_{\text{max}}$	+GW170817	+ NICER 1,2	+ NICER 3	+ NICER 5
$\rho_{0}$	$\mathrm{fm}^{-3}$	$0.153_{-0.007}^{+0.007}$	$0.153_{-0.007}^{+0.007}$	$0.153_{-0.007}^{+0.007}$	$0.153_{-0.007}^{+0.007}$	$0.154_{-0.006}^{+0.007}$
$E_{0}$	$\mathrm{MeV}$	$-16.157_{-0.597}^{+0.595}$	$-16.148_{-0.603}^{+0.586}$	$-16.167_{-0.595}^{+0.589}$	$-16.169_{-0.591}^{+0.596}$	$-16.172_{-0.586}^{+0.596}$
$M^{*}/M$	$-$	$0.705_{-0.080}^{+0.058}$	$0.725_{-0.068}^{+0.044}$	$0.717_{-0.059}^{+0.041}$	$0.724_{-0.053}^{+0.037}$	$0.732_{-0.044}^{+0.032}$
$K$	$\mathrm{MeV}$	$260.480_{-35.364}^{+54.002}$	$252.967_{-31.690}^{+57.774}$	$253.869_{-29.762}^{+57.498}$	$250.626_{-29.024}^{+59.213}$	$243.619_{-25.036}^{+61.378}$
$E_{\mathrm{sym}}$	$\mathrm{MeV}$	$28.929_{-2.231}^{+1.955}$	$28.875_{-2.190}^{+1.853}$	$28.928_{-2.196}^{+1.804}$	$28.879_{-2.170}^{+1.936}$	$28.775_{-2.168}^{+1.878}$
$L$	$\mathrm{MeV}$	$55.963_{-7.221}^{+7.361}$	$54.600_{-6.820}^{+7.367}$	$54.639_{-6.533}^{+7.197}$	$53.909_{-6.489}^{+7.112}$	$52.572_{-6.231}^{+7.160}$

The 90% posterior ranges of the nuclear saturation parameters are given in Table 3. We also show the 68% credible interval for each parameter on top of each diagonal panel of Figure 1. From both sets of data, it is evident that the effective mass $M^{*}$, the symmetry energy $E_{\text{sym}}$, and its slope $L$ at saturation are remarkably constrained. Although the constraint on incompressibility $K$ at saturation is broader, it still shows a significant improvement relative to its chosen prior range. Overall, we can say that we have very tight constraints on all of our parameters.

Let us first analyze the marginal distributions of Figure 1. The posterior median for $\rho_{0}$ sits at approximately 0.153 $\mathrm{fm}^{-3}$ across all datasets till NICER 3. The addition of NICER 5 shifts it towards a little high value but it still remains under the baseline uncertainty. The very stability of $\rho_{0}$ demonstrates that astrophysical data to date-GW170817 and NICER-act as consistency checks on saturation density rather than refining it appreciably beyond what is known from nuclear theory and experiments.

The median values of $E_{0}$ remain remarkably stable across all cases ranging narrowly from -16.148 $\mathrm{MeV}$ (+ GW170817) to -16.172 $\mathrm{MeV}$ (+ NICER-5) indicating that the underlying nuclear physics, particularly the chiral effective field theory ($\chi$-EFT) input combined with the observed maximum neutron star mass, already constrains $E_{0}$ tightly. The addition of +NICER 5 gives the highest median (-16.172 $\mathrm{MeV}$). However, the shift is only $\sim 0.015$ $\mathrm{MeV}$ upward from the baseline, which is again tiny compared to the baseline uncertainty. It does hint that the most recent NICER release might be pushing the posterior upward, but not yet at a statistically significant level. The $68\%$ credible intervals remain nearly constant in width ($\pm\,0.36$ $\mathrm{MeV}$), even after incorporating constraints from GW170817 and multiple NICER data releases. This stability suggests that the astrophysical measurements are fully consistent with the nuclear-theory predictions but are not yet precise enough to significantly reduce the uncertainty in $E_{0}$.

The posterior median of the nucleon effective mass shifts to higher values when astrophysical data from GW170817 and NICER (1, 2, 3, and 5) are incorporated, increasing from 0.705 to 0.732. Since a higher effective mass generally leads to a softer equation of state (EoS), this trend suggests a preference for softer EoSs under most constraint sets. Among these, only NICER 1,2 slightly favors a lower effective mass, consistent with a stiffer EoS. Notably, NICER 5 yields the highest median effective mass, indicating the strongest preference for a soft EoS. In addition to this systematic shift in the median, the associated uncertainties are also reduced as astrophysical observations are applied.

Starting from a baseline of $K=260.48\,\mathrm{MeV}$ ($\chi$-EFT+M_max), incorporating GW170817 shifts the median downward to 252.97 $\mathrm{MeV}$. Since a lower incompressibility generally leads to a softer equation of state (EoS) it signals that tidal‐deformability measurements favor a softer EoS. Also successive NICER radius measurements further pull the central value to 250.63 $\mathrm{MeV}$ (NICER-3). Among these as well, only NICER 1,2 favors a slightly higher incompressibility, consistent with a stiffer EoS. The most recent NICER-5 release has the strongest impact, driving the median to 243.62 $\mathrm{MeV}$ and contracting the 68% interval by roughly $\approx 17.6\%$ as compared to the baseline, underscoring the progressive power of NICER mass-radius measurements in softening the EoS and sharpening our estimate of nuclear incompressibility.

In our analysis the symmetry energy $E_{\rm sym}$ is driven from a broad, flat prior of $25–40\ \mathrm{MeV}$ to a tightly peaked posterior at around $29\ \mathrm{MeV}$, with very little change as we impose the constraints such as GW170817 and NICER in succession. Concretely, $\chi$-EFT+M_max already yields $E_{\rm sym}\approx 28.93$ $\mathrm{MeV}$, and adding GW or any of the NICER datasets shifts the median by no more than ±0.15 $\mathrm{MeV}$ and leaves the credible bounds essentially intact. Still it’s important to note that the GW and NICER measurements favour a slightly lower symmetry energy at saturation (except NICER 1,2) especially NICER 5. The lack of significant shift or tightening underscores that current astrophysical observations namely GW170817 and NICER do not yet surpass nuclear theory and maximum mass constraints on $E_{\mathrm{sym}}$. To put simply, $\chi$-EFT and maximum mass alone impose very tight constraints on symmetry energy.

The slope of the symmetry energy, $L$, is driven from a completely uninformative uniform prior over 30-130 $\mathrm{MeV}$ to a sharply peaked posterior of roughly 56 $\mathrm{MeV}$ with a $\pm 4-5$ $\mathrm{MeV}$ 68% credible interval as soon as we impose the $\chi$-EFT+M_max constraint. Adding the GW170817 tidal-deformability constraint shifts the median downward by $\approx 1.36\ \mathrm{MeV}$, indicating it favors a slightly softer density dependence. Incorporating NICER-1,2 hardly moves the median ($L\approx 54.64\ \mathrm{MeV}$) but trims the lower bound by $\approx 0.16\ \mathrm{MeV}$, showing consistent support without new tightening. NICER-3 then nudges $L$ down to 53.91 $\mathrm{MeV}$ and shrinks the error bars again. Finally NICER-5 exerts the strongest pull, reducing the median by $\approx 1.34\ \mathrm{MeV}$ to 52.57 $\mathrm{MeV}$ and contracting the 68% interval by $\sim 0.6\ \mathrm{MeV}$ on the low side highlighting the constraining power of NICER mass-radius measurements.

Having examined the marginal distributions in detail, we now turn our attention to the correlations between parameters. In the two‐dimensional panels of the corner plot, we show the 90% credible contours to capture the parameter uncertainties and correlations. It is evident from $90\%$ regions of the joint posterior between $\rho_{0}$ and $E_{0}$ that these two quantities are already well constrained within a very small region of the 2D plot and also have almost zero correlation.

In contrast, $\rho_{0}$ vs. $M^{*}/M$ develops a modest anti‐correlation as each constraint is added, especially NICER‐3 and 5 , indicating that higher saturation density favors lower effective mass under tighter mass-radius measurements. Similarly, there is a slight positive correlation between $\rho_{0}$ and $K$ which grows with NICER 3 and NICER 5 data the most. For $\rho_{0}$ and $E_{\rm sym}$ there remains a slight positive correlation with almost no change across all datasets.

The slope $L$ shows a moderate negative correlation with $\rho_{0}$, strongest after NICER‐5, implying that stiffer symmetry energy slope pair with slightly lower saturation density. The $E_{0}$-$E_{\rm sym}$ anti‐correlation set by nuclear theory remains unchanged by astrophysical data, and the mild positive $K$-$E_{\rm sym}$ correlation becomes somewhat more pronounced under successive constraints. Also, if we look at the joint distributions we observe that, as additional constraints are added sequentially, some credible regions significantly decrease in size. This sequential reduction is a result of the contraction of the posterior distributions, i.e., improved parameter estimation and reduced uncertainty as more information is incorporated into the analysis.

IV.2 Couplings

Figure 2: Same as Figure 1 but now for coupling constants.

Table 4: 90% Posterior Credible Intervals for the Coupling Constants

Coupling	$\chi$-EFT+ $\text{M}_{\text{max}}$	+GW170817	+ NICER 1,2	+ NICER 3	+ NICER 5
$g_{\sigma}$	$9.650_{-0.798}^{+0.895}$	$9.440_{-0.680}^{+0.767}$	$9.517_{-0.618}^{+0.652}$	$9.437_{-0.571}^{+0.593}$	$9.349_{-0.525}^{+0.494}$
$g_{\omega}$	$10.511_{-1.322}^{+1.595}$	$10.084_{-1.036}^{+1.395}$	$10.233_{-0.920}^{+1.194}$	$10.066_{-0.820}^{+1.091}$	$9.867_{-0.720}^{+0.924}$
$g_{\rho}$	$9.026_{-1.273}^{+1.010}$	$9.035_{-1.194}^{+0.969}$	$9.044_{-1.192}^{+0.979}$	$9.065_{-1.167}^{+0.985}$	$9.084_{-1.157}^{+0.952}$
$b$	$0.004_{-0.002}^{+0.003}$	$0.005_{-0.002}^{+0.003}$	$0.004_{-0.002}^{+0.002}$	$0.005_{-0.002}^{+0.002}$	$0.005_{-0.002}^{+0.002}$
$c$	$-0.003_{-0.003}^{+0.006}$	$-0.004_{-0.003}^{+0.008}$	$-0.004_{-0.003}^{+0.006}$	$-0.004_{-0.003}^{+0.007}$	$-0.004_{-0.003}^{+0.007}$
$\Lambda_{\omega\rho}$	$0.043_{-0.013}^{+0.022}$	$0.048_{-0.015}^{+0.023}$	$0.047_{-0.014}^{+0.018}$	$0.049_{-0.014}^{+0.019}$	$0.054_{-0.016}^{+0.020}$

The coupling constants $g_{\sigma}$, $g_{\omega}$, $g_{\rho}$, $b$, $c$, and $\Lambda_{\omega\rho}$ in the RMF theory represent the strengths of interactions between nucleons via meson exchanges and self-interactions. Figure 2 shows the posterior distributions of the coupling constants and their correlations and Table 4 records the corresponding 90% credible intervals. We begin by examining the one‑dimensional marginalized distributions.

The trends observed in the coupling constants reflect how the RMF model adjusts to fit both nuclear theory and astrophysical observations. The scalar coupling $g_{\sigma}$ shows a noticeable decrease in median value from 9.650 $\chi$-EFT+M_max to 9.349 (NICER 5), with a steady reduction in uncertainty. The $g_{\omega}$ coupling, which provides the dominant repulsive interaction, also decreases slightly from 10.511 $\chi$-EFT+M_max to 9.867(NICER 5). After a drop with GW170817, NICER 1,2 briefly increases them, then they drop again with NICER 3 and NICER 5. NICER 5 plays the most significant role in this suppression.

The stability of $g_{\rho}$ indicates that the isospin asymmetry, important for neutron stars with large neutron-proton differences, is well-captured by the baseline and does not require significant modification but still there is a gradual increase as constraints are added on top. This points to a slightly stronger isospin force being preferred as more precise constraints are applied. The small, stable values of $b$ and $c$ highlight the role of nonlinear self-interactions in fine-tuning the equation of state, particularly the incompressibility, which affects how nuclear matter responds to compression. These are already well constrained by nuclear theory and are not strongly impacted by astrophysical measurements.

The increase in $\Lambda_{\omega\rho}$ is particularly noteworthy, as it suggests that the interaction between the $\omega$-meson and $\rho$-meson which controls the density dependence of the symmetry energy becomes more significant when astrophysical data are included. It shows a clear upward trend, from 0.043 at baseline to 0.054 with NICER 5, along with steadily narrowing uncertainties. The sensitivity of $\Lambda_{\omega\rho}$ to these constraints highlights the growing role of NICER data in shaping the isovector channel’s behavior at supranuclear densities.

The joint posterior distributions of the Lagrangian couplings reveal several meaningful parameter correlations, reflecting how different sectors of the model adjust to simultaneously satisfy nuclear theory and astrophysical observations. Most prominent is the very strong positive correlation between the scalar and vector couplings, $g_{\sigma}$ and $g_{\omega}$. The correlation between $g_{\sigma}$ and $g_{\omega}$ comes from the criterion of saturation of nuclear forces for symmetric matter. This correlation is already strong at baseline and becomes more compact, particularly under GW170817 and NICER-5–indicating that increases in scalar attraction must be compensated by enhanced vector repulsion to maintain consistency with the saturation properties and maximum mass constraints.

Within the scalar sector, $g_{\sigma}$ shows a moderate negative correlation with the nonlinear scalar coupling $b$ at the baseline, which gradually weakens as observational constraints are added. A similar trend is observed between $g_{\omega}$ and $b$. But for these a slight negative correlation remains at the end. In contrast, neither $g_{\sigma}$ nor $g_{\omega}$ exhibits meaningful correlation with the quartic scalar coupling $c$ across datasets, though a slight negative trend begins to emerge with NICER-5.

Both $g_{\sigma}$ and $g_{\omega}$ also display weak but growing negative correlations with the mixed isovector coupling $\Lambda_{\omega\rho}$, becoming moderately anti-correlated by NICER-5. This reflects a trade-off between the strength of isoscalar interactions and the density dependence of the symmetry energy, especially under tight mass-radius constraints of NICER 5. Among the nonlinear couplings, $b$ and $c$ are consistently negatively correlated, but this structure remains stable across constraints, with only the sampled area shifting. Meanwhile, $b$ and $\Lambda_{\omega\rho}$ maintain a moderate positive correlation that also persists under constraint tightening.

Together, these evolving correlations among the couplings reveal how the underlying structure of the theory becomes increasingly constrained as more precise constraints are applied. In particular, the narrowing of $g_{\sigma}$–$g_{\omega}$ correlation highlights the key roles of isoscalar and isovector balance in controlling the stiffness of the equation of state under astrophysical observations.

IV.3 Equation of state

Next, we present plots for various physical quantities and observables. To generate these plots, we follow the procedure outlined below. First, we use the entire set of posterior samples to generate the corresponding EoS and mass–radius (MR) relations for each parameter set. Then we fix a discrete set of mass values and, for each mass point, collect the corresponding radii from all the generated MR curves. This results in a distribution of radii at each mass. From these distributions, we compute the desired credible interval—in this case, the 90% interval. Repeating this process across all mass points yields a band representing the 90% credible region in the MR plane. The same method is applied across all successive constraints to generate the final diagram shown in Figure 4. We adopt a similar procedure to generate other similar plots.

Figure 3: (a) Pressure vs baryon number density of matter in $\beta-$equilibrium. The different colored pairs of curves denote the 90 % credible‐interval boundaries of the posterior samples as constraints are successively imposed in the order as shown in the legend.

Let’s now analyze the individual figures one by one. A detailed examination of Figure 3 reveals the progressive impact of successive astrophysical constraints on the pressure–number density relation. The inclusion of GW170817 significantly narrows the 90% credible band, notably shifting its upper boundary to lower pressure values compared to the $\chi$-EFT+$M_{\rm max}$ baseline. It indicates that GW170817 favors softer equations of state (EoSs). The addition of NICER 1 and 2 data further tightens the band, and the inclusion of NICER 3 reduces the area enclosed by the 90% credible contours by approximately one-third relative to the baseline. The most substantial constraint comes from NICER 5, which shifts the upper boundary even closer to the lower one, effectively eliminating nearly half of the previously allowed region. Overall, the sequential application of astrophysical observations systematically narrows the 90% credible region by pushing the upper pressure boundary downward, increasingly favoring softer EoSs. In particular, the latest NICER measurement provides the most stringent constraint, strongly preferring the softest EoSs.

IV.4 Neutron Star Properties

Figure 4: The left panel (a) shows the 90% mass-radius contours and the right panel (b) mass vs tidal deformability contours for neutron stars corresponding to the posterior samples.

Next, we turn to Figure Figure 4, which presents the MR diagram. The widest contour (black) spans radii from approximately 12.3 to 13.5 km at $M=1.4\,M_{\odot}$, reflecting the fact that, in the absence of astrophysical data such as gravitational waves or NICER measurements, the sole requirement to support the most massive known pulsars allows a wide range of EoS stiffness. The inclusion of tidal deformability constraints from GW170817 removes a substantial portion of the MR space, particularly shifting the right boundary toward smaller radii and disfavoring larger radius values. This indicates that GW170817, when combined with $\chi$-EFT and the maximum mass constraint, favors softer EoSs. Adding the first two NICER mass–radius posteriors further trims the allowed region. The inclusion of NICER 3 tightens the contour even more and shifts it slightly leftward again, reinforcing the preference for softer EoSs. Finally, with the addition of NICER 5, the entire 90% credible region shifts further to the left, significantly narrowing the band and strongly favoring the softest EoS models considered. The dramatic reduction in the final contour relative to the baseline $\chi$-EFT+$M_{\rm max}$ region underscores the power of multimessenger observations in determining neutron-star properties.

In Figure Figure 4, we present the mass–tidal deformability ($M$–$\Lambda$) relationship for neutron stars. The applied constraints follow the same sequential order as in previous figures. As before, we observe a progressive narrowing of the tidal deformability bands with each added constraint. Initially, the deformability at each mass exhibits considerable uncertainty. The inclusion of the GW170817 measurement significantly reduces the 90% credible interval, reflecting a strong preference for lower tidal deformability values—indicative of a softer equation of state (EoS). Subsequent incorporation of the NICER 1 and 2 observations further shrinks the 90% region. The addition of NICER 3 and 5 has a comparable effect to that of GW170817, tightening the bands even more and reinforcing the trend toward a softer EoS. At higher masses, the tidal deformability bands begin to widen slightly. It suggests that current nuclear and astrophysical results are less effective at constraining $\Lambda$ for massive stars. Nonetheless, the overall uncertainty is still significantly reduced compared to the initial, unconstrained scenario.

IV.5 Composition of neutron star matter

Figure 5: Same as Fig. 3 but for (a) Proton fraction vs Number density, and (b) Sound Speed squared vs number density.

We show the $90\%$ credible intervals for proton fraction $x_{p}$ as a function of density in Figure 5. Across all datasets, the $x_{p}$ increases monotonically from nearly zero at $n\approx 0.05\,\mathrm{fm}^{-3}$, with the upper boundary of the credible interval band reaching just below 0.15 and the lower boundary slightly above 0.10 by $n=1.0\,\mathrm{fm}^{-3}$. Up to $n\approx 0.2\,\mathrm{fm}^{-3}$, all five $90\%$ credible bands are nearly indistinguishable. Above $n\approx 0.2$ $\mathrm{fm}^{-3}$, all the bands begin to widen, with the baseline $\chi$-EFT+M_max band exhibiting the greatest broadening. Incorporating GW constraints narrows this baseline band by $\sim 10\%$ and shifts it slightly downward at around $n=0.4$ $\mathrm{fm}^{-3}$. While the inclusion of NICER 1 and 2 data has little impact, NICER 3 and 5 significantly tighten the credible bands, highlighting the complementary role of mass-radius measurements in constraining proton fraction at these densities. This information is crucial for determining the threshold of the direct Urca process, which governs the cooling of their neutron stars and their bulk viscosity. At higher densities ($n\sim 0.9$ $\mathrm{fm}^{-3}$), the bands begin to converge again, with little variation between the different datasets.

Since we are working within a relativistic mean-field framework, the theory is inherently causal. It is confirmed in Figure 5, where the sound speed never approaches the speed of light. At low densities, the sound speed remains small and does not vary much between EoSs, as indicated by the narrow 90% credible interval. However, beyond $\sim 0.2$ fm^-3, the sound speed increases rapidly, and the uncertainty band broadens, reaching its maximum width near 0.5 fm^-3. Beyond this point, the sound speed grows more slowly, and the credible band narrows slightly. The impact of successive observational constraints follows a similar trend to that seen in previous figures. Notably, the NICER 5 data impose the most stringent constraint, significantly narrowing the baseline band (defined by $\chi$-EFT and the maximum mass) across the relevant density range.

Figure 6: Corner plot depicting the posterior distributions of core-crust transition density and transition pressure. The diagonal panels display the one-dimensional marginalized distributions for each, while the off-diagonal panels show the 90% contours for the two-dimensional joint probability density functions. The colored contours and distributions represent successive constraints applied in the order indicated by the legend.

Turning to Figure 6, we present a corner plot depicting the core–crust transition density ($\rho_{cc}$) and pressure ($P_{cc}$) of the neutron star EoS, with successive observational constraints applied in the same manner as in previous figures. We observe that both the transition density and pressure exhibit a decrease in their median values following the inclusion of the GW170817 constraint. These values increase again upon the application of NICER 1 and 2 constraints, before decreasing once more with the addition of NICER 3 and NICER 5 data. This trend is particularly evident in the one-dimensional marginal distribution of the transition pressure, where the median shows a substantial drop with the inclusion of GW170817 and NICER 5. Additionally, the 2D plot shows a clear positive correlation between the transition density and transition pressure for all the cases, as expected.

Examining the trends across all the plots studied so far, we observe a consistent pattern in how each constraint influences the EoS with respect to the previous one. The GW170817 constraint favors softer EoSs, NICER 1 and 2 slightly favor stiffer ones, while NICER 3 and 5 again prefer softer EoSs. This trend directly impacts the core–crust transition point in a neutron star. In particular, constraints that favor softer EoSs tend to result in lower transition densities and pressures. Since we use the same crust for all the EoSs, all neutron stars governed by them have the same density and pressure at the surface. As a result, GW170817 supports stars with thinner crusts, NICER 1 and 2 indicate a modestly thicker crust, and NICER 3 and 5 again point to thinner crusts.

IV.6 The 1.4 M_⊙ and the Maximum Mass Neutron Star

Figure 7: The three panels above show the posterior probability distributions of radius (R), dimensionless tidal deformability $(\Lambda)$, and central pressure $(P_{c})$ for a canonical 1.4 $M_{\odot}$ neutron star. The constraints are applied in the order as shown in the legend.

Table 5: 90% Posterior Credible Intervals for Key Physical Properties of a 1.4 $M_{\odot}$ Neutron Star

Physical Quantity	Unit	$\chi$-EFT+ $M_{\mathrm{max}}$	+GW170817	+ NICER 1,2	+ NICER 3	+ NICER 5
$R$	km	$12.889_{-0.383}^{+0.416}$	$12.696_{-0.317}^{+0.337}$	$12.740_{-0.269}^{+0.287}$	$12.652_{-0.259}^{+0.271}$	$12.508_{-0.241}^{+0.257}$
$\Lambda$	–	$579.048_{-114.860}^{+154.847}$	$517.686_{-86.204}^{+110.062}$	$532.906_{-75.304}^{+94.378}$	$506.908_{-69.035}^{+83.901}$	$468.558_{-60.619}^{+72.046}$
$P_{c}$	MeV/fm3	$52.007_{-9.086}^{+9.819}$	$56.799_{-8.089}^{+8.633}$	$55.404_{-6.703}^{+6.926}$	$57.580_{-6.523}^{+6.862}$	$61.080_{-6.360}^{+6.748}$

In Figure 7, we show the posterior distribution of radius, tidal deformability, and central pressure for 1.4$M_{\odot}$ neutron stars. The median radius, based on the baseline analysis incorporating $\chi$-EFT and the maximum mass constraint, is found to be approximately 12.89 km. Upon including the GW170817 constraint, this value experiences the most significant reduction, dropping to around 12.70 km, indicating that GW170817 favors a softer EoS. A softer EoS generates less pressure at a given density. As a result, it leads to more compact stars (i.e., smaller radii for a given mass, $1.4\,M_{\odot}$ here), where the higher central densities can generate sufficient pressure to balance gravity. The inclusion of NICER data sets 1, 2, and 3 does not significantly alter the radius of a 1.4 $M_{\odot}$ neutron star. Nevertheless, they maintain the overall trend observed earlier: the median radius first shows a slight increase to 12.74 km, followed by a mild decrease to 12.65 km. However, the addition of NICER 5 on top of the previous constraints results in a more pronounced reduction in the median radius to approximately 12.50 km, once again indicating a strong preference for a softer EoS. Furthermore, the associated uncertainties are appreciably reduced compared to the baseline scenario, consistent with the trends seen in earlier analyses.

Other physical quantities, such as tidal deformability and central pressure, consistently reflect the same behavior. For instance, if GW170817 favors a softer EoS based on the radius distribution, it similarly supports a softer EoS when examined through the lens of tidal deformability and other related quantities. The median and credible intervals of the three quantities in Figure 7 are shown in Table 5.

Figure 8: One‐dimensional posterior probability densities for (a) the maximum neutron-star mass and (b) their associated central pressure. The constraints are applied in the order as shown in the legend.

Table 6: 90% Posterior Credible Intervals for Key Physical Properties of Maximum Mass Neutron Star

Physical Quantity	Unit	$\chi$-EFT+ $M_{\mathrm{max}}$	+GW170817	+ NICER 1,2	+ NICER 3	+ NICER 5
$M$	$M_{\odot}$	$2.299_{-0.234}^{+0.309}$	$2.217_{-0.179}^{+0.267}$	$2.245_{-0.159}^{+0.229}$	$2.214_{-0.142}^{+0.206}$	$2.174_{-0.123}^{+0.174}$
$R$	km	$11.386_{-0.767}^{+1.095}$	$11.110_{-0.598}^{+0.898}$	$11.194_{-0.521}^{+0.772}$	$11.079_{-0.476}^{+0.691}$	$10.922_{-0.432}^{+0.586}$
$\mathcal{E}_{c}$	MeV/fm3	$1210.880_{-233.366}^{+216.565}$	$1283.168_{-216.267}^{+176.826}$	$1259.228_{-184.379}^{+149.198}$	$1289.506_{-173.916}^{+138.793}$	$1330.866_{-155.701}^{+129.207}$
$P_{c}$	MeV/fm3	$561.085_{-71.098}^{+78.915}$	$581.967_{-59.319}^{+72.306}$	$578.755_{-53.221}^{+61.949}$	$588.550_{-51.671}^{+62.121}$	$606.106_{-54.353}^{+58.972}$

In Figure 8, we show the posterior distributions of mass and central pressure corresponding to the maximum mass neutron star. The inference about EoS stiffness drawn from the 1.4 $M_{\odot}$ neutron star is consistently supported by the maximum mass neutron star plots as well. After incorporating all the constraints, we find the median maximum mass of the neutron star to be $2.17\,M_{\odot}$. In Table 6 we show the posterior 90% credible intervals of the same quantities which we have in Figure 8, as well as radius and central energy density of the maximum mass neutron star.

IV.7 Neutron Skin Thickness

Figure 9: Neutron skin distributions for ${}^{48}\mathrm{Ca}$ and ${}^{208}\mathrm{Pb}$ under different constraint sets.

The posterior distributions of neutron skin thickness ($r_{\text{skin}}$) for Calcium-48 (${}^{48}\mathrm{Ca}$) and Lead-208 (${}^{208}\mathrm{Pb}$), calculated within the RMF model, provide insights into the nuclear EoS and its connection to neutron star properties. The neutron skin thickness, defined as the difference between the neutron and proton root-mean-square radii ($r_{\text{skin}}=<r_{n}>-<r_{p}>$), is sensitive to the density dependence of symmetry energy especially its slope $L$, which governs the behavior of neutron-rich matter in both finite nuclei and neutron stars. The plots (Figure 9) present the 90% credible intervals for $r_{\text{skin}}$ under two constraint sets: one with all constraints, including chiral effective field theory ($\chi$-EFT), maximum neutron star mass ($M_{\text{max}}$), GW170817, and NICER observations (NICER 1, 2, 3, 5), and another with all constraints except $\chi$-EFT. These distributions are compared to experimental results from the Calcium Radius Experiment (CREX) for ${}^{48}\mathrm{Ca}$ ($0.121\pm 0.026\pm 0.024$ fm) and the Lead Radius Experiment (PREX) for ${}^{208}\mathrm{Pb}$ ($0.283\pm 0.071$ fm). This section discusses the trends in these distributions, the effect of $\chi$-EFT constraints, and their implications for nuclear and astrophysics.

IV.7.1 Neutron Skin Thickness of ${}^{48}\mathrm{Ca}$

The posterior distributions for the neutron skin thickness of ${}^{48}\mathrm{Ca}$ are shown in the top panels of Figure 9. In the fully constrained model, including $\chi$-EFT, $M_{\text{max}}$, GW170817, and NICER observations (denoted as +NICER 5), the distribution peaks sharply at approximately 0.17 fm, with a narrow 68% interval (CI). This peak is higher than the CREX experimental mean of approximately 0.12 fm, which has a $1\,\sigma$ range of 0.07–0.17 fm. More than half of the lower half of +NICER 5 distribution falls within the CREX $1\,\sigma$ range, indicating consistency with the experimental result, though with a slight upward shift in the predicted skin thickness. In the model excluding $\chi$-EFT constraints (i.e., including only $M_{\text{max}}$, GW170817, and NICER observations), the distribution peaks at approximately 0.16 fm, with a comparatively broader 68% CI. This peak position is less than that of the fully constrained model, suggesting that the exclusion of $\chi$-EFT has a significant effect on the predicted neutron skin thickness for ${}^{48}\mathrm{Ca}$. Also, now almost complete 68% CI falls in the $1\sigma$ range of CREX but the uncertainty has increased. The slight upward shift relative to the CREX mean in both cases may be attributed to the influence of astrophysical constraints. These constraints favor a moderately soft EoS, which correlates with a slightly larger median neutron skin than CREX in lighter nuclei like ${}^{48}\mathrm{Ca}$.

IV.7.2 Neutron Skin Thickness of ${}^{208}\mathrm{Pb}$

The posterior distributions for the neutron skin thickness of ${}^{208}\mathrm{Pb}$ are shown in the bottom panels of Fig. 11. In the fully constrained model (+NICER 5, including $\chi$-EFT), the distribution peaks sharply at approximately 0.16 fm, with a narrow 68% CI. This value is notably lower than the PREX experimental mean, which is centered around 0.28 fm with a broad $1\,\sigma$ range of 0.21–0.35 fm. The +NICER 5 peak lies outside of the PREX $1\,\sigma$ range, suggesting a potential tension between the model predictions and the PREX result. In the model excluding $\chi$-EFT constraints, the distribution peaks at approximately 0.13 fm, with a comparatively broader 68% CI. This difference indicates that $\chi$-EFT has significant impact on the predicted neutron skin thickness for ${}^{208}\mathrm{Pb}$. The peak position is significantly below the PREX mean in both cases, suggests that the astrophysical constraints are the primary drivers of the reduced skin thickness. The NICER observations, which constrain the EoS through mass-radius measurements, favor a softer EoS that correlates with a smaller neutron skin thickness in heavy nuclei like ${}^{208}\mathrm{Pb}$.

The discrepancy between the model predictions and the PREX result for ${}^{208}\mathrm{Pb}$ is noteworthy. The PREX experiment, which measures the neutron skin via parity-violating electron scattering, reported a relatively large skin thickness, implying a stiff symmetry energy slope ($L$) at saturation density. In contrast, the constrained RMF model, influenced heavily by GW and NICER data, predicts a smaller skin thickness, consistent with a softer EoS and a lower $L$. This tension may indicate limitations in the RMF model.

Figure 10: Neutron skin thickness of ²⁰⁸Pb as a function of $L$.

Finally, in Figure 10, we present the neutron skin thickness of ²⁰⁸Pb as a function of $L$, using the posterior samples obtained from our Bayesian analysis that incorporates both the $\chi$-EFT inputs and all astrophysical constraints. Interestingly, the results do not exhibit a clear correlation between these two quantities. This finding contrasts with earlier studies, which reported a strong correlation between the neutron skin thickness and $L$ Viñas et al. (2014); Essick et al. (2021).

V Conclusion

In this study, we employed a relativistic mean-field model to investigate and constrain the EoS of neutron stars. By sequentially applying constraints from chiral effective field theory ($\chi$-EFT), the maximum observed neutron star mass (M${}_{\text{max}}$), GW170817, and NICER measurements, we have obtained tight posterior distributions for both nuclear and astrophysical observables, providing a comprehensive picture of the dense nuclear matter.

The sequential application of constraints reveals a clear trend. Within the baseline interval set by $\chi$-EFT and M${}_{\text{max}}$, GW170817 favors a softer EoS, NICER 1 and 2 suggest a slightly stiffer EoS, and NICER 3 and 5 reinforce the preference for a softer EoS, with NICER 5 providing the most stringent constraints.

The RMF model parameters, particularly those governing high-density behavior, are well-constrained. The nuclear saturation parameters exhibit tight posterior distributions, as shown in Table 3. Since the coupling constants are derived from these saturation properties, they are similarly well-constrained (Fig. 2). A significant finding is the growing importance of the $\omega$-$\rho$ coupling ($\Lambda_{\omega\rho}$), which increases from 0.043 to 0.054 as constraints are added, with steadily narrowing uncertainties (Fig. 2). This parameter, which governs the density dependence of the symmetry energy, becomes increasingly critical under astrophysical constraints, particularly NICER 5, highlighting its role in describing neutron-rich matter in neutron star cores.

The crust-core transition density and pressure, as shown in Fig. 6, provide insights into neutron star crust thickness. Constraints favoring a softer EoS, such as GW170817 and NICER 5, lead to lower transition density and pressure, resulting in a thinner crust. In contrast, NICER 1 and 2, which favor a slightly stiffer EoS, support a slightly thicker crust. NICER 3 and 5 further reduce the transition parameters, reinforcing the preference for a thinner crust, consistent with the overall softening of the EoS.

A highlight of this study is the prediction of neutron skin thickness for ${}^{48}\mathrm{Ca}$ and ${}^{208}\mathrm{Pb}$, which probes the symmetry energy and its slope ($L$). For ⁴⁸Ca, the model predicts a neutron skin thickness consistent with the CREX measurement in both cases, whether $\chi$-EFT constraints are included or not, reinforcing the model’s robustness in reproducing this experimental result. For ${}^{208}\mathrm{Pb}$, the predicted skin thickness is significantly lower than the PREX mean of 0.28 fm. When $\chi$-EFT is excluded, the peak shifts to 0.13 fm, with the tail of the distribution falling within the PREX 1$\sigma$ range. This discrepancy in the neutron skin suggests potential tension between nuclear (PREX) and astrophysical (NICER, GW170817) constraints, possibly due to limitations in the RMF model. Furthermore, we did not find any clear correlation between the neutron skin thickness of ²⁰⁸Pb and the symmetry energy slope $L$, unlike previous studies.

The pronounced softness of the EoS, driven by GW170817 and NICER 5, may indicate the presence of exotic degrees of freedom such as hyperons and quarks Nandi and Char (2018); Nandi et al. (2019) in neutron star cores, which we plan to explore in future work within our model. The tension between the ${}^{208}\mathrm{Pb}$ neutron skin thickness and PREX results warrants further investigation. One possibility is to include more terms in the expansion of symmetry energy, which we plan to study next.

In conclusion, this study reestablishes the power of combining nuclear theory with multi-messenger astrophysics to constrain the neutron star EoS and nuclear properties. The tight constraints on model parameters, the softening of the EoS, and the successful reproduction of CREX highlight the robustness of our approach, while the discrepancy with PREX opens new avenues for research. Future gravitational wave observations and more precise NICER data could further tighten these constraints, enhancing our understanding of dense matter physics.