# On the Moment of Inertia of PSR J0737-3039 A from LIGO/Virgo and NICER

Zhiqiang Miao1; Ang Li1; Zi-Gao Dai2,3
1Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, China;
<EMAIL_ADDRESS>
2Department of Astronomy, School of Physical Sciences, University of Science
and Technology of China, Hefei 230026, China
3School of Astronomy and Space Science, Nanjing University, Nanjing 210023,
China

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We perform a Bayesian analysis of neutrons star moment of inertia by utilizing
the available gravitational-wave data from LIGO/Virgo (GW170817 and GW190425)
and mass-radius measurements from the Neutron Star Interior Composition
Explorer (PSR J0030+0415 and PSR J0740+6620), incorporating the possible phase
transition in the pulsar inner core. We find that the moment of inertia of
pulsar A in the double pulsar binary J0737-3039 is $\sim 1.30\times
10^{45}\,{\rm g\,cm^{2}}$, which only slightly depends on the employed
hadronic equation of states. We also demonstrate how a moment of inertia
measurement would improve our knowledge of the equation of state and the mass-
radius relation for neutron stars and discuss whether a quark deconfinement
phase transition is supported by the available data and forthcoming data that
could be consistent with this hypothesis. We find that if pulsar A is a quark
star, that its moment of inertia is a large value of $\sim 1.55\times
10^{45}\,{\rm g\,cm^{2}}$ suggesting the possibility of distinguishing it from
(hybrid-)neutron stars with measurements of PSR J0737-3039A moment of inertia.
We finally demonstrate the moment-of-inertia–compactness universal relations
and provide analytical fits for both (hybrid-)neutron star and quark star
results based on our analysis.

###### keywords:

dense matter – equation of state – stars: neutron – gravitational waves –
pulsars: individual (PSR J0737-3039)

††pagerange: On the Moment of Inertia of PSR J0737-3039 A from LIGO/Virgo and
NICER–LABEL:lastpage

##  1 Introduction

Neutron stars are the densest and smallest stars observed in the Universe.
Observations of neutron stars provide valuable information on the internal
structure of these objects and, therefore, on the nuclear equation of state
(EOS) of dense matter in their interior. Various observables are tightly
related with microphysics inside the star and can help to set reliable
constraints on its EOS. Among them, there are the mass (Demorest et al., 2010;
Antoniadis et al., 2013; Cromartie et al., 2020; Fonseca et al., 2021), the
radius (Riley et al., 2019a, 2021a; Miller et al., 2019, 2021), and the tidal
deformability (Abbott et al., 2017, 2018, 2020). In addition, measurements are
planned for the moment of inertia (MOI) of PSR J0737-3039 A (hereafter A)
(Lyne et al., 2004), the $1.338\,M_{\odot}$ primary component of the first
double pulsar system PSR J0737-3039 (Burgay et al., 2003), based on the long-
term pulsar timing to determine the periastron advance due to relativistic
spin-orbit coupling (Damour & Schafer, 1988). The MOI of A is expected to be
measured with $\sim 10\%$ accuracy within the next decade, which will provide
complementary constraints on the EOS (Bejger, Bulik, & Haensel, 2005; Lattimer
& Schutz, 2005; Greif et al., 2020; Hu et al., 2020).

The neutron star MOI has been theoretically calculated using a wide range of
EOS models (e.g., Yuan & Zhang, 1999; Morrison et al., 2004; Lattimer &
Schutz, 2005; Bejger, Bulik, & Haensel, 2005; Rikovska Stone et al., 2007;
Klähn et al., 2007; Read et al., 2009; Urbanec, Miller, & Stuchlík, 2013; Li,
2015; Li et al., 2016; Motta et al., 2019; Sun, Liu, & Xing, 2019; Wei et al.,
2019; Lourenço et al., 2020; Margaritis, Koliogiannis, & Moustakidis, 2020).
Several studies have been performed using EOS models, including those
incorporating non-nucleon baryonic components (such as $\Delta$-isobars,
hyperons and antikaon condensates, etc.) (e.g., Klähn et al., 2007;
Bandyopadhyay et al., 2018; Sen & Jha, 2019; Ribes et al., 2019; Fortin et
al., 2020; Raduta, Oertel, & Sedrakian, 2020).

In several recent studies, the MOI of A has been predicted based on EOSs
constrained using physical stability arguments (Raithel, Özel, & Psaltis,
2016), neutron star observations (Steiner et al., 2015; Landry & Kumar, 2018;
Kumar & Landry, 2019; Jiang et al., 2020; Silva et al., 2021), or the results
of nuclear experiments (Worley, Krastev, & Li, 2008; Lim, Holt, & Stahulak,
2019). However, a phase transition was not explicitly taken into account in
these previous studies. In the inner stellar core, there may be a possibility
that the hadron pressure is so high that nucleons should reasonably melt down
to the deconfined phase, and EOSs with an early phase transition at $\sim
2n_{0}$ (where $n_{0}=0.16\,{\rm~{}fm^{-3}}$ denotes the nuclear saturation
density) are compatible with current neutron star observations (Al-Mamun et
al., 2021; Li et al., 2020, 2021a; Miao et al., 2020; Xie & Li, 2021),
therefore a description of the high-density core matter in terms of only
nucleons might be inadequate for the the MOI study of A.

Furthermore, in the extreme case of quark deconfinement phase transition,
neutron stars may convert completely to quark stars (e.g., Olinto, 1987;
Olesen & Madsen, 1994; Ouyed, Dey, & Dey, 2002; Herzog & Röpke, 2011),
resulting in various astrophysical implications, such as neutrino and
gravitational waves emission (e.g., Pagliara, Herzog, & Röpke, 2013; Bombaci
et al., 2016) as well as short gamma-ray bursts (e.g., Cheng & Dai, 1996;
Ouyed, Dey, & Dey, 2002; Berezhiani et al., 2003). There are also discussions
on the trigger of such transition due to self-annihilation of dark matter
(Pérez-García Silk, & Stone, 2010) and the characteristic delay due to the
transition that could lead to observational signatures helping distinguish
different central engine models of short gamma-ray bursts (Pérez-García,
Daigne, & Silk, 2013). Therefore, it is interesting and meaningful to consider
the possibility of A as a quark star and discuss all three different MOI
results when assuming A is a pure neutron star, a hybrid star, and a pure
quark star. By doing so, we are in the right place to investigate whether a
future MOI measurement to $10\%$ precision is sufficient to distinguish among
the three qualitatively different models. What particularly interesting is
that, since the quark matter in hybrid stars and quark stars are described
under different stability conditions, and both of them can’t coexist stably
(Weber, 2005), a confirmation (or exclusion) of either of them can shed light
on the remaining-unsolved QCD at finite baryon density (Gorda et al., 2021).

This paper is organized as follows. In Section 2, we introduce the EOSs used
to model the hadronic phase of neutron stars and a general parameterization
for describing the phase transition and the quark matter phase of neutron
stars. The neutron star observations used and the Bayesian analysis method are
detailed in Section 3. We present our results and a discussion in Section 4
and summarize our paper in Section 5.

##  2 Neutron stars with or without a quark core

Following our previous study (Li et al., 2021a), we describe neutron stars
with and without a quark core below.

For neutron stars without a quark core, we employ the unified QMF (Zhu, Zhou,
& Li, 2018) and DD2 (Fortin et al., 2016) EOS models for soft and stiff
hadronic matter, while using the standard BPS EOS for the neutron star outer
crust (Baym, Pethick, & Sutherland, 1971). Solving the Tolman-Oppenheimer-
Volkoff (TOV) equations yields the neutron star maximum mass as $M_{\rm
TOV}^{\rm QMF}=2.07\,M_{\odot}$ and $M_{\rm TOV}^{\rm DD2}=2.42\,M_{\odot}$,
with the radius of a typical $1.4\,M_{\odot}$ star as $R_{1.4}^{\rm
QMF}=11.77\,{\rm km}$ and $R_{1.4}^{\rm DD2}=13.17\,{\rm km}$, respectively.
Furthermore, as the rotation frequency of A, $\Omega_{\rm A}=276.8\,{\rm Hz}$
(Burgay et al., 2003) is far smaller than the Kepler frequency, one can employ
the slow rotation approximation (Hartle, 1967) to calculate the MOI of the
$1.338$-solar-mass A, which is $I_{A}^{\rm QMF}=1.33\times 10^{45}\,{\rm
g\,cm^{2}}$ and $I_{A}^{\rm DD2}=1.60\times 10^{45}\,{\rm g\,cm^{2}}$ using
the two EOSs.

For neutron stars with a quark core, we implement a sharp first-order phase
transition between the hadronic and quark phases by employing the “constant-
speed-of-sound” (CSS) parameterization (Alford, Han, & Prakash, 2013), while
maintaining the hadronic EOS as QMF or DD2. The full EOS can be expressed as
follows (we use units where $\hbar=c=1$):

$\varepsilon(p)=\left\\{\\!\begin{aligned} &\varepsilon_{\rm HM}(p)\
,&p<p_{\rm trans}\\\ &\varepsilon_{\rm HM}(p_{\rm
trans})+\Delta\varepsilon+c_{\rm QM}^{-2}(p-p_{\rm trans})\ ,&p>p_{\rm
trans}\end{aligned}\right.$ (1)

where $p_{\rm trans}$ is the pressure of the transition, $\Delta\varepsilon$
is the discontinuity in the energy density, and $c_{\rm QM}$ is the sound
speed in quark matter. Correspondingly, three dimensionless parameters are
chosen as the EOS parameters of hybrid stars, namely, the transition density
$n_{\rm trans}/n_{0}$, the transition strength
$\Delta\varepsilon/\varepsilon_{\rm trans}$, and the sound speed squared in
quark matter $c_{\rm QM}^{2}$, where $n_{\rm trans}\equiv n_{\rm HM}(p_{\rm
trans})$ and $\varepsilon_{\rm trans}\equiv\varepsilon_{\rm HM}(p_{\rm
trans})$.

In the following, within the context of a first-order phase transition from
hadronic matter to quark matter in the neutron star interior, we intend to
make a direct estimate of the MOI of A based on a Bayesian analysis of
observational data for the tidal deformability of the GW170817 (Abbott et al.,
2017) and GW190425 (Abbott et al., 2020) binary neutron star mergers, as
detected by LIGO/Virgo, and the masses and radii of PSR J0030+0451 (Lommen et
al., 2000; Becker et al., 2000) and PSR J0740+6620 (Stovall et al., 2014), as
detected by NICER (Miller et al., 2019; Riley et al., 2019a; Miller et al.,
2021; Riley et al., 2021a). The EOS parameters are varied to obtain a 3D
parameter space for the EOSs, and a Bayesian analysis is performed to select
preferred EOSs based on the neutron star observations from LIGO/Virgo and
NICER. The corresponding study on quark stars are then followed (see below in
Sec. 4.5). The Bayesian analysis framework is introduced in the next section.

##  3 Observational constraints and Bayesian analysis

In the present study, our observational data set $\boldsymbol{d}$ contains two
gravitational wave events (GW170817 and GW190425) and two NICER mass-radius
measurements of pulsars (PSR J0030+0451 and PSR J0740+6620). Assuming each
dataset in $\boldsymbol{d}$ is independent of one another, the total
likelihood can be expressed as:

$\begin{split}\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta},\boldsymbol{p}_{c})=\prod_{i}\mathcal{L}_{{\rm
GW},i}\times\prod_{j}\mathcal{L}_{{\rm NICER},j}\ ,\end{split}$ (2)

where the priors of the parameters are
$p(\boldsymbol{\theta},\boldsymbol{p}_{c})=p(\boldsymbol{\theta})p(\boldsymbol{p}_{c})$,
with $\boldsymbol{p}_{c}$ being the central pressure of the NICER’s pulsars.
We employ the same priors as in our previous study (Li et al., 2021a) for EOS
parameters $\boldsymbol{\theta}=\\{n_{\rm
trans}/n_{0},\Delta\varepsilon/\varepsilon_{\rm trans},c_{\rm QM}^{2}\\}$.

GW170817 and GW190425. We calculate the likelihood of each gravitational event
through a high-precision interpolation of the likelihood developed in
Hernandez Vivanco et al. (2020) from fitting the strain data released by
LIGO/Virgo, which is encapsulated in the python package toast, i.e.,

$\mathcal{L}_{{\rm GW},i}=F_{i}(\Lambda_{1},\Lambda_{2},\mathcal{M}_{c},q)\ ,$
(3)

where $F_{i}(\cdot)$ is the interpolation function of the gravitational wave
event $i$. The chirp mass
$\mathcal{M}_{c}\equiv(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$ and the mass
ratio $q\equiv M_{2}/M_{1}$. The tidal deformabilities can be derived from the
component masses, i.e., $\Lambda_{1}=\Lambda_{1}(\boldsymbol{\theta};M_{1})$
and $\Lambda_{2}=\Lambda_{2}(\boldsymbol{\theta};M_{2})$. We mention here that
since the likelihood involves no EOS assumption about the phase state of the
merging neutron stars, we can adopt it for the study of both (hybrid) neutron
stars and quark stars.

NICER mass-radius measurements. Recently, the NICER collaboration reported
simultaneous mass-radius measurements of PSR J0030+0451 and PSR J0740+6620.
The results for PSR J0030+0451 were $M=1.34_{-0.16}^{+0.15}\,M_{\odot}$,
$R=12.71_{-1.19}^{+1.14}\,{\rm km}$ (Riley et al., 2019a) or
$M=1.44_{-0.14}^{+0.15}\,M_{\odot}$, $R=13.02_{-1.06}^{+1.24}\,{\rm km}$
(Miller et al., 2019). The results for PSR J0740+6620 were
$M=2.072_{-0.066}^{+0.067}\,M_{\odot}$, $R=12.39_{-0.98}^{+1.30}\,{\rm km}$
(Riley et al., 2021a) or $M=2.062_{-0.091}^{+0.090}\,M_{\odot}$,
$R=13.71_{-1.50}^{+2.61}\,{\rm km}$ (Miller et al., 2021). To incorporate each
of these measurements, we write the likelihood as:

$\mathcal{L}_{{\rm NICER},j}={\rm KDE}(M,R|\boldsymbol{S}_{j})\ ,$ (4)

where the right-hand side is a kernel density estimation of the posterior
samples $\boldsymbol{S}_{j}$ of the mass-radius measurement $j$. Similar to Li
et al. (2021a), we use the ST+PST model samples from Riley et al. (2019b) for
PSR J0030+0451. For PSR J0740+6620, we use the NICER x XMM samples from Riley
et al. (2021b). The mass and radius $(M,R)$ are mapped from the EOS parameters
$\boldsymbol{\theta}$ and the central pressure $\boldsymbol{p}_{c}$ through
the TOV equations. For PSR J0030+0451, we set $p_{c}$ (in unit ${\rm
dyn/cm^{2}}$) to uniformly distribute in logarithmic space between
$\log(p_{c})=34.5$ and $\log(p_{c})=35.5$, while for PSR J0740+6620, $p_{c}$
uniformly distributes in logarithmic space from $\log(p_{c})=35$ to
$\log(p_{c})=36$. In addition, we consider parameter sets
$\\{\boldsymbol{\theta},\boldsymbol{p}_{c}\\}$ that correspond to stable
neutron stars, i.e., those $\\{\boldsymbol{\theta},\boldsymbol{p}_{c}\\}$
correspond to unstable neutron stars should be rejected in sampling. See more
discussions in our previous study (Li et al., 2021a).

The posterior distribution is proportional to the product of the prior and the
likelihood function, i.e.,

$p(\boldsymbol{\theta},\boldsymbol{p}_{c}|\boldsymbol{d})\propto
p(\boldsymbol{\theta},\boldsymbol{p}_{c})\mathcal{L}(\boldsymbol{d}|\boldsymbol{\theta},\boldsymbol{p}_{c})\
.$ (5)

In practice, we sample from the posterior distribution by using the nested
sampling software PyMultiNest (Buchner et al., 2014). We finally obtain a
total of approximately 3,000 posterior samples of the CSS model parameters for
both QMF+CSS and DD2+CSS. These EOS parameter samples are used in the study of
(hybrid)neutron star properties.

##  4 Results and discussion

Figure 1: Posterior distribution of $\bar{I}\equiv I/M^{3}$ vs. $\Lambda$ for
a $1.338\,M_{\odot}$ star. The dots are calculated from the corresponding EOS
parameter posterior samples. The marginalized posterior distributions of
$\Lambda$ and $\bar{I}$ are shown in the upper and right panel, respectively.
The black dashed line shows the fitting formula from Landry & Kumar (2018),
updated from the original relation of Yagi & Yunes (2013) with a larger set of
EOSs.

###  4.1 MOI of PSR J0737-3039 A

Fig. 1 presents the posterior distributions of the dimensionless MOI
$\bar{I}\equiv I/M^{3}$ vs. the tidal deformability $\Lambda$. The analysis is
performed for a fixed star mass $M=1.338\,M_{\odot}$, corresponding to A. Our
results can be effectively reconciled with the universal relation of neutron
stars (shown as a black line in Fig. 1) (Yagi & Yunes, 2013; Landry & Kumar,
2018), with an error less than 0.5%, indicating that the I-Love relation
remains valid within the present hybrid star scenario (see also in, e.g., Miao
et al., 2020). Such universality is expected since hybrid stars and neutron
stars differ only in the high-density core part and it is the low-density
crust part that mainly determines the $\bar{I}$-$\Lambda$ curves with
different EOSs (Yagi & Yunes, 2013).

Figure 2: Theoretical bounds on the MOI of J0737-3039 A as a function of its radius, within both QMF+CSS and DD2+CSS (90% credible interval). The marginalized posterior distributions of $R$ and $I$ are shown in the upper and right panel, respectively. The black curve corresponds to an extreme bound from Raithel, Özel, & Psaltis (2016), assuming the AP4 EOS up to nuclear saturation density $n_{0}$ and a constant-density core inside the neutron stars.  Table 1: Theoretical predictions on the MOI of J0737-3039 A ($M=1.338\,M_{\odot}$), given in $10^{45}\,{\rm g\,cm^{2}}$. Except the values of Worley, Krastev, & Li (2008), the quoted error bars refer to symmetric 90% credible intervals about the median. We mention here that the current $90\%$ upper limit from recent double pulsar timing measurements is $3\times 10^{45}\,\rm g\,cm^{2}$ (Kramer et al., 2021).  Reference | $I_{A}$
---|---
Worley, Krastev, & Li (2008) | $[1.30,1.63]$
Landry & Kumar (2018) | $1.15_{-0.24}^{+0.38}$
Lim, Holt, & Stahulak (2019) | $1.36_{-0.26}^{+0.12}$
Jiang et al. (2020) | $1.35_{-0.14}^{+0.26}$
Silva et al. (2021) | $1.68_{-0.48}^{+0.53}$
This work (QMF+CSS) | $1.27_{-0.14}^{+0.18}$
This work (DD2+CSS) | $1.29_{-0.15}^{+0.26}$

In Fig. 2, we show the 90% bounds on the MOI for A as a function of the A
radius. Also shown is the extreme bound from Raithel, Özel, & Psaltis (2016),
which is derived by placing one or two constant-density cores inside the
neutron stars based on stability arguments, whereas the outer layer (with a
density lower than the nuclear saturation density) is described by the AP4
EOS. The EOS modelling results based on currently available neutron star
observations fall within a fairly narrow region of the $I-R$ plane.

In Table 1, we compare our MOI prediction for A with those of previous
studies. Using EOS-insensitive relations, Landry & Kumar (2018) resulted in an
MOI of $\sim 1.15\times 10^{45}\,{\rm g\,cm^{2}}$ based on the GW170817 tidal
measurement. In Silva et al. (2021), a larger MOI of $\sim 1.68\times
10^{45}\,{\rm g\,cm^{2}}$ was obtained based on the NICER mass-radius
measurement of PSR J0030+0451 at 90% credible interval. In Jiang et al.
(2020), both data sets were used with parameterized EOSs to obtain an MOI
between $\sim 1.35\times 10^{45}\,{\rm g\,cm^{2}}$. In the present study, we
not only explicitly include the possibility of a phase transition in the EOS
prior but also incorporate two additional sources, GW190425 and PSR
J0740+6620, into the analysis. Our results indicate that
$I_{A}=1.27_{-0.14}^{+0.18}\times 10^{45}\,{\rm g\,cm^{2}}$ (or $I/M_{\rm
A}^{3/2}=41.3_{-4.5}^{+5.8}{\rm km^{2}}/\,M_{\odot}^{1/2}$) within QMF+CSS or
$1.29_{-0.15}^{+0.26}\times 10^{45}\,{\rm g\,cm^{2}}$ (or $I/M_{\rm
A}^{3/2}=41.9_{-4.9}^{+8.4}{\rm km^{2}}/\,M_{\odot}^{1/2}$) within DD2+CSS at
$90\%$ credible interval. Predictions have also been made based on EOSs
constrained by relativistic heavy-ion collisions (Worley, Krastev, & Li, 2008)
or microscopic EOSs consistent with empirical data for finite nuclei (Lim,
Holt, & Stahulak, 2019), as shown in Table 1.

Figure 3: The $90\%$ credible regions of $I/M^{3/2}$ as a function of $M$,
transformed from the EOS parameter posteriors. Also shown in solid grey lines
are the results of two neutron star EOSs (QMF and DD2). The vertical line
corresponds to $M=1.338\,M_{\odot}$, while the error bar represents a
hypothetical $10\%$ precision measurement centered on $I/M_{\rm
A}^{3/2}=40\,{\rm km^{2}}/\,M_{\odot}^{1/2}$.

For subsequent calculations, Fig 3 shows the posterior distributions of
$I/M^{3/2}$ as a function of the mass $M$. The scaled MOI is found to lie
within the range $I/M^{3/2}\sim 35-50\,{\rm km^{2}}/\,M_{\odot}^{1/2}$. Note
that in Steiner et al. (2015), a similar and somewhat narrower range of $\sim
35-45\,{\rm km^{2}}/\,M_{\odot}^{1/2}$ was calculated based on two sets of
parameterized EOSs (which were denoted as Models A and C) obtained from
analyzing several mass-radius measurements from photospheric radius expansion
(PRE) bursts or quiescent low-mass X-ray binaries.

Below, we use the implications of the LIGO/Virgo and NICER data for the MOI of
A (as shown in Fig. 3) to analyze the available MOI data and the mocked data
consisting of a future MOI measurement of A with an accuracy of $10\%$. In
Sec. 4.2, we demonstrate how the available MOI measurements add no new
information on the EOS of neutron star matter due to its huge uncertainties.
In Sec. 4.3, we present updated constraints on the EOS and the corresponding
mass-radius relations of the stars from a mocked MOI data with an accuracy of
$10\%$. In Sec. 4.4, we explore further whether a phase with deconfined quarks
is supported by current and future observational data. In Sec. 4.5, we further
discuss the possibility of distinguishing quark stars from (hybrid-)neutron
stars with the MOI measurements of A.

###  4.2 Constraints with available MOI measurements

Recently, an MOI measurement of pulsar A, with a 90% upper limit
$I_{A}<3\times 10^{45}\,{\rm g\ cm^{2}}$ is reported in Kramer et al. (2021),
by fitting the post-Keplerian parameters in the binary system. It is mentioned
in Kramer et al. (2021) that such a measurement gives a loose upper limit for
A’ radius of 22 km (with 90% confidence) and is regarded as not constraining
to the EOS, in comparison with the accurate observations from, e.g.,
LIGO/Virgo, where $R_{\rm 1.4}\leqslant 13.5$ km is obtained (Abbott et al.,
2018). We here report a similar conclusion in Fig. 4, where we compare the
prior distribution and the posterior distribution of the MOI of pulsar A.

The prior of $I_{A}$ is calculated from the QMF+CSS model with the prior of
the CSS model parameters in Sec. 3. The posterior is obtained from the prior
and the likelihood of the MOI measurement reported in Kramer et al. (2021). In
practice, we equal the likelihood function to the posterior distribution of
the MOI of A reported in Kramer et al. (2021), which can be well fitted by a
Gaussian distribution with mean $\mu_{I}=0.283\times 10^{45}\,{\rm g\ cm^{2}}$
and standard deviation $\sigma_{I}=1.705\times 10^{45}\,{\rm g\ cm^{2}}$,
i.e.,

$\mathcal{L}_{I}=\frac{1}{\sqrt{2\pi}\sigma_{I}}\exp\left\\{-\frac{(I_{A}-\mu_{I})^{2}}{2\sigma_{I}^{2}}\right\\}\
.$ (6)

Again, we sample from posterior distribution and the result is shown in Fig.
4. We also calculate the Kullback–Leibler divergence from prior to posterior
in bits for $I_{A}$, $D_{\rm KL}=0.066$. The small value of $D_{\rm KL}$ and
the visual comparison of the prior and posterior indicate that the available
MOI measurements do provide no helpful constraints on the EOS, with respect to
the LIGO/Virgo and NICER data.

Figure 4: Prior and posterior distributions of the MOI of pulsar A. The prior of $I_{A}$ is obtained from the QMF+CSS model with the prior of the CSS model parameters in Sec. 3, while the posterior is conditioned on the the MOI result from Kramer et al. (2021). The Kullback–Leibler divergence, $D_{\rm KL}$, from prior to posterior in bits for $I_{A}$ is also shown.  Table 2: The resulting sound speed squared in quark matter $c_{\rm QM}^{2}$, the transition pressure $p_{\rm trans}$ and discontinuity on the energy density $\Delta\varepsilon$, the maximum mass $M_{\rm TOV}$ and the radius of a $1.4\,M_{\odot}$ star $R_{1.4}$ are collected for both cases of QMF+CSS and DD2+CSS, from our previous work (Li et al., 2021a) (in the 2nd column) and the present work (in the 3rd column). In the last four columns, we show the results of four tests when a mocked MOI measurement is added, being $I/M_{\rm A}^{3/2}=35,40,45,50\,{\rm km^{2}}/\,M_{\odot}^{1/2}$, respectively. The quoted error bars refer to symmetric 90% credible intervals about the median. The Bayes factors of the hybrid star against the normal neutron star, $\mathcal{B}_{\rm NS}^{\rm HS}$, are also shown for each case. See text for details. | GW170817+ PSR J0030 | GW170817+GW190425 | +Mocked MOI measurement
---|---|---|---
(Li et al. 2021) | +PSR J0030+ PSR J0740 | Test i | Test ii | Test iii | Test iv
QMF+CSS
$c_{\rm QM}^{2}$ | $0.81_{-0.28}^{+0.17}$ | $0.73_{-0.30}^{+0.24}$ | $0.77_{-0.27}^{+0.21}$ | $0.74_{-0.27}^{+0.24}$ | $0.71_{-0.28}^{+0.26}$ | $0.70_{-0.30}^{+0.27}$
$p_{\rm trans}/{\rm MeV\,fm^{-3}}$ | $17.6_{-14.0}^{+37.9}$ | $16.8_{-13.9}^{+42.5}$ | $17.1_{-14.1}^{+31.2}$ | $19.3_{-16.4}^{+38.1}$ | $17.9_{-15.0}^{+41.9}$ | $12.8_{-10.3}^{+45.8}$
$\Delta\varepsilon/{\rm MeV\,fm^{-3}}$ | $87.0_{-72.8}^{+84.2}$ | $94.9_{-84.0}^{+87.8}$ | $134.5_{-78.4}^{+72.8}$ | $101.9_{-82.4}^{+81.8}$ | $78.5_{-68.9}^{+86.4}$ | $63.6_{-55.8}^{+76.0}$
$M_{\rm TOV}/\,M_{\odot}$ | $2.36_{-0.26}^{+0.49}$ | $2.23_{-0.21}^{+0.42}$ | $2.22_{-0.20}^{+0.29}$ | $2.24_{-0.22}^{+0.35}$ | $2.23_{-0.22}^{+0.44}$ | $2.25_{-0.24}^{+0.54}$
$R_{1.4}/{\rm km}$ | $11.70_{-0.74}^{+0.85}$ | $11.61_{-0.76}^{+0.75}$ | $11.21_{-0.60}^{+0.60}$ | $11.52_{-0.67}^{+0.48}$ | $11.74_{-0.68}^{+0.60}$ | $11.88_{-0.61}^{+0.75}$
$\mathcal{B}_{\rm NS}^{\rm HS}$ | 1.49 | 0.63 | 1.77 | 0.66 | 0.43 | 0.44
DD2+CSS
$c_{\rm QM}^{2}$ | $0.80_{-0.29}^{+0.18}$ | $0.71_{-0.29}^{+0.26}$ | $0.78_{-0.26}^{+0.20}$ | $0.73_{-0.27}^{+0.25}$ | $0.69_{-0.28}^{+0.27}$ | $0.68_{-0.28}^{+0.29}$
$p_{\rm trans}/{\rm MeV\,fm^{-3}}$ | $14.1_{-9.6}^{+18.2}$ | $15.2_{-11.3}^{+21.1}$ | $12.8_{-8.9}^{+15.6}$ | $12.9_{-9.2}^{+18.9}$ | $15.2_{-11.4}^{+20.5}$ | $20.8_{-16.7}^{+18.0}$
$\Delta\varepsilon/{\rm MeV\,fm^{-3}}$ | $141.2_{-89.4}^{+91.4}$ | $149.9_{-96.6}^{+83.3}$ | $195.7_{-53.0}^{+68.1}$ | $168.0_{-70.9}^{+67.8}$ | $134.0_{-74.8}^{+77.9}$ | $106.1_{-73.3}^{+83.1}$
$M_{\rm TOV}/\,M_{\odot}$ | $2.39_{-0.28}^{+0.42}$ | $2.23_{-0.21}^{+0.39}$ | $2.23_{-0.21}^{+0.27}$ | $2.25_{-0.22}^{+0.35}$ | $2.25_{-0.22}^{+0.41}$ | $2.22_{-0.20}^{+0.45}$
$R_{1.4}/{\rm km}$ | $11.95_{-0.94}^{+1.04}$ | $11.89_{-0.93}^{+1.15}$ | $11.28_{-0.62}^{+0.59}$ | $11.60_{-0.68}^{+0.81}$ | $12.06_{-0.84}^{+0.89}$ | $12.51_{-0.99}^{+0.65}$
$\mathcal{B}_{\rm NS}^{\rm HS}$ | 1.78 | 2.80 | $4.48\times 10^{4}$ | 95.20 | 3.72 | 1.40
Figure 5: The 90% credible regions of the EOSs for QMF+CSS. The dashed lines
show the results from joint analysis of the present data, including two
gravitational wave data and two NICER’s mass-radius measurements. The red
shaded regions show the updated posterior distributions from four tests
(i-iv), in which we add a mocked MOI measurement for PSR J0737-3039 A to be
$I/M_{\rm A}^{3/2}=35,40,45,50\,{\rm km^{2}}/\,M_{\odot}^{1/2}$, respectively.
Figure 6: Same as Fig.5, but for the mass-radius relations.

###  4.3 Updated constraints with mocked MOI measurements of PSR J0737-3039 A

The MOI of A is expected to be measured with 11 percent precision at a 68%
confidence level in 2030 (Hu et al., 2020), and the results will provide
further constraints on the EOS and neutron star properties. To determine how
such a precise MOI measurement would improve our present knowledge of the
mass-radius relation, we carry out four tests (Tests i-iv) incorporating
mocked data for future measurements of the MOI for A into our Bayesian
analysis framework.

For Test i to Test iv, we add MOI measurements scaled by $M_{\rm A}^{3/2}$ as
$I/M_{\rm A}^{3/2}=35,40,45,50\,{\rm km^{2}}/\,M_{\odot}^{1/2}$ (or
equivalently, $I=1.08,1.23,1.39,1.54\times 10^{45}\,{\rm g\ cm^{2}}$), with a
10% accuracy for all four cases. We repeat the sampling process for the
neutron stars and report the results in Table 2, that is, the sound speed
squared in quark matter $c_{\rm QM}^{2}$, the maximum mass $M_{\rm TOV}$ and
the radius of a typical $1.4\,M_{\odot}$ star $R_{1.4}$ for hybrid stars,
along with results from Li et al. (2021a). The 90% credible regions of the
EOSs are shown in Fig. 5 with and without the mocked MOI measurement for the
representative QMF+CSS case. Fig. 6 shows the corresponding results for the
mass-radius relations.

We found in Li et al. (2021a) that the maximum mass of neutron stars is
$M_{\rm TOV}=2.36_{-0.26}^{+0.49}\,M_{\odot}$
($2.39_{-0.28}^{+0.42}\,M_{\odot}$) and the radius of a $1.4\,M_{\odot}$ star
is $R_{1.4}=11.70_{-0.74}^{+0.85}\,{\rm km}$ ($11.95_{-0.94}^{+1.04}\,{\rm
km}$) at 90% credible interval for QMF+CSS (DD2+CSS) by incorporating both the
GW170817 data and the NICER mass-radius measurement of PSR J0030+0451, as well
as the $2.14\,M_{\odot}$ mass measurement of PSR J0740+6620 from radio timing.
In the present study, the peak value of the maximum mass shifts from $\sim
2.4\,M_{\odot}$ to $\sim 2.2\,M_{\odot}$ by updating the PSR J0740+6620 data
with a more restrictive mass-radius measurement from the NICER x XMM data
111Apart from updating the PSR J0740+6620 data, there are two other
differences between this study and Li et al. (2021a): 1) the likelihood of
encapsulating gravitational wave data obtained using the interpolation
function and 2) additional GW190425 data. For the former, we checked that the
results obtained using the interpolation function are consistent with those
obtained using gravitational wave data, whereas for the latter, we proved that
the addition of GW190425 data had a small effect on the results. Therefore, we
consider the main difference between this study and Li et al. (2021a) to be
the PSR J0740+6620 data update., indicating that a softer EOS is required.
This result can also be directly inferred from the posterior distribution of
the sound speed squared in quark matter $c_{\rm QM}^{2}$, whose median shifts
from $\sim 0.8$ to $\sim 0.7$. See e.g., Xia et al. (2021) for more
discussions on the sound speed. As the resulting maximum mass is insensitive
to the hadronic EOS (Li et al., 2021a), we use the two posterior sample sets
(QMF+CSS and DD2+CSS) to re-generate two equal size sample sets and combined
them as the weighted-average posterior samples. We finally find $M_{\rm
TOV}=2.23_{-0.21}^{+0.41}\,M_{\odot}$ at 90% credible interval. Such results
suggest that the O3b compact object GW200210_092254 (Abbott et al., 2021),
with a mass around $2.8\,M_{\odot}$, should be a black hole rather than a
neutron star. There is also a slight decrease in the radius due to this
update. The inferred radius changes from $11.70_{-0.74}^{+0.85}\,{\rm km}$
($11.95_{-0.94}^{+1.04}\,{\rm km}$) to $11.61_{-0.76}^{+0.75}\,{\rm km}$
($11.89_{-0.93}^{+1.15}\,{\rm km}$) for QMF+CSS (DD2+CSS) at 90% credible
interval.

Next, we determine the effects of incorporating a mocked MOI measurement into
the data. Fig. 5 shows that the additional MOI measurement provides further
constraints on the EOS. These constraints on the EOS mainly occur at
intermediate densities $\sim 300-600\,{\rm MeV/fm^{3}}$ (or equivalently,
$\sim 2-4~{}\varepsilon_{0}$, where $\varepsilon_{0}\approx 157\,{\rm
MeV/fm^{3}}$ is the nuclear saturation mass density). Previously, Miller,
Chirenti, & Lamb (2020) varied the uncertainty of the mocked MOI measurements
(from 5% to 50%) but fixed the central value to the currently most favoured
value of MOI, and found a 10% precision MOI measurement could further improve
the estimation on the radii for neutron stars with masses of
$M=1.8-2\,M_{\odot}$ but had a negligible effect on the corresponding
estimates for masses of approximately $M=1\,M_{\odot}$. Their result is
consistent with our present result for Test ii (shown in the upper right panel
of Fig. 6), since the mocked MOI in the Test ii case is close to the most
probable value obtained from the present analysis, $I/M_{\rm A}^{3/2}\sim
41~{}{\rm km}^{2}/\,M_{\odot}^{1/2}$ (as reported earlier in Sec. 4.1).
Presently, we varied the central value of the MOI measurement but fixed the
uncertainty to 10%. The results of the other three tests based on different
MOI measurements show that it is possible to improve our estimates of neutron
star radii over the entire mass span.

Additionally, as detailed in Table 2, as the measured MOI increases, the
inferred sound speed in quark matter decreases. Taking the QMF+CSS case as an
example, $c_{\rm QM}^{2}$ drops from $0.77_{-0.27}^{+0.21}$ at $I/M_{\rm
A}^{3/2}=35\,{\rm km^{2}/\,M_{\odot}^{1/2}}$ (Test i) to
$0.70_{-0.30}^{+0.27}$ for a measured $I/M_{\rm A}^{3/2}=50\,{\rm
km^{2}}/\,M_{\odot}^{1/2}$ (Test iv) because larger stellar radii have softer
quark matter EOS (i.e., smaller sound speeds). Finally, it is important to
note that the MOI measurement has a limited effect on the resulting maximum
mass. The inferred maximum mass is found to be insensitive to the MOI
measurement, with the most likely value of the maximum mass remaining at $\sim
2.2\,M_{\odot}$.

###  4.4 With or without deconfined quarks?

In this section, we analyze whether the current data favour a strong first-
order phase transition inside neutron star inner cores. For this purpose, we
compute the Bayes factor between the EOS models with and without a phase
transition, $\mathcal{B}_{\rm NS}^{\rm HS}$. The results are reported in the
ninth and sixteenth rows of Table 2.

In general, the Bayes evidences are comparable for both soft QMF and stiff DD2
cases, indicating that the current data are compatible with the possibility of
a phase transition, namely the compact objects are likely hybrid stars,
besides the pure hadronic nature as commonly assumed. The results are
nevertheless slightly dependent on the two hadronic EOSs employed. Following
the quantifying of the significance in Kasen & Raftery (1995), when only
considering available observational data, the possibility of a phase
transition is possible for the stiff DD2 case (with $x\equiv
2\ln(\mathcal{B}_{\rm NS}^{\rm HS})=2.05$), but not supported in the soft QMF
case (with $x<0$). The evidence of a phase transition for the DD2 case is
strengthened if we include the MOI measurement in the data: Table 2 shows the
large Bayes factors, especially for relatively low MOI measurements with $x>2$
for Test i, ii and iii. Note that the recent PREX-II measurement of the
neutron skin of 208Pb indicates a rather stiff symmetry energy (Reed et al.,
2021), resulting in a large predicted radius for a $1.4\,M_{\odot}$ star,
which may conflict with the tidal deformability constraint from GW170817
(Abbott et al., 2018). An early phase transition, e.g., at $\sim 2n_{0}$ (Li
et al., 2021a), could resolve this conflict, because the transition could lead
to a lower radius from the softening of the EOS at high densities. Supporting
evidence for a phase transition also justifies our introduction of a quark
core into the present analysis.

###  4.5 (Hybrid) neutron stars or quark stars?

There is a theoretical possibility that three-flavor quark matter may be even
strongly bound with respect to the ${}^{56}\rm Fe$, namely the energy per
baryon $(E/A)_{\rm uds}$ is smaller than $930\;\text{MeV}$, as proposed by
(Bodmer, 1971; Terazawa, 1979; Witten, 1984). 222See Yuan et al. (2022) for
recent discussions on the stability of both strange and non-strange quark
matter. Since it has so far proved very difficult to unambiguously distinguish
quark stars from normal neutron stars from astrophysical observations
(Lattimer & Prakash, 2007), all pulsar-like objects should be quark stars, if
not all of them. Following similar procedure for the (hybrid)neutron star
analysis in Sec. 3, we consider different structural model for A in the
present section and redo the bayesian analysis from quark matter EOSs based on
the available gravitational-wave and x-ray observational data. Also, the
crustal effects on the mass-radius relations of quark stars are found to be at
most 2-3% (Li et al., 2021b), we, therefore, consider the quark stars to be
bare ones in the present analysis. Nevertheless, it is not clear at all
whether quark stars have a crust or not. See discussion in Weber (2005).

The quark star matter is composed of u, d, s quarks with the charge neutrality
maintained by the inclusion of leptons. Since the QCD equation of motion can
not be solved on the lattice at finite baryon number density, we use
perturbative QCD (which is applicable above $40n_{0}$) and an MIT-type bag
model for the EOS of quark star matter, with the possible existence of
superfluidity for quarks being included (Alford et al., 2001; Rajagopal &
Wilczek, 2001; Lugones & Horvath, 2002). Then three EOS parameters appear in
the likelihood of Eq. (2) in the quark star case as
$\boldsymbol{\theta}=\\{B_{\rm eff},a_{4},\Delta\\}$: the bag parameter
$B_{\rm eff}$, the QCD correction coefficient $a_{4}$ and the quark pairing
gap $\Delta$. Following our previous studies (Zhou, Zhou, & Li, 2018; Miao et
al., 2021; Li et al., 2021b), we assign wide boundaries to them as $B_{\rm
eff}^{1/4}\in[125,150]\,{\rm MeV}$, $\Delta\in[0,100]\,{\rm MeV}$,
$a_{4}\in[0.4,1]$ as theoretically estimated. We mention here that additional
stability requirements should be imposed on the EOS parameters of quark stars:
First, the energy per baryon for non-strange quark matter should satisfy
$(E/A)_{\rm ud}\geq 934\,\;\text{MeV}$ to guarantee the observed stability of
atomic nuclei; Second, $(E/A)_{\rm uds}\leq 930\,\;\text{MeV}$ is required,
according to the hypothesis that strange quark matter is absolutely stable.
The most probable values of the EOS parameters (with their $90\%$ confidence
boundaries) are reported in Table 3, based on the joint analysis of two
gravitational wave data and two NICER’s mass-radius measurements.

Table 3: Most probable intervals of the EOS parameters ($90\%$ confidence level) for quark stars. See text for details.  Parameters | GW170817+GW190425+PSR J0030+ PSR J0740
---|---
$B_{\rm eff}^{1/4}/\rm MeV$ | $144.4_{-7.6}^{+4.6}$
$a_{4}$ | $0.73_{-0.21}^{+0.22}$
$\Delta/\rm MeV$ | $71.2_{-52.4}^{+24.5}$
Figure 7: Posterior distributions of $I$ vs. $R$ for a $M=1.338\,M_{\odot}$ star, transformed from the EOS parameters posteriors. The probability density of $R$ and $I$ are shown in the upper and right panel, respectively, with the 68% credible interval in dashed lines and the 95% credible interval in shade regions. Two markers denote the two results for pure neutron stars.  Table 4: Fit parameters of the $\bar{I}-C$ relation for hybrid stars and quark stars in the functional form of $\bar{I}=a_{1}C^{-1}+a_{2}C^{-2}+a_{3}C^{-3}+a_{4}C^{-4}$ (Breu & Rezzolla, 2016). The coefficients of determination are also listed in the last column.  | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $R^{2}$
---|---|---|---|---|---
Hybrid star | $6.759\times 10^{-1}$ | $3.356\times 10^{-1}$ | $-1.820\times 10^{-2}$ | $5.514\times 10^{-4}$ | 0.99757
Quark star | $3.461\times 10^{-1}$ | $4.013\times 10^{-1}$ | $-1.198\times 10^{-3}$ | $6.499\times 10^{-5}$ | 0.99999
Figure 8: Posterior distribution of $\bar{I}\equiv I/M^{3}$ shown as a
function of stellar compactness $C\equiv M/R$, transformed from the EOS
parameters posteriors. The solid and dashed lines represent the fitting
functions for the (hybrid-)neutron star (with both QMF+CSS and DD2+CSS cases)
and quark star results, respectively.

In Fig. 7, we show the posterior distributions of $I$ as functions of $R$ for
both quark stars and (hybrid) neutron stars. One observes that the two $I$
vs.$R$ distributions are located in separate regions and become more and more
separated with each other when increasing $R$ (equivalently decreasing the
compactness $M/R$), suggesting that it can be possible to distinguish them
from simultaneous measurement of the stellar radius and the MOI. In
particular, we calculate the MOI of A as a quark star to be
$I_{A}=1.55_{-0.11}^{+0.13}\times 10^{45}\,{\rm g\,cm^{2}}$, at 90% credible
interval. Such predicted results based on quark star EOSs are obviously larger
than those based on neutrons star EOSs. Consequently, If the MOI is measured
large for A, for example, $I_{A}\gtrsim 1.4\times 10^{45}\,{\rm g\,cm^{2}}$,
it is most likely a quark star rather than a neutron star with or without a
quark core (see Table 1), provided that the accuracy of the radius measurement
is at least $\sim 1\,\hbox{$\,{\rm km}$}$. Also, there is a good correlation
between the radius and the MOI of quark stars, which can be fitted by
${I}/({10^{45}\,{\rm g\,cm^{2}}})=aR^{b}$ with $a=1.69\times 10^{-2}$,
$b=1.87$ and the coefficient of determination is $R^{2}=0.997$. The fitted
formula allows us to translate accurately from a MOI measurement to its radius
and vice versa, under the condition of A being a quark star.

Before the end of the section, we come back to the universal relations
mentioned at the beginning of Sec. 4. It is expected that the dimensionless
MOI $\bar{I}$ of quark stars follow the similar universal relations with
neutron stars and hybrid stars when expressed as a function of $\Lambda$, as
shown in, e.g., Yagi & Yunes (2017). Nevertheless, the $\bar{I}$ values of
quark stars should behave quite differently from those of neutron stars and
hybrid stars when expressed as a function of the stellar compactness $C\equiv
M/R$ (Ravenhall & Pethick, 1994; Bejger & Haensel, 2002; Lattimer & Schutz,
2005). We show in Fig. 8 the dimensionless MOI $\bar{I}$ as a function of $C$
for both (hybird-)neutron stars and quark stars. We see that the neutron star
and quark star results are both well distributed in universal
$\bar{I}$-versus-$C$ curves in the form of
$\bar{I}=a_{1}C^{-1}+a_{2}C^{-2}+a_{3}C^{-3}+a_{4}C^{-4}$ (Breu & Rezzolla,
2016) (the fit parameters are collected in Table 4). The two universal lines
reconcile in the high-compactness regime, but it is evident that the
(hybrid-)neuron star line deviates from the quark star one and moves downwards
in the low-compactness regime. This is because the matter in quark stars is
relatively evenly distributed; The surface density is actually close to the
central density, resulting in a larger MOI than in the (hybrid-)neutron star
case when decreasing the stellar mass or compactness (Yagi & Yunes, 2013).

##  5 Summary and conclusions

There is a high expectation of an imminent precise measurement of the neutron
star MOI for the very first time for the first double pulsar system PSR
J0737-3039. Continuing efforts are being made to theoretically predict the MOI
of PSR J0737-3039 A based on many modern nuclear EOSs within the constraints
provided by the results of nuclear experiments and astrophysical observations.
Considering that the pulsar interior may have a quark matter core, it is
meaningful to update the MOI analysis of previous studies by incorporating
strangeness phase transitions. Furthermore, such an MOI measurement would
place new constraints on the neutron star EOS. A question that naturally
follows is what information this measurement will contribute to our
understanding of dense matter.

In this study, a hadron-quark phase transition in the EOS prior modelled by
the generic CSS parameterization is explicitly considered and a Bayesian
analysis is performed on the MOI of PSR J0737-3039 A using available data from
LIGO/Virgo and NICER. Both soft QMF and stiff DD2 are applied for the low-
density hadronic phase to test the model-dependence of the MOI prediction. A
MOI of $1.27_{-0.14}^{+0.18}\times 10^{45}\,{\rm g\,cm^{2}}$ is obtained using
QMF+CSS, and a similar result of $1.29_{-0.15}^{+0.26}\times 10^{45}\,{\rm
g\,cm^{2}}$ is obtained using DD2+CSS. We also demonstrated in detail that an
accurate MOI measurement, when available, can significantly reduce the EOS
uncertainties, especially when the MOI value is a few times the nuclear
saturation density. In addition, we evaluated the relative probability of a
phase transition occurring using the pure hadronic star as a reference state,
based on the available data and a potential MOI constraint. The results show
that the hybrid star scenario is justified as one of the many possibilities of
the nature of pulsars, and the maximum mass is predicted to be $M_{\rm
TOV}\sim 2.2\,M_{\odot}$. On the other hand, a relatively larger MOI result is
found, $I_{A}\sim 1.55\times 10^{45}\,{\rm g\,cm^{2}}$, if assuming PSR
J0737-3039 A is a quark star, compared to that of (hybrid-)neutron-star
modelling, suggesting it is still possible to distinguish between the two if
more precise radius measurement is possible (at least $\sim 1\,\hbox{$\,{\rm
km}$}$ accuracy), possible in the near future with joint efforts of the
multimessenger missions, like aLIGO/Virgo, NICER, eXTP, SKA. We finally
examine the MOI–compactness universal relations and provide easy-to-use
analytic fits on the results for both (hybrid-)neutron-star and quark star
cases.

## ACKNOWLEDGMENTS

We are thankful to Prof. B.-A. Li and Tong Liu for helpful discussions. The
work is supported by National SKA Program of China (No. 2020SKA0120300), the
National Key Research and Development Program of China (No. 2017YFA0402600),
the National Natural Science Foundation of China (Nos. 11873040 and 11833003),
the science research grants from the China Manned Space Project (No. CMS-
CSST-2021-B11) and the Youth Innovation Fund of Xiamen (No. 3502Z20206061).

## Data Availability

The data underlying this article are available in the article.

## References

  * Abbott et al. (2017) Abbott B. P., Abbott R., Abbott T. D., Acernese F., Ackley K., Adams C., Adams T., et al., 2017, PhRvL, 119, 161101. doi:10.1103/PhysRevLett.119.161101
  * Abbott et al. (2018) Abbott B. P., Abbott R., Abbott T. D., Acernese F., Ackley K., Adams C., Adams T., et al., 2018, PhRvL, 121, 161101. doi:10.1103/PhysRevLett.121.161101
  * Abbott et al. (2020) Abbott B. P., Abbott R., Abbott T. D., Abraham S., Acernese F., Ackley K., Adams C., et al., 2020, ApJL, 892, L3. doi:10.3847/2041-8213/ab75f5
  * Abbott et al. (2021) Abbott R., Abbott T. D., Acernese F., Ackley K., et al., 2021, arXiv, arXiv:2111.03606
  * Al-Mamun et al. (2021) Al-Mamun M., Steiner A. W., Nättilä J., Lange J., O’Shaughnessy R., Tews I., Gandolfi S., et al., 2021, PhRvL, 126, 061101. doi:10.1103/PhysRevLett.126.061101
  * Alford, Han, & Prakash (2013) Alford M. G., Han S., Prakash M., 2013, PhRvD, 88, 083013. doi:10.1103/PhysRevD.88.083013
  * Alford et al. (2001) Alford M., Rajagopal K., Reddy S., Wilczek F., 2001, PhRvD, 64, 074017. doi:10.1103/PhysRevD.64.074017
  * Antoniadis et al. (2013) Antoniadis J., Freire P. C. C., Wex N., Tauris T. M., Lynch R. S., van Kerkwijk M. H., Kramer M., et al., 2013, Sci, 340, 448. doi:10.1126/science.1233232
  * Bandyopadhyay et al. (2018) Bandyopadhyay D., Bhat S. A., Char P., Chatterjee D., 2018, EPJA, 54, 26. doi:10.1140/epja/i2018-12456-y
  * Baym, Pethick, & Sutherland (1971) Baym G., Pethick C., Sutherland P., 1971, ApJ, 170, 299. doi:10.1086/151216
  * Becker et al. (2000) Becker W., Trümper J., Lommen A. N., Backer D. C., 2000, ApJ, 545, 1015. doi:10.1086/317839
  * Bejger, Bulik, & Haensel (2005) Bejger M., Bulik T., Haensel P., 2005, MNRAS, 364, 635. doi:10.1111/j.1365-2966.2005.09575.x
  * Bejger & Haensel (2002) Bejger M., Haensel P., 2002, A&A, 396, 917. doi:10.1051/0004-6361:20021241
  * Berezhiani et al. (2003) Berezhiani Z., Bombaci I., Drago A., Frontera F., Lavagno A., 2003, ApJ, 586, 1250. doi:10.1086/367756
  * Bodmer (1971) Bodmer A. R., 1971, PhRvD, 4, 1601. doi:10.1103/PhysRevD.4.1601
  * Bombaci et al. (2016) Bombaci I., Logoteta D., Vidaña I., Providência C., 2016, EPJA, 52, 58. doi:10.1140/epja/i2016-16058-5
  * Breu & Rezzolla (2016) Breu C., Rezzolla L., 2016, MNRAS, 459, 646. doi:10.1093/mnras/stw575
  * Buchner et al. (2014) Buchner J., Georgakakis A., Nandra K., Hsu L., Rangel C., Brightman M., Merloni A., et al., 2014, A&A, 564, A125. doi:10.1051/0004-6361/201322971
  * Burgay et al. (2003) Burgay M., D’Amico N., Possenti A., Manchester R. N., Lyne A. G., Joshi B. C., McLaughlin M. A., et al., 2003, Natur, 426, 531. doi:10.1038/nature02124
  * Cheng & Dai (1996) Cheng K. S., Dai Z. G., 1996, PhRvL, 77, 1210. doi:10.1103/PhysRevLett.77.1210
  * Cromartie et al. (2020) Cromartie H. T., Fonseca E., Ransom S. M., Demorest P. B., Arzoumanian Z., Blumer H., Brook P. R., et al., 2020, NatAs, 4, 72. doi:10.1038/s41550-019-0880-2
  * Damour & Schafer (1988) Damour T., Schafer G., 1988, NCimB, 101B, 127. doi:10.1007/BF02828697
  * Demorest et al. (2010) Demorest P. B., Pennucci T., Ransom S. M., Roberts M. S. E., Hessels J. W. T., 2010, Natur, 467, 1081. doi:10.1038/nature09466
  * Fonseca et al. (2021) Fonseca E., Cromartie H. T., Pennucci T. T., Ray P. S., Kirichenko A. Y., Ransom S. M., Demorest P. B., et al., 2021, ApJL, 915, L12. doi:10.3847/2041-8213/ac03b8
  * Fortin et al. (2016) Fortin M., Providência C., Raduta A. R., Gulminelli F., Zdunik J. L., Haensel P., Bejger M., 2016, PhRvC, 94, 035804. doi:10.1103/PhysRevC.94.035804
  * Fortin et al. (2020) Fortin M., Raduta A. R., Avancini S., Providência C., 2020, PhRvD, 101, 034017. doi:10.1103/PhysRevD.101.034017
  * Gorda et al. (2021) Gorda T., Kurkela A., Paatelainen R., Säppi S., Vuorinen A., 2021, PhRvL, 127, 162003. doi:10.1103/PhysRevLett.127.162003
  * Greif et al. (2020) Greif S. K., Hebeler K., Lattimer J. M., Pethick C. J., Schwenk A., 2020, ApJ, 901, 155. doi:10.3847/1538-4357/abaf55
  * Hartle (1967) Hartle J. B., 1967, ApJ, 150, 1005. doi:10.1086/149400
  * Hernandez Vivanco et al. (2020) Hernandez Vivanco F., Smith R., Thrane E., Lasky P. D., 2020, MNRAS, 499, 5972. doi:10.1093/mnras/staa3243
  * Herzog & Röpke (2011) Herzog M., Röpke F. K., 2011, PhRvD, 84, 083002. doi:10.1103/PhysRevD.84.08300
  * Hu et al. (2020) Hu H., Kramer M., Wex N., Champion D. J., Kehl M. S., 2020, MNRAS, 497, 3118. doi:10.1093/mnras/staa2107
  * Jiang et al. (2020) Jiang J.-L., Tang S.-P., Wang Y.-Z., Fan Y.-Z., Wei D.-M., 2020, ApJ, 892, 55. doi:10.3847/1538-4357/ab77cf
  * Klähn et al. (2007) Klähn T., Blaschke D., Sandin F., Fuchs C., Faessler A., Grigorian H., Röpke G., et al., 2007, PhLB, 654, 170. doi:10.1016/j.physletb.2007.08.048
  * Kramer et al. (2021) Kramer M., Stairs I. H., Manchester R. N., Wex N., Deller A. T., Coles W. A., Ali M., et al., 2021, PhRvX, 11, 041050. doi:10.1103/PhysRevX.11.041050
  * Kumar & Landry (2019) Kumar B., Landry P., 2019, PhRvD, 99, 123026. doi:10.1103/PhysRevD.99.123026
  * Landry & Kumar (2018) Landry P., Kumar B., 2018, ApJL, 868, L22. doi:10.3847/2041-8213/aaee76
  * Lattimer & Prakash (2007) Lattimer J. M., Prakash M., 2007, PhR, 442, 109. doi:10.1016/j.physrep.2007.02.003
  * Lattimer & Schutz (2005) Lattimer J. M., Schutz B. F., 2005, ApJ, 629, 979. doi:10.1086/431543
  * Li (2015) Li A., 2015, ChPhL, 32, 079701. doi:10.1088/0256-307X/32/7/079701
  * Li et al. (2016) Li A., Dong J. M., Wang J. B., Xu R. X., 2016, ApJS, 223, 16. doi:10.3847/0067-0049/223/1/16
  * Li et al. (2021a) Li A., Miao Z., Han S., Zhang B., 2021, ApJ, 913, 27. doi:10.3847/1538-4357/abf355
  * Li et al. (2021b) Li A., Miao Z.-Q., Jiang J.-L., Tang S.-P., Xu R.-X., 2021, MNRAS, 506, 5916. doi:10.1093/mnras/stab2029
  * Li et al. (2020) Li A., Zhu Z.-Y., Zhou E.-P., Dong J.-M., Hu J.-N., Xia C.-J., 2020, JHEAp, 28, 19. doi:10.1016/j.jheap.2020.07.001
  * Lim, Holt, & Stahulak (2019) Lim Y., Holt J. W., Stahulak R. J., 2019, PhRvC, 100, 035802. doi:10.1103/PhysRevC.100.035802
  * Lommen et al. (2000) Lommen A. N., Zepka A., Backer D. C., McLaughlin M., Cordes J. M., Arzoumanian Z., Xilouris K., 2000, ApJ, 545, 1007. doi:10.1086/317841
  * Lourenço et al. (2020) Lourenço O., Dutra M., Lenzi C. H., Biswal S. K., Bhuyan M., Menezes D. P., 2020, EPJA, 56, 32. doi:10.1140/epja/s10050-020-00040-z
  * Lugones & Horvath (2002) Lugones G., Horvath J. E., 2002, PhRvD, 66, 074017. doi:10.1103/PhysRevD.66.074017
  * Lyne et al. (2004) Lyne A. G., Burgay M., Kramer M., Possenti A., Manchester R. N., Camilo F., McLaughlin M. A., et al., 2004, Sci, 303, 1153. doi:10.1126/science.1094645
  * Margaritis, Koliogiannis, & Moustakidis (2020) Margaritis C., Koliogiannis P. S., Moustakidis C. C., 2020, PhRvD, 101, 043023. doi:10.1103/PhysRevD.101.043023
  * Miao et al. (2021) Miao Z., Jiang J.-L., Li A., Chen L.-W., 2021, ApJL, 917, L22. doi:10.3847/2041-8213/ac194d
  * Miao et al. (2020) Miao Z., Li A., Zhu Z., Han S., 2020, ApJ, 904, 103. doi:10.3847/1538-4357/abbd41
  * Miller, Chirenti, & Lamb (2020) Miller M. C., Chirenti C., Lamb F. K., 2020, ApJ, 888, 12. doi:10.3847/1538-4357/ab4ef9
  * Miller et al. (2019) Miller M. C., Lamb F. K., Dittmann A. J., Bogdanov S., Arzoumanian Z., Gendreau K. C., Guillot S., et al., 2019, ApJL, 887, L24. doi:10.3847/2041-8213/ab50c5
  * Miller et al. (2021) Miller M. C., Lamb F. K., Dittmann A. J., Bogdanov S., Arzoumanian Z., Gendreau K. C., Guillot S., et al., 2021, ApJL, 918, L28. doi:10.3847/2041-8213/ac089b
  * Morrison et al. (2004) Morrison I. A., Baumgarte T. W., Shapiro S. L., Pandharipande V. R., 2004, ApJL, 617, L135. doi:10.1086/427235
  * Motta et al. (2019) Motta T. F., Kalaitzis A. M., Antić S., Guichon P. A. M., Stone J. R., Thomas A. W., 2019, ApJ, 878, 159. doi:10.3847/1538-4357/ab218e
  * Olesen & Madsen (1994) Olesen M. L., Madsen J., 1994, PhRvD, 49, 2698. doi:10.1103/PhysRevD.49.2698
  * Olinto (1987) Olinto A. V., 1987, PhLB, 192, 71. doi:10.1016/0370-2693(87)91144-0
  * Ouyed, Dey, & Dey (2002) Ouyed R., Dey J., Dey M., 2002, A&A, 390, L39. doi:10.1051/0004-6361:20020982
  * Pagliara, Herzog, & Röpke (2013) Pagliara G., Herzog M., Röpke F. K., 2013, PhRvD, 87, 103007. doi:10.1103/PhysRevD.87.103007
  * Pérez-García, Daigne, & Silk (2013) Pérez-García M. A., Daigne F., Silk J., 2013, ApJ, 768, 145. doi:10.1088/0004-637X/768/2/145
  * Pérez-García Silk, & Stone (2010) Pérez-García M. A., Silk J., Stone J. R., 2010, PhRvL, 105, 141101. doi:10.1103/PhysRevLett.105.141101
  * Raduta, Oertel, & Sedrakian (2020) Raduta A. R., Oertel M., Sedrakian A., 2020, MNRAS, 499, 914. doi:10.1093/mnras/staa2491
  * Rajagopal & Wilczek (2001) Rajagopal K., Wilczek F., 2001, PhRvL, 86, 3492. doi:10.1103/PhysRevLett.86.3492
  * Raithel, Özel, & Psaltis (2016) Raithel C. A., Özel F., Psaltis D., 2016, PhRvC, 93, 032801. doi:10.1103/PhysRevC.93.032801
  * Ravenhall & Pethick (1994) Ravenhall D. G., Pethick C. J., 1994, ApJ, 424, 846. doi:10.1086/173935
  * Reed et al. (2021) Reed B. T., Fattoyev F. J., Horowitz C. J., Piekarewicz J., 2021, PhRvL, 126, 172503. doi:10.1103/PhysRevLett.126.172503
  * Read et al. (2009) Read J. S., Lackey B. D., Owen B. J., Friedman J. L., 2009, PhRvD, 79, 124032. doi:10.1103/PhysRevD.79.124032
  * Ribes et al. (2019) Ribes P., Ramos A., Tolos L., Gonzalez-Boquera C., Centelles M., 2019, ApJ, 883, 168. doi:10.3847/1538-4357/ab3a93
  * Rikovska Stone et al. (2007) Rikovska Stone J., Guichon P. A. M., Matevosyan H. H., Thomas A. W., 2007, NuPhA, 792, 341. doi:10.1016/j.nuclphysa.2007.05.011
  * Riley et al. (2019a) Riley T. E., Watts A. L., Bogdanov S., Ray P. S., Ludlam R. M., Guillot S., Arzoumanian Z., et al., 2019, ApJL, 887, L21. doi:10.3847/2041-8213/ab481c
  * Riley et al. (2019b) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, v1.0.0, Zenodo, doi:10.5281/zenodo.3386449
  * Riley et al. (2021a) Riley T. E., Watts A. L., Ray P. S., Bogdanov S., Guillot S., Morsink S. M., Bilous A. V., et al., 2021, ApJL, 918, L27. doi:10.3847/2041-8213/ac0a81
  * Riley et al. (2021b) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, v1.0.0, Zenodo, doi:10.5281/zenodo.4697625
  * Sen & Jha (2019) Sen D., Jha T. K., 2019, JPhG, 46, 015202. doi:10.1088/1361-6471/aaf0b0
  * Silva et al. (2021) Silva H. O., Holgado A. M., Cárdenas-Avendaño A., Yunes N., 2021, PhRvL, 126, 181101. doi:10.1103/PhysRevLett.126.181101
  * Steiner et al. (2015) Steiner A. W., Gandolfi S., Fattoyev F. J., Newton W. G., 2015, PhRvC, 91, 015804. doi:10.1103/PhysRevC.91.015804
  * Stovall et al. (2014) Stovall K., Lynch R. S., Ransom S. M., Archibald A. M., Banaszak S., Biwer C. M., Boyles J., et al., 2014, ApJ, 791, 67. doi:10.1088/0004-637X/791/1/67
  * Sun, Liu, & Xing (2019) Sun B. Y., Liu Z. W., Xing R. Y., 2019, AIPC, 2127, 020020. doi:10.1063/1.5117810
  * Terazawa (1979) Terazawa, H., INS, University of Tokyo Report No. INSReport-336, 1979
  * Urbanec, Miller, & Stuchlík (2013) Urbanec M., Miller J. C., Stuchlík Z., 2013, MNRAS, 433, 1903. doi:10.1093/mnras/stt858
  * Weber (2005) Weber F., 2005, PrPNP, 54, 193. doi:10.1016/j.ppnp.2004.07.001
  * Wei et al. (2019) Wei J.-B., Figura A., Burgio G. F., Chen H., Schulze H.-J., 2019, JPhG, 46, 034001. doi:10.1088/1361-6471/aaf95c
  * Witten (1984) Witten E., 1984, PhRvD, 30, 272. doi:10.1103/PhysRevD.30.272
  * Worley, Krastev, & Li (2008) Worley A., Krastev P. G., Li B.-A., 2008, ApJ, 685, 390. doi:10.1086/589823
  * Xia et al. (2021) Xia C., Zhu Z., Zhou X., Li A., 2021, ChPhC, 45, 055104. doi:10.1088/1674-1137/abea0d
  * Xie & Li (2021) Xie W.-J., Li B.-A., 2021, PhRvC, 103, 035802. doi:10.1103/PhysRevC.103.035802
  * Yagi & Yunes (2013) Yagi K., Yunes N., 2013, PhRvD, 88, 023009. doi:10.1103/PhysRevD.88.023009
  * Yagi & Yunes (2017) Yagi K., Yunes N., 2017, PhR, 681, 1. doi:10.1016/j.physrep.2017.03.002
  * Yuan et al. (2022) Yuan W.-L., Li A., Miao Z., Zuo B., Bai Z., 2022, PhRvD, 105, 123004. doi:10.1103/PhysRevD.105.123004
  * Yuan & Zhang (1999) Yuan Y. F., Zhang J. L., 1999, ApJ, 525, 950. doi:10.1086/307921
  * Zhou, Zhou, & Li (2018) Zhou E.-P., Zhou X., Li A., 2018, PhRvD, 97, 083015. doi:10.1103/PhysRevD.97.083015
  * Zhu, Zhou, & Li (2018) Zhu Z.-Y., Zhou E.-P., Li A., 2018, ApJ, 862, 98. doi:10.3847/1538-4357/aacc28
  * Kasen & Raftery (1995) Kass, R.E., & Raftery, A.E. 1995, Journal of the American Statistical Society, 90, 773