Bayesian Analysis of Neutron-star Properties with Parameterized Equations of State: The Role of the Likelihood Functions

The EOS models we construct are a patchwork of several different components, which we briefly summarize below and that follow a procedure that is similar to the one discussed in Altiparmak et al. (2022), to which we refer the reader for details. At densities below 0.5 n_s , we use the Baym–Pethick–Sutherland (BPS) prescription (Baym et al. 1971). We then describe the EOS in the range n/n_s ∈ [0.5–1.1] with a single polytrope of the form p = Kn^Γ (Rezzolla & Zanotti 2013), where Γ is uniformly sampled in [1.77, 3.23], such that the pressure resides entirely between the soft and stiff model of Hebeler et al. (2013), while the polytropic constant K is fixed by matching to the BPS EOS at 0.5 n_s .

For number densities between 1.1 and ≈40 n_s we use a series of piecewise-linear segments for the sound speed as a function of the chemical potential as starting point to construct thermodynamic quantities. Throughout this work we use N = 11 segments of the following sound speed form:

$$\begin{eqnarray}&&{c}_{s}^{2}(\mu )=\displaystyle \frac{\left({\mu }_{i+1}-\mu \right){c}_{s,i}^{2}+\left(\mu -{\mu }_{i}\right){c}_{s,i+1}^{2}}{{\mu }_{i+1}-{\mu }_{i}},\end{eqnarray} \tag{ 1 }$$

where the values of the chemical potential μ₀ and the sound speed at the first matching point, ${c}_{s,0}^{2}$, are determined by the corresponding polytrope; the values of μ_N = 2.6 GeV, as well as ${c}_{s,N}^{2}$ are determined by the perturbative QCD boundary conditions discussed below. In our previous work (Altiparmak et al. 2022) we randomly sample the remaining values of ${c}_{s,i}^{2}$ and μ_i . But in this work we fix the chemical potentials that are not determined by boundary conditions to be log-equidistantly separated at

$$\begin{eqnarray}&&{\mu }_{i}=1.02{\left(\displaystyle \frac{2.6}{1.02}\right)}^{i/N}\,\mathrm{GeV}\,\,\,\,\mathrm{for}\,\,\,\,i=1,\ldots ,N-1,\end{eqnarray} \tag{ 2 }$$

and only sample the corresponding ${c}_{s,i}^{2}$ uniformly in the range [0, 1] set by thermodynamic stability and causality. We note that within our approach for the construction of the low-density EOS, the maximum chemical potential that can be reached by the polytrope at n = 1.1 n_s when varying the Γ parameter is μ₀ = 1.02 GeV, which is then used as the lower boundary of the log-equidistant chemical potentials. Using a relatively large number of N = 11 segments with log-equidistant chemical potentials allows us to ensure sufficient resolution at densities relevant to the NS interior. Keeping the matching chemical potentials at constant values allows us to reduce the number of free parameters in the Bayesian analysis and make it numerically feasible for the large number of segments we use.

Finally, at chemical potential μ = 2.6 GeV, corresponding to a number density of ≈40 n_s , we impose the parameterized next-to-next-to-leading order (2NLO) perturbative QCD result of Fraga et al. (2014) for the pressure of cold quark matter

$$\begin{eqnarray}&&{p}_{\mathrm{QCD}}(\mu ,X)=\displaystyle \frac{3}{4{\pi }^{2}}\displaystyle \frac{{\mu }^{4}}{{3}^{4}}\left[{c}_{1}-\displaystyle \frac{{d}_{1}{X}^{-{\nu }_{1}}}{(\mu /\mathrm{GeV})-{d}_{2}{X}^{-{\nu }_{2}}}\right],\end{eqnarray} \tag{ 3 }$$

where c₁ = 0.9008, d₁ = 0.5034, d₂ = 1.452, ν₁ = 0.3553, ν₂ = 0.9101, and the effective renormalization scale parameter X is chosen in the range [1, 4].

There has been recent progress in extending Equation (3) to partial 3NLO (Gorda et al. 2021a, 2021b) resulting in a small correction to which our analysis is however not sensitive. Naively one could expect that imposing the constraint (3) at ≲40 n_s has little impact on the EOS at NS densities that are ≲7 n_s . However, a number of recent studies have shown that high-density QCD is indeed constraining, not only for the EOS (Gorda et al. 2022; Komoltsev & Kurkela 2022; Somasundaram et al. 2022), but also for related quantities, such as the sound speed and the conformal anomaly inside massive NSs (Ecker & Rezzolla 2023).

In addition to theoretical constraints, there exists a growing amount of astrophysical data that further constrain the dense matter EOS and in consequence the predictions for NS properties. Arguably the largest impact has direct mass measurements from Shapiro delay of binary systems that have massive radio pulsars like PSR J0740+6620 ($M={2.08}_{-0.07}^{+0.07}\,{M}_{\odot }$) by Fonseca et al. (2021), PSR J0348+0432 ($M={2.01}_{-0.04}^{+0.04}\,{M}_{\odot }$) by Antoniadis et al. (2013), and PSR J1614-2230 ($M={1.908}_{-0.016}^{+0.016}\,{M}_{\odot }$) by Arzoumanian et al. (2018). These measurements set lower limits on the maximum mass M_TOV of static NSs. We will therefore refer to them collectively as TOV constraints from here on.

In addition, the NICER experiment has provided combined mass and radius measurements from the X-ray pulse-profile modeling of PSR J0030+0451 by Miller et al. (2019) ($M={1.44}_{-0.14}^{+0.15}\,{M}_{\odot },\,R={13.02}_{-1.06}^{+1.24}\,\mathrm{km}$) and by Riley et al. (2019) ($M={1.34}_{-0.16}^{+0.15}\,{M}_{\odot },\,R={12.71}_{-1.19}^{+1.14}\,\mathrm{km}$), and of PSR J0740+6620 by Miller et al. (2021) ($M={2.08}_{-0.07}^{+0.07}\,{M}_{\odot },\,R={13.7}_{-1.5}^{+2.6}\,\mathrm{km}$) and by Riley et al. (2021) ($M={2.072}_{-0.066}^{+0.067}\,{M}_{\odot },\,R={12.39}_{-0.98}^{+1.30}\,\mathrm{km}$). These measurements offer bounds for the minimum and maximum radii of typical (M ≈ 1.4 M_⊙) and massive (M ≈ 2 M_⊙) NSs. We refer to these collectively as NICER constraints.

There also exist constraints on the tidal deformability of NSs in binary systems that have been obtained from the direct gravitational-wave (GW) detection by the LIGO/Virgo Collaboration from the merger event GW170817 (Abbott et al. 2017). From this measurement, upper bounds on the tidal deformability Λ_1.4 ≤ 580 (Abbott et al. 2018) of individual NS with M = 1.4 M_⊙, as well as for the binary tidal deformability parameter ${\tilde{{\rm{\Lambda }}}}_{1.186}\leqslant 720$ (Abbott et al. 2019) and for the chirp mass

$$\begin{eqnarray}&&{{ \mathcal M }}_{\mathrm{chirp}}:={({M}_{1}{M}_{2})}^{3/5}{({M}_{1}+{M}_{2})}^{-1/5}={1.186}_{-0.001}^{+0.001}\,{M}_{\odot }\end{eqnarray} \tag{ 4 }$$

have been derived. We then refer to these constraints obtained from GW170817 as the GW constraints.

Finally, there are other constraints that could be used but will not be imposed here, either because they are theoretical constraints and hence with model-dependent uncertainties (see, e.g., Bauswein et al. 2017; Koeppel et al. 2019; Tootle et al. 2021, for some constraints on the minimum radii derived from prompt collapse of the binary to a black hole), or because the measurement uncertainties are excessively large. In particular, there exist predictions for the maximum mass by a number of groups on the basis of the GW170817 event and the corresponding gamma-ray burst event GRB170817A, namely, ${M}_{\mathrm{TOV}}\lesssim {2.16}_{-0.15}^{+0.17}\,{M}_{\odot }$ (Margalit & Metzger 2017; Rezzolla et al. 2018; Ruiz et al. 2018; Shibata et al. 2019). However, such upper bounds on M_TOV also require a number of assumptions about the collapse to a black hole of the merged object and about the ejected mass leading to the kilonova signal of AT 2017gfo (see, e.g., Nathanail et al. 2021).

We also note several recent NS-mass measurements have triggered considerable interest. For example, observations of the black-widow pulsar PSR J0952-0607 have led to an estimated mass of M = 2.35 ± 0.17 M_⊙ (Romani et al. 2022), which exceeds any previous measurements, including that of PSR J2215+5135 ($M={2.27}_{-0.15}^{+0.17}\,{M}_{\odot }$) by Linares et al. (2018). Similarly, the recent discovery of a light central compact object with a mass of $M={0.77}_{-0.17}^{+0.20}\,{M}_{\odot }$ within the supernova remnant HESS J1731-347 (Doroshenko et al. 2022) has been interpreted as the lightest NS known or a strange star; the associated radius has been estimated to be $R={10.4}_{-0.78}^{+0.86}\,\,\mathrm{km}$. In Appendix A we contrast this estimate with the one that can be inferred from our Bayesian analysis.

In the following, we discuss how we implement the TOV, NICER, and GW constraints in our Bayesian analysis, starting with the traditional way that uses variable likelihood functions, followed by the simpler sharp cutoff method that uses constant step functions for the likelihoods. We will refer to the former as the variable-likelihood (VL) analysis and to the latter as the constant-likelihood (CL) analysis. In each case, we employ more than 1.5 × 10⁵ causal and thermodynamically consistent posterior EOSs that satisfy the nuclear theory and perturbative QCD boundary conditions, as well as the observational constraints from pulsars and GW measurements.

Bayesian analysis has been extensively used to infer the probability distribution of unknown parameters by exploiting the knowledge of some measured variables. This is accomplished by constructing likelihood functions that correlate unknown parameters and measured variables. In this approach, the probability of a certain set of parameters is given by the posterior probability

$$\begin{eqnarray}&&P({\boldsymbol{\theta }}| {\boldsymbol{d}})=\displaystyle \frac{L({\boldsymbol{d}}| {\boldsymbol{\theta }})\pi ({\boldsymbol{\theta }})}{Z},\end{eqnarray} \tag{ 5 }$$

where the vector θ collects all the free parameters, d denotes the imposed measurement data, π( θ ) is the so-called prior, and the normalization constant Z is also called the evidence. In our case, the likelihood L can be written as the product of the likelihoods of the three constraints we impose

$$\begin{eqnarray}&&L({\boldsymbol{d}}| {\boldsymbol{\theta }})={{ \mathcal L }}^{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}}){{ \mathcal L }}^{\mathrm{NICER}}({\boldsymbol{\theta }})\,{{ \mathcal L }}^{\mathrm{GW}}({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 6 }$$

where the individual likelihood functions ${{ \mathcal L }}^{\mathrm{TOV}},{{ \mathcal L }}^{\mathrm{NICER}}$, and ${{ \mathcal L }}^{\mathrm{GW}}$ are defined in detail below and the free parameters vectors are given by

$$\begin{eqnarray}&&{\boldsymbol{\theta }}:={{\boldsymbol{\theta }}}_{\mathrm{EOS}}\cup \{{\mu }_{c}^{{\ell }}| \,{\ell }=1,2,3,4\},\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}_{\mathrm{EOS}}:=\{{c}_{{\rm{s}},i}^{2}| \,i=1,2,\ldots ,N-1\}\cup {\rm{\Gamma }}.\end{eqnarray} \tag{ 8 }$$

The four chemical potentials ${\mu }_{c}^{{\ell }}$ refer to the values at the center of the NSs, two of which refer to pulsars with simultaneous mass–radius measurements, while the other two to the two NSs composing the GW170817 event. The three types of parameters in this construction, i.e., ${c}_{{\rm{s}}}^{2},{\rm{\Gamma }}$ and ${\mu }_{c}^{{\ell }}$, are uniformly distributed in the following ranges:

$$\begin{eqnarray}&&{c}_{{\rm{s}},i}^{2}\in [0,1],\,\,\,\,{\rm{\Gamma }}\in [1.77,3.23],\,\,\,\,{\mu }_{c}^{{\ell }}\in [1,2.1]\,\mathrm{GeV}.\end{eqnarray} \tag{ 9 }$$

For the parameter sampling, we use the Pymultinest (Buchner 2016) algorithm implemented in Bilby (Ashton et al. 2019).

For simplicity, to approximate the mass measurements of PSR J0740+6620, PSR J0348+0432, and PSR J1614-2230, we use Gaussian distributions and then the cumulative density function (CDF) of each Gaussian is used to build the likelihood contribution of the corresponding pulsar, namely,

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\rm{i}}}^{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})=\displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(\displaystyle \frac{{M}_{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-{\bar{M}}_{i}}{\sqrt{2}{\sigma }_{i}}\right)\right],\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\mathrm{erf}(x):=\displaystyle \frac{2}{\sqrt{\pi }}{\int }_{0}^{x}{e}^{-{t}^{2}}{dt},\end{eqnarray} \tag{ 11 }$$

is the error function, while ${\bar{M}}_{i}$ and σ_i (i = 1, 2, 3) are the mean and the standard deviation of the ith mass measurement, respectively. The total likelihood for the TOV constraints is then given by the product of the individual likelihoods of these three pulsars:

$$\begin{eqnarray}&&{{ \mathcal L }}^{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}}):=\displaystyle \prod _{i=1}^{3}{{ \mathcal L }}_{{\rm{i}}}^{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}}).\end{eqnarray} \tag{ 12 }$$

The probability distribution function of mass and radius of PSR J0030+0451 ⁴ and PSR J0740+6620 ⁵ are all estimated by the kernel density estimation (KDE; Scott 1992) with a Gaussian kernel using the released posterior (M,R) samples S _k . The likelihood of NICER constraints can then be represented as the product of the KDE of each pulsar, namely,

$$\begin{eqnarray}&&{{ \mathcal L }}^{\mathrm{NICER}}({\boldsymbol{\theta }}):=\displaystyle \prod _{k=1}^{2}{\mathrm{KDE}}_{k}\left(M({{\boldsymbol{\theta }}}_{\mathrm{EOS}},{\mu }_{c}^{k}),R({{\boldsymbol{\theta }}}_{\mathrm{EOS}},{\mu }_{c}^{k})| {{\boldsymbol{S}}}_{k}\right),\end{eqnarray} \tag{ 13 }$$

where ${{\boldsymbol{\theta }}}_{\mathrm{EOS}}\cup \{{\mu }_{c}^{k}| \,k=1,2\}$ are parameters needed to construct this likelihood and they constitute a subset of θ , while the mass M and the radius R of the kth pulsar are functions of the sampled EOS parameters θ _EOS and of the central chemical potential parameters ${\mu }_{c}^{k}$.

When imposing a measurement of a single GW event, the corresponding likelihood can be expressed as (Abbott et al. 2020)

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal L }}^{\mathrm{GW}}({\boldsymbol{d}}| {{\boldsymbol{\theta }}}_{\mathrm{GW}},{\boldsymbol{W}})\\ \,\propto \,\displaystyle \prod _{i}^{{N}_{d}}\exp \left[-2{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{\left|{\tilde{d}}_{i}(f)-{\tilde{h}}_{i}(f;{{\boldsymbol{\theta }}}_{\mathrm{GW}},{\boldsymbol{W}})\right|}^{2}}{{S}_{n}^{i}(f)}\,{df}\right],\end{array}\end{eqnarray} \tag{ 14 }$$

where W represents the waveform model, i and N_d denote the ith and total number of GW detectors that have detected this event, respectively, while Sⁱ _n , ${\tilde{d}}_{i}$, and ${\tilde{h}}_{i}$ represent the power spectral density of the detector noise, the detected strain signal, and the expected strain from the waveform model, respectively.

The vector θ _GW includes parameters that are particularly useful to infer EOS properties, namely,

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}_{\mathrm{GW}}^{\mathrm{EOS}}:=\{{M}_{d}^{1},{M}_{d}^{2},{{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2}\},\end{eqnarray} \tag{ 15 }$$

where M^j _d (j = 1, 2) are the individual masses in the detector frame, while the tidal deformabilities of the two NSs are computed as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{j}:=\displaystyle \frac{2}{3}{k}_{2}{\left(\displaystyle \frac{{R}_{j}}{{M}_{s}^{j}}\right)}^{5}={{\rm{\Lambda }}}_{j}({\mu }_{c}^{j},{{\boldsymbol{\theta }}}_{\mathrm{EOS}}).\end{eqnarray} \tag{ 16 }$$

In the expression above, R_j is the radius of the jth NS in the binary, k₂ is the second tidal Love number, M^j _s is the mass in the source frame, and ${\mu }_{c}^{j}$ denotes the central chemical potential. The relation between the masses in the two frames is simply mediated by the cosmological redshift z, namely,

$$\begin{eqnarray}&&{M}_{d}^{j}=(1+z){M}_{s}^{j}({\mu }_{c}^{j},{{\boldsymbol{\theta }}}_{\mathrm{EOS}}),\end{eqnarray} \tag{ 17 }$$

where z = 0.0099 for the GW170817 event (Abbott et al. 2019). To this scope, we use the interpolated likelihood function for GW170817 that is implemented in Toast (Hernandez Vivanco et al. 2020) and that is marginalized over all the nuisance parameters ${{\boldsymbol{\theta }}}_{\mathrm{GW}}^{\mathrm{nui}}:={{\boldsymbol{\theta }}}_{\mathrm{GW}}\setminus {{\boldsymbol{\theta }}}_{\mathrm{GW}}^{\mathrm{EOS}}$ in Equation (14) to construct the likelihood for the GW constraint. These nuisance parameters consist of a dozen of additional parameters that are useful when performing a Bayesian analysis on GW-emitting binaries but that are not relevant for our EOS inference. Including these parameters would significantly slow down the sampling process and hence we use the likelihood of Toast that marginalizes over these parameters.

In the CL analysis, we use the same set of astrophysical data as in the VL analysis, but the corresponding likelihood functions are constructed differently. More specifically, the likelihood functions are implemented as Heaviside step functions of the form

$$\begin{eqnarray}\begin{array}{rcl}{H}_{\pm } & = & {H}_{\pm }\left({x}_{i}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-{x}_{i}^{\mathrm{cut}}\right)\\ & = & \left\{\begin{array}{ll}1, & \,\mathrm{for}\,\pm \left[{x}_{i}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-{x}_{i}^{\mathrm{cut}}\right]\geqslant 0,\\ 0, & \,\mathrm{else},\end{array}\right.\end{array}\end{eqnarray} \tag{ 18 }$$

where x_i ( θ _EOS) is a quantity constructed from the prior of EOS parameters θ _EOS and is either constrained from above (−) or from below (+) by some sharp cutoff value ${x}_{i}^{\mathrm{cut}}$ obtained from the ith measurement of this quantity. In the following, in order to construct the CL functions, we will use the same sharp cutoff values ${x}_{i}^{\mathrm{cut}}$ for the astrophysical constraints employed by Altiparmak et al. (2022), which we discuss in detail below. On the other hand, for the EOS parameterization, the priors of the EOS parameters and the theoretical constraints from pQCD and CET will be the same as those in the VL analysis discussed in the previous section.

For the TOV constraints, we choose a single lower bound of M_TOV ≥ 2.01 M_⊙, which is motivated by the mass measurements of PSR J0740+6620 and PSR J0348+0432. The likelihood function that implements this constraint is then simply given by

$$\begin{eqnarray}&&{{ \mathcal L }}^{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})={H}_{+}\left({M}_{\mathrm{TOV}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-2.01\,{M}_{\odot }\right).\end{eqnarray} \tag{ 19 }$$

For the NICER constraints, most relevant are the lower bounds on the NS radii since the upper bounds are typically set by the binary tidal deformability for the GW170817 event

$$\begin{eqnarray}&&\tilde{{\rm{\Lambda }}}:=\displaystyle \frac{16}{13}\displaystyle \frac{\left(12{M}_{2}+{M}_{1}\right){M}_{1}^{4}{{\rm{\Lambda }}}_{1}+\left(12{M}_{1}+{M}_{2}\right){M}_{2}^{4}{{\rm{\Lambda }}}_{2}}{{\left({M}_{1}+{M}_{2}\right)}^{5}}\end{eqnarray} \tag{ 20 }$$

as these are effectively more constraining. It is also important to note that measurement of NS masses ≳2.2 M_⊙ make the minimum-radius bounds obtained by NICER essentially ineffective, as recently shown in Ecker & Rezzolla (2023). In practice, we choose a lower bound of R > 10.8 km at M = 1.1 M_⊙ and a lower bound of R > 10.75 km at M = 2.0 M_⊙ to approximate the radius measurements of PSR J0030+0451 and PSR J0740+6620, respectively. The corresponding likelihood function can then be written as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal L }}^{\mathrm{NICER}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}}) & = & {H}_{+}\left({R}_{2.0}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-10.75\,\mathrm{km}\right)\\ & & \times \,{H}_{+}\left({R}_{1.1}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-10.80\,\mathrm{km}\right).\end{array}\end{eqnarray} \tag{ 21 }$$

Finally, as a GW constraint, we impose the upper bound $\tilde{{\rm{\Lambda }}}\leqslant 720$ derived from GW170817 only for a low-spin prior. ⁶ The likelihood can be expressed as

$$\begin{eqnarray}&&{{ \mathcal L }}^{\mathrm{GW}}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})={H}_{-}\left({\tilde{{\rm{\Lambda }}}}_{1.186}({{\boldsymbol{\theta }}}_{\mathrm{EOS}})-720\right),\end{eqnarray} \tag{ 22 }$$

which has to be evaluated at a fixed chirp mass ${{ \mathcal M }}_{\mathrm{chirp}}=1.186\,{M}_{\odot }$ and mass ratios q:=M₁/M₂ ∈ [0.7, 1], as imposed by the GW170817 detection. We note that the results will vary if different values for the maximum mass are considered. These changes have been explored by Ecker & Rezzolla (2023), who have shown that the changes introduced by larger values of M_TOV are mostly quantitative and that no different qualitative behavior emerges. This has been verified by taking into account the mass measurement of PSR J0952-0607 (Romani et al. 2022) for both the VL and the CL approaches. As discussed in more detail in Appendix C, we find that the difference between the results from the CL and VL methods is larger when considering a larger maximum mass, but also that this is simply due to the larger error bar in the mass measurement of PSR J0952-0607, which inevitably will affect differently the VL, where the functional form of the likelihood is fundamental, and the CL approach, which instead is not sensitive to the error bar. Finally, strictly speaking the results presented here will be modified if different cutoff values are chosen for the confidence limits (e.g., changing from 90% to 95%); these changes, however, will be proportional to the difference in the cutoff limits and hence of the order of few percent at most. For this reason, and again because they will not introduce qualitative differences, we do not consider them in our analysis.

In this section, we first analyze the EOS properties that can be inferred from our two Bayesian setups. We first show in the left panel of Figure 1 the PDFs of the sound speed that we computed using the VL approach. The PDFs show a structure very similar to the ones computed by Altiparmak et al. (2022) and Marczenko et al. (2023), which were obtained making use of the CL approach. In particular, it is possible to note that the sound speed rises rapidly until energy densities e ≈ 600 MeV fm⁻³, where it significantly exceeds the conformal limit ${c}_{\mathrm{CFT}}^{2}=1/3$. This is a well-known feature that is necessary to explain maximum-mass measurements M ≈ 2 M_⊙ (see discussion by Bedaque & Steiner 2015; Hoyos et al. 2016; Moustakidis et al. 2017; Kanakis-Pegios et al. 2020; Brandes et al. 2023; Gorda et al. 2022) that require a large sound speed at moderate densities, and is further enhanced if larger maximum-mass constraints are imposed, namely, if M_TOV > 2 M_⊙ (see Ecker & Rezzolla 2023).

Zoom In Zoom Out Reset image size

The decrease of the sound speed after the maximum at higher densities is a consequence of the pQCD constraint imposed at very large energy densities, as well as of the tidal-deformability constraint from GW170817, that disfavors stiff models with too large sound speed. Overall, when comparing the 90% confidence level contours for the VL and CL approaches (red-solid and orange-dashed lines, respectively), it is apparent that not only the qualitative behavior is remarkably similar, but also the quantitative ones. These differences become more pronounced when comparing the 100% confidence level contours (red-dotted and orange-dotted lines, respectively) as it is in the tails of the likelihood functions that the largest differences are imprinted in the two methods. More importantly, the statistical relevance of these stellar models is commensurate to the very low probability with which they appear.

We should note that there are also two artificial features in the PDF that are due to the parameterization method. One is the local minimum in the upper bound of the 90% confidence level at ${c}_{s}^{2}\approx 0.7$ and e ≈ 350 MeV fm⁻³, which is an artifact of fixing the chemical potentials at constant values and which vanishes when sampling the chemical potentials randomly. The second one is the pronounced maximum at large densities, close to where the pQCD boundary conditions are imposed. This maximum is caused by cases in which the pressure at large densities is relatively low, but where it is still possible to reach the allowed interval for p_QCD in a causal way with strongly varying sound speed slightly below μ = 2.6 GeV (see also Altiparmak et al. 2022, where this is discussed in more detail).

Also shown in the left panel of Figure 1, respectively, with solid and dashed black (purple) lines are the 90% confidence levels of the PDFs relative to the centers of NSs with M = 1.4 M_⊙ (or maximally massive, with M = M_TOV) for the VL and the CL approaches, respectively. Also in this case, it is possible to note that the differences between the two evaluations of the likelihood functions are very small. More specifically, we find that in the VL (CL) approach the median value of the sound speed at the center of the NSs is $({c}_{s,c}{)}_{1.4}^{2}=0.60\,(0.63)$ for the typical mass M = 1.4 M_⊙ and that 97% of the EOSs have central sound speed ${c}_{s,c}^{2}\gt 1/3$. Similarly, we find that at the center of maximally massive NSs with M = M_TOV, the sound speed are $({c}_{s,c}{)}_{\mathrm{TOV}}^{2}=0.23\,(0.20)$, in agreement with the results of Ecker & Rezzolla (2023), while only 28 (24)% of the EOSs have central sound speed ${c}_{s,c}^{2}\gt 1/3;$ a detailed description of the results in Figure 1 can be found in Table 1.

Table 1. NS and EOS Properties (90% Confidence Levels) from the VL and the CL Methods

Method	R1.4	R2.0	RTOV	Λ1.4	ΛTOV	nc, 1.4	nc, TOV	pc, 1.4	pc, TOV	${({c}_{s,c}^{2})}_{1.4}$	${({c}_{s,c}^{2})}_{\mathrm{TOV}}$
	(km)	(km)	(km)			(ns )	(ns )	(MeV fm−3)	(MeV fm−3)
VL	${12.2}_{-0.9}^{+0.9}$	${12.4}_{-1.1}^{+1.0}$	${11.8}_{-1.0}^{+1.3}$	${450}_{-190}^{+230}$	${11}_{-5}^{+18}$	${2.4}_{-0.5}^{+0.7}$	${5.5}_{-1.2}^{+1.2}$	${60}_{-20}^{+30}$	${350}_{-170}^{+210}$	${0.60}_{-0.22}^{+0.23}$	${0.23}_{-0.17}^{+0.38}$
CL	${12.4}_{-1.0}^{+0.5}$	${12.6}_{-1.4}^{+0.7}$	${12.1}_{-1.3}^{+1.1}$	${490}_{-220}^{+130}$	${13}_{-7}^{+26}$	${2.3}_{-0.4}^{+0.8}$	${5.2}_{-1.1}^{+1.5}$	${50}_{-10}^{+30}$	${310}_{-150}^{+220}$	${0.63}_{-0.25}^{+0.23}$	${0.20}_{-0.16}^{+0.38}$

Note. Listed are the radii R of NSs with mass M = 1.4, 2.0 M_⊙, with maximum mass M = M_TOV, the corresponding tidal deformabilities Λ_1.4 and Λ_TOV, as well as the number density n_c, the pressure p_c and the sound speed ${c}_{s,c}^{2}$ at the center of typical (M = 1.4 M_⊙) and maximally massive (M = M_TOV) NSs.

Download table as:  ASCIITypeset image

In full similarity, we show in the right panel of Figure 1 the two PDFs relative to the EOSs. Also in this case, the features of the PDF of the VL approach are very similar to that obtained previously via the CL method (Altiparmak et al. 2022). As remarked in several previous studies, a most interesting feature is the clear change of slope (also referred to as a kink) of the PDF at energy densities e ≈ 600 MeV fm⁻³. Annala et al. (2020) have interpreted this feature as a sign of the appearance of a deconfinement phase transition from dense baryonic matter to quark matter (see also Motta et al. 2021 for an alternative interpretation in terms of hyperons). The 90% confidence intervals of the two approaches are again shown by red-solid and orange-dashed lines and they agree remarkably well. Again, deviations can be appreciated when comparing the 100% confidence intervals, marked by red and orange-dotted lines, which exhibit differences at energy densities e ≈ 300 MeV fm⁻³. The origin of these differences is reasonably well understood. In particular, the difference in the upper bound of these contours is due to the diverse treatment of the GW constraint: the VL approach allows for larger variations than the CL method since in the latter the binary tidal deformability is assumed to be strictly $\tilde{{\rm{\Lambda }}}\lt 720$. As a result, the VL approach allows for EOSs that can have radii for stars M ≈ 1.4 M_⊙ that are larger than those in the CL method (this point will be discussed also in the next section).

Finally, the right panel of Figure 1 also reported by black (purple) solid and dashed lines the values of the pressure and energy density reached at the center of typical (maximally massive) NSs with M = 1.4 M_⊙ (M_TOV) for the VL and the CL methods, respectively. Also in this case, the two sets of curves are almost identical. These contours are useful to determine the median values of the energy density and pressure at the center of representative stars. More specifically, within the VL (CL) approach we find e_c = 390 (380)MeV fm⁻³ and p_c = 60 (50)MeV fm⁻³ for typical mass (M = 1.4 M_⊙) NSs, while e_c = 1100 (1030)MeV fm⁻³ and p_c = 350 (310)MeV fm⁻³ inside maximally massive NSs (M = M_TOV) (see also Table 1).

We next move on to the NS properties obtained from the VL and the CL approaches. As in the previous section we show for clarity only the PDFs for VL, but the confidence intervals for both VL and CL. Let us start with a discussion of the mass–radius relation shown in Figure 2. Overall, the PDFs from the VL and CL analyses are very similar and centered around R ≈ 12 km for NSs with masses M ≲ 2.0 M_⊙, while they start to differ for more massive stars, where the VL approach allows for NSs with R ∼ 14 km although with rather small probability. The slightly wider PDF in the VL approach can be explained by the fact that the NICER and the GW constraints are imposed with distributions that have some support also beyond the sharp cutoff values that are employed in the CL approach. This can be seen more clearly by comparing the outer contours (100% intervals) shown as red (VL) and orange (CL) dotted lines. As a result, the VL PDF yields upper and lower bounds for the NS radii that are less restrictive than those of the CL approach. On the other hand, when comparing the 90% confidence intervals of the two approaches (solid and dashed lines) it is possible to conclude that they are essentially identical, thus underlying that the methodological differences in the use of the likelihood functions two approaches do not result in significant changes in the PDFs. Furthermore, both 90% confidence levels are in good agreement with the analytic lower bound $R/\mathrm{km}\gtrsim 8.91+2.66{(M/{M}_{\odot })-0.88(M/{M}_{\odot })}^{2}$ (black-dashed line) derived when using the detection of GW170817 and the estimates on the threshold mass to prompt collapse (Koeppel et al. 2019; Tootle et al. 2021).

Zoom In Zoom Out Reset image size

Also shown in the right panel of Figure 2 are the one-dimensional PDF slices obtained in the two approaches for NSs with fixed masses M = 1.4 and 2.0 M_⊙, together with the corresponding median estimates (vertical solid lines) for their radii. These cuts show more explicitly that the VL approach results in distributions that are wider and less skewed than those for the CL case, but also that the resulting median values differ only by ≈200 m (see Table 1 for details); overall, the differences between the two approaches are clearly much smaller than the current observational uncertainties.

Next, we discuss the PDFs for the tidal deformabilities of isolated and binary NSs that are plotted in Figure 3. More specifically, in the left panel, we show how the PDFs of the tidal deformability Λ of isolated NSs depend on the gravitational mass M in these two Bayesian analyses. There is a simple explanation for the overall trend of the PDF: for generic EOSs that are viable, the relative change in the NS radius is ≲6% in the relevant mass range of M ≈ 1–2 M_⊙, so that the tidal deformability is (see Equation (16)) Λ ∼ k₂ M⁻⁵ ∼ M^−p , with p ≳ 5. Comparing the confidence intervals of the VL and CL approaches leads to conclusions that are very similar to those drawn for the mass–radius relation in Figure 2. More specifically, while the 100% confidence interval of the VL approach is significantly wider than that of the CL method, the 90% contours are again almost identical, so that the differences between the two methods affect mostly the less probable and therefore less important regions of the PDFs.

Zoom In Zoom Out Reset image size

In the right panel of Figure 3, we show instead the dependence of the binary tidal deformability parameter $\tilde{{\rm{\Lambda }}}$ on the chirp mass ${{ \mathcal M }}_{\mathrm{chirp}}$. The PDF is based on more than 9.0 × 10⁶ binary models with uniformly distributed mass ratios q ∈ [0.4, 1.0], i.e., a range that is reasonable according to the mass ratios in binary NS systems detected so far. Analogous results from the CL approach have been shown by Altiparmak et al. (2022) and were later generalized to large-mass constraints by Ecker & Rezzolla (2023) (see also Zhao & Lattimer 2018, for similar studies but where no PDF is computed). Note that, in contrast to previous plots, the red-solid and orange-dashed lines show here the 99% confidence levels of the PDFs. As discussed by Ecker & Rezzolla (2022), these bounds—and in particular the lower one—are of great importance to constrain the tidal deformability (and hence the EOS) using accurate measurements of the chirp mass of future GW detections of binary NS mergers. ⁷ Interestingly, the lower bounds of these intervals are again almost identical, but the upper bounds differ. The reason for this is that the upper bound is directly set by the GW constraint of GW170817, while the lower bound mostly depends on the TOV constraint (Ecker & Rezzolla 2023).

In particular, the sharp cutoff on $\tilde{{\rm{\Lambda }}}$ at ${{ \mathcal M }}_{\mathrm{chirp}}=1.186\,{M}_{\odot }$ employed in the CL approach strongly skews the distribution toward the larger values of $\tilde{{\rm{\Lambda }}}$. This is clearly visible in the bottom part of the right panel in Figure 3, where we report cuts of the binary tidal deformability at ${{ \mathcal M }}_{\mathrm{chirp}}=1.186{M}_{\odot }$ for the two approaches; clearly, the CL cut has a rapid falloff for $\tilde{{\rm{\Lambda }}}\simeq 750$, while the VL cut is almost Gaussian with a longer tail up to $\tilde{{\rm{\Lambda }}}\simeq 1000$. As a result, the CL approach slightly underestimates the upper bound on $\tilde{{\rm{\Lambda }}}$, while the VL approach gives a more conservative estimate. Following Altiparmak et al. (2022), we have calculated an analytic fitting function that approximates these bounds and is given by ⁸

$$\begin{eqnarray}&&{\tilde{{\rm{\Lambda }}}}_{\min (\max )}=a+b{{ \mathcal M }}_{\mathrm{chirp}}^{c},\end{eqnarray} \tag{ 23 }$$

where a = −16, b = 560, and c = −5.1 for the lower bound in both VL and CL analyses, while a = −19 (0.4), b = 2200(1900), and c = −5.1 ( −5.5) for the upper bound in the VL (CL) analysis.

Notwithstanding these small differences, the median values (thin vertical lines) for the case of the GW170817 event, whose chirp mass was ${{ \mathcal M }}_{\mathrm{chirp}}=1.186{M}_{\odot }$, are quite similar: ${\tilde{{\rm{\Lambda }}}}_{1.186}={480}_{-190}^{+260}$ (VL) versus ${510}_{-210}^{+180}$ (CL). Furthermore, the $\tilde{{\rm{\Lambda }}}({{ \mathcal M }}_{\mathrm{chirp}})$ relation can also be used the other way around. When in the future the EOS will be better constrained, it will be then possible to constrain the source-frame chirp mass. Our results can be combined with the detector frame chirp mass from the waveform-matched filter to find the redshift of the luminosity distance and thus constrain cosmological models (see, e.g., Messenger et al. 2014).

We next switch our attention to the PDFs of the sound speed ${c}_{s}^{2}(r/R)$ in the two approaches. This is presented in the left and right panels of Figure 4, where we report the spatial dependence (i.e., as a function of the normalized radial coordinate r/R) of the PDFs of the sound speed inside typical (M = 1.4 M_⊙) and maximally massive (M = M_TOV) NSs.

Zoom In Zoom Out Reset image size

The overall behavior of the PDF for the VL approach shown here is similar to the one from the CL method presented by Ecker & Rezzolla (2022, 2023). In particular, it is interesting to note that in typical NSs the sound speed rises monotonically from the surface toward the center, where it reaches values ${({c}_{s,c}^{2})}_{1.4}\approx 0.6$ that are significantly larger than the conformal limit. On the other hand, in maximally massive NSs the sound speed is nonmonotonic and has a local maximum with ${c}_{s,\max }^{2}\gt 1/3$ in the outer layers (i.e., r/R ∼ 0.6) and a local minimum at the center with ${({c}_{s,c}^{2})}_{\mathrm{TOV}}\lesssim 1/3$ (see Table 1 for details). The confidence intervals and median values from VL and CL approaches are remarkably similar, even quantitatively. More specifically, for NSs with mass M = 1.4 M_⊙ (left panel) the VL prediction for the central sound speed, ${({c}_{s,c}^{2})}_{1.4}={0.60}_{-0.22}^{+0.23}$, is slightly lower than the one by CL, ${({c}_{s,c}^{2})}_{1.4}={0.63}_{-0.25}^{+0.23}$. Note that also in this case, the major differences are in the upper 90% confidence levels, where the sharp cutoff on $\tilde{{\rm{\Lambda }}}$ employed in the CL analysis indirectly causes a slight overestimate of the maximum possible values of the sound speed. Instead, for maximally massive NSs with M = M_TOV (right panel), an opposite behavior is observed: the central sound speed ${({c}_{s,c}^{2})}_{\mathrm{TOV}}={0.23}_{-0.17}^{+0.38}$, from the VL approach is slightly larger than the one from the CL method, ${({c}_{s,c}^{2})}_{\mathrm{TOV}}={0.20}_{-0.16}^{+0.38}$. This is due to the pQCD constraint, which introduces an anticorrelation between the sound speed parameters at low and high densities, namely, the large values of c_s ² in the low-density regions need to be compensated by lower values in the high-density region.

Finally, we show in Figure 5 the correlation between the radius R_TOV and the central number density n_{c, TOV} for an NS at the maximum mass M_TOV. Here again, red-solid and orange-dashed ellipses enclose the 90% confidence levels obtained from the VL and the CL approaches, respectively. These intervals are essentially identical; hence, the relation between R_TOV and n_{c, TOV} is largely independent of the choice of the likelihood functions. In addition, we show results for a number of microphysical multipurpose EOS models of pure nuclear matter (nucleons), with hyperons (hyperons), with quark matter (quark matter) and models derived from holographic QCD. All EOSs are taken from the public database CompOSE ⁹ (Typel et al. 2015) and correspond to beta-equilibrium slices at the lowest available temperature with electrons. Some of the corresponding properties are listed in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Properties of the EOS Models Used in Figure 5

Type	EOS	R1.4	R2.0	RTOV	MTOV	${\tilde{{\rm{\Lambda }}}}_{1.186}^{{\rm{q}}=1\,(0.7)}$	Ref.
		(km)	(km)	(km)	(M⊙)
Nucleons	DD2	13.76	13.45	12.04	2.42	816 (784)	Hempel & Schaffner-Bielich (2010)
	LS220	12.66	11.34	10.62	2.04	574 (569)	Schneider et al. (2017)
	SFHx	12.4	11.78	10.9	2.13	464 (451)	Steiner et al. (2013)
	SLy4	11.73	10.67	10.0	2.05	313 (309)	Schneider et al. (2017)
	KOST2	11.58	11.21	10.2	2.22	358 (352)	Togashi et al. (2017)
Hyperons	DD2 Λϕ	13.75	12.84	11.77	2.1	813 (772)	Banik et al. (2014)
	DD2Y	13.75	12.33	11.53	2.04	814 (767)	Marques et al. (2017)
	DD2YΔ 1.1–1.1	13.47	12.22	10.63	2.17	698 (687)	Raduta et al. (2020)
Quark matter	RDF1.6	12.46	10.31	10.04	2.01	509 (501)	Bastian (2021)
	RDF1.7	12.46	11.44	10.76	2.11	508 (506)	Bastian (2021)
	TM1B145	12.92	12.25	10.57	2.27	625 (612)	Sagert et al. (2010)
	V-QCD soft	12.42	12.01	11.91	2.02	537 (517)	Demircik et al. (2022)
	V-QCD interm	12.51	12.4	11.86	2.14	565 (543)	Demircik et al. (2022)
	V-QCD stiff	12.65	12.85	11.89	2.34	617 (591)	Demircik et al. (2022)

Note. The properties listed are the corresponding NS radii R_1.4, R_2.0 and R_TOV at M = 1.4, 2.0 M_⊙ and the maximum mass M_TOV, respectively, as well as the binary tidal deformability parameter ${\tilde{{\rm{\Lambda }}}}_{1.186}$ of GW170817-like binaries with chirp mass ${{ \mathcal M }}_{c}=1.186\,{M}_{\odot }$ for two different mass ratios q = 1 and 0.7.

Download table as:  ASCIITypeset image

As a concluding remark we note that the correlation between R_TOV and a normalized number density n_{c, TOV}/n_s is described by a second-order polynomial:

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{c,\mathrm{TOV}}}{{n}_{s}}={d}_{0}\,\left[1-\left(\displaystyle \frac{{R}_{\mathrm{TOV}}}{10\,\mathrm{km}}\right)\right]+{d}_{1}\,{\left(\displaystyle \frac{{R}_{\mathrm{TOV}}}{10\,\mathrm{km}}\right)}^{2},\end{eqnarray} \tag{ 24 }$$

where the coefficients are given by d₀ = 27.6 and d₁ = 7.5. The relative residuals are well described by a Gaussian distribution centered on zero and have a standard deviation of 3.7%. The importance of Equation (24) is that it allows estimating the maximum central density from future radius measurements of very massive NSs. We have found that a relation similar to Equation (24) holds also when considering the radii of NSs with M = 1.4 M_⊙. In this case, the coefficients are given by d₀ = 26.6 and d₁ = 7.6, but the scatter is 4 times larger, reaching deviations of ≃30%; although larger than in the case of maximally massive stars, these deviations are much smaller than the uncertainty in n_c,TOV/n_s ≃ 3.7–7.0.