# Effect of geometrically thin discs on precessing, thick flows: Relevance to
type-C QPOs

D. A. Bollimpalli1,2,3, P. C. Fragile3,4, W. Kluźniak2
1Max-Planck-Institut für Astrophysik, 85741 Garching b. München, Germany
2Nicolaus Copernicus Astronomical Center, ul. Bartycka 18, PL 00-716 Warsaw,
Poland
3Department of Physics and Astronomy, College of Charleston, Charleston, SC
29424, USA
4Kavli Institute for Theoretical Physics, Kohn Hall, University of California,
Santa Barbara, CA 93107, USA Contact e-mail<EMAIL_ADDRESS>
###### Abstract

Type-C quasi-periodic oscillations (QPOs) are the low-frequency QPOs most
commonly observed during the hard spectral state of X-ray binary systems. The
leading model for these QPOs is the Lense-Thirring precession of a hot,
geometrically thick accretion flow that is misaligned with respect to the
black hole spin axis. However, none of the work done to date has accounted for
the effects of a surrounding, geometrically thin disc on this precession, as
would be the case in the truncated disc picture of the hard state. To address
this, we perform a set of GRMHD simulations of truncated discs misaligned with
the spin axes of their central black holes. Our results confirm that the
inner-hot flow still undergoes precession, though at a rate that is only 5
percent of what is predicted for an isolated, precessing torus. We find that
the exchange of angular momentum between the outer, thin and the inner, thick
disc causes this slow-down in the precession rate and discuss its relevance to
type-C QPOs.

###### keywords:

accretion, accretion discs — magnetohydrodynamics — methods: numerical —
relativistic processes – stars: black holes — X-rays: binaries

††pubyear: 2020††pagerange: Effect of geometrically thin discs on precessing,
thick flows: Relevance to type-C QPOs–A

##  1 Introduction

Many accreting black holes and neutron stars exhibit rapid variability in
their X-ray light curves, which is often quasi-periodic in nature (indicated
by broad peaks of power in the frequency domain). Further, they occasionally
undergo outbursts during which the source transits through several different
spectral states, as it simultaneously varies in luminosity and spectral
hardness, tracing a “q” shape in a hardness-intensity diagram, (Fender et al.,
2004; McClintock & Remillard, 2006). The current best models that provide a
framework for understanding these outbursts relate the changing spectra to a
changing geometry of the accretion flow and the presence or absence of
different physical components (e.g., a corona or jet) (Esin et al., 1997;
Begelman & Armitage, 2014). The X-ray spectra observed in these systems is
mainly composed of a soft, thermal component (below approximately 3 keV) and a
hard power-law component (extending roughly between 10 - 100 keV). The soft
component is thought to be due to thermal, blackbody-like radiation from a
relatively cold, optically thick and geometrically thin disc. The non-thermal,
power-law component is thought to be generated from the inverse-Compton
scattering of soft seed photons emitted from the disc by a radiatively
inefficient cloud of hot electrons, often referred to as the “corona.” The
location and geometry of the corona are still under much debate, with some
models placing it at the base of the jet (Martocchia & Matt, 1996; Fabian et
al., 2009; Parker et al., 2015), while others consider it to sandwich the disc
(Galeev et al., 1979; Haardt & Maraschi, 1993; Svensson, 1996; Beloborodov,
1999). The truncated-disc model envisions a third possible geometry with a
standard-thin disc truncated well outside the last stable orbit, and the
corona as a geometrically thick, radiatively inefficient accretion flow
filling in the inner gap (Eardley et al., 1975; Esin et al., 1997; Liu et al.,
2007). Although all the above mentioned models can reproduce the observed
X-ray spectra, the truncated disc model has the advantage of being able to
simultaneously explain some of the observed variability properties.
Irrespective of the model, a better understanding of the physical origins of
the variability observed in each spectral state can serve as a powerful probe
of the inner regions of accretion discs.

The nature of quasi-periodic oscillations (QPOs) varies from one spectral
state to another. Broadly speaking, the QPOs observed in black hole systems
are classified as high-frequency ($\gtrsim 60$ Hz) and low-frequency
($\lesssim 30$ Hz) (Remillard & McClintock, 2006; Belloni, 2010). Low-
frequency QPOs (LFQPOs) are further classified into type-A, B and C, of which
the first two are observed in the intermediate states, during the transition
from hard to soft spectral states (Motta, 2016). Type-C QPOs are commonly
observed in most accretion states (Motta et al., 2012; Motta, 2016), although
they are particularly prominent in the hard state. For these QPOs, both the
fractional QPO amplitude and the phase-lag between the hard and soft photons
measured at the QPO frequency are found to have an inclination dependence
(Motta et al., 2015; van den Eijnden et al., 2017). Such a dependence strongly
suggests a geometric origin for this QPO. For a detailed review of type-C
QPOs, we refer the reader to Ingram & Motta (2019).

Most of the models invoked to explain type-C QPOs are based either on the
geometry or the intrinsic properties of the disc. Possible mechanisms that
have been invoked include trapped corrugation modes ($c$-modes) in the inner
regions (Kato & Fukue, 1980; Wagoner, 1999; Silbergleit et al., 2001) or
spiral density waves propagating back and forth between the disc inner radius
and the Lindblad radius at some frequency (Tagger & Pellat, 1999; Rodriguez et
al., 2002). Along the same lines, oscillations of a standing shock in the
accretion flow could also produce a type-C QPO; however, to produce the
required frequencies, such a shock would need to be located quite far from the
black hole (Chakrabarti & Titarchuk, 1995). Another popular model is the
relativistic precession model, which associates the type-C QPO with the nodal
precession frequency of a geodesic particle orbit in relativity (Stella &
Vietri, 1998; Stella et al., 1999).

A variation of this model, which seems to match observations of the type-C QPO
particularly well, is, instead of having a particle or infinitesimal ring
precess, have a finite inner region of the accretion flow precess (Ingram et
al., 2009). This model ties in with the truncated disc picture of the hard
accretion state, in that it is the hot, thick flow inside the truncation
radius (i.e., the corona) that is precessing (Ingram & Done, 2011). This has
the advantages of lowering the precession frequency into the proper range,
naturally explaining the rise in the QPO frequency during the early phases of
an outburst cycle (when the truncation radius is moving inward), and easily
reproducing the inclination dependence seen in QPO observations (since
precession is more apparent the closer to the orbital plane one observes).
This model also naturally explains the observed correlation between the QPO
phase and the energy of the iron emission lines seen in phase-resolved
spectroscopic studies (Ingram et al., 2016).

Previous numerical work involving magnetohydrodynamical simulations of
isolated, tilted, thick discs have demonstrated global Lense-Thirring
precession at frequencies matching observations (Fragile et al., 2007; Liska
et al., 2018; White et al., 2019). However, none of those studies consider the
full truncated-disc geometry and especially the effects of the surrounding
geometrically thin disc on the type-C QPO frequencies. The main concern for
QPO models based on the Lense-Thirring precession of the hot, thick flow is
its coupling to the outer thin disc. Questions that need to be answered are:
1) How does the surrounding thin disc alter the precession, both in terms of
its frequency and amplitude; 2) Will the coupling slow down the precession and
eventually stop it or perhaps boost the strength of the QPOs? The main
objective of this paper is to address some of these questions.

For this purpose, we perform general relativistic magnetohydrodynamics (GRMHD)
simulations of tilted, truncated discs around rotating black holes. We evolve
the simulations long enough to establish their steady-state behavior which
amounts to a significant fraction of the disc’s precession period. The
remainder of this paper is organized as follows: In Section 2, we describe the
initial setup of our simulations. In Section 3, we discuss the main result of
our paper, that the outer, thin disc slows down the precession of the inner,
thick disc. We conclude with a discussion in Section 4.

##  2 Initial Setup

###  2.1 Physical setup

We performed three, fully 3D, GRMHD simulations of truncated accretion discs,
all with the angular momentum axes of the disc and the black hole misaligned
by an angle of $\beta_{0}=15^{\circ}$, as listed in Table 1; one is a high-
resolution, 4-level simulation (with dimensionless spin parameter
$a_{*}=a/M=0.9$), and two are 3-level, low-resolution simulations with
$a_{*}=0.9$ and 0.5.

Simulation | Spin | Tilt | Resolution | $t_{\mathrm{end}}$
---|---|---|---|---
| ($a_{*}$) | ($\beta_{0}$) |  | [$r_{\rm g}$$/c$]
a9b15r15L4 | 0.9 | $15^{\circ}$ | 384x256x256 | 25,000
a9b15r15L3 | 0.9 | $15^{\circ}$ | 192x128x128 | 14,850
a5b15r15L3 | 0.5 | $15^{\circ}$ | 192x128x128 | 16,700
Table 1: Simulations presented in this paper.

The initial setup of all the simulations consists of a finite torus,
surrounded by a thin disc. The torus is initialized following the procedure of
Chakrabarti (1985) assuming a constant specific angular momentum, with its
inner edge at $r_{\rm in}=6.5\,$$r_{\rm g}$, pressure maximum at $r_{\rm
cen}=9\,$$r_{\rm g}$, and polytropic index of $\gamma=5/3$, where $r_{\rm
g}=GM/c^{2}$ is the gravitational radius. In all the simulations, a thin slab
with fixed height, $H=0.4$$r_{\rm g}$, surrounds the torus extending from
$15$$r_{\rm g}$to the outer boundary of the simulations domain at $\approx
250$$r_{\rm g}$, with a density distribution

$\rho(R,z)\propto\frac{\exp(-z^{2}/2H^{2})}{\\{1+\exp[(15r_{\rm
g}-R)/H]\\}\\{1+\exp[(R-(40r_{\rm g})^{1.5})/H]\\}}\leavevmode\nobreak\ ,$ (1)

where $R=r\sin\theta$ is the cylindrical radius. The angular momentum in the
thin disc is Keplerian. In order to maintain the desired thin structure, we
use an artificial cooling function, as described in Fragile et al. (2012a),
but only applied to radii beyond $15$$r_{\rm g}$.

Both the torus and the slab are threaded with numerous small poloidal loops of
magnetic field with alternating polarity, starting from a magnetic vector
potential of the form

$A_{\phi}\propto\frac{1}{1+\exp(\delta)}\sqrt{P_{\rm g}}\sin(2\pi
R/5H)\leavevmode\nobreak\ ,$ (2)

where $P_{\rm g}$ is the local gas pressure and
$\delta=10[(z/H)^{2}+H^{2}/(R-r_{\rm ms})^{2}+H^{2}/(40r_{\rm
g}^{1.5}-R)^{2}-1]$, with $H=0.6$ and $0.4$$r_{\rm g}$for the torus and thin
disc regions, respectively, and $r_{\rm ms}$ is the location of the marginally
stable orbit. The resulting magnetic field is normalized such that the ratio
of the gas pressure to the magnetic pressure is initially $\geq 10$ throughout
the torus and thin disc.This field configuration prevents the accumulation of
strong net flux in the inner regions, so we do not expect a magnetically
arrested disc (MAD) to form. The background region surrounding the thin disc
and torus is initialized with a non-magnetic gas with low-density,
$\rho=10^{-5}\rho_{\rm max}r^{-1.5}$ and internal energy density,
$e=10^{-7}e_{\rm max}r^{-2.5}$, where $\rho_{\rm max}$ and $e_{\rm max}$ are
the maximum gas density and internal energy density in the disc/torus.

###  2.2 Numerical setup

All simulations are performed using Cosmos++ (Anninos et al., 2005), with its
five-stage, strong-stability-preserving Runge-Kutta (SSPRK) time-stepping
scheme (Spiteri & Ruuth, 2002) and third-order, piecewise-parabolic spatial
reconstruction (Colella & Woodward, 1984). The fluxes are calculated at zone
interfaces using a two-wave, HLL, approximate Riemann solver (Harten et al.,
1983).

The simulations are performed in modified Kerr-Schild spherical-polar
coordinates, which allow us to place the inner domain boundary inside the
black hole event horizon. The grid is logarithmically spaced in $r$, extending
over $r\in[1.4\,r_{\rm g},\,40^{1.5}\,r_{\rm g}]$, and uniformly spaced in
$\phi$, covering the full $[0,2\pi]$ domain. The polar angle $\theta$ is
uniformly spaced in the coordinate $x_{2}$, ranging from 0 to 1, with

$\theta=\pi x_{2}+\frac{1-h}{2}\sin(2\pi x_{2})\leavevmode\nobreak\ ,$ (3)

where $h=0.5$ is used to concentrate the grid cells close to the midplane
(McKinney, 2006). We cut out a small cone of opening angle $10^{-15}\pi$ near
the poles to avoid calculating metric terms on the pole.

All simulations use a base resolution of $48\times 32\times 32$ with an
additional two or three levels of static mesh refinement added, with adjacent
refinement levels differing in resolution by a factor of $2$ in each
dimension, to bring most of the grid up to a fiducial resolution of $192\times
128\times 128$ or $384\times 256\times 256$ for “3-level" and “4-level"
simulations, respectively. The initial two levels of refinement are used to
resolve both the torus and the thin slab, which is common for both the 3-level
and 4-level simulations. The torus region is refined with two refinement
levels in the regions within $1.4<r/r_{\rm g}<24$ and
$0.14\pi<\theta<0.86\pi$, and $1.4<r/r_{\rm g}<21$ and $0.2\pi<\theta<0.8\pi$,
respectively. The thin disc region within $12.6<r/r_{\rm g}<253$ and
$0.4\pi<\theta<0.6\pi$ is refined twice. For the 4-level simulation, an
additional refinement layer is used to better resolve the thin disc over
$r\in[10.5\,r_{\rm g},60\,r_{\rm g}]$.

Outflow boundary conditions are used at both the inner and the outer radial
boundaries. We use pole-axis boundary conditions near the poles, where
information from the corresponding zone across the pole is used for
calculating gradients, although fluxes, magnetic fields, and electric fields
(emf’s) are zeroed out on the faces and edges that touch the pole. Lastly,
periodic boundary conditions are implemented along the azimuthal direction.
The zero-divergence of the magnetic fields is ensured using a constrained
transport scheme (Fragile et al., 2012b).

As expected, once the simulations are set to evolve, the initial magnetic
fields trigger the magnetorotational instability, causing turbulence to
transport the angular momentum and thus allow accretion. The heat generated in
this process inflates the torus, thus forming a hot thick flow in the inner
regions. The cooling function, which is set to work in the outer thin slab
region ($r\geq 15$$r_{\rm g}$), acts as a heat sink and helps to maintain the
target thickness, $H/r\approx 0.05$. The simulations have been run for more
than $25,000$$r_{\rm g}$$/c$ for the 4-level simulation and around
$15,000$$r_{\rm g}$$/c$ for the 3-level simulations. All of them have
successfully maintained a two-component flow, as anticipated.

##  3 Results

###  3.1 Precession angle

We compute the precession angle, $\gamma$, by adopting the same formalism as
in Fragile et al. (2007); this angle quantifies how much the disc angular
momentum vector has twisted, or precessed, around the black hole spin axis, as

$\displaystyle\gamma(r)=\cos^{-1}\left[\frac{\bm{J}_{\rm BH}\times\bm{J}_{\rm
disc}(r)}{\lvert\bm{J}_{\rm BH}\times\bm{J}_{\rm
disc}(r)\rvert}\cdot\hat{\textbf{y}}\right]\leavevmode\nobreak\ ,$ (4)

where $\bm{J}_{\rm
BH}=\left(-aM\sin\beta_{0}\hat{x},\,0,\,aM\cos\beta_{0}\hat{z}\right)$ is the
angular momentum vector of the black hole, and

$\displaystyle\left(\bm{J}_{\rm
disc}\right)_{\eta}(r)=\frac{\epsilon_{\mu\nu\sigma\eta}L^{\mu\nu}S^{\sigma}}{2\sqrt{-S^{\alpha}S_{\alpha}}}$
(5)

for $\eta=1,2,3$ corresponds to the Cartesian vector components of the angular
momentum vector of the disc as measured in asymptotically flat space. Here

$\displaystyle L^{\mathrm{\mu\nu}}=\int\left(x^{\mu}T^{\mathrm{\nu
0}}-x^{\nu}T^{\mathrm{\mu 0}}\right){\rm d}^{\mathrm{3}}x$ (6)

and $S^{\sigma}=\int T^{\sigma 0}{\rm d}^{\mathrm{3}}x$, where $T^{\mu\nu}$ is
the MHD stress-energy tensor of the fluid [see Appendix A]. As a reminder, the
angular momentum vector of the disc at $t=0$ is aligned along the direction of
the $\hat{z}$ unit vector. To avoid the degeneracy in cosine for angles
greater than $180^{\circ}$, we also look at the projection of $\bm{J}_{\rm
BH}\times\bm{J}_{\rm disc}(r)$ onto $\hat{\textbf{x}}$ .

Fig. 1 shows the evolution of the precession angle for simulation a9b15r15L4.
There is a clear signature of a precession front starting from the inner
thick-disc region, which propagates to the outer radii (from lower left to
upper right) over time. The propagation happens on the bending wave timescale
shown by the black, dashed curve 111Bending waves propagate at roughly half
the sound speed, $c_{\rm s}$, so the dashed curve is given by $\int 2{\rm
d}r/\langle c_{\rm s}(r,t)\rangle_{\theta,\phi}$.. Radii very close to the
black hole undergo rapid precession to angles up to $90^{\circ}$, holding
there for a period until the surrounding regions begin to catch up. Beyond 20
$r_{\rm g}$, the initial precession is much more modest ($\gamma\lesssim
30^{\circ}$). This value also freezes in after a very short period of time
($<2000$ $r_{\rm g}$$/c$ after the bending wave passes); this frozen-in twist
is noted by the vertical stripes of color in Fig. 1, such as the vertical, red
patch at $\approx 40$ $r_{\rm g}$. Inside of 20 $r_{\rm g}$, however, the
behavior is much different. Rather than vertical stripes of constant color, we
see nearly simultaneous variation of color between 5 and 20 $r_{\rm g}$. This
simultaneous variation indicates solid-body precession; all radii are
experiencing the same degree of precession at roughly the same time. An
example is the nearly horizontal, dark green band around a time of $10^{4}$
$r_{\rm g}$$/c$, indicating that the inner torus has precessed $\approx
60^{\circ}$ by this time.

Figure 1: Space-time plot of the precession angle, $\gamma$ (measured in
degrees) for simulation a9b15r15L4. The dashed curve represents the bending
wave propagating at half the sound speed. The consistent precession angle seen
between 5 and 20$r_{\rm g}$ at any given time suggest that this region is
undergoing rigid-body precession.

###  3.2 Effect of the outer thin disc on the precessing inner thick disc

The Lense-Thirring precession frequency of a test-particle in an orbit
misaligned with the spin-axis of the black hole is given by

$\Omega_{\rm
LT}=\Omega_{\phi}\left[1-\sqrt{1-\frac{4a_{*}}{r^{3/2}}+\frac{3a_{*}^{2}}{r^{2}}}\right]\leavevmode\nobreak\
.$ (7)

For an isolated, torus-like structure between an inner and outer radius,
$r_{\rm i}$ and $r_{\rm o}$, precessing as a solid body, the angular frequency
of precession is computed from weighting $\Omega_{\rm LT}$ by the angular
momentum of the torus:

$\Omega_{\rm prec}=\frac{\int_{r_{\rm i}}^{r_{\rm o}}\Omega_{\rm LT}\Sigma
r^{3}\Omega_{\phi}{\rm d}r}{\int_{r_{\rm i}}^{r_{\rm o}}\Sigma
r^{3}\Omega_{\phi}{\rm d}r}\leavevmode\nobreak\ ,$ (8)

where $\Omega_{\phi}(r)$ is the local angular velocity. Assuming a power-law
profile for the surface density, $\Sigma(r)\sim\Sigma(r_{\rm i})(r/r_{\rm
i})^{-\zeta}$, one can arrive at an analytical estimate of $\Omega_{\rm prec}$
given by eq.(43) in Fragile et al. (2007).

Fig. 2 shows the cumulative precession angle for matter between 5 and
15$r_{\rm g}$ in our three simulations. We exclude the region inside
$r=5$$r_{\rm g}$, because tilted discs have effective inner radii well outside
of the marginally stable limit (Fragile, 2009). The expected precession rates
[eq. (8)] for tori extending over this range of radii are represented by the
grey dashed lines, with $\zeta=-6.1$ for the $a_{*}=0.9$ simulation and $-5.2$
for the $a_{*}=0.5$ one. We see that the predicted precession rates agree
reasonably well with the early ($t\lesssim 1000r_{\rm g}/c$) behavior of the
simulations. However, after $t\approx 2000r_{\rm g}/c$, which is roughly when
the bending wave reaches the transition region ($15\lesssim r\lesssim
20$$r_{\rm g}$, see Fig. 1), the slopes of the cumulative precession curves
change dramatically. The precession rate decreases by more than 95 percent to
a rate comparable to what eq. (8) would give if $r_{\rm o}$ were $40r_{\rm
g}$. This establishes the key results of our paper: the outer thin disc
clearly affects the precession of the inner thick disc by reducing its rate.
It is encouraging that we find nearly identical precession rates in the low-
and high-resolution, $a_{*}=0.9$ simulations, suggesting that this conclusion
is robust. To ensure that the decrease in the precession rate is due to the
outer thin disc, we performed another simulation identical to a9b15r15L4,
except that the outer thin disc is absent, and found that it precesses
continuously at the rate predicted by eq. 8.

Figure 2: Time evolution of the cumulative precession angle over the range
$5$-$15$$r_{\rm g}$. The slopes of these curves give the precession rates. The
grey, dashed lines represent the precession rates estimated from eq. (8). The
coloured, dotted curves represent the precession estimated using eq. (12) for
the simulations a9b15r15L4 and a5b15r15L3.

The decrease in the precession rate in the truncated disc simulations can be
attributed to the exchange of angular momentum between the thick and thin
discs. We discuss how to estimate this effect in Appendix A. In particular, if
we include both the torque and flux terms in eq. (12), then we find remarkable
agreement between eq. (9) and the simulation results, as shown by the dotted
curves in Fig. 2. It is worth pointing out that the if we exclude the flux
terms in eq. (12), we recover the expected precession rates from eq. (8).

##  4 Discussion

Over the last decade, the Lense-Thirring precession model by Ingram et al.
(2009) has gained wide popularity for explaining the type-C QPOs observed in
X-ray binary systems. The dependence of the precession frequency on the outer
radius of the torus in this model, facilitates an independent measurement of
the truncation radius when compared to observations. According to Ingram et
al. (2009), the estimated frequencies for an isolated torus with $r_{\rm i}$
between $5-10$$r_{\rm g}$ and $r_{\rm o}$ varying from $10-50$$r_{\rm g}$ can
produce frequencies in the observed range of $0.01-10$ Hz. However, as we
showed in the previous section, the presence of a thin disc outside the hot,
thick flow decreases the precession frequency by about 95 percent. This has
direct implications for the truncation radius estimation.

For a typical black hole mass of $10\,M_{\odot}$, the precession frequencies
resulting from our simulations are roughly 0.14 Hz and 0.05 Hz for $a_{*}=0.9$
and 0.5, respectively. These frequencies fall near the lower end of the type-C
QPO range (Remillard & McClintock, 2006). To reach the upper end of the
frequency range requires decreasing the size of the precessing torus. In this
way, the slower precession frequencies noted in this work actually help
relieve some of the tensions between the Lense-Thirring precession model and
some type-C QPO observations. One such example is the 2.2 Hz QPO observed in
the black hole binary system GRS 1915+105 (Nathan et al., 2022). According to
eq. (8), the outer edge of an isolated torus precessing around a
$10.1\,M_{\odot}$ black hole with $a_{*}=0.998$ would need to be around
$41$$r_{\rm g}$ to produce a frequency this low. However, the inferred
truncation radius from phase-resolved spectroscopy of this particular QPO is
$1.43^{+0.01}_{-0.02}$$r_{\rm g}$ (Nathan et al., 2022). One way this conflict
could be resolved is if the outer, thin disc decreases the Lense-Thirring
precession frequency of the inner torus as found in our work. Assuming that
the percentage decrease in the precession frequency does not vary with the
truncation radius, we can match the 2.2 Hz frequency with an outer edge now
only at $\approx 11$$r_{\rm g}$. This is still in conflict with the phase-
resolved spectroscopy results, but moves things considerably in the right
direction. Of course, we do not know yet how the precession rate will vary as
the truncation radius changes, something that will need to be explored in
future studies.

We conclude that in a truncated-disc geometry, the rate of Lense-Thirring
precession of the hot, thick flow around a black hole can be significantly
reduced from the isolated case by the advected angular momentum fluxes from
the surrounding thin disc. In principle, this result could potentially be
relevant to other astrophysical phenomena based on Lense-Thirring precession,
such as tilted discs in tidal disruption events (TDEs), where the exchange of
angular momentum between the incoming debris stream from the disrupted star
and the precessing disc around the supermassive black hole could potentially
slow down the precession or even halt it if the fluxes are strong enough; this
is something that will be studied with future simulations.

## Acknowledgements

We want to thank Pavel Ivanov for helpful discussions regarding this work.
This research was partly supported by the Polish NCN grants
2018/29/N/ST9/02139, 2019/33/B/ST9/01564 and Polish National Agency for
Academic Exchange grant PPN/IWA/2018/1/00099/U/0001. Simulations in this work
were performed using the Cyfronet Prometheus cluster, part of the PL-Grid
computing network and the Max Planck Computing & Data Facility. This work also
used the Extreme Science and Engineering Discovery Environment (XSEDE), which
is supported by National Science Foundation grant (NSF) number ACI-1548562.
PCF gratefully acknowledges the support of the NSF through grants AST1616185,
PHY-1748958, and AST-1907850. This work was performed in part at the Aspen
Center for Physics, which is supported by NSF grant PHY-1607611.

## Data Availability

The data underlying this paper will be shared upon reasonable request to the
corresponding author.

## References

  * Anninos et al. (2005) Anninos P., Fragile P. C., Salmonson J. D., 2005, ApJ, 635, 723
  * Begelman & Armitage (2014) Begelman M. C., Armitage P. J., 2014, ApJ, 782, L18
  * Belloni (2010) Belloni T. M., 2010, in Belloni T., ed., , Vol. 794, Lecture Notes in Physics, Berlin Springer Verlag. p. 53
  * Beloborodov (1999) Beloborodov A. M., 1999, ApJ, 510, L123
  * Chakrabarti (1985) Chakrabarti S. K., 1985, ApJ, 288, 1
  * Chakrabarti & Titarchuk (1995) Chakrabarti S., Titarchuk L. G., 1995, ApJ, 455, 623
  * Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
  * Eardley et al. (1975) Eardley D. M., Lightman A. P., Shapiro S. L., 1975, ApJ, 199, L153
  * Esin et al. (1997) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 865
  * Fabian et al. (2009) Fabian A. C., et al., 2009, Nature, 459, 540
  * Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105
  * Fragile (2009) Fragile P. C., 2009, ApJ, 706, L246
  * Fragile et al. (2007) Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668, 417
  * Fragile et al. (2012a) Fragile P. C., Wilson J., Rodriguez M., 2012a, MNRAS, 424, 524
  * Fragile et al. (2012b) Fragile P. C., Wilson J., Rodriguez M., 2012b, MNRAS, 424, 524
  * Galeev et al. (1979) Galeev A. A., Rosner R., Vaiana G. S., 1979, ApJ, 229, 318
  * Haardt & Maraschi (1993) Haardt F., Maraschi L., 1993, ApJ, 413, 507
  * Harten et al. (1983) Harten A., Lax P. D., Leer B. V., 1983, SIAM Review, 25, 35
  * Ingram & Done (2011) Ingram A., Done C., 2011, MNRAS, 415, 2323
  * Ingram & Motta (2019) Ingram A. R., Motta S. E., 2019, New Astron. Rev., 85, 101524
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Ingram et al. (2016) Ingram A., van der Klis M., Middleton M., Done C., Altamirano D., Heil L., Uttley P., Axelsson M., 2016, MNRAS, 461, 1967
  * Kato & Fukue (1980) Kato S., Fukue J., 1980, PASJ, 32, 377
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liu et al. (2007) Liu B. F., Taam R. E., Meyer-Hofmeister E., Meyer F., 2007, ApJ, 671, 695
  * Martocchia & Matt (1996) Martocchia A., Matt G., 1996, MNRAS, 282, L53
  * McClintock & Remillard (2006) McClintock J. E., Remillard R. A., 2006, in , Vol. 39, Compact stellar X-ray sources. pp 157–213
  * McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
  * Motta (2016) Motta S. E., 2016, Astronomische Nachrichten, 337, 398
  * Motta et al. (2012) Motta S., Homan J., Muñoz Darias T., Casella P., Belloni T. M., Hiemstra B., Méndez M., 2012, MNRAS, 427, 595
  * Motta et al. (2015) Motta S. E., Casella P., Henze M., Muñoz-Darias T., Sanna A., Fender R., Belloni T., 2015, MNRAS, 447, 2059
  * Nathan et al. (2022) Nathan E., et al., 2022, MNRAS, 511, 255
  * Parker et al. (2015) Parker M. L., et al., 2015, ApJ, 808, 9
  * Pringle (1992) Pringle J. E., 1992, MNRAS, 258, 811
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Rodriguez et al. (2002) Rodriguez J., Varnière P., Tagger M., Durouchoux P., 2002, A&A, 387, 487
  * Silbergleit et al. (2001) Silbergleit A. S., Wagoner R. V., Ortega-Rodríguez M., 2001, ApJ, 548, 335
  * Spiteri & Ruuth (2002) Spiteri R. J., Ruuth S. J., 2002, SIAM Journal on Numerical Analysis, 40, 469
  * Stella & Vietri (1998) Stella L., Vietri M., 1998, ApJ, 492, L59
  * Stella et al. (1999) Stella L., Vietri M., Morsink S. M., 1999, ApJ, 524, L63
  * Svensson (1996) Svensson R., 1996, A&AS, 120, 475
  * Tagger & Pellat (1999) Tagger M., Pellat R., 1999, A&A, 349, 1003
  * Wagoner (1999) Wagoner R. V., 1999, Phys. Rep., 311, 259
  * White et al. (2019) White C. J., Quataert E., Blaes O., 2019, ApJ, 878, 51
  * van den Eijnden et al. (2017) van den Eijnden J., Ingram A., Uttley P., Motta S. E., Belloni T. M., Gardenier D. W., 2017, MNRAS, 464, 2643

##  Appendix A Calculation of precession rate

Let $J_{1}$ and $J_{2}$ be orthogonal vector components of the angular
momentum of the precessing torus projected onto the equatorial plane of the
black hole. We can then define the precession angle, $\gamma$, of the torus
extending from $r_{\rm i}$ to $r_{\rm o}$ as $\tan\gamma=J_{2}/J_{1}$.
Differentiating with respect to time gives us the precession rate,

$\displaystyle\frac{\partial\gamma}{\partial
t}=\frac{1}{J_{1}^{2}+J_{2}^{2}}\left(J_{1}\frac{\partial J_{2}}{\partial
t}-J_{2}\frac{\partial J_{1}}{\partial t}\right)\leavevmode\nobreak\ .$ (9)

All we need to do is compute the terms $J_{i}$ and $\partial J_{i}/\partial t$
(for $i=1,2$) from the simulation data.

For $J_{i}$, we take

$\displaystyle J_{i}=\int_{r_{\rm i}}^{r_{\rm o}}{\rm d}r{\rm d}\theta{\rm
d}\phi\sqrt{-g}l_{i}\leavevmode\nobreak\ ,$ (10)

where $l_{1}=l^{x}\cos\beta_{0}+l^{z}\sin\beta_{0}$ and $l_{2}=l^{y}$; the
vector components of the angular momentum density ($\bm{l}$) in Cartesian
coordinates are $l^{x}=yT^{tz}-zT^{ty}$, $l^{y}=zT^{tx}-xT^{tz}$, and
$l^{z}=xT^{ty}-yT^{tx}$; and $T^{\mu\nu}=\left(\rho h+2P_{\rm
m}\right)u^{\mu}u^{\nu}+\left(P_{\rm g}+P_{\rm
m}\right)g^{\mu\nu}-b^{\mu}b^{\nu}$ is the stress-energy tensor of the fluid
in the MHD limit, where $h=1+\epsilon+P_{\rm g}/\rho$ is the specific
enthalpy, $\epsilon$ is the specific internal energy density, $P_{\rm m}$ is
the magnetic pressure, $u^{\mu}$ is the fluid four-velocity, and $b^{\mu}$ is
the four-vector magnetic field measured by an observer comoving with the
fluid.

To predict the time-derivatives of the angular momentum in a simple model, we
use the angular momentum conservation equation (Pringle, 1992)

$\displaystyle\frac{\partial\bm{L}}{\partial
t}+\frac{1}{r}\frac{\partial}{\partial
r}\left(rV^{r}\bm{L}\right)=\bm{\Omega}_{\rm
LT}\times\bm{L}\leavevmode\nobreak\ ,$ (11)

where $\bm{L}$ is the angular momentum vector per unit surface area and
$V^{r}$ is the radial velocity. Technically, this equation was derived in the
diffusive regime, but as we show, it works well even for our truncated discs.
We have ignored the viscous terms in the above conservation equation, as we
find them to be negligible when compared to the Lense-Thirring torque and
angular momentum fluxes in our simulations. Integrating the vector components
of the above equation over the surface area, we get

$\displaystyle\frac{\partial J_{i}}{\partial t}=\mp\int_{r_{\rm i}}^{r_{\rm
o}}{\rm d}r{\rm d}\theta{\rm d}\phi\sqrt{-g}\Omega_{\rm
LT}l_{j}-\left[F_{i}|_{r_{\rm o}}-F_{i}|_{r_{\rm i}}\right],$ (12)

where the minus sign and $j=2$ are taken for $i=1$, and the plus sign and
$j=1$ are taken for $i=2$. The first term on the right-hand-side is the usual
Lense-Thirring torque term, and the terms within braces represent the angular
momentum flux through the boundaries, $r_{\rm i}$ and $r_{\rm o}$, given by

$\displaystyle F_{i}=\int{\rm d}\theta{\rm
d}\phi\sqrt{-g}l_{i}V^{r}\leavevmode\nobreak\ .$ (13)


# X-ray studies of neutron stars and the equation of state

Zhiqiang Miao ID 1,2 Liqiang Qi3 Juan Zhang3 Ang Li ID 2<EMAIL_ADDRESS>Mingyu Ge3<EMAIL_ADDRESS>1 Tsung-Dao Lee Institute, Shanghai Jiao Tong
University, Shanghai 201210, China;
2 Department of Astronomy, Xiamen University, Xiamen 361005, China;
3 Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, China

(August 28, 2024)

###### Abstract

The understanding of neutron star equation of state hinges on a comprehensive
analysis of multi-messenger, multi-wavelength data. The recent scrutiny of PSR
J0030+0451 data by NICER introduces complexities, unveiling a tension with
another X-ray observation of the central compact object in HESS J1731-347,
specifically concerning the mass-radius constraint of low-mass neutron stars.
This tension persists when integrating NICER’s updated data with LIGO/Virgo’s
gravitational-wave data from the GW170817 binary neutron star merger. Despite
attempts to reconcile these disparate observations, the current combined data
still can not distinguish different types of neutron stars – whether they are
pure neutron stars or hybrid stars. Bayesian inference indicates only modest
changes in the posterior ranges of parameters related to the nuclear matter
and deconfinement phase transition. This ongoing exploration underscores the
intricate challenges in precisely characterizing neutron stars. It also points
out that it is possible to probe the equation of state at different density
regimes from future more accurate radii of neutron stars with various masses.

##  I Introduction

The determination of the equation of state (EoS) for neutron stars (NSs)
through quantum many-body computations proves to be a formidable task due to
the complicated nonperturbative problems in physics (Huber, 2020). Effective
discrimination between different EoSs based on nuclear experiments and
astrophysical observations is further hindered by a deficit and limited
precision of measurements (Li, 2022). Consequently, the study of EoS, namely
how the pressure $P$ of NS matter system changes with the stellar (energy)
density $\varepsilon$, relies on an iterative interplay between experiment and
theory, with their reciprocal feedback proving crucial for progress (Baym  _et
al._ , 2018; Li  _et al._ , 2020; Burgio  _et al._ , 2021; Lattimer, 2021; Li
_et al._ , 2021a). Assuming that NSs obey general relativity, the EoS
establishes a unique sequence of stars in hydrostatic equilibrium using the
Tolman-Oppenheimer-Volkoff (TOV) equation. Accurate measurements of mass and
radius of NSs can then be employed to constrain physical quantities in the
EoS, describing the properties of dense nuclear matter (i.e. the
incompressibility, the symmetry energy) and strangeness phase transitions
(i.e. the critical chemical potential of the phase transitions, the sound
velocity).

The direct measurement of the mass and radius of NSs is extremely challenging.
An indirect measurement has been proposed to constrain the mass and radius of
NSs based on the pulse-profile modeling of the thermal emission from hot spots
on the stellar surface of NSs (Watts  _et al._ , 2016; Bogdanov  _et al._ ,
2019; Watts, 2012). The characteristics of the observed pulse profile of a
millisecond pulsar (MSP) are correlated to the mass, radius, and emission
configurations of hot regions on the stellar surface via general relativity
and special relativity effects, including the light bending, gravitational
redshift, Doppler shifts, relativistic aberration, and propagation time
differences. Bayesian inference allows the estimation of the mass and radius
of NSs conditional on the simultaneous spectral and timing measurements of
soft X-rays. It has been successfully applied to measure the mass and radius
of two MSPs, PSR J0030+0451 and PSR J0740+6620, as demonstrated by the NICER
collaboration (Riley  _et al._ , 2019; Miller  _et al._ , 2019, 2021; Riley
_et al._ , 2021). Conversely, the constraints on mass and radius in X-ray flux
oscillations in accreting pulsars are highly dependent on the adopted models
(Watts, 2012; Miller and Lamb, 2015; Watts  _et al._ , 2016).

Despite the sensitivity of the NICER likelihood function to the geometry and
temperature configurations of the surface hot regions, the background
estimation is completely free when analyzing PSR J0030+0451 Riley  _et al._
(2019); Miller  _et al._ (2019). It can cause under-estimation or over-
estimation of the background when performing the NICER-only analysis, which
consists of instrumental noise, cosmic background, and contamination sources
around the target. To address this, Ref. Riley  _et al._ (2021) conducted the
joint analysis of the NICER and XMM-Newton datasets of PSR J0740+6620, where
the phase-averaged X-ray data from XMM-Newton EPIC was utilized to constrain
the unpulsed portion of the pulse profile. Consequently, the XMM-Newton
likelihood function significantly reduced the NICER posterior volume,
impacting the inferred mass-radius and corresponding radiation geometries. The
inclusion of XMM-Newton data led to adjustments in the posteriors, favoring
less compact stars. It highlighted the importance of background estimation
through instrument cross-calibration. The overall shift in the posterior
probability distribution function (PDF) of the radius for PSR J0740+6620 was
found to be much smaller than the measurement uncertainty, given the
informative mass prior Miller (2016). This observed shift is anticipated to
have an insignificant effect on EoS inference with typical EoS priors Greif
_et al._ (2019); Raaijmakers  _et al._ (2021).

An interesting question arises regarding the effect of a joint analysis of the
NICER and XMM-Newton data for PSR J0030+0451 on mass-radius inference. In a
recent study Vinciguerra  _et al._ (2024), the radius and mass of PSR
J0030+0451 were reanalyzed with an updated NICER response matrix and an
upgraded analysis framework. To better estimate the background, i.e. the
phase-invariant components that do not come from the hot regions on the
stellar surface, the spectral data of the XMM-Newton EPIC MOS1 and MOS2 were
considered. The inclusion of the XMM-Newton data is expected to yield a more
robust constraint on the background, attributed to its imaging capability.
Substantial shifts in solution(s) were identified compared to those reported
in Ref. Riley  _et al._ (2019). Importantly, the newly introduced models infer
configurations with one hot spot on the same hemisphere as the observer, a
feature absent in Ref. Riley  _et al._ (2019). The introduction of the XMM-
Newton data led to a reduction in background estimates from NICER-only data,
resulting in higher compactness values Riley  _et al._ (2021); Salmi  _et al._
(2022). In the absence of an informative mass prior from radio timing
measurements of PSR 0030+0451, the overall shift of both mass and radius is
evident and is anticipated to significantly impact the underlying EoS of NSs.
In the present work, we are primarily interested in how these improved mass
and radius constraints affect the EoS study of NSs with and without
strangeness phase transitions111The present work will focus on the study of
the EoS describing the NS core and does not discuss the properties of the NS
crust. In our analysis, the BPS (Baym  _et al._ , 1971) and NV (Negele and
Vautherin, 1973) EoSs are employed for the outer and inner crust,
respectively. .

On the other hand, the challenge arises in reconciling the large radii of PSR
J0740+6620 and PSR J0030+0451 with other observations, such as the tidal
deformability inferred from the inspiring gravitational waves of GW170817
Abbott  _et al._ (2018, 2019). In particular, a recent analysis of the
supernova remnant HESS J1731-347 suggests a very low mass and a small radius
for the central compact object (CCO) within it (Doroshenko  _et al._ , 2022).
Various theories and models have been proposed to interpret the observed mass
and radius, often incorporating exotic degrees of freedom beyond nucleons (see
discussions in e.g., (Yuan and Li, 2023)), but none have yielded conclusive
results thus far. Our interest extends to the exploration of multi-messenger
observations and their potential to distinguish different phase states of the
EoS modeling.

The paper is organized as follows: In Sec. II, we provide an overview of the
relevant data utilized in this study. Sec. III offers a brief introduction to
the fundamentals of EoS modeling. In Sec. IV we present the results of our
inference for the global NS properties and microscopic EoS parameters, for
both nuclear matter and phase transitions. We conclude in Sec. V.

##  II Multimessenger mass–radius constraints

###  II.1 Millisecond pulsars PSR J0030+0451 and PSR J0740+6620

The mass-radius measurements of PSR J0030+0451 and PSR J0740+6620, as reported
by the NICER collaboration, are considered in the present work. The
contribution of the unpulsed component to the thermal emission of hot spots is
constrained by jointly analyzing the NICER and XMM-Newton data.

The Bayesian inference of the radius and mass of the rotation-powered MSP PSR
J0030+0451 was initially performed independently by Miller et al. Miller  _et
al._ (2019) and Riley et al. Riley  _et al._ (2019), conditional on the pulse-
profile data sets from the NICER observations. Different configurations of hot
regions were proposed by Miller et al. Miller  _et al._ (2019) and Riley et
al. Riley  _et al._ (2019), respectively. The two studies proposed different
configurations of hot regions, with Miller et al. Miller  _et al._ (2019)
suggesting three oval, uniform-temperature emitting spots while Riley et al.
Riley  _et al._ (2019) proposed one small angular spot and the other being an
azimuthally extended narrow crescent. Despite these differences in inferred
hot region properties, the estimation of the equatorial radius and
gravitational mass are consistent. Miller et al. Miller  _et al._ (2019)
reported $R=13.02_{-1.06}^{+1.24}\,\;\text{km}$ and
$M=1.44_{-0.14}^{+0.15}\,\,M_{\odot}$, while Riley et al.Riley  _et al._
(2019) reported $R=12.71_{-1.19}^{+1.14}\;\text{km}$ and
$M=1.34_{-0.16}^{+0.15}\,\,M_{\odot}$, at the 68% confidence level.

As introduced earlier, Ref. Vinciguerra  _et al._ (2024) recently employed an
improved pipeline to conduct a joint analysis of the NICER and XMM-Newton
datasets for PSR J0030+0451. They employed four models to describe the surface
heating distribution, with increasing complexity, ST-U (two single-temperature
spherical caps), ST+PST (a single-temperature spherical cap and the other one
protruding single-temperature), ST+PDT (a single-temperature spherical cap and
the other one protruding dual-temperature) and PDT-U (two protruding dual-
temperature spherical caps). The inferred radius and mass of the ST-U and
ST+PST models are consistent with previous studies Miller  _et al._ (2019);
Riley  _et al._ (2019), but are strongly disfavoured with smaller evidence
values compared to those of ST+PDT and PDT-U. The newly adopted models (ST+PDT
and PDT-U) exhibit a higher preference in terms of Bayesian evidence.
Conditioned on the most favored PDT-U model, NICER data alone yields
$M=1.41_{-0.19}^{+0.20}\,M_{\odot}$ and $R=13.12_{-1.21}^{+1.35}\;\text{km}$.
When incorporating data from both NICER and XMM-Newton, the derived values are
$M=1.70_{-0.19}^{+0.18}\,M_{\odot}$ and $R=14.44_{-1.05}^{+0.88}\;\text{km}$.
In this study, we employ a bivariate Gaussian distribution to emulate the two
recent results:

$p(\bm{x})=\frac{\exp\left[-\frac{1}{2}(\bm{x}-\bm{\mu})^{T}\Sigma^{-1}(\bm{x}-\bm{\mu})\right]}{2\pi\sqrt{|\bm{\Sigma}|}}$
(1)

where $\bm{x}=(M,R)^{T}$. The mean vector $\bm{\mu}$ and the covariance matrix
$\bm{\Sigma}$ are

$\bm{\mu}=\begin{pmatrix}\mu_{M}\\\ \mu_{R}\\\
\end{pmatrix},\quad\bm{\Sigma}=\begin{pmatrix}\sigma_{M}^{2}&\rho\sigma_{M}\sigma_{R}\\\
\rho\sigma_{M}\sigma_{R}&\sigma_{R}^{2}\\\ \end{pmatrix}$ (2)

For results using only NICER data, we use $\mu_{M}=1.41\,M_{\odot}$,
$\mu_{R}=13.12\;\text{km}$, $\sigma_{M}=0.20\,M_{\odot}$,
$\sigma_{R}=1.3\;\text{km}$ and $\rho=0.9$. While for results using both NICER
and XMM-Newton data, we use $\mu_{M}=1.70\,M_{\odot}$,
$\mu_{R}=14.44\;\text{km}$, $\sigma_{M}=0.19\,M_{\odot}$,
$\sigma_{R}=0.9\;\text{km}$ and $\rho=0.9$.

Similarly, we utilize the posteriors derived from the ST-U model for PSR
J0740+6620, incorporating both NICER data alone and the joint analysis of
NICER and XMM-Newton data. In this case, the mass constraint of $2.08\pm
0.07\,M_{\odot}$ obtained from radio timing measurement (Fonseca  _et al._ ,
2021) has been incorporated in the Bayesian inference, leading to a more
precise estimation of the mass-radius of NSs. When considering NICER data
exclusively, the estimated equatorial radius and gravitational mass are are
reported as $R=12.39_{-0.98}^{+1.30}\;\text{km}$ and
$M=2.072_{-0.066}^{+0.067}\,M_{\odot}$ from (Riley  _et al._ , 2021) and
$R=13.71_{-1.50}^{+2.61}\;\text{km}$ and
$M=2.062_{-0.091}^{+0.090}\,M_{\odot}$ from (Miller  _et al._ , 2021). With
the inclusion of both NICER and XMM-Newton data, the inferred radius and mass
are presented as $R=11.29_{-0.81}^{+1.20}\;\text{km}$ and
$M=2.078_{-0.063}^{+0.066}\,M_{\odot}$ from (Riley  _et al._ , 2021) and
$R=11.51_{-1.13}^{+1.87}\;\text{km}$ and
$M=2.072_{-0.094}^{+0.087}\,M_{\odot}$ from (Miller  _et al._ , 2021). For
consistency, we adopt the results from (Riley  _et al._ , 2021) in the present
work.

Figure 1: Gravitational mass as a function of the stellar radius based on some
representative EoS models with or without phase transition taken from Refs.
(Burgio  _et al._ , 2021; Zhu  _et al._ , 2018) (see text for details). The
mass-radius measurements from the NICER mission for the pulsars PSR J0030+0451
and PSR J0740+6620 are displayed, with the contours representing the results
obtained solely from NICER data and the shaded areas reflecting the results
obtained by utilizing both NICER and XMM-Newton data. Additionally, the mass-
radius inferred from the GW170817 tidal deformability measurement and the
mass-radius measurement of the central compact object (CCO) in the supernova
remnant HESS J1731-347 are included. All these measurements are presented at
the 90% confidence level. The grey line denotes the lower boundary of mass-
radius, corresponding to the currently known maximum frequency $\nu=716\,{\rm
Hz}$. The range for the minimum spin period of MSPs from Ref. (Liu  _et al._ ,
2023) ($P_{\rm min}=1.43_{-0.09}^{+0.03}\,{\rm ms}$) is also presented at the
68% confidence level.

Figure 2: Posteriors of mass-radius relations for NSs (left two panels) and
HSs (right two panels), all at the 90% confidence level. The upper two panels
present the results without GW170817 constraints, while the lower two panels
incorporate both NICER and GW170817 constraints. The results with XMM-Newton
data as depicted within shaded regions show a noticeable overall shift towards
the right when compared to those obtained without. Other shaded regions
correspond identically to those depicted in Fig. 1.

Figure 3: Posteriors of EoSs and the sound speed squared in the matter of NSs
(left two panels) and HSs (right two panels), all at the 90% confidence level.
The additional upper horizontal axis in each panel indicates the corresponding
baryon number density within the most favored EoS parameter sets considering
both GW170817 and NICER with XMM-Newton data..  Figure 4: Posterior
distributions of the stellar properties: the maximum mass $M_{\rm TOV}$, the
radius for a $1.4\,M_{\odot}$ star $R_{1.4}$ and the tidal deformability of a
$1.4\,M_{\odot}$ star $\Lambda_{1.4}$, conditioned by different constraints as
outlined in Sec. II. The priors are also shown in dotted curves.

###  II.2 Supernova remnant HESS J1731-347

CCOs, located at the centres of supernova remnants, are radio-quiet isolated
NSs emitting steady thermal X-ray. The majority of confirmed CCOs exhibit an
absence of detectable pulsations (Alford and Halpern, 2023), indicative of a
relatively uniform temperature distribution across the NS surface. Their
energy spectra can be well-fitted with thermal emissions (Alford and Halpern,
2023). It is a challenge to generate such a simple thermal spectrum,
considering the possible energetic physical processes occurring in the NS
shell, atmosphere, or magnetosphere. The simplicity of CCO energy spectra
implies that its emission is uncorrelated to the geometrical uncertainties of
the emission region and the typical accreting aspects, such as accretion
physics and the radiative transfer in strong magnetic fields (Doroshenko  _et
al._ , 2022). These specific characteristics position CCOs as excellent
sources for studying the EoS of dense matter (Doroshenko  _et al._ , 2022).

The CCO in the supernova remnant HESS J1731-347 stands out as the brightest
one in its class and attracts extensive observations. Analyzing XMM-Newton
data, Ref. Klochkov  _et al._ (2015) found the blackbody spectrum model to be
unfavorable, and hydrogen atmosphere models led to an unreasonably large
distance. Instead, the carbon atmosphere model is compatible with the source
location around 3 kpc Klochkov  _et al._ (2015). A reanalysis by Ref.
Doroshenko  _et al._ (2022) using XMM-Newton and Suzaku observations, with
similar conditions as used in Klochkov et al. Klochkov  _et al._ (2015) (an
energy range of $1-10\,\rm keV$, a uniform-temperature carbon atmosphere
model, a distance of $3.2\rm\,kpc$ and a wabs absorption model) resulted in
the same outcomes with slightly tighter parameter ranges. Including the low-
energy band $<1\,\rm keV$ and substituting the outdated wabs absorption model
with tbabs, the NS mass was lowered from $\sim 1.4\,M_{\odot}$ to $\sim
1.0\,M_{\odot}$. However, the lack of a reliable distance introduces
significant uncertainties in the NS parameters. Ref. Doroshenko  _et al._
(2016) provided evidence of an embedded optical star in the same shell as the
CCO within the supernova remnant HESS J1731-347, suggesting that the optical
star and the CCO are originally in a binary system. Estimating the CCO’s
distance based on that of the optical star, available from Gaia parallax
measurements, yields $2.5\,\rm kpc$ (see Doroshenko  _et al._ (2022) and
references therein). With this revised distance, Ref. Doroshenko  _et al._
(2022) further constrains the mass to $0.77^{+0.20}_{-0.17}M_{\odot}$ and the
radius to $10.4^{+0.86}_{-0.78}\;\text{km}$. When treating the distance as a
free parameter with informative priors, the mass and radius parameters of the
CCO are $M=0.83^{+0.17}_{-0.13}\,M_{\odot}$ and
$R=11.25^{+0.53}_{-0.37}\;\text{km}$.

###  II.3 GW170817

The EoS governs not only the stable configuration of a single star but also
the dynamics of NS mergers. During the inspiral phase, the influence of the
EoS is evident in the tidal polarizability $\Lambda$, which can be computed by
treating the tidal field from the companion as perturbations to general-
relativistic hydrodynamic equation. Apart from the mass and radius
measurement, gravitational wave observations from binary NS mergers uniquely
provide independent constraints on the EoS (Miao  _et al._ , 2023). In this
study, we calculate the likelihood of gravitational wave event GW170817
(Abbott  _et al._ , 2017) through a high-precision interpolation of the
likelihood, which is developed in (Hernandez Vivanco  _et al._ , 2020) from
fitting the strain data released by LIGO/Virgo. Encapsulated in the PYTHON
package toast222https://git.ligo.org/francisco.hernandez/toast, the likelihood
function is expressed as:

$\mathcal{L}_{\rm GW}=F(M_{\rm ch},q,\Lambda_{1},\Lambda_{2})$ (3)

where $F(\cdot)$ represents the interpolation function. $M_{\rm ch}$ and $q$
denote chirp mass and mass ratio respectively, while tidal deformability
$\Lambda_{i}$ is evaluated from the component mass and the EoS, expressed as
$\Lambda_{i}\equiv\Lambda_{i}(M_{i};{\rm EoS})$. Note $M_{1}=M_{\rm
ch}(1+q)^{1/5}/q^{3/5}$ and $M_{2}=M_{1}q$.

###  II.4 Critical spin frequency

The Keplerian frequency is one of the most studied physical quantities for
pulsars (Zhou  _et al._ , 2021). An EoS that predicts Kepler frequencies that
are smaller than the observed rotational frequencies is to be rejected, as it
is not compatible with observation. The Kepler frequency, namely the maximum
spin frequency, for an NS, is approximately given by $\nu_{\rm
max}=1045\left(\frac{M}{M_{\odot}}\right)^{1/2}\left(\frac{R}{10{\rm
km}}\right)^{-3/2}\,{\rm Hz}$, marking the threshold beyond which mass
shedding initiates. This frequency depends significantly on the EoS governing
the internal composition of pulsars (Haensel  _et al._ , 2009). The fastest
known pulsar is PSR J1748-2446ad, with a frequency of $\nu=716\,{\rm Hz}$
(period $P=1.396$ ms). Recently, we studied the period distribution of MSPs
(Liu  _et al._ , 2023) utilizing the largest available samples from known
pulsar surveys and found no statistical support for a minimum spin period
cutoff. Bayesian evidence favoring certain period distribution models suggests
a minimum spin period of $P_{\rm min}=1.43^{+0.03}_{-0.09}\,\rm ms$ is
indicated. This constraint serves as an upper limit in the mass-radius
relationship (see Figs. 1 and 2), with no capacity yet to distinguish between
different models of EoSs.

Table 1: Most probable intervals of the stellar properties (quoted as median+68% credible interval): the maximum mass $M_{\rm TOV}$ in $\,M_{\odot}$, the radius of a $1.4\,M_{\odot}$ star $R_{1.4}$ in km and the dimensionless tidal deformability of a $1.4\,M_{\odot}$ star $\Lambda_{1.4}$.  | $R_{1.4}$ | $\Lambda_{1.4}$ | $M_{\rm TOV}$
---|---|---|---
NS (w.o. XMM) | $12.24_{-0.36}^{+0.49}$ | $406.10_{-64.56}^{+88.31}$ | $2.09_{-0.07}^{+0.08}$
NS (with XMM) | $12.82_{-0.46}^{+0.36}$ | $506.88_{-88.41}^{+97.11}$ | $2.13_{-0.08}^{+0.10}$
NS (plus HESS) | $12.11_{-0.26}^{+0.29}$ | $398.44_{-57.14}^{+63.66}$ | $2.13_{-0.09}^{+0.10}$
HS (w.o. XMM) | $11.69_{-0.63}^{+0.75}$ | $277.69_{-86.50}^{+161.95}$ | $2.15_{-0.09}^{+0.16}$
HS (with XMM) | $12.85_{-0.51}^{+0.45}$ | $534.71_{-120.78}^{+136.63}$ | $2.21_{-0.14}^{+1.15}$
HS (plus HESS) | $12.20_{-0.47}^{+0.35}$ | $437.62_{-118.18}^{+104.86}$ | $2.35_{-0.25}^{+0.94}$

##  III Theoretical modelling of EoS

Before delving into our inferences on the EoS, we observe the key constraints
detailed above in Sec. II, which are illustrated in Fig. 1. These constraints
are juxtaposed with several representative microscopic EoS models (as taken
from the review article (Burgio  _et al._ , 2021)), starting from realistic
nucleon-nucleon potentials (Argonne V18 (Wiringa  _et al._ , 1995), Bonn B
(Brockmann and Machleidt, 1990)) explicitly including many-body forces; BOB
(N+Y) and V18 (N+Y) are EoSs containing hyperons. Two phenomenological
relativistic-field EoSs are QMF40 and QMF60 from Ref. Zhu  _et al._ (2018),
where the isospin-asymmetric contribution of the nuclear EoS, characterized by
the symmetry energy $E_{\rm sym}(n)$, is purposely changed to be studied
later:

$\displaystyle E/A(n,\beta)=E/A(n,\beta=0)+E_{\rm sym}(n)\beta^{2}+...$ (4)

where $n=n_{n}+n_{p}$ is the total nucleon number density and
$\beta=(n_{n}-n_{p})/n$ is the isospin asymmetry of nuclear matter.
$E/A=\varepsilon/n-m_{N}$ is the energy per nucleon, $m_{N}=939\,\text{MeV}$
being the nucleon mass in free space.

The minimum of the $E/A(n)$ curve for symmetric nuclear matter ($\beta=0$)
corresponds to a bound equilibrium state at zero pressure. The values of $E/A$
and $n$ at this minimum will be denoted by $E_{0}/A$ and $n_{0}$, where
$n_{0}$ is the so-called saturation density $n_{0}=0.16\;\text{fm}^{-3}$
resulting from the interplay of the short-distance repulsion and the long-
distance attraction in the nucleon-nucleon interaction. Around the saturation
density $n_{0}$, both $E/A(n,\beta=0)$ and $E_{\rm sym}(n)$ can be expanded in
isospin asymmetry $\beta$ as follows:

$\displaystyle E/A(n,0)=E_{0}/A+\frac{1}{18}K_{0}\frac{n-n_{0}}{n_{0}}+...$
(5) $\displaystyle E_{\rm sym}(n)=E_{\rm
sym}^{0}+\frac{1}{3}L_{0}\frac{n-n_{0}}{n_{0}}+...$ (6)

where $K_{0}$, $E_{\rm sym}^{0}$ and $L_{0}=3n_{0}({dE_{\rm
sym}}/{dn})_{n_{0}}$ are the incompressibility, the symmetry energy and its
slope at the saturation point, respectively. $K_{0}$ gives the curvature of
$E_{0}/A$ at $n=n_{0}$ and the associated increase of the energy per nucleon
of symmetric nuclear matter due to a small compression or rarefaction, whereas
$E_{\rm sym}^{0}$ determines the increase in the energy per nucleon due to a
small asymmetry $\beta$ in asymmetric nuclear matter.

Strong interactions play a dominating role and are the basic ingredient of
calculations in the NS EoS, i.e.,
$P(\varepsilon)=n(d\varepsilon/dn)-\varepsilon$. The theoretical underpinning
of strong interactions between baryons is the quantum chromodynamics (QCD),
where the fundamental fields are those of quarks and gluons. A comprehensive
exploration of constructing baryonic interactions based on lattice QCD is
available in reviews such as (Beane  _et al._ , 2011; Aoki  _et al._ , 2012).
However, lattice QCD calculations are currently extremely expensive and often
limited to scenarios with large quark masses. An alternative approach, rooted
in the realm of quark degrees of freedom, is formulated within the framework
of effective field theory (Weinberg, 1990, 1991), developed as chiral
effective field theory ($\chi$EFT) approach (Hebeler  _et al._ , 2011). This
approach organizes terms based on their dependence on the physical parameter
$q/m_{N}$, where $m_{N}$ denotes the nucleon mass and $q$ represents a generic
momentum in the relevant Feynman diagrams. This systematic ordering
facilitates the calculation of various terms constituting the baryon force,
yielding nucleon-nucleon interactions that reasonably reproduce two-body data.
Extending chiral perturbation theory to nuclear matter calculations has been a
subject of considerable effort. However, the assumption of a small $q/m_{N}$
parameter inherently limits the applicability of these forces, with the safe
maximum density typically around the saturation value, $n_{0}$. Extrapolation
becomes essential for densities beyond $n_{0}$, and theoretical uncertainty
grows rapidly. Several calculations have extended up to $2\,n_{0}$ (e.g.,
(Drischler  _et al._ , 2020a, b)), with uncertainties estimated by analyzing
the order-by-order convergence in the chiral expansion and many-body
perturbation theory. See more discussions in e.g., (Burgio  _et al._ , 2021;
Li  _et al._ , 2019). In the ensuing sections, we provide illustrative
examples of their application to the study of NS EoS, as depicted in Fig. 5.

Figure 5: Posteriors of baryon number density vs. pressure for NSs (left
panel) and HSs (right panel). Also shown are the 68% and 95% credible regions
derived in Ref. (Drischler  _et al._ , 2021) using N3LO $\chi$EFT calculations
(Drischler  _et al._ , 2020a, b). The data points in both panels with error
bars depict the inferred pressure at $\sim 2\,n_{0}$ from Ref. (Legred  _et
al._ , 2021) without including the updated radius measurement of J0030+0451
incorporating the XMM-Newton data.

Currently, it remains a formidable challenge to connect QCD with the low-
energy nuclear physics phenomena since it is in the non-perturbative regime.
Consequently, one has to resort to phenomenological approaches. Starting in
the 1950s, the development of meson-exchange models for nucleon interaction
marked a significant step forward, grounded in Yukawa’s idea. Namely, the
strong interaction among the different baryons originates from the coupling of
baryons to various meson fields, where quarks and gluons do not appear
explicitly. Correspondingly, such meson-exchange interaction models operate
with the baryon and meson fields, and the meson couplings are described by
corresponding Lagrangian densities, depending on the symmetry behavior of a
meson field under rotations and reflections. The inclusion of mesons with rest
mass below $1\,\rm GeV$ suffices for these models. In the meson-exchange
picture, one can construct the bare/realistic interaction to describe nucleon-
nucleon scattering data and phase shifts at laboratory energies $\leq
350\,\text{MeV}$, like the above-mentioned Argonne V18, Bonn B interactions,
serving as fundamental inputs for microscopic many-body frameworks; One can
also construct in-medium/effective interactions, wherein phenomenological
parameters reproduce various properties of nuclear systems like finite-nuclei
and NS observations. By employing the latter effective interactions, one can
compute the energy of hadronic systems through mean-field approximation. It is
the framework of relativistic mean-field model (Walecka, 1974; Serot, 1992)
that is widely adopted in NS physics. Its Lorentz invariance guarantees that
the sound speed in dense matter does not exceed the speed of light at any
density. The mean-field approximation, applied to a spatially uniform ground
state of hadronic matter, entails hadrons occupying momentum states within
Fermi spheres. The eigenfunctions of hadrons manifest as plane-wave solutions
to the Dirac equations. The self-consistent solution involves solving the
Klein-Gordon equations for meson fields alongside the Dirac equations for the
hadrons. The many-body state is subsequently constructed as an independent
particle, derived from the single-particle wave functions. Hadrons are
represented by four-component Dirac spinors, and meson-field operators are
substituted with their expectation values in the ground state.

Figure 6: Same as Fig. 4, but for various properties for nuclear matter
saturation (the incompressibility $K_{0}$ (in MeV), the symmetry energy
$E_{\rm sym}^{0}$ (in MeV), the symmetry energy slope $L_{0}$ (in MeV) and the
ratio between the effective mass and free nucleon mass $m^{\ast}_{0}/m_{N}$)
and phase transition (the critical density ratio $n_{\rm t}/n_{0}$, the
critical chemical potential $\mu_{\rm crit}$ (in MeV), the transition strength
$\Delta\varepsilon/\varepsilon_{\rm t}$, and the squared sound velocity
$c_{s}^{2}$) in high-density (quark) matter. See text in Sec. III.

In the computation of the NS EoS, the selected theoretical approaches also
depend on the relevant degrees of freedom of the problem, from nuclei and
nucleons at lower densities to additional particles, such as hyperons and
quarks, at high densities (Li  _et al._ , 2008; Zhu  _et al._ , 2016). In the
simplest and conservative picture, the NS core is modeled as an electrically
neutral uniform fluid of neutrons, protons, electrons and muons in equilibrium
with respect to the weak interaction ($\beta$-stable nuclear matter). The
isoscalar scalar $\sigma$ meson and the isoscalar vector $\omega$ meson
mediate the long and short-range part of the interaction, respectively, in
symmetric nuclear matter, while isovector mesons (like the isovector vector
$\rho$ meson) need to be included as well to treat isospin-asymmetric matter.
One Lagrangian density based on such effective interaction can be written as:

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\overline{\psi}\\{(i\gamma_{\mu}\partial^{\mu}-m^{\ast}-\gamma_{\mu}(g_{\omega_{N}}\omega^{\mu}+\frac{1}{2}g_{\rho_{N}}\rho_{\mu}\vec{\tau})\\}\psi$
(7)
$\displaystyle+\frac{1}{2}\partial^{\mu}\sigma\partial_{\mu}\sigma-\frac{1}{2}m_{\sigma}^{2}\sigma^{2}-\frac{1}{3}g_{2}\sigma^{3}-\frac{1}{4}g_{3}\sigma^{4}$
$\displaystyle-\frac{1}{4}W^{\mu\nu}W_{\mu\nu}+\frac{1}{2}m_{\omega}^{2}\omega^{\mu}\omega_{\mu}$
$\displaystyle-\frac{1}{4}\vec{R}^{\mu\nu}\vec{R}_{\mu\nu}+\frac{1}{2}m_{\rho}^{2}\vec{\rho}^{\mu}\vec{\rho}_{\mu}$
$\displaystyle+\Lambda_{v}(g_{\omega N}^{2}\omega^{\mu}\omega_{\mu})(g_{\rho
N}^{2}\vec{\rho}^{\mu}\vec{\rho}_{\mu})\ ,$

where $\psi$ represents the wave function of nucleons. $\sigma$,
$\omega_{\mu}$, $\vec{\rho}_{\mu}$ denote the fields of $\sigma$, $\omega$,
and $\rho$ mesons, respectively. $W^{\mu\nu}$ and $\vec{R}^{\mu\nu}$ are the
anti-symmetry tensor fields of $\omega$ and $\rho$ mesons. Non-linear meson
self-interaction terms in the Lagrangian are introduced to account for medium
dependence of the effective mean-field interactions (see e.g., (Boguta and
Bodmer, 1977)). The cross-coupling from the $\omega$ meson and $\rho$ meson is
introduced to achieve a reasonable slope of symmetry energy (Horowitz and
Piekarewicz, 2001). $m_{\sigma}=510\leavevmode\nobreak\ \rm{MeV}$,
$m_{\omega}=783\leavevmode\nobreak\ \rm{MeV}$, and
$m_{\rho}=770\leavevmode\nobreak\ \rm{MeV}$ are the meson masses. $g_{\omega
N}$ and $g_{\rho N}$ are the nucleon coupling constants for $\omega$ and
$\rho$ mesons. From the simple quark counting rule, we obtain $g_{\omega
N}=3g_{\omega q}$ and $g_{\rho N}=g_{\rho q}$ from the corresponding coupling
constants with quarks. $m^{\ast}=m_{N}-g_{\sigma_{N}}\sigma$ is the (Dirac)
effective nucleon mass (see discussions in e.g., (Li  _et al._ , 2016)).
Simply, $m^{\ast}$ can be evaluated from the calculation of the confined
quarks as a function of the $\sigma$ field, $g_{\sigma N}=-\partial
m^{\ast}/\partial\sigma$, which defines the $\sigma$ coupling with nucleons
(depending on the parameter $g_{\sigma q}$). For example, by constructing a
nucleon from confined quarks with a two-body confining potential in a quark
mean-field (QMF) model (Zhu and Li, 2018), the mass of the nucleon in the
nuclear medium is expressed as the binding energy of three quarks
$E_{0}=\sum_{q}\epsilon_{q}^{\ast}$, defined by the zeroth-order term after
solving the Dirac equation of each quark inside a nucleon in a (harmonic
oscillator) confinement potential Zhu  _et al._ (2018):

$\displaystyle m^{\ast}=E_{0}-\epsilon_{\rm c.m.}+\delta m^{\pi}+(\Delta
E)_{g}\ ,$ (8)

where the effective single quark energy is given by
$\epsilon_{q}^{*}=\epsilon_{q}-g_{q\omega}\omega-\tau_{3q}g_{q\rho}\rho$, and
the effective quark mass is given by $m_{q}^{\ast}=m_{q}-g_{\sigma q}\sigma$.
In Eq. (8), we shall consider three corrections in the zeroth-order nucleon
mass in the nuclear medium, including the contribution of the center-of-mass
(c.m.) correction $\epsilon_{\rm c.m.}$, pionic correction $\delta m^{\pi}$
and gluonic correction $(\Delta E)_{g}$. The pion correction is generated by
the chiral symmetry of QCD theory and the gluon correction by the short-range
exchange interaction of quarks. See details in (Li  _et al._ , 2020; Zhu  _et
al._ , 2023a).

Following standard quantum field theory, with the Euler-Lagrange equation, the
equations of motion of the nucleon and mesons can be obtained and subsequently
solved self-consistently in terms of the mean-field approximation. Note that
the space components of the vector mesons are removed in the parity
conservation system. In addition, the spatial derivatives of nucleons and
mesons are neglected in the infinite nuclear matter due to the transformation
invariance. The energy density and pressure are then generated by the energy-
momentum tensor $T_{\mu\nu}$ and are then calculated as their ground-state
expectation values,

$\displaystyle
T_{\mu\nu}=-g_{\mu\nu}\mathcal{L}+\frac{\partial\phi_{i}}{\partial
x^{\nu}}\frac{\mathcal{L}}{\partial(\partial\phi_{i})/\partial x_{\mu}}\ ,$
(9)

where $\phi_{i}$ denotes the nucleon and various mesons.

The relativistic-field EoS model detailed above is characterized by six
parameters ($g_{\sigma q},g_{\omega q},g_{\rho q},g_{2},g_{3},\Lambda_{v}$).
To determine these parameters, we reproduce the six experimental parameters of
nuclear matter: the saturation density $n_{0}$ and the corresponding values at
the saturation point of the binding energy $E_{0}/A$, the incompressibility
$K_{0}$, the symmetry energy $E_{\rm sym}^{0}$, the symmetry energy slope
$L_{0}$ and the nucleon effective mass $m^{\ast}_{0}$. In the subsequent
section, following a methodology akin to our prior work in (Zhu  _et al._ ,
2023a, b), we vary four parameters ($K_{0},E_{\rm
sym}^{0},L_{0},m^{\ast}_{0}$), while keeping the saturation density $n_{0}$
and energy per nucleon $E_{0}/A$ are constants. Uniform prior distributions
within the empirical ranges, as presented in Table 2, are employed for these
parameters. For instance, for the symmetry energy slope $L_{0}$, besides the
two typical values of 40 and 60 MeV (see discussions in (Lattimer, 2023)) as
illustrated in Fig. 1, a broad range from 20 to 120 MeV is considered to
respect the somewhat contradictory results from PREX-II experiment (Adhikari
_et al._ , 2021) and $\chi$EFT calculation (Hu  _et al._ , 2022).

Table 2: Most probable intervals of various properties for nuclear matter saturation and phase transition (quoted as median+68% credible interval). The priors of these parameters are also listed in the second row.  | $K_{0}$ | $E_{\rm sym}^{0}$ | $L_{0}$ | $m^{*}_{0}/m_{N}$ | $n_{\rm t}/n_{0}$ | $\mu_{\rm crit}$ | $\Delta\varepsilon/\varepsilon_{\rm t}$ | $c_{s}^{2}$
---|---|---|---|---|---|---|---|---
prior | U$(220,280)$ | U$(30,35)$ | U$(20,120)$ | U$(0.65,0.8)$ | U$(1,4)$ | / | U$(0,2)$ | U$(1/3,1)$
NS (w.o. XMM) | $259.63_{-15.53}^{+13.92}$ | $32.64_{-1.77}^{+1.58}$ | $53.56_{-21.02}^{+23.05}$ | $0.77_{-0.02}^{+0.01}$ | / | / | / | /
NS (with XMM) | $265.06_{-14.85}^{+10.70}$ | $32.62_{-1.65}^{+1.61}$ | $77.50_{-21.34}^{+10.90}$ | $0.76_{-0.02}^{+0.02}$ | / | / | / | /
NS (plus HESS) | $265.11_{-15.81}^{+10.23}$ | $32.76_{-1.74}^{+1.51}$ | $38.39_{-11.41}^{+17.68}$ | $0.76_{-0.02}^{+0.02}$ | / | / | / | /
HS (w.o. XMM) | $259.85_{-18.83}^{+14.01}$ | $32.66_{-1.79}^{+1.61}$ | $62.13_{-25.00}^{+20.44}$ | $0.77_{-0.03}^{+0.02}$ | $2.28_{-0.72}^{+0.92}$ | $1078.82_{-69.80}^{+125.22}$ | $0.38_{-0.24}^{+0.35}$ | $0.77_{-0.18}^{+0.16}$
HS (with XMM) | $265.86_{-17.86}^{+10.25}$ | $32.73_{-1.77}^{+1.52}$ | $76.60_{-23.76}^{+12.79}$ | $0.75_{-0.02}^{+0.03}$ | $3.11_{-1.05}^{+0.69}$ | $1311.67_{-239.94}^{+1063.26}$ | $0.31_{-0.22}^{+0.95}$ | $0.68_{-0.21}^{+0.21}$
HS (plus HESS) | $260.99_{-19.61}^{+13.15}$ | $32.73_{-1.63}^{+1.51}$ | $41.74_{-13.27}^{+20.48}$ | $0.76_{-0.02}^{+0.03}$ | $2.38_{-1.11}^{+1.21}$ | $1101.85_{-116.53}^{+1198.91}$ | $0.37_{-0.26}^{+0.65}$ | $0.68_{-0.21}^{+0.20}$

As for NSs with a transition to deconfined quark phase in their interior (at
some transition pressure $P_{\rm t}$), same with our previous works Miao  _et
al._ (2020); Li  _et al._ (2021b), we employ the general constant-speed-of-
sound parametrization Alford  _et al._ (2013) from the transition onset up to
the maximum central pressure of a star (see detailed discussions in e.g., (Xia
_et al._ , 2021)), while keeping using QMF for the hadronic matter:

$\varepsilon(P)=\left\\{\\!\begin{array}[]{ll}\varepsilon_{\rm t}(P)&P<P_{\rm
t}\\\ \varepsilon_{\rm t}(P_{\rm t})+\Delta\varepsilon+c_{\rm s}^{-2}(P-P_{\rm
t})&P>P_{\rm t}\end{array}\right.\ \\\ $

In addition to the four EoS parameters used in the QMF model, we introduce
three additional parameters for the quark matter phase, i.e., the transition
density $n_{\rm t}/n_{0}$, the transition strength
$\Delta\varepsilon/\varepsilon_{\rm t}$ and the sound speed squared in quark
matter $c_{\rm s}^{2}$. The priors for these parameters are chosen based on
Ref. (Li  _et al._ , 2021b) and are presented in Table 2. Note that the value
of the sound speed can be varied from $1/\sqrt{3}$ (the perturbative QCD
value) to 1 (the casual limit).

From Fig. 1, it is evident that the radius of a $1.4\,M_{\odot}$ star, denoted
as $R_{1.4}$, is subject to nuclear matter properties, specifically $L_{0}$ at
the saturation point; The EoS models QMF40 and QMF60, discussed previously,
belong to a set of NS EoSs constructed in Ref. Zhu  _et al._ (2018), enabling
fine-tuning of the symmetry energy slope $L_{0}$ and a systematic exploration
of the $R$ (or $\Lambda)$ versus $L$ relationship. with well-reproduced robust
observables from laboratory nucleons, nuclear saturation, heavy-ion collisions
(HIC), and heavy pulsars. Rather, the maximum mass and the possibility of a
strangeness phase transition (e.g., involving hyperons or deconfined quarks)
are associated with the EoS at supersaturation densities. The XMM-Newton data
of PSR J0030+0451, in conjunction with the NICER data, can shed light on the
EoS at low density (manifested by $L_{0}$) at $\sim n_{0}$, intermediate
density (manifested by $R_{\rm 1.4}$) at $\sim 2\,n_{0}$ and high density
(manifested by the maximum mass $M_{\rm TOV}$) at $\sim 5\,n_{0}$. On the
other hand, the mass and radius data of HESS J1731-347 should primarily
complement other EoS constraints at low and intermediate densities at $\sim
1$-$2\,n_{0}$. It is expected that the stiffening effect introduced by the
XMM-Newton data and the softening effect introduced by the HESS J1731-347 data
will lead to a compromise, resulting in a flattening of the credible ranges
for both global NS properties and microscopic EoS parameters.

##  IV Bayesian inference of EoS

In this section, we present the outcomes of our Bayesian inference regarding
NS andEoS properties. Our analysis incorporates astronomical observations
within the Bayesian framework detailed in (Zhu  _et al._ , 2023a) and (Li  _et
al._ , 2021b). We report stellar properties in Figs. 2, 4 and Table 1.
VariousEoS properties are detailed in Figs. 3, 5, 6, and Table 2. We discuss
primarily the impact of the updated mass and radius analysis for PSR
J0030+0451 data and explore as well how different measurements of NSs with low
and high masses jointly narrow theEoS parameter space. The current quantitive
results (at some confidence level) of theEoS at supersaturation density are
also confronted with the extrapolation of the microscopic calculations of the
EoS e.g., in the $\chi$EFT method, and general improvement is achieved (see
below in Fig. 5).

Analyzing Fig. 2, it becomes apparent that a tension exists for $R_{1.4}$
between the constraints derived from the NICER sources and HESS J1731-347.
This tension appears more pronounced than the one observed between the NICER
sources and GW17817, a notable point frequently discussed in existing
literature. The incorporation of XMM-Newton data for the masses and radii of
two NICER sources further intensifies this tension. As a consequence, we
anticipate that the uncertain ranges of $R_{1.4}$ (as well as $L_{0}$) will
experience effective narrowing through the combination of multimessenger data,
incorporating radio, gravitational wave, and x-ray observations.

In Fig. 3, we present the EoS and sound speed for both NSs (left panels) and
HSs (right panels). The updated radius measurements, incorporating XMM-Newton
data, enhance the precision of the EoS inference by tightening the constraint
on the pressure at densities around $\sim 1$-$2\,n_{0}$ ($\sim 1$-$5\,n_{0}$)
for NSs (HSs). EoS parameter spaces are much enlarged when considering
possible strangeness phase transitions and large uncertainties reside
especially at high densities above $\sim 2\,n_{0}$, leading to unresolved
predictions for the maximum mass ($M_{\rm TOV}$) and the sound speed at high-
density cores of NSs. For instance, the inclusion of XMM-Newton data for the
masses and radii of two NICER sources raises $M_{\rm TOV}$ from
$2.09\,M_{\odot}$ ($2.15\,M_{\odot}$) to $2.13\,M_{\odot}$ ($2.21\,M_{\odot}$)
in the case of pure NS (HS), as depicted in Fig. 2. If we further incorporate
HESS data, we observe an enlargement of the EoS parameter space for both NSs
and HSs due to the compromising effects mentioned earlier. In Fig. 4, the
compromising effects introduced by the XMM data (comparing dash-dotted curves
with dashed ones) and the HESS data (comparing dashed curves with solid ones)
become more evident. Consequently, $R_{1.4}$ is refined to a much narrower
range around $12.11\;\text{km}$ ($12.20\;\text{km}$) for pure NSs (HSs). A
similar trend is observed for $\Lambda_{1.4}$. Nevertheless, the posterior PDF
of $M_{\rm TOV}$ exhibits flattening toward higher values. For instance, the
$1\sigma$ upper bound increases from $2.31\,M_{\odot}$ with only the NICER
data to over $\sim 3.3\,M_{\odot}$ when adding the XMM and/or HESS data (see
Table 1).

From Fig. 5, it is seen that the updated mass and radius measurements,
incorporating XMM-Newton data, do not inform the EoS within the applicable
regime of $\chi\text{EFT}$ (i.e., below $n_{0}$), and effectively lift the
lower bounds of the pressure above $n_{0}$ in both cases of pure NSs and HSs.
Consequently, by combining the information of $\chi\text{EFT}$ at
subsaturation densities with multiple observations of radio, gravitational
wave, and X-ray data at high densities, the constraint on the pressure at
densities $\sim 1$-$3\,n_{0}$ is effectively tightened. The posterior
+GW170817+NICER (w.o.XMM) is comparable to the seventh column of Table II in
Ref. (Legred  _et al._ , 2021), where the mass and radius constraints for PSR
J0030+0451 and J0740+6620 from NICER are utilized, but without the updated
radius measurement of J0030+0451 incorporating the XMM-Newton data. The
current analysis with NICER$\times$XMN for light PSR J0030+0451 sets a much
narrower limit on the stiffness at $\sim 2\,n_{0}$, evident when comparing
$P_{2n_{0}}=13.4$-$30.4\,\text{MeV/fm}^{3}$ from Ref. (Legred  _et al._ ,
2021) with the present $P_{2n_{0}}=22.0$-$27.1\,\text{MeV/fm}^{3}$, both at
the $1\sigma$ credible level, in the left panel of the NS case.

As shown in Table 1, the $1\sigma$ credible intervals of HSs are much larger
than those of pure NSs. Similar results for $R_{1.4}$ (around $12\;\text{km}$)
are found for both pure NSs and HSs, as revealed by a consistent framework
(Zhu  _et al._ , 2018) for dealing with single nucleons and nuclear many-body
systems from the quark level. Many previous combined analyses of
multimessenger observations also support this consistency (e.g., in (Capano
_et al._ , 2020)). However, there remains relatively large uncertainty for
$M_{\rm TOV}$ in the range of $\sim 2.0$-$3.3\,M_{\odot}$.

From Fig. 6, we learn that $K_{0}$ is more sensitive to the XMM-Newton data,
while $L_{0}$ is more sensitive to both the XMM-Newton data and the HESS data.
The HESS data has a significant modulation effect on $L_{0}$ (correspondingly
to $R_{1.4}$ and $\Lambda_{1.4}$ as shown in Fig. 4) in both cases of pure NSs
and HSs, as well as on $n_{\rm t}$ for HSs. While the saturation property
$L_{0}$ benefits from such compensation, yielding a value of $41.7$
($38.4\,\text{MeV}$) in the case of pure NS (HS), a conflict arises for
$n_{\rm t}$ in the range of $1.08$-$3.98\,n_{0}$ at the $1\sigma$ credible
level; Similarly, the predicted $1\sigma$ credible interval for the sound
speed squared is possibly in the range of $0.39$-$0.97$ at high-density cores
of NSs (See Table 2). The influence of individual mass and radius constraints
on the phase transition parameters exhibits similarities to our previous
investigations (Miao  _et al._ , 2020; Li  _et al._ , 2021b). See discussions
also in Ref. (Drischler  _et al._ , 2022). In particular, consistent
predictions are obtained for the saturation properties $K_{0}$, $E_{\rm
sym}^{0}$, and $m^{*}_{0}$ for the cases with or without strangeness phase
transition, with slightly larger ranges for those in the case of HSs, as
detailed in Table 2. When accounting for the possibility of a phase
transition, the resulting $L_{0}$ is around $41.7\,\text{MeV}$, close to the
EoS (Zhu  _et al._ , 2018) from the quark level previously constructed in the
QMF framework.

##  V Summary and perspective

NSs provide a distinctive environment for exploring fundamental nuclear
physics. It offers a unique opportunity to study EoS of cold matters in
equilibrium, whose density is a few orders of magnitude higher than that of
typical nuclear matters on Earth. The pulse-profile modeling method has been
successfully applied to infer the mass and radius of NSs conditional on the
thermal emission of hot spots on the stellar surface. The latest mass and
radius inferences of PSR J0030+0451 and PSR J0740+6620 are obtained from the
joint analysis of the NICER and XMM-Newton data. It provides a preliminary
constraint on the EoS together with the observations of the supernova remnant
HESS J1731-347 and gravitational waves of GW170817. As expected, when
incorporating the updated radius measurement of light PSR J0030+0451
integrating XMM-Newton data, the present analysis establishes a refined
constraint on the stiffness at $\sim 2\,n_{0}$. In particular, to compromise
the different results regarding the radius of light NSs around
$1.4\,M_{\odot}$ between the HESS J1731-347 data and the updated PSR
J0030+0451 data, the previously unknown isospin-asymmetric contribution of the
EoS can be better determined. For example, the symmetry energy slope tends to
be a low value of around $40\,\text{MeV}$, even with the inclusion of possible
non-nucleon degrees of freedom. It emphasizes the need for collaborative data
analysis and methodological refinements.

In the current work, seven parameters collaboratively determine the EoS prior,
ach parameter carrying specific physical significance associated with the
composition and phase state of NS matter. This framework allows for a direct
connection between NS observations and microphysics-driven investigations of
the EoS across various density regimes and facilitates a connection with the
ongoing research efforts in the field of relativistic HIC (Li  _et al._ ,
2023). Moreover, the versatility of this analysis allows for straightforward
extensions to higher-dimensional models, accommodating the inclusion of
additional particles within NSs, such as hyperons (Sun  _et al._ , 2023) and
dark matter (Miao  _et al._ , 2022a), to contribute a more comprehensive
understanding of the interplay between different phase states of NS matter.
While establishing a direct connection between QCD and NS EoS remains
challenging, the current framework offers a pragmatic approach to translating
NS observations into the EoS microphysics. This approach supports quantitative
studies of the EoS at different density regimes, inferring the values of
unresolved properties of nuclear saturation ($K_{0},E_{\rm
sym}^{0},L_{0},m_{0}^{*}$) and deconfinement phase transition ($n_{\rm
t}/n_{0},\Delta\epsilon/\epsilon_{\rm t},c_{s}^{2},\mu_{\rm crit}$),
particularly with the increasing accumulation of important NS observables,
such as mass and radius.

Nevertheless, the uncertainties of the background estimation and the emission
configurations of hot regions limit the accuracy of the inferred mass and
radius of NSs. There is a possibility that the current radius measurement,
contingent on the XMM data set, which provides an indirect NICER background
constraint, may face potential advancements in the future (Riley  _et al._ ,
2021). Addressing non-thermal emission components, exemplified by cases like
PSR J0437-4715, is crucial (Bogdanov, 2013; Guillot  _et al._ , 2016;
González-Caniulef  _et al._ , 2019). Consequently, achieving more accurate
measurements requires the future X-ray telescopes with large effective area,
high energy resolution, and high time resolution in the soft X-ray band, as
well as adequate imaging capability to discriminate the background, e.g. the
planned enhanced X-ray Timing and Polarimetry Mission (eXTP) Watts  _et al._
(2019) and the Spectroscopic Time-Resolving Observatory for Broadband Energy
X-rays (STROBE-X) (Ray  _et al._ , 2019). The follow-up X-ray Telescope (FXT),
one of the payloads aboard the recently launched Einstein Probe (EP)
mission333https://ep.bao.ac.cn/ep/ has high space and time resolution (better
than 100 $\mu$s (Cui  _et al._ , 2023)), of which the former property is
capable of providing accurate background estimation and the latter could help
recognize the timing evolution of the CCO and obtain a better pulse profile if
any. The long-term monitoring of EP on CCOs and PSRs will afford more
observation constraints on the EoS. In addition, given the thermal evolution
of young isolated CCOs, such as HESS J1731-347, is primarily governed by
neutrino emission, the monitoring of their X-ray emission serves as a valuable
avenue for acquiring complementary insights into the stellar interior
(Beznogov  _et al._ , 2018). With expectations of improved measurements in the
near future, ongoing endeavors aim to enhance our understanding of the inner
structure of NSs Miao  _et al._ (2022b); Ascenzi  _et al._ (2024).

###### Acknowledgements.

We are thankful to Profs. J. B. Wei, X. Zhou and the XMU neutron star group
for helpful discussions. The work is supported by National SKA Program of
China (No. 2020SKA0120300), National Natural Science Foundation of China
(grant Nos. 12273028, 12373051 and 12333007), and International Partnership
Program of Chinese Academy of Sciences (Grant No.113111KYSB20190020).

## References

  * Huber (2020) M. Q. Huber, Physics Reports 879, 1 (2020).
  * Li (2022) A. Li, in  _European Physical Journal Web of Conferences_, European Physical Journal Web of Conferences, Vol. 260 (2022) p. 04001.
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, _et al._ , Reports on Progress in Physics 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Li  _et al._ (2020) A. Li, Z. Y. Zhu, E. P. Zhou, _et al._ , Journal of High Energy Astrophysics 28, 19 (2020), arXiv:2007.05116 [nucl-th] .
  * Burgio  _et al._ (2021) G. F. Burgio, H. J. Schulze, I. Vidaña, and J. B. Wei, Progress in Particle and Nuclear Physics 120, 103879 (2021), arXiv:2105.03747 [nucl-th] .
  * Lattimer (2021) J. M. Lattimer, Annual Review of Nuclear and Particle Science 71, 433 (2021).
  * Li  _et al._ (2021a) A. Li, Z. Q. Miao, J. L. Jiang, _et al._ , Mon. Not. R. Astron. Soc. 506, 5916 (2021a), arXiv:2009.12571 [astro-ph.HE] .
  * Watts  _et al._ (2016) A. L. Watts, N. Andersson, D. Chakrabarty, _et al._ , Reviews of Modern Physics 88, 021001 (2016), arXiv:1602.01081 [astro-ph.HE] .
  * Bogdanov  _et al._ (2019) S. Bogdanov, F. K. Lamb, S. Mahmoodifar, _et al._ , Astrophys. J. Lett. 887, L26 (2019), arXiv:1912.05707 [astro-ph.HE] .
  * Watts (2012) A. L. Watts, Annual Review of Astronomy & Astrophysics 50, 609 (2012), arXiv:1203.2065 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley, A. L. Watts, S. Bogdanov, _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller, F. K. Lamb, A. J. Dittmann, _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley, A. L. Watts, P. S. Ray, _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Miller and Lamb (2015) M. C. Miller and F. K. Lamb, Astrophys. J. 808, 31 (2015), arXiv:1407.2579 [astro-ph.HE] .
  * Miller (2016) M. C. Miller, Astrophys. J. 822, 27 (2016), arXiv:1602.00312 [astro-ph.HE] .
  * Greif  _et al._ (2019) S. K. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, and A. L. Watts, Mon. Not. R. Astron. Soc. 485, 5363 (2019), arXiv:1812.08188 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2021) G. Raaijmakers, S. K. Greif, K. Hebeler, _et al._ , Astrophys. J. Lett. 918, L29 (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Vinciguerra  _et al._ (2024) S. Vinciguerra, T. Salmi, A. L. Watts, _et al._ , Astrophys. J. 961, 62 (2024), arXiv:2308.09469 [astro-ph.HE] .
  * Salmi  _et al._ (2022) T. Salmi, S. Vinciguerra, D. Choudhury, _et al._ , Astrophys. J. 941, 150 (2022), arXiv:2209.12840 [astro-ph.HE] .
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Negele and Vautherin (1973) J. W. Negele and D. Vautherin, Nuclear Physics A 207, 298 (1973).
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ , Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ , Physical Review X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Doroshenko  _et al._ (2022) V. Doroshenko, V. Suleimanov, G. Pühlhofer, and A. Santangelo, Nature Astronomy 6, 1444 (2022).
  * Yuan and Li (2023) W.-L. Yuan and A. Li, arXiv e-prints , arXiv:2312.17102 (2023), arXiv:2312.17102 [nucl-th] .
  * Fonseca  _et al._ (2021) E. Fonseca, H. T. Cromartie, T. T. Pennucci, P. S. Ray, _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Zhu  _et al._ (2018) Z.-Y. Zhu, E.-P. Zhou, and A. Li, Astrophys. J. 862, 98 (2018), arXiv:1802.05510 [nucl-th] .
  * Liu  _et al._ (2023) X.-J. Liu, Z.-Q. You, Z.-C. Chen, _et al._ , arXiv e-prints , arXiv:2312.12185 (2023), arXiv:2312.12185 [astro-ph.HE] .
  * Alford and Halpern (2023) J. A. J. Alford and J. P. Halpern, Astrophys. J. 944, 36 (2023), arXiv:2302.05893 [astro-ph.HE] .
  * Klochkov  _et al._ (2015) D. Klochkov, V. Suleimanov, G. Pühlhofer, _et al._ , Astron. Astrophys. 573, A53 (2015), arXiv:1410.1055 [astro-ph.HE] .
  * Doroshenko  _et al._ (2016) V. Doroshenko, G. Pühlhofer, P. Kavanagh, _et al._ , Mon. Not. R. Astron. Soc. 458, 2565 (2016), arXiv:1508.03557 [astro-ph.HE] .
  * Miao  _et al._ (2023) Z. Miao, E. Zhou, and A. Li, arXiv e-prints , arXiv:2305.08401 (2023), arXiv:2305.08401 [nucl-th] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ , Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Hernandez Vivanco  _et al._ (2020) F. Hernandez Vivanco, R. Smith, E. Thrane, and P. D. Lasky, Mon. Not. R. Astron. Soc. 499, 5972 (2020), arXiv:2008.05627 [astro-ph.HE] .
  * Zhou  _et al._ (2021) X. Zhou, A. Li, and B.-A. Li, Astrophys. J. 910, 62 (2021), arXiv:2011.11934 [astro-ph.HE] .
  * Haensel  _et al._ (2009) P. Haensel, J. L. Zdunik, M. Bejger, and J. M. Lattimer, Astron. Astrophys. 502, 605 (2009), arXiv:0901.1268 [astro-ph.SR] .
  * Wiringa  _et al._ (1995) R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, Phys. Rev. C 51, 38 (1995), arXiv:nucl-th/9408016 [nucl-th] .
  * Brockmann and Machleidt (1990) R. Brockmann and R. Machleidt, Phys. Rev. C 42, 1965 (1990).
  * Beane  _et al._ (2011) S. R. Beane, W. Detmold, K. Orginos, and M. J. Savage, Progress in Particle and Nuclear Physics 66, 1 (2011), arXiv:1004.2935 [hep-lat] .
  * Aoki  _et al._ (2012) S. Aoki, T. Doi, T. Hatsuda, _et al._ , Progress of Theoretical and Experimental Physics 2012, 01A105 (2012), arXiv:1206.5088 [hep-lat] .
  * Weinberg (1990) S. Weinberg, Physics Letters B 251, 288 (1990).
  * Weinberg (1991) S. Weinberg, Nuclear Physics B 363, 3 (1991).
  * Hebeler  _et al._ (2011) K. Hebeler, S. K. Bogner, R. J. Furnstahl, _et al._ , Phys. Rev. C 83, 031301 (2011), arXiv:1012.3381 [nucl-th] .
  * Drischler  _et al._ (2020a) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020a), arXiv:2004.07232 [nucl-th] .
  * Drischler  _et al._ (2020b) C. Drischler, J. A. Melendez, R. J. Furnstahl, and D. R. Phillips, Phys. Rev. C 102, 054315 (2020b), arXiv:2004.07805 [nucl-th] .
  * Li  _et al._ (2019) A. Li, J. Hu, S. Bao, _et al._ , Nuclear Physics Review 36, 1 (2019).
  * Drischler  _et al._ (2021) C. Drischler, S. Han, J. M. Lattimer, _et al._ , Phys. Rev. C 103, 045808 (2021), arXiv:2009.06441 [nucl-th] .
  * Legred  _et al._ (2021) I. Legred, K. Chatziioannou, R. Essick, _et al._ , Phys. Rev. D 104, 063003 (2021), arXiv:2106.05313 [astro-ph.HE] .
  * Walecka (1974) J. D. Walecka, Annals of Physics 83, 491 (1974).
  * Serot (1992) B. D. Serot, Reports on Progress in Physics 55, 1855 (1992).
  * Li  _et al._ (2008) A. Li, G. F. Burgio, U. Lombardo, and G. X. Peng, International Journal of Modern Physics E 17, 1635 (2008).
  * Zhu  _et al._ (2016) Z.-Y. Zhu, A. Li, J.-N. Hu, and H. Sagawa, Phys. Rev. C 94, 045803 (2016), arXiv:1607.04007 [nucl-th] .
  * Boguta and Bodmer (1977) J. Boguta and A. R. Bodmer, Nuclear Physics A 292, 413 (1977).
  * Horowitz and Piekarewicz (2001) C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001), arXiv:astro-ph/0010227 [astro-ph] .
  * Li  _et al._ (2016) A. Li, J. N. Hu, X. L. Shang, and W. Zuo, Phys. Rev. C 93, 015803 (2016), arXiv:1511.07979 [nucl-th] .
  * Zhu and Li (2018) Z.-Y. Zhu and A. Li, Phys. Rev. C 97, 035805 (2018), arXiv:1802.07441 [nucl-th] .
  * Zhu  _et al._ (2023a) Z. Zhu, A. Li, J. Hu, and H. Shen, Phys. Rev. C 108, 025809 (2023a), arXiv:2305.16058 [nucl-th] .
  * Zhu  _et al._ (2023b) Z. Zhu, A. Li, and T. Liu, Astrophys. J. 943, 163 (2023b), arXiv:2211.02007 [astro-ph.HE] .
  * Lattimer (2023) J. M. Lattimer, Particles 6, 30 (2023), arXiv:2301.03666 [nucl-th] .
  * Adhikari  _et al._ (2021) D. Adhikari, H. Albataineh, D. Androic, K. Aniol, _et al._ , Phys. Rev. Lett. 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
  * Hu  _et al._ (2022) B. Hu, W. Jiang, T. Miyagi, _et al._ , Nature Physics 18, 1196 (2022), arXiv:2112.01125 [nucl-th] .
  * Miao  _et al._ (2020) Z. Miao, A. Li, Z. Zhu, and S. Han, Astrophys. J. 904, 103 (2020), arXiv:2006.00839 [nucl-th] .
  * Li  _et al._ (2021b) A. Li, Z. Miao, S. Han, and B. Zhang, Astrophys. J. 913, 27 (2021b), arXiv:2103.15119 [astro-ph.HE] .
  * Alford  _et al._ (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Xia  _et al._ (2021) C. Xia, Z. Zhu, X. Zhou, and A. Li, Chinese Physics C 45, 055104 (2021), arXiv:1906.00826 [nucl-th] .
  * Capano  _et al._ (2020) C. D. Capano, I. Tews, S. M. Brown, _et al._ , Nature Astronomy 4, 625 (2020), arXiv:1908.10352 [astro-ph.HE] .
  * Drischler  _et al._ (2022) C. Drischler, S. Han, and S. Reddy, Phys. Rev. C 105, 035808 (2022), arXiv:2110.14896 [nucl-th] .
  * Li  _et al._ (2023) A. Li, G.-C. Yong, and Y.-X. Zhang, Phys. Rev. D 107, 043005 (2023), arXiv:2211.04978 [nucl-th] .
  * Sun  _et al._ (2023) X. Sun, Z. Miao, B. Sun, and A. Li, Astrophys. J. 942, 55 (2023), arXiv:2205.10631 [astro-ph.HE] .
  * Miao  _et al._ (2022a) Z. Miao, Y. Zhu, A. Li, and F. Huang, Astrophys. J. 936, 69 (2022a), arXiv:2204.05560 [astro-ph.HE] .
  * Bogdanov (2013) S. Bogdanov, Astrophys. J. 762, 96 (2013), arXiv:1211.6113 [astro-ph.HE] .
  * Guillot  _et al._ (2016) S. Guillot, V. M. Kaspi, R. F. Archibald, _et al._ , Mon. Not. R. Astron. Soc. 463, 2612 (2016), arXiv:1512.03957 [astro-ph.HE] .
  * González-Caniulef  _et al._ (2019) D. González-Caniulef, S. Guillot, and A. Reisenegger, Mon. Not. R. Astron. Soc. 490, 5848 (2019), arXiv:1904.12114 [astro-ph.HE] .
  * Watts  _et al._ (2019) A. L. Watts, W. Yu, J. Poutanen, _et al._ , Science China Physics, Mechanics, and Astronomy 62, 29503 (2019), arXiv:1812.04021 [astro-ph.HE] .
  * Ray  _et al._ (2019) P. S. Ray, Z. Arzoumanian, D. Ballantyne, _et al._ , arXiv e-prints , arXiv:1903.03035 (2019), arXiv:1903.03035 [astro-ph.IM] .
  * Cui  _et al._ (2023) W. Cui, H. Wang, X. Zhao, _et al._ , Experimental Astronomy 55, 603 (2023).
  * Beznogov  _et al._ (2018) M. V. Beznogov, E. Rrapaj, D. Page, and S. Reddy, Phys. Rev. C 98, 035802 (2018), arXiv:1806.07991 [astro-ph.HE] .
  * Miao  _et al._ (2022b) Z. Miao, A. Li, and Z.-G. Dai, Mon. Not. R. Astron. Soc. 515, 5071 (2022b), arXiv:2107.07979 [astro-ph.HE] .
  * Ascenzi  _et al._ (2024) S. Ascenzi, V. Graber, and N. Rea, arXiv e-prints , arXiv:2401.14930 (2024), arXiv:2401.14930 [astro-ph.HE] .