††thanks: Corresponding author

# Neutron Star vs Quark Star in the Multimessenger Era

Zheng Cao  School of Physics and Astronomy, Shanghai Key Laboratory for
Particle Physics and Cosmology, and Key Laboratory for Particle Astrophysics
and Cosmology (MOE), Shanghai Jiao Tong University, Shanghai 200240, China
Lie-Wen Chen<EMAIL_ADDRESS>School of Physics and Astronomy, Shanghai
Key Laboratory for Particle Physics and Cosmology, and Key Laboratory for
Particle Astrophysics and Cosmology (MOE), Shanghai Jiao Tong University,
Shanghai 200240, China

###### Abstract

Neutron stars (NSs) which could contain exotic degrees of freedom in the core
and the self-bound quark stars (QSs) made purely of absolutely stable
deconfined quark matter are still two main candidates for the compact objects
observed in pulsars and gravitational wave (GW) events in binary star mergers.
We perform a Bayesian model-agnostic inference of the properties of NSs and
QSs by combining multi-messenger data of GW170817, GW190425, PSR J0030+0451,
PSR J0740+6620, PSR J1614-2230, PSR J0348+0432 as well as ab initio
calculations from perturbative quantum chromodynamics and chiral effective
field theory. We find the NS scenario is strongly favored against the QS
scenario with a Bayes factor of NS over QS
$\mathcal{B}^{\text{NS}}_{\\_}\text{QS}=11.5$. In addition, the peak of the
squared sound velocity $c_{\\_}s^{2}\sim 0.5c^{2}$ around $3.5$ times nuclear
saturation density $n_{\\_}0$ observed in the NS case disappears in the QS
case which suggests that the $c_{\\_}s^{2}$ first increases and then saturates
at $c_{\\_}s^{2}\sim 0.5c^{2}$ above $\sim 4n_{\\_}0$. The sound velocity and
trace anomaly are found to approach the conformal limit in the core of heavy
NSs with mass $M\gtrsim 2M_{\\_}{\odot}$, but not in the core of QSs.

Introduction.— Understanding the nature of compact stars (CSs) observed in
pulsars and gravitational wave (GW) events in binary star mergers is one of
fundamental questions in contemporary nuclear physics, astrophysics and
cosmology. The baryon number density in the core of CSs can reach several
times nuclear saturation density ($n_{\\_}0=0.16\text{ fm}^{-3}$), making the
CSs be ideal laboratories to study the properties of dense nuclear matter and
QCD phase diagram at extreme high densities and low temperatures Lattimer and
Prakash (2004); Weber (2005); Alford  _et al._ (2008); Fukushima and Hatsuda
(2011); Baym  _et al._ (2018); Blaschke and Chamel (2018); Burgio  _et al._
(2021); Li  _et al._ (2021a); Sedrakian  _et al._ (2023), which is
unaccessible in terrestrial labs. Theoretically, it is still a big challenge
to determine the properties of dense nuclear matter at several times of
$n_{\\_}0$ from ab initio QCD calculations due to the complicated
nonperturbative feature of QCD Brambilla  _et al._ (2014), and thus the
composition inside the CSs is largely unknown. As discussed in Ref. Weber
(2005), the CSs could be neutron stars (NSs) for which besides the
conventional neutrons and protons, some exotic degrees of freedom such as
hyperons, meson condensates and even quark matter may appear in the core. A
popular alternative for CSs is the self-bound quark stars (QSs) made purely of
absolutely stable deconfined quark matter (QM) composed of $u$, $d$ and $s$
Bodmer (1971); Witten (1984); Terazawa (1989); Farhi and Jaffe (1984) or $u$
and $d$ Holdom  _et al._ (2018); Zhao  _et al._ (2019); Cao  _et al._ (2022)
quarks with some leptons. The CSs could be even self-bound strangeon stars in
a solid state comprised of strangeons (quark-clusters with three-light-flavor
symmetry) Xu (2003, 2018).

Thanks to the fast development in astrophysical observation facilities,
significant progress has been made in the last decades for the measurement of
CSs. For example, the mass of several heavy pulsars with mass $M\sim
2M_{\\_}{\odot}$ was measured precisely by Shapiro delay Fonseca  _et al._
(2021). The mass and radius of PSR J0030+0451 with $M\sim 1.4M_{\\_}{\odot}$
and PSR J0740+6620 with $M\sim 2M_{\\_}{\odot}$ were determined simultaneously
by NICER via pulse-profile modeling Miller  _et al._ (2019, 2021); Riley  _et
al._ (2019, 2021). Especially, in recent years, two gravitational wave (GW)
events GW170817 Abbott  _et al._ (2017, 2019) and GW190425 Abbott  _et al._
(2020) from binary star mergers were reported by the LIGO Scientific and Virgo
Collaborations (LVC), which inaugurates a new era of multimessenger astronomy.

Theoretically, ab initio calculations of dense matter have also made
significant progress in recent years. At the low density limit, chiral
effective field theory (ChEFT) Holt  _et al._ (2016); Drischler  _et al._
(2021), which is the low-energy realization of QCD, provides a satisfactory
constraint on the equation of state (EOS) of the NS matter up to densities
$n\approx 1\sim 2n_{\\_}0$ with controllable uncertainties Hebeler  _et al._
(2013); Wellenhofer  _et al._ (2015). At asymptotically high densities with
baryon chemical potential of multi GeV, perturbative QCD (pQCD) computations
become feasible Gorda  _et al._ (2021) and provide potential constraints on
the EOS of dense matter at intermediate densities inside CSs by combining the
results at low densities Annala  _et al._ (2018, 2020, 2022); Gorda  _et al._
(2023). Based on these state-of-the-art ab initio calculations together with
the multi-messenger data, it is extremely interesting to perform a comparative
study on NSs and QSs, which may provide valuable information on the nature of
CSs and the properties of dense matter.

We perform here a Bayesian model-agnostic inference of the properties of NSs
and QSs by combining the data on the mass of heavy pulsars with $M\sim
2M_{\\_}{\odot}$ determined by Shapiro delay, the mass and radius of PSR
J0030+0451 and PSR J0740+6620 from NICER, the tidal deformabilities of CSs
from GW170817 and GW190425 together with ab initio calculations from pQCD and
ChEFT. We find the current multi-messenger data and constraints from pQCD and
ChEFT strongly favor the NS scenario against the QS scenario. In addition, our
analyses on the sound velocity and trace anomaly suggest that the conformal
limit is violated inside QSs, but reached in the core of heavy NSs with
$M\gtrsim 2M_{\\_}{\odot}$.

Model-agnostic EOS— To construct model-agnostic EOSs for NS matter, we adopt
the EOS derived from ChEFT at low densities and extrapolate it to high density
to match the pQCD constraints by speed of sound extension approach Tews  _et
al._ (2018). In particular, following Ref. Annala  _et al._ (2020), in the
density region of $n\in[0.58n_{\\_}0,1.1n_{\\_}0]$, we choose “Soft”,
“Intermediate” and “Stiff” EOSs of Ref. Hebeler  _et al._ (2013), and match
them to BPS EOS Baym  _et al._ (1971) below $0.58n_{\\_}0$. The speed of sound
extension Tews  _et al._ (2018) is then utilized to obtain EOS of NS matter
from $1.1n_{\\_}0$ to $12n_{\\_}0$, in which we uniformly sample a sequence of
stitching points $\\{(n_{\\_}i,c_{\\_}{s,i}^{2})\\}_{\\_}{i=1}^{N}$
($n_{\\_}j>n_{\\_}k$ for $j>k$). These matching points are then connected
using piecewise-linear function to obtain $c_{\\_}{\mathrm{s}}^{2}(n)$ as

$c_{\\_}{\mathrm{s}}^{2}(n)=\frac{\left(n_{\\_}{i+1}-n\right)c_{\\_}{\mathrm{s},i}^{2}+\left(n-n_{\\_}i\right)c_{\\_}{\mathrm{s},i+1}^{2}}{n_{\\_}{i+1}-n_{\\_}i},$
(1)

and the EOS of NS matter from $1.1n_{\\_}0$ to $12n_{\\_}0$ can then be
obtained using the fundamental thermodynamic relation (see, e.g., Ref. Gorda
_et al._ (2023))

$\displaystyle\mu(n)$
$\displaystyle=\mu_{\\_}1\exp\left[\int_{\\_}{n_{\\_}1}^{n}dn^{\prime}\frac{c_{\\_}s^{2}\left(n^{\prime}\right)}{n^{\prime}}\right],$
(2) $\displaystyle\varepsilon(n)$
$\displaystyle=\varepsilon_{\\_}1+\int_{\\_}{n_{\\_}1}^{n}dn^{\prime}\mu\left(n^{\prime}\right),$
(3) $\displaystyle p(n)$ $\displaystyle=-\varepsilon(n)+\mu(n)n.$ (4)

We set $n_{\\_}1=1.1n_{\\_}0$ and $n_{\\_}N=12n_{\\_}0$, and
$c_{\\_}{s,1}^{2}$ (also the corresponding chemical potential $\mu_{\\_}1$ and
energy density $\varepsilon_{\\_}1$) is fixed at the corresponding value from
ChEFT, $n_{\\_}i$ ($i=2$, $\cdots$, $N-1$) and $c_{\\_}{s,i}^{2}$ ($i=2$,
$\cdots$, $N$) is uniformly sampled in $[1.1n_{\\_}0,12n_{\\_}0]$ and [0, 1]
respectively. In the following we use $N=6$ and we note that the results just
change slightly when $N$ vary from 5 to 10. The nested sampler pymultinest
Buchner  _et al._ (2014) (which is installed in bilby Ashton  _et al._ (2019))
is then used to sample over parameter
$\\{(n_{\\_}i,c_{\\_}{s,i}^{2})\\}_{\\_}{i=1}^{N}$ and generate posterior
distribution.

For QS matter (also for strangeon matter), its EOS is unknown even at low
densities, and so we only consider its basic self-bound property with minimum
assumption that the pressure becomes to zero at finite baryon number density
$n_{\\_}1$ corresponding to the QS surface. Assuming $n_{\\_}1>n_{\\_}0$ and
$n_{\\_}N=12n_{\\_}0$, we uniformly sample $n_{\\_}i$ ($i=1$, $\cdots$, $N-1$)
($n_{\\_}j>n_{\\_}k$ for $j>k$) and $c_{\\_}{s,i}^{2}$ ($i=1$, $\cdots$, $N$)
in $[n_{\\_}0,12n_{\\_}0]$ and [0, 1], respectively. At the same time, we also
uniformly sample the chemical potential $\mu_{\\_}1$ at $n_{\\_}1$ in
$[\mu_{\\_}{1,\rm min},930]$ MeV with $\mu_{\\_}{1,\rm min}=500$ MeV,
considering the fact that $\mu_{\\_}1$ should be less than the binding energy
per baryon of the observed stable nuclei (i.e., 930 MeV) to satisfy the
absolutely stable condition Bodmer (1971); Witten (1984); Terazawa (1989);
Farhi and Jaffe (1984). The energy density at $n_{\\_}1$ can then be obtained
as $\varepsilon_{\\_}1=\mu_{\\_}1n_{\\_}1$. The full EOS of QS matter can be
then obtained similarly as in the case of NS.

Bayesian analysis— We use Bayesian hierarchical model to combine constraints
from multiple observations with uncertaintes and then make parameter estimate.
According to Bayes’ theorem, as discussed in Ref Hernandez Vivanco  _et al._
(2020); Thrane and Talbot (2019), for the given data set $\vec{d}$ and
hypothesis $\mathcal{H}$, the posterior distribution of EOS parameters
$\theta$ can be written as

$p(\theta|\vec{d},\mathcal{H})=\frac{\prod_{\\_}i\mathcal{L}(d_{\\_}i|\theta,\mathcal{H})\pi(\theta|\mathcal{H})}{\mathcal{Z}_{\\_}{\mathcal{H}}(\vec{d})},$
(5)

where $i$ runs over individual constraints and each constraint is independent
of each other, $\pi(\theta|\mathcal{H})$ is the hyper-prior distribution for
$\theta$ and here is chosen as uniform distribution,
$\mathcal{L}(d_{\\_}i|\theta,\mathcal{H})$ is the likelihood of the EOS
parameters under the assumption of $\mathcal{H}$ for data $d_{\\_}i$, and the
$\mathcal{Z}_{\\_}{\mathcal{H}}(\vec{d})\equiv\int\prod_{\\_}i\mathcal{L}(d_{\\_}i|\theta,\mathcal{H})\pi(\theta|\mathcal{H})\mathrm{d}\theta$
is a normalization factor called evidence which quantifies how much the
hypothesis is preferred by the data. Based on the data set $\vec{d}$, the
Bayes factor $\mathcal{B}^{\text{NS}}_{\\_}\text{QS}$ for CSs as NSs against
QSs can be obtained as Jeffreys (1961)

$\mathcal{B}^{\text{NS}}_{\\_}\text{QS}=\mathcal{Z}_{\\_}\text{NS}(\vec{d})/\mathcal{Z}_{\\_}\text{QS}(\vec{d}).$
(6)

In the present Bayesian analyses, we combine the Shapiro delay mass
measurements of heavy CSs, the NICER mass-radius of PSR J0030+0451 and PSR
J0740+6620 analyzed by Miller et al. Miller  _et al._ (2019, 2021) (Similar to
Ref. Huth  _et al._ (2022)), the tidal information of GW170817 Abbott  _et
al._ (2017, 2019) and GW190425 Abbott  _et al._ (2020) together with the pQCD
constraints at ultra high densities Gorda  _et al._ (2021); Komoltsev and
Kurkela (2022); Gorda  _et al._ (2023); LHQ  (and ChEFT results for NSs) as
our default data set $\vec{d}_{\\_}\text{def}$.

For the mass measurements, the prior distribution of CS mass for a given EOS
parameter $\theta$ can be written as Landry  _et al._ (2020); Biswas (2022)
$\pi(m|\theta)=\frac{\mathbf{1}_{\\_}{[M_{\\_}\text{low},M_{\\_}{\text{TOV}}(\theta)]}}{M_{\\_}{\text{TOV}}(\theta)-M_{\\_}\text{low}}$,
where $M_{\\_}\text{low}=0.1M_{\\_}\odot$ is the assumed lower bound of the
mass of CSs and $M_{\\_}{\text{TOV}}(\theta)$ is the maximum mass of static
CSs determined by the EOS. For a given mass measurement data $d_{\\_}M$, the
likelihood is
$\mathcal{L}(d_{\\_}M|\theta)=\int\mathrm{d}m\mathcal{L}(d_{\\_}M|m)\pi(m|\theta)$.
We consider here the precise Shapiro delay mass measurement of PSR J1614-2230
Demorest  _et al._ (2010); Fonseca  _et al._ (2016); Arzoumanian  _et al._
(2018), PSR J0348+0432 Antoniadis  _et al._ (2013), PSR J0740+6620 Fonseca
_et al._ (2021) and use Gaussian function $\mathcal{N}(1.908,0.016^{2})$,
$\mathcal{N}(2.01,0.04^{2})$, $\mathcal{N}(2.08,0.07^{2})$ to approximate
these mass measurements $\mathcal{L}(d_{\\_}M|m)$, respectively.

Similarly, for the NICER mass-radius measurement data $d_{\\_}\text{R}$ of PSR
J0030+0451 or PSR J0740+6620, the likelihood can be written as
$\mathcal{L}(d_{\\_}R|\theta)=\int\mathrm{d}m\mathcal{L}[d_{\\_}R|m,R(m,\theta)]\pi(m|\theta)$
and we use Kernel Density Estimation (KDE) to approximate the posterior mass-
radius distribution Mil ; web . To avoid double counting, for PSR J0740+6620
Fonseca  _et al._ (2021), we do not use its Shapiro delay mass data if we
include its NICER mass-radius data. For PSR J0437-4715, we use 2D Gaussian
distribution $\mathcal{N}([13.6,1.44]^{\top},\text{diag}(0.85^{2},0.07^{2}))$
to approximate its mass-radius joint distribution.

For the measurements of GW events, the likelihood can be written as
$\mathcal{L}(d_{\\_}{\text{GW}}|\theta)=\int\mathrm{d}\omega\mathcal{L}(d_{\\_}{\text{GW}}|\omega)\pi(\omega|\theta)$
Hernandez Vivanco  _et al._ (2019), where
$\omega=\\{\mathcal{M}_{\\_}c,q,\Lambda_{\\_}1,\Lambda_{\\_}2\\}$ and
$\mathcal{L}(d_{\\_}\text{GW}|\omega)$ is the nuisance-marginalized likelihood
Hernandez Vivanco  _et al._ (2020) which has marginalized over extrinsic
parameters of the source. With the convention $m_{\\_}1\geq m_{\\_}2$, the
tidal parameter $\Lambda_{\\_}i$ of each compact object could be uniquely
determined by the chirp mass $\mathcal{M}_{\\_}c$, mass ratio $q$, and EOS
parameter $\theta$, i.e. $\Lambda_{\\_}i(\mathcal{M}_{\\_}c,q,\theta)$.

Based on thermodynamic stability and causality, the results from pQCD can be
utilized to constrain the EOS at intermediate density region by fully taking
advantage of thermodynamic potentials Komoltsev and Kurkela (2022). At a high
chemcial potential $\mu_{\\_}H=2.6$ GeV (i.e. $n\approx 40n_{\\_}0$), the
uncertainties of thermodynamics quantity could be parameterized by a
dimensionless parameter $X$ and could be expressed by a set
$\vec{\beta}_{\\_}{\mathrm{pQCD}}(X)=\left\\{p_{\\_}{\mathrm{pQCD}}\left(\mu_{\\_}H,X\right),n_{\\_}{\mathrm{pQCD}}\left(\mu_{\\_}H,X\right),\mu_{\\_}H\right\\}$.
By integrating the uncertainties, one can obtain the corresponding likelihood
$\mathcal{L}(\mathrm{pQCD}|\theta)=\int
d\vec{\beta}_{\\_}HP(\vec{\beta}_{\\_}H)\mathbf{1}_{\\_}{\left[\Delta
p_{\\_}{\min},\Delta p_{\\_}{\max}\right]}(\Delta p)$ Gorda  _et al._ (2023),
with $\Delta p=p_{\\_}\text{pQCD}-p_{\\_}L$ and $p_{\\_}L$ is the pressure of
the last point (i.e., $n_{\\_}N=12n_{\\_}0$ here) of interpolated EOS. It
should be noted that the $p_{\\_}L$ value depends on the EOS at the low
density $n_{\\_}1$ as well as the sequence of stitching points
$\\{(n_{\\_}i,c_{\\_}{s,i}^{2})\\}_{\\_}{i=1}^{N}$ in the speed of sound
extension.

Results and discussions— Using the default data set $\vec{d}_{\\_}\text{def}$,
we perform Bayesian model-agnostic inference of the properties of NSs and QSs.
Firstly, for $M_{\\_}{\rm TOV}$, our present analyses indicate that it is
$M_{\\_}{\rm TOV,NS}=2.17_{\\_}{-0.15}^{+0.26}M_{\\_}\odot$ for NSs and
$M_{\\_}{\rm TOV,QS}=2.49_{\\_}{-0.35}^{+0.47}M_{\\_}\odot$ for QSs in 90%
credible interval (CI), indicating the QSs would have a significantly larger
$M_{\\_}{\rm TOV}$ than the NSs. The $M_{\\_}{\rm
TOV,NS}=2.17_{\\_}{-0.15}^{+0.26}M_{\\_}\odot$ is in nice agreement with the
value of $2.18_{\\_}{-0.13}^{+0.27}M_{\\_}\odot$ estimated recently by taking
advantage of the various structures sampling by a single-layer feed-forward
neural network model embedded in the Bayesian nonparametric inference Han  _et
al._ (2023), implying our present result is independent of the detailed
realization of the model-agnostic EOSs.

Secondly, for the radius $R_{\\_}\text{1.4}$ of CSs with canonical mass of
$1.4M_{\\_}\odot$, its value is estimated to be $R_{\\_}\text{1.4,
NS}=12.44_{\\_}{-0.71}^{+0.74}$ km for NSs and $R_{\\_}\text{1.4,
QS}=11.41_{\\_}{-0.61}^{+0.64}$ km for QSs in 90% CI, suggesting NSs have a
larger $R_{\\_}\text{1.4}$ than QSs. $R_{\\_}\text{1.4,
NS}=12.44_{\\_}{-0.71}^{+0.74}$ km well agrees with
$R_{\\_}{1.4}=12.42_{\\_}{-0.99}^{+0.52}\mathrm{~{}km}$ (95% CI) reported in
Ref. Altiparmak  _et al._ (2022) where an ensemble of EOSs is generated in
advance and weight each EOS according to the likelihood. Our result on
$R_{\\_}\text{1.4, NS}$ is also consistent with
$R_{\\_}{1.4}=11.98_{\\_}{-0.40}^{+0.35}$ (90% CI) Pang  _et al._ (2022)
obtained recently by combining all the EOS-sensitive observations, including
data of the kilonovae and the GRB afterglow. In addition, though different
methods are adopted to construct EOS, our result on $R_{\\_}\text{1.4, NS}$ is
in agreement with previous work Huth  _et al._ (2022); Pang  _et al._ (2021);
Legred  _et al._ (2021); Biswas  _et al._ (2021); Li  _et al._ (2021b); Han
_et al._ (2021); Dietrich  _et al._ (2020); Jiang  _et al._ (2020); Radice and
Dai (2019) within the uncertainty. For the case of QSs, our result also agrees
with $R_{\\_}{1.4}=11.50_{\\_}{-0.55}^{+0.52}$ km obtained by Miao et al. Miao
_et al._ (2021) within the MIT bag model.

Thirdly, the tidal deformability $\Lambda_{\\_}\text{1.4}$ of a
$1.4M_{\\_}\odot$ CS is estimated to be
$\Lambda_{\\_}\text{1.4,NS}=504_{\\_}{-174}^{+223}$ for NSs and
$\Lambda_{\\_}\text{1.4,QS}=642_{\\_}{-204}^{+260}$ for QSs in 90% CI, and
thus QSs would have a significantly larger $\Lambda_{\\_}\text{1.4}$ than NSs.
The obtained $\Lambda_{\\_}\text{1.4,NS}$ nicely agrees with the result
$507_{\\_}{-242}^{+234}$ in Ref. Legred  _et al._ (2021) where the Gaussian
processes are applied to construct the model-independent EOSs. Our result on
$\Lambda_{\\_}\text{1.4,QS}$ is also consistent with
$\Lambda_{\\_}{1.4}=650_{\\_}{-190}^{+230}$ obtained in Ref. Miao  _et al._
(2021) within the MIT bag model.

Figure 1:  Posterior distribution of radius of NSs and QSs at different masses
(90% CI). The NICER mass-radius posterior distributions (90% CI) of PSR
J0030+0451 Miller  _et al._ (2019) and PSR J0740+6620 Miller  _et al._ (2021)
are also shown for comparison.

The above discussions indicate that the NS and QS scenarios lead to different
predictions of $M_{\\_}{\rm TOV}$, $R_{\\_}\text{1.4}$ and
$\Lambda_{\\_}\text{1.4}$ for CSs. To assess the preference of the NS and QS
hypotheses in the description of the current multi-messenger data under the
constraints from pQCD and ChEFT, we evaluate the Bayes factor of NS over QS
and find $\mathcal{B}^{\text{NS}}_{\\_}\text{QS}=11.5$. This value of
$\mathcal{B}^{\text{NS}}_{\\_}\text{QS}$ means that the NS hypothesis for CSs
is strongly preferred against the QS hypothesis according to the
interpretation of Bayes factor, i.e., the
$\mathcal{B}_{\\_}{H_{\\_}0}^{H_{\\_}1}\in[10,30]$ indicates strong evidence
for hypothesis $H_{\\_}1$ Lee and Wagenmakers (2014). To see the main reason
leading to the large value of $\mathcal{B}^{\text{NS}}_{\\_}\text{QS}=11.5$,
we calculate the $\mathcal{B}^{\text{NS}}_{\\_}\text{QS}$ by removing
individually the data from the default data set $\vec{d}_{\\_}\text{def}$, and
we find the value of $\mathcal{B}^{\text{NS}}_{\\_}\text{QS}$ changes to 16.1,
5.1, 10.6, 12.4 and 1.3 by removing the data/constraints of pQCD, GW170817,
GW190425, PSR J0740+6620 and PSR J0030+0451, respectively. The large value of
$\mathcal{B}^{\text{NS}}_{\\_}\text{QS}=11.5$ is thus mainly due to the
constraint from PSR J0030+0451, and next from GW170817. To see more clearly
the influence of PSR J0030+0451, we show in Fig. 1 the 90% CI of radius at
different masses for NSs and QSs. It is seen that the NS hypothesis can indeed
describe the NICER measurements of PSR J0030+0451 much better than the QS
hypothesis as the latter just marginally overlaps with the NICER mass-radius
of PSR J0030+0451 where the probability density is relatively low. Therefore,
our results suggest that the multi-messenger data prefer CSs as NSs over QSs,
implying the conjecture of QM (as well as the strangeon matter) as the true
ground state of QCD matter Bodmer (1971); Witten (1984); Terazawa (1989);
Farhi and Jaffe (1984) is disfavored. This may provide a natural explanation
on the fact that there is so far no definite evidence for the existence of
strangelet-like exotic objects after decades experimental and observational
searching (See, e.g., Refs. Abelev  _et al._ (2007); Han (2009)).

Figure 2: Posterior distributions (68% CI) and the corresponding median values
for squared speed of sound $c_{\\_}s^{2}$ (a), pressure $p$ (b) and trace
anomaly $\Delta$ (c) of NS matter and QS matter as functions of baryon number
density.

The sound velocity $c_{\\_}s$ is an important quantity to feature the EOS of
density matter. Shown in Fig. 2(a) is the $68\%$ CI of $c_{\\_}s^{2}$ as a
function of baryon density for NS matter and QS matter. One sees the squared
speed of sound $c_{\\_}\text{s,NS}^{2}$ for NS matter first increases with
baryon density and reaches a peak value of $c_{\\_}\text{s,NS,max}^{2}\sim
0.5c^{2}$ ($c$ is the speed of light in vacuum) around $n\approx 3.5n_{\\_}0$
(i.e., $n_{\\_}\text{pk,NS}=0.55_{\\_}{-0.14}^{+0.19}$ fm-3), then decreases
and approaches the conformal limit $c^{2}/3$ above $n\approx 4.5n_{\\_}0$.
This peak structure may be related to the quarkyonic matter McLerran and Reddy
(2019) or the high density behavior of the symmetry energy Zhang and Li
(2023). It is interesting to mention that according to percolation theory, the
critical density that nucleons begin overlap with each other is estimated to
be $0.57_{\\_}{-0.09}^{+0.12}$ fm -3 Marczenko  _et al._ (2023), very close to
the $n_{\\_}\text{pk,NS}$. We have checked that the peak structure will
disappear if the pQCD constraint is removed from $\vec{d}_{\\_}\text{def}$.
Furthermore, if the constraints of heavy CSs with $M\approx 2M_{\\_}\odot$ are
excluded from $\vec{d}_{\\_}\text{def}$, the $c_{\\_}\text{s,NS}^{2}$ will
increase monotonously until $12n_{\\_}0$. Therefore, the constraints from pQCD
and heavy CSs with $M\approx 2M_{\\_}\odot$ are necessary conditions for the
sound velocity peak structure in NS matter.

In contrast to $c_{\\_}\text{s,NS}^{2}$, it is interesting to see from Fig.
2(a) that the squared speed of sound $c_{\\_}\text{s,QS}^{2}$ for QS matter
first increases with baryon density and then essentially saturates at about
$0.5c^{2}$ above $n\sim 4n_{\\_}0$. Therefore, the peak structure is not
present in $c_{\\_}\text{s,QS}^{2}$ although the constraints of pQCD and heavy
CSs with $M\approx 2M_{\\_}\odot$ are both considered. This feature clearly
shows that the pQCD limits on the EOS and $c_{\\_}\text{s}^{2}$ of dense
matter at intermediate densities inside CSs significantly depends on the input
EOS at low densities. Indeed, the low density EOS is very different for NS and
QS matters, as illustrated in Fig. 2(b) where the $68\%$ CI of pressure as a
function of baryon density is dislayed for NS and QS matters. One sees from
Fig. 2(b) that the pressure of NS matter is well constrained by BPS EOS and
ChEFT below $1.1n_{\\_}0$ but the pressure of QS matter around $n_{\\_}0$
rapidly drops to zero due to absolutely stable condition although with large
uncertainties. Quantitatively, the energy per baryon $\mu_{\\_}1$ and the
number density $n_{\\_}1$ at zero pressure point of QS matter are estimated to
be $[654,812]$ MeV and $[0.24,0.35]$ fm-3 ($68\%$ CI), respectively.

Very recently, the trace anomaly normalized by the energy density, i.e.,
$\Delta=1/3-p/\epsilon$, is proposed as a new measure of conformality Fujimoto
_et al._ (2022). Shown in Fig. 2(c) is the $68\%$ CI of the (normalized) trace
anomaly $\Delta$ as a function of baryon density for NS and QS matters. One
sees the $\Delta$ for NS matter first decreases with baryon density and then
essentially approaches the conformal limit $\Delta=0$ above $n\approx
4.5n_{\\_}0$. On the other hand, for QS matter, the $\Delta$ decreases
monotonously and becomes negative above $n\approx 5n_{\\_}0$. As pointed out
in Ref. Fujimoto  _et al._ (2022), the sound velocity can be decomposed into
the derivative and the nonderivative terms in terms of $\Delta$, i.e.,
$c_{\\_}\text{s}^{2}/c^{2}=1/3-\Delta-\epsilon d\Delta/d\epsilon$, and the
sound velocity peak observed in NS matter can be attributed to the derivative
term from $\Delta$. In addition, our results indicate while the NS matter
seems to obey the conjecture Fujimoto  _et al._ (2022) that the matter part of
the trace anomaly is positive definite, the QS matter violates it.

Furthermore, we note that the central density in $2M_{\\_}\odot$ and maximum
mass ($M=2.17_{\\_}{-0.15}^{+0.26}M_{\\_}\odot$) NS are estimated to be
$n_{\\_}{\text{c,NS,}2M_{\\_}\odot}=0.56_{\\_}{-0.10}^{+0.14}$ fm -3 and
$n_{\\_}\text{c,NS,max}=0.90_{\\_}{-0.12}^{+0.11}$ fm -3, respectively. The
corresponding values for QS are
$n_{\\_}{\text{c,QS,}2M_{\\_}\odot}=0.54_{\\_}{-0.09}^{+0.11}$ fm -3 and
$n_{\\_}\text{c,QS,max}=0.97_{\\_}{-0.14}^{+0.16}$ fm -3, respectively. These
results imply that the conformal symmetry may be restored in the core of heavy
NSs with $M\gtrsim 2M_{\\_}\odot$, consistent with the conclusion obtained
recently from other groups Fujimoto  _et al._ (2022); Marczenko  _et al._
(2023); Annala  _et al._ (2023).

Conclusion.— Based on Bayesian model-agnostic inference of the properties of
NSs and QSs by combining the multi-messenger data and ab initio calculations
from pQCD and ChEFT, we find that the NS scenario is strongly favored against
the QS scenario for the CSs, and the NS and QS matters display rather
different density behaviors of sound velocity and trace anomaly. Our finding
sheds light on the nature of CSs observed in pulsars and gravitational wave
events in binary star mergers and provides valuable information on the
properties of dense matter inside CSs.

Acknowledgments.— The authors would like to thank Tyler Gorda, Sohpia Han,
Aleksi Kurkela, Ang Li, Yifeng Sun, Renxin Xu, Zhen Zhang and Zhenyu Zhu for
useful discussions. This work was supported by the National SKA Program of
China No. 2020SKA0120300 and the National Natural Science Foundation of China
under Grant Nos. 12235010 and 11625521. The computations in this paper were
run on the Siyuan-1 cluster supported by the Center for High Performance
Computing at Shanghai Jiao Tong University.

## References

  * Lattimer and Prakash (2004) J. M. Lattimer and M. Prakash, Science 304, 536 (2004), arXiv:astro-ph/0405262 .
  * Weber (2005) F. Weber, Progress in Particle and Nuclear Physics 54, 193 (2005).
  * Alford  _et al._ (2008) M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schäfer, Rev. Mod. Phys. 80, 1455 (2008), arXiv:0709.4635 [hep-ph] .
  * Fukushima and Hatsuda (2011) K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Blaschke and Chamel (2018) D. Blaschke and N. Chamel, Astrophys. Space Sci. Libr. 457, 337 (2018), arXiv:1803.01836 [nucl-th] .
  * Burgio  _et al._ (2021) G. F. Burgio, H. J. Schulze, I. Vidana, and J. B. Wei, Prog. Part. Nucl. Phys. 120, 103879 (2021), arXiv:2105.03747 [nucl-th] .
  * Li  _et al._ (2021a) B.-A. Li, B.-J. Cai, W.-J. Xie, and N.-B. Zhang, Universe 7, 182 (2021a), arXiv:2105.04629 [nucl-th] .
  * Sedrakian  _et al._ (2023) A. Sedrakian, J.-J. Li, and F. Weber, Prog. Part. Nucl. Phys. 131, 104041 (2023), arXiv:2212.01086 [nucl-th] .
  * Brambilla  _et al._ (2014) N. Brambilla _et al._ , Eur. Phys. J. C 74, 2981 (2014), arXiv:1404.3723 [hep-ph] .
  * Bodmer (1971) A. R. Bodmer, Phys. Rev. D 4, 1601 (1971).
  * Witten (1984) E. Witten, Phys. Rev. D 30, 272 (1984).
  * Terazawa (1989) H. Terazawa, J. Phys. Soc. Jap. 58, 3555 (1989).
  * Farhi and Jaffe (1984) E. Farhi and R. L. Jaffe, Phys. Rev. D 30, 2379 (1984).
  * Holdom  _et al._ (2018) B. Holdom, J. Ren, and C. Zhang, Phys. Rev. Lett. 120, 222001 (2018), arXiv:1707.06610 [hep-ph] .
  * Zhao  _et al._ (2019) T. Zhao, W. Zheng, F. Wang, C.-M. Li, Y. Yan, Y.-F. Huang, and H.-S. Zong, Phys. Rev. D 100, 043018 (2019), arXiv:1904.09744 [nucl-th] .
  * Cao  _et al._ (2022) Z. Cao, L.-W. Chen, P.-C. Chu, and Y. Zhou, Phys. Rev. D 106, 083007 (2022), arXiv:2009.00942 [astro-ph.HE] .
  * Xu (2003) R.-X. Xu, Astrophys. J. Lett. 596, L59 (2003), arXiv:astro-ph/0302165 .
  * Xu (2018) R. Xu, Sci. China Phys. Mech. Astron. 61, 109531 (2018), arXiv:1802.04465 [astro-ph.HE] .
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Holt  _et al._ (2016) J. W. Holt, M. Rho, and W. Weise, Phys. Rept. 621, 2 (2016), arXiv:1411.6681 [nucl-th] .
  * Drischler  _et al._ (2021) C. Drischler, J. W. Holt, and C. Wellenhofer, Ann. Rev. Nucl. Part. Sci. 71, 403 (2021), arXiv:2101.01709 [nucl-th] .
  * Hebeler  _et al._ (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Wellenhofer  _et al._ (2015) C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev. C 92, 015801 (2015), arXiv:1504.00177 [nucl-th] .
  * Gorda  _et al._ (2021) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Phys. Rev. Lett. 127, 162003 (2021), arXiv:2103.05658 [hep-ph] .
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Annala  _et al._ (2022) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Phys. Rev. X 12, 011058 (2022), arXiv:2105.05132 [astro-ph.HE] .
  * Gorda  _et al._ (2023) T. Gorda, O. Komoltsev, and A. Kurkela, Astrophys. J. 950, 107 (2023), arXiv:2204.11877 [nucl-th] .
  * Tews  _et al._ (2018) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018), arXiv:1801.01923 [nucl-th] .
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Buchner  _et al._ (2014) J. Buchner, A. Georgakakis, K. Nandra, L. Hsu, C. Rangel, M. Brightman, A. Merloni, M. Salvato, J. Donley, and D. Kocevski, Astron. Astrophys. 564, A125 (2014), arXiv:1402.0004 [astro-ph.HE] .
  * Ashton  _et al._ (2019) G. Ashton _et al._ , Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Hernandez Vivanco  _et al._ (2020) F. Hernandez Vivanco, R. Smith, E. Thrane, and P. D. Lasky, Mon. Not. Roy. Astron. Soc. 499, 5972 (2020), arXiv:2008.05627 [astro-ph.HE] .
  * Thrane and Talbot (2019) E. Thrane and C. Talbot, Publ. Astron. Soc. Austral. 36, e010 (2019), [Erratum: Publ.Astron.Soc.Austral. 37, e036 (2020)], arXiv:1809.02293 [astro-ph.IM] .
  * Jeffreys (1961) H. Jeffreys,  _Theory of Probability_ (Oxford University Press, New York, 1961).
  * Huth  _et al._ (2022) S. Huth _et al._ , Nature 606, 276 (2022), arXiv:2107.06229 [nucl-th] .
  * Komoltsev and Kurkela (2022) O. Komoltsev and A. Kurkela, Phys. Rev. Lett. 128, 202701 (2022), arXiv:2111.05350 [nucl-th] .
  * (46) Likelihood function of QCD, https://github.com/OKomoltsev/QCD-likelihood-function.
  * Landry  _et al._ (2020) P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
  * Biswas (2022) B. Biswas, Astrophys. J. 926, 75 (2022), arXiv:2106.02644 [astro-ph.HE] .
  * Demorest  _et al._ (2010) P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Fonseca  _et al._ (2016) E. Fonseca _et al._ , Astrophys. J. 832, 167 (2016), arXiv:1603.00545 [astro-ph.HE] .
  * Arzoumanian  _et al._ (2018) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Suppl. 235, 37 (2018), arXiv:1801.01837 [astro-ph.HE] .
  * Antoniadis  _et al._ (2013) J. Antoniadis _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * (53) NICER PSR J0030+0451 Illinois-Maryland MCMC Samples, J00303spotRM.txt, https://zenodo.org/record/3473466.
  * (54) NICER PSR J0740+6620 Illinois-Maryland MCMC Samples, NICER+XMMJ0740RM.txt, https://zenodo.org/record/4670689.
  * Hernandez Vivanco  _et al._ (2019) F. Hernandez Vivanco, R. Smith, E. Thrane, P. D. Lasky, C. Talbot, and V. Raymond, Phys. Rev. D 100, 103009 (2019), arXiv:1909.02698 [gr-qc] .
  * Han  _et al._ (2023) M.-Z. Han, Y.-J. Huang, S.-P. Tang, and Y.-Z. Fan, Sci. Bull. 68, 913 (2023), arXiv:2207.13613 [astro-ph.HE] .
  * Altiparmak  _et al._ (2022) S. Altiparmak, C. Ecker, and L. Rezzolla, Astrophys. J. Lett. 939, L34 (2022), arXiv:2203.14974 [astro-ph.HE] .
  * Pang  _et al._ (2022) P. T. H. Pang _et al._ , (2022), arXiv:2205.08513 [astro-ph.HE] .
  * Pang  _et al._ (2021) P. T. H. Pang, I. Tews, M. W. Coughlin, M. Bulla, C. Van Den Broeck, and T. Dietrich, Astrophys. J. 922, 14 (2021), arXiv:2105.08688 [astro-ph.HE] .
  * Legred  _et al._ (2021) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Phys. Rev. D 104, 063003 (2021), arXiv:2106.05313 [astro-ph.HE] .
  * Biswas  _et al._ (2021) B. Biswas, P. Char, R. Nandi, and S. Bose, Phys. Rev. D 103, 103015 (2021), arXiv:2008.01582 [astro-ph.HE] .
  * Li  _et al._ (2021b) Y. Li, H. Chen, D. Wen, and J. Zhang, Eur. Phys. J. A 57, 31 (2021b), arXiv:2008.02955 [nucl-th] .
  * Han  _et al._ (2021) M.-Z. Han, J.-L. Jiang, S.-P. Tang, and Y.-Z. Fan, Astrophys. J. 919, 11 (2021), arXiv:2103.05408 [hep-ph] .
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Jiang  _et al._ (2020) J.-L. Jiang, S.-P. Tang, Y.-Z. Wang, Y.-Z. Fan, and D.-M. Wei, Astrophys. J. 892, 1 (2020), arXiv:1912.07467 [astro-ph.HE] .
  * Radice and Dai (2019) D. Radice and L. Dai, Eur. Phys. J. A 55, 50 (2019), arXiv:1810.12917 [astro-ph.HE] .
  * Miao  _et al._ (2021) Z. Miao, J.-L. Jiang, A. Li, and L.-W. Chen, Astrophys. J. Lett. 917, L22 (2021), arXiv:2107.13997 [astro-ph.HE] .
  * Lee and Wagenmakers (2014) M. D. Lee and E.-J. Wagenmakers,  _Bayesian Cognitive Modeling_ (Cambridge University Press, Cambridge, England, 2014).
  * Abelev  _et al._ (2007) B. I. Abelev _et al._ (STAR), Phys. Rev. C 76, 011901 (2007), arXiv:nucl-ex/0511047 .
  * Han (2009) K. Han (LSSS), J. Phys. G 36, 064048 (2009), arXiv:0812.3885 [nucl-ex] .
  * McLerran and Reddy (2019) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019), arXiv:1811.12503 [nucl-th] .
  * Zhang and Li (2023) N.-B. Zhang and B.-A. Li, Eur. Phys. J. A 59, 86 (2023), arXiv:2208.00321 [nucl-th] .
  * Marczenko  _et al._ (2023) M. Marczenko, L. McLerran, K. Redlich, and C. Sasaki, Phys. Rev. C 107, 025802 (2023), arXiv:2207.13059 [nucl-th] .
  * Fujimoto  _et al._ (2022) Y. Fujimoto, K. Fukushima, L. D. McLerran, and M. Praszalowicz, Phys. Rev. Lett. 129, 252702 (2022), arXiv:2207.06753 [nucl-th] .
  * Annala  _et al._ (2023) E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, J. Nättilä, and A. Vuorinen, (2023), arXiv:2303.11356 [astro-ph.HE] .