# Hybrid star models in the light of new multi-messenger data

Jia Jie Li1, Armen Sedrakian2,3, Mark Alford4 1School of Physical Science and
Technology, Southwest University, Chongqing 400715, China<EMAIL_ADDRESS>
2Frankfurt Institute for Advanced Studies, D-60438 Frankfurt am Main, Germany;
<EMAIL_ADDRESS>
3Institute of Theoretical Physics, University of Wroclaw, 50-204 Wroclaw,
Poland
4Department of Physics, Washington University, St. Louis, MO 63130, USA;
<EMAIL_ADDRESS>

###### Abstract

Recent astrophysical mass inferences of compact stars (CSs) HESS J1731-347 and
PSR J0952-0607, with extremely small and large masses respectively, as well as
the measurement of the neutron skin of Ca in the CREX experiment challenge and
constrain the models of dense matter. In this work, we examine the concept of
hybrid stars - objects containing quark cores surrounded by nucleonic
envelopes - as models that account for these new data along with other
(multimessenger) inferences. We employ a family of 81 nucleonic equations of
state (EoSs) based on covariant density functional with variable skewness and
slope of symmetry energy at saturation density and a constant speed-of-sound
EoS for quark matter. For each nucleonic EoS, a family of hybrid star EoS is
generated by varying the transition density from nucleonic to quark matter,
the density jump at the transition, and the speed-of-sound. These models are
tested against the data from GW170817 and J1731-347, which favor low-density
soft EoS and PSR J0592-0607 and J0740+6620, which require high-density stiff
EoS. We then examine the occurrence of twin configurations and quantify the
ranges of masses and radii that they can possess. It is shown that including
J1731-347 data favors EoS models which predict low-mass twin stars with masses
$M\lesssim 1.3\,M_{\odot}$ that can be realized if the deconfinement
transition density is low. If combined with a large speed-of-sound in quark
matter such models allow for maximum masses of hybrid stars in the range
$2.0$-$2.6\,M_{\odot}$.

###### Subject headings:

Compact objects (228); Neutron stars (1108); Nuclear astrophysics (1129); High
energy astrophysics (739)

##  1\. Introduction

Although the concepts of quark stars and hybrid stars containing a quark core
enclosed in a nucleonic envelope were proposed long ago and have been
intensively studied over the following decades (for reviews, see Alford et al.
(2008); Anglani et al. (2014); Baym et al. (2018); Sedrakian (2023), they
remain at the forefront of the exploration of superdense matter and compact
stars (CSs). They may arise through the onset of a first-order phase
transition between the nucleonic and quark phases, in which case a new branch
of stellar equilibria populated by stable hybrid stars arise. One intriguing
aspect of hybrid stars is the existence of twin configurations, i.e, two stars
that possess the same mass but different radii, whereby the larger star
consists entirely of nucleonic matter, while the more compact star is a hybrid
star (Glendenning & Kettner, 2000; Zdunik & Haensel, 2013; Benic et al., 2015;
Alford & Sedrakian, 2017; Alvarez-Castillo et al., 2019; Blaschke et al.,
2020; Christian & Schaffner-Bielich, 2022; Li et al., 2023a). The study of
hybrid stars and their twin configurations continues to be an active area of
research, with ongoing efforts to explore their properties and to investigate
potential observational signatures that could distinguish them from hadronic
CSs (see Alford et al., 2005; Bonanno & Sedrakian, 2012; Masuda et al., 2013;
Klähn et al., 2013; Christian et al., 2019; Malfatti et al., 2020; Li et al.,
2020; Tan et al., 2020, 2022; Brodie & Haber, 2023; Tsaloukidis et al., 2023,
and references therein).

Figure 1.—  Mass-radius relations for hybrid EoSs with $\rho_{\rm
tr}=1.5\,\rho_{0}$ and $2.0\,\rho_{0}$, for $c^{2}_{s}=1$ and $2/3$ in the
quark phase. In each panel, the hybrid EoSs are built upon nucleonic models
grouped by $L_{\rm sym}=20$, 40, 60, 80 and 100 MeV; for each, $Q_{\rm sat}$
is set to -600 and 1000 MeV. The circle on each curve denotes the
configuration with central density $\rho_{\rm tr}$; lines ending with a filled
circle are for the isoscalar-soft nucleonic model with $Q_{\rm sat}=-600$ MeV
and lines ending with an empty circle are for the isoscalar-stiff one with
$Q_{\rm sat}=1000$ MeV. Observational constraints from multi-messenger
astronomy are shown for comparison (note there are two overlapping regions for
HESS J1731-347 (Doroshenko et al., 2022)). For each nucleonic EoS, thick lines
show the hybrid models with a maximum energy jump that yields twin
configurations and a branch passing through the 95% confidence region for the
mass-radius constraints for PSR J0740+6620 and J0030+0451. The thin lines show
the model with a critical value of $\Delta\varepsilon$ for which a higher
value leads a disconnected mass-radius curve. Solid lines represent models
that are satisfying all constraints, while dash-dotted lines show those
failing for J1731-347.

In recent years, observations of CSs have reached unprecedented levels of
precision, providing valuable constraints on the EoS of stellar matter and its
composition. Among these constraints, the most interesting ones are depicted
in Figure 1. These include:
(a) The most massive pulsar PSR J0740+6620 (Fonseca et al., 2021), which has
pushed the maximum mass of CSs to $2.0\,M_{\odot}$.
(b) The simultaneous mass and radius measurements for PSR J0030+0451 (Riley et
al., 2019; Miller et al., 2019) and J0740+6620 (Riley et al., 2021; Miller et
al., 2021) extracted from X-ray observations using the NICER instrument.
(c) The detection of gravitational waves from the binary CS merger event
GW170817 (Abbott et al., 2019), which provides us with an estimate of the
tidal deformability of the binary that can further yield constraints on the
radii of the individual components.

In addition to these well-established constraints, some more speculative
results have been reported very recently:
(d) The companion of the “black widow” pulsar PSR J0952-0607 (the second
fastest known spinning CS with a frequency of 707 Hz) has been detected in the
Milky Way Galaxy, with a mass estimate of $M=2.35^{+0.17}_{-0.17}\,M_{\odot}$
(68% credible interval) (Romani et al., 2022).
(e) The X-ray spectrum of the central compact object within the supernova
remnant HESS J1731-347 was modeled and in combination with the distance
estimates from Gaia observations yielded a very low mass of
$M=0.77^{+0.20}_{-0.17}\,M_{\odot}$ and a relatively small radius of
$R=10.4^{+0.86}_{-0.78}$ km (68% credible interval) (Doroshenko et al., 2022).

Given these new (albeit somewhat speculative) observations, it is worthwhile
to take a fresh look at the concept of a hybrid star by allowing for a broad
range of admissible parameter space within models that have been tested on a
more limited set of observations. In the present work, we implement such a
program for the construction of hybrid EoSs and apply filtering to extract the
range of EoS that is consistent with the current data. Our further focus is on
twin stars and their compatibility with the aforementioned multi-messenger
information after the inclusion of extreme domains of mass mentioned above in
the analysis.

##  2\. Hybrid models

To describe low-density matter, we use a covariant density functional of
nuclear matter in which baryons are coupled to mesons with density-dependent
couplings (Typel & Wolter, 1999; Oertel et al., 2017; Sedrakian et al., 2023).
The nuclear matter properties are conveniently assessed using the coefficients
of Taylor expansion of energy density with respect to the baryon density and
isospin asymmetry at the saturation density $\rho_{0}$. The low-order
coefficients, the saturation energy $E_{\rm sat}$, compressibility $K_{\rm
sat}$, and symmetry energy $J_{\rm sym}$ are either strongly constrained or
have no noticeable impact on the gross properties (e.g., mass, radius,
deformability, etc.) of CSs (Margueron et al., 2018; Li & Sedrakian, 2019).

Figure 2.—  Mass-radius relations for hybrid EoSs with $\rho_{\rm tr}\simeq
2.5\,\rho_{0}$ and $3.0\,\rho_{0}$, for $c^{2}_{s}=1$ and $2/3$ in the quark
phase. In each panel, the hybrid EoSs are constructed from nucleonic models
grouped by $L_{\rm sym}=20$, 40, 60, 80 and 100 MeV; for each, $Q_{\rm sat}$
is chosen such that the resultant mass-radius curve passes through the 95%
confidence region for HESS J1731-347 (solid lines) or GW170817 (dash-dotted
lines), respectively. The circle on each curve denotes the configurations with
central density $\rho_{\rm tr}$. Thick lines show the hybrid models with a
maximum $\Delta\varepsilon$ that yield twin configurations and a branch
passing through the 95% confidence region for PSR J0740+6620. Thin lines show
the models with a connected mass-radius curve with a critical value of
$\Delta\varepsilon$ set by the constraint for J0740+6620 or appearance of twin
configurations. The inset illustrates hybrid EoSs built upon an isoscalar-
stiff nucleonic model for the case $M^{\rm Q}_{\rm max}\leq M^{\rm H}_{\rm
max}$ which yields ultra-compact configurations for the case $c^{2}_{s}=1$.
Solid lines represent models satisfying all constraints, while dash-dotted
lines show those failing for J1731-347.

Here we adopt a family of nucleonic EoS consisting of 81 EoS which share the
same low-order coefficients as the original DDME2 parameterization (Lalazissis
et al., 2005), but have alternative density-dependence of the couplings, thus,
resulting in different high-order coefficients (Li & Sedrakian, 2023).
Specifically, the skewness coefficient for symmetric matter (characterizing
the isoscalar sector) varies in the range $-600\leq Q_{\rm sat}\leq 1000$ MeV
and the slope parameter of the symmetry energy (characterizing the isovector
sector) varies in the range $20\leq L_{\rm sym}\leq 100$ MeV. Note that values
of $Q_{\rm sat}$ and $L_{\rm sym}$ control, respectively, the high- and
intermediate-density behaviors of the nucleonic EoS, thereby influencing the
maximum mass of static nucleonic stars and the radius of the canonical-mass
star (Margueron et al., 2018; Li & Sedrakian, 2019). The value of $L_{\rm
sym}$, commonly considered to be close to 60 MeV, has been challenged by the
recent analyses of the neutral weak form factor of 48Ca in the CREX
experiment, which suggests rather low values of $L_{\rm sym}=20\pm 30$ MeV
(Adhikari et al., 2022; Lattimer, 2023). This can be contrasted to the values
of $L_{\rm sym}=106\pm 37$ MeV obtained by some groups via the analysis of the
neutron skin thickness of 208Pb in the PREX-II experiment (Adhikari et al.,
2021).

Figure 3.—  EoS collection utilized in the work. Solid thick lines represent
models that satisfy all constraints, while dashed thin lines show those that
fail for HESS J1731-347. The blue regions display the allowed region of EoS
from Annala et al. (2022), which includes EoS mimicking first-order phase
transitions. Insets show EoS that give rise to ultra-compact stars.

The quark phase is described by the constant sound speed (CSS) parametrization
(Zdunik & Haensel, 2013; Alford et al., 2013)

$\displaystyle p(\varepsilon)=\left\\{\begin{array}[]{ll}p_{\rm
tr},&\varepsilon_{\rm tr}<\varepsilon<\varepsilon_{\rm
tr}\\!+\\!\Delta\varepsilon,\\\\[2.15277pt] p_{\rm
tr}+c^{2}_{s}\,\bigl{[}\varepsilon-(\varepsilon_{\rm
tr}\\!+\\!\Delta\varepsilon)\bigr{]},&\varepsilon_{\rm
tr}\\!+\\!\Delta\varepsilon<\varepsilon,\end{array}\right.$ (3)

where $p_{\rm tr}$ is the pressure and $\varepsilon_{\rm tr}$ is the energy
density at the phase transition from nucleonic to quark matter at zero
temperature and nucleonic density $\rho_{\rm tr}$; $\Delta\varepsilon$ is the
energy discontinuity and $c_{s}$ is the sound speed in the quark phase. The
first-order phase transition described by this EoS leads to a new stable
branch populated by hybrid stars which may or may not be separated from the
nucleonic branch by an instability region. Furthermore, for a range of
parameters, this leads to the phenomenon of twin configurations — two stable
CSs that have identical masses but different radii and consequently, tidal
deformabilities (Alford et al., 2013; Benic et al., 2015; Paschalidis et al.,
2018; Alvarez-Castillo et al., 2019; Christian et al., 2019; Montana et al.,
2019; Christian & Schaffner-Bielich, 2021). In this work, we adopt two values
of $c_{s}^{2}$: the maximally stiff EoS with $c^{2}_{s}=1$ (in natural units)
allows us to explore the limits of ranges of masses and radii accessible to
models of the type considered in this work; see also Tsaloukidis et al.
(2023). The alternative value $c^{2}_{s}=2/3$ is intermediate between the
maximally stiff and the asymptotic high-density value $c^{2}_{s}=1/3$
corresponding to a conformally symmetric quark phase. Other values of
$c^{2}_{s}$ have been considered very recently by several authors (Brodie &
Haber, 2023; Sagun et al., 2023) in the context of HESS J1731-347 for
alternate nucleonic EoSs.

##  3\. Computational Scheme

To cover the parameter space of our model, we choose from the family of 81
nucleonic EoS models given in Li & Sedrakian (2023) a subset which are
identified by the values of $Q_{\rm sat}$ and $L_{\rm sym}$. We then choose
the values for $\rho_{\rm tr}=1.5,\,2.0,\,2.5,\,3.0\,\rho_{0}$ to match these
nucleonic EoS to quark matter EoS, for a magnitude of the density jump
$\Delta\varepsilon$ and sound speed $c_{s}^{2}=2/3$ or 1.

For each EoS in our collection, we solve the Tolman-Oppenheimer-Volkoff
equations for static and spherically symmetric stellar objects. If all the
mass and radius ranges of astronomical data are intersected by a sequence
corresponding to a particular EoS it is accepted, otherwise it is rejected.
More specifically, we consider the set of constraints (a)-(c) listed in the
Introduction, i.e., those based on NICER inferences of mass and radius for
J0030+0451 and J0740+6620, GW170817 inference of radii and masses in the
binary CS merger, and radio-wave measurements of the mass of J0740+6620. We
then supplement these by the mass and radius of J1731-347 and the mass of
J0952-0607. Note that, if we assume that uniform rotation, the increase of
mass of the maximum-mass configuration due to the rotation of 707 Hz is about
of $0.05\,M_{\odot}$. Therefore, the lower limit on the mass of J0952-0607 is
about the same as that of J0740+6620 at a 95% confidence level. The parameters
for quark EoS, $p_{\rm tr}$ and $\Delta\varepsilon$, can be transformed into
observational properties of CSs, i.e., the maximum mass of the nucleonic
branch $M^{\rm H}_{\rm max}$ and that of the hybrid branch $M^{\rm Q}_{\rm
max}$, respectively, in the $M$-$R$ diagram (Li et al., 2023b). We further
introduce the minimal mass in the hybrid branch, denoted as $M^{\rm Q}_{\rm
min}$, which we shall use for assessing the range of mass where twin
configurations exist.

Figure 4.—  Ranges of mass and radius that characterize twin configurations
for the case A, i.e., $M^{\rm Q}_{\rm max}\geq M^{\rm H}_{\rm max}$. In panels
(a, b) the hybrid EoS models ($c^{2}_{s}=1$) are built using an isovector-soft
nucleonic model ($Q_{\rm sat}=800,\,L_{\rm sym}=20$ MeV), while in panels (c,
d) - using an isovector-stiff nucleonic model ($Q_{\rm sat}=800,\,L_{\rm
sym}=80$ MeV). Panels (a, c) show the mass ranges of twins $\Delta M_{\rm
twin}$ using color coding. Panels (b, d) do the same but for the radii of twin
configurations $\Delta R_{\rm twin}$. By open squares, we show the models that
are ruled out by combined observational constraints, whereas by open circles
we show those that are excluded by HESS J1731-347 data.

##  4\. Analysis of the results

Figures 1 and 2 show the $M$-$R$ diagrams for hybrid EoS models where the
nucleonic EoS is specified by values of $Q_{\rm sat}$ and $L_{\rm sym}$ and
the quark matter EoS is specified by the sound speed value $c^{2}_{s}=1$ and
2/3. Figure 1 is for lower transition densities $\rho_{\rm tr}=1.5$ or
$2.0\,\rho_{0}$; Figure 2 is for higher transition densities
$\rho_{\text{tr}}=2.5$ or $3.0\,\rho_{0}$. The isovector properties of
nucleonic EoS are measured by the value of $L_{\rm sym}$. Consider first the
isovector-soft ($L_{\rm sym}=20$ MeV) EoSs, which contain members that are
consistent with astrophysical constraints at the 95% confidence level for any
of the chosen parameters shown in panels of Figures 1 and 2. It is seen that
as the transition density increases in the considered range
$1.5$-$3.0\,\rho_{0}$ the masses of twin stars increase as well. A larger
speed-of-sound allows the hybrid branch to stretch to the region of smaller
radii, and more generally to cover a broader range, as can be seen from the
ultra-compact members shown in the insets in panels (b, d) of Figure 2. These
models can be contrasted with the isovector-stiff ($L_{\rm sym}=80$ MeV) EoSs.
In this case, the transition to quark matter needs to occur at a sufficiently
low density $\rho_{\rm tr}\lesssim 2.0\,\rho_{0}$, to satisfy the constraints
imposed by HESS J1731-347 and GW170817, as displayed in Figure 1. Furthermore,
compatibility with the large radius of PSR J0740+6620 is achieved for
relatively small energy density discontinuity at the phase transition. The
resulting sequences satisfy the two-solar mass constraint because the maximum
of the hybrid branch $M^{\rm Q}_{\rm max}\sim 2.0\,M_{\odot}$, which is also
the maximum mass for the entire sequences. In this case, only low-mass twin
configurations appear that are consistent with astrophysical constraints.

Figure 3 shows the high-density portion of all the EoS used in this work,
i.e., including those that fail to account for HESS J1731-347. The EoS of
ultra-compact stars, with radii in the range $6\leq R\leq 10$ km and much
higher central densities, is separated in the insets. It is seen that the
addition of J1731-347 does not change the range of allowed EoS and that the
high-density range covered by our collection is broadly compatible with the
one deduced in Annala et al. (2022) by imposing current astrophysical
constraints on a very large randomly generated ensemble of physics-agnostic
EoS covering the intermediate density range. We also note that the low-density
range of the EoS, which can be found in Li & Sedrakian (2023), is consistent
with the low-density microscopic (chiral) EoS, as presented in Annala et al.
(2022), only for $40\leq L_{\rm sym}\leq 60$.

To facilitate the following discussion, we divided the stellar configurations
into those where the maximum-mass star is on the hybrid branch (case A) and
those where the maximum-mass star is on the nucleonic branch (case B). This
type of classification has been considered previously; case A was studied in
Li et al. (2021) and Christian & Schaffner-Bielich (2022) and case B in Li et
al. (2023a).

Case A: Figure 4 summarizes the contours representing the ranges of mass and
radius that characterize twin configurations for hybrid EoS models with
$M^{\rm H}_{\rm max}\leqslant M^{\rm Q}_{\rm max}$ for $c_{s}^{2}=1$. Filled
symbols indicate models that satisfy astrophysical constraints, while empty
symbols represent models that do not. The color of the square shows the range
of mass or radius over which twins can be found. The case $c_{s}^{2}=2/3$
shows the same trends but the ranges of masses and radii are significantly
narrower. The shape of the sequences in Case A dictates that the range of
masses where twin configurations exist is defined as $\Delta M_{\rm
twin}=M^{\rm H}_{\rm max}-M^{\rm Q}_{\rm min}$. Then, the differences in radii
$\Delta R_{\rm twin}$ for twin configurations is determined by comparing the
radius of $M^{\rm H}_{\rm max}$ nucleonic star to that of the hybrid
counterpart with an identical mass. It can be observed from Figures 1 and 2
that this is also the maximal radius difference for a hybrid EoS predicting
twin configurations. A lesson that we learn from Figure 4 is that the possible
combinations $M^{\rm Q}_{\rm max}$ and $M^{\rm H}_{\rm max}$ for which twin
configurations exist is large for the isovector-soft nucleonic model [(panel
(a)] and is small for the isovector-stiff model [(panel (c)]. Also these two
cases predict different $M^{\rm Q}_{\rm max}$ and $M^{\rm H}_{\rm max}$ for
which the largest range of masses $\leq 0.2M_{\odot}$ (shown by dark squares)
is achieved. The corresponding values of $\Delta R_{\rm twin}$ shown in panels
(b, d) indicate that the largest difference between the twin radii is achieved
for sequences that have both maxima around the two-solar mass. The twins with
$\Delta R_{\rm twin}\sim 4$-5 km correspond to ultracompact hybrid stars
discussed earlier by us (Li et al., 2023a).

Figure 5.—  Same as in Figure 4 but for the case B corresponding to $M^{\rm
Q}_{\rm max}\leq M^{\rm H}_{\rm max}$. The hybrid EoS models ($c^{2}_{s}=1$)
are built upon an isovector-soft nucleonic model ($Q_{\rm sat}=800,\,L_{\rm
sym}=20$ MeV). Ultra-compact configurations having large $\Delta R_{\rm twin}$
are present only for models with $M^{\rm H}_{\rm max}\simeq
2.0$-$2.15\,M_{\odot}$.

Case B: Figure 5 summarizes the ranges of masses and radii that characterize
twin configurations for hybrid EoS models in the case $M^{\rm Q}_{\rm
max}\leqslant M^{\rm H}_{\rm max}$, for $c_{s}^{2}=1$. Here the range of
masses is defined $\Delta M_{\rm twin}=M^{\rm Q}_{\rm max}-M^{\rm Q}_{\rm
min}$ and the differences in radii $\Delta R_{\rm twin}$ is extracted from the
radius of $M^{\rm Q}_{\rm max}$ hybrid star and its nucleonic twin. For such
cases, as illustrated in the insets of Figure 2, we require the nucleonic
branch passing through all the 95% confidence regions of constraints, this can
be achieved only for isovector-soft models. It is seen that ultra-compact
configurations with radii $6\leq R\leq 10$ km, associated with large values of
$\Delta R_{\rm twin}$, are present only for models with $M^{\rm H}_{\rm
max}\simeq 2.0$-$2.15\,M_{\odot}$. These ultra-compact members studied in more
detail in (Li et al., 2023a) together with those appearing in the case A offer
an alternative to strange stars, which up until recently were considered as
the only possible objects with small radii (of the order of 6–9 km) (see
Bombaci et al., 2021, and references therein). It is worthwhile also to
mention that the mass values found for twins in this case are well within the
range of mass obtained for PSR J0592-0607 at a 95% confidence level (Romani et
al., 2022).

##  5\. Conclusions

In this work, we investigated the consistency of the hybrid star models with
the current astrophysical and laboratory data and explored the range of
parameters that lead to twin stars. To do so, we employed a recently
constructed family of nucleonic EoSs parametrized by the $Q_{\rm sat}$ and
$L_{\rm sym}$ parameters of the Taylor expansion of the energy density of
nucleonic matter. This was then combined with a CSS parameterization of the
quark matter EoS for a range of values of the transition density, energy
density jump, and sound speed. We have identified two classes of models with
low and high transition densities (as shown in Figures 1 and 2, respectively)
that are consistent with the current data. This is achieved through a
compromise between the requirement for a soft EoS at low densities, driven by
the data from GW170817, and the need for stiff EoS at high densities to
account for the radio-astronomy data for PSR J0592-0607 and J0740+6620. The
ranges of parameters of twin configurations that are allowed by our hybrid
models were found as shown in Figures 4 and 5. Notably, the data from HESS
J1731-347 supports the existence of low-mass twin configurations with masses
$M\lesssim 1.3\,M_{\odot}$, regardless of the chosen nucleonic EoS, but it
requires a low value for the transition density. These hybrid models also
allow for maximum masses $M^{\rm Q}_{\rm max}$ ranging from
$2.0$-$2.6\,M_{\odot}$ depending on the value of the sound speed in quark
matter. In addition, ultra-compact configurations, which lie outside the
constraints imposed by GW170817 and J1731-347, emerge for isoscalar-stiff
nucleonic EoS models and for large speed-of-sound $c_{s}=1$ in the case where
the nucleonic branch extends up to $2.0$-$2.15\,M_{\odot}$.

In conclusion, the current astrophysical and nuclear physics data do not rule
out hybrid stars, and leave room for hybrid-nucleonic mass twins if a strong
first-order phase transition occurs in dense matter for reasonable values of
nucleonic isoscalar and isovector characteristics. It is even plausible to
argue that all the stars with observational information we have considered in
this study are hybrid configurations. But this does not exclude the
possibility that current observational data is accounted for by the nucleonic
branch of hybrid sequences, whereas the hybrid branch lies outside of ranges
currently covered, as would be the case for ultra-compact members of
sequences.

## Acknowledgements

J. L. acknowledges the support of the National Natural Science Foundation of
China (Grant No. 12105232), the Fundamental Research Funds for the Central
Universities (Grant No. SWU-020021), and by the Venture & Innovation Support
Program for Chongqing Overseas Returnees (Grant No. CX2021007). A. S.
acknowledges the DFG No. SE 1836/5-2 and the Polish NCN Grant No.
2020/37/B/ST9/01937 at Wrocław University. M. A. is partly supported by the
U.S. Department of Energy, Office of Science, Office of Nuclear Physics under
Award No. DE-FG02-05ER41375.

## References

  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbot, T. D., et al. 2019, PhRvX, 9, 011001, doi: 10.1103/PhysRevX.9.011001
  * Adhikari et al. (2021) Adhikari, D., Albataineh, H., Androic, D., et al. 2021, PhRvL, 126, 172502, doi: 10.1103/PhysRevLett.126.172502
  * Adhikari et al. (2022) —. 2022, PhRvL, 129, 042501, doi: 10.1103/PhysRevLett.129.042501
  * Alford et al. (2005) Alford, M., Braby, M., Paris, M. W., & Reddy, S. 2005, ApJ, 629, 969, doi: 10.1086/430902
  * Alford & Sedrakian (2017) Alford, M., & Sedrakian, A. 2017, PhRvL, 119, 161104, doi: 10.1103/PhysRevLett.119.161104
  * Alford et al. (2013) Alford, M. G., Han, S., & Prakash, M. 2013, PhRvD, 88, 083013, doi: 10.1103/PhysRevD.88.083013
  * Alford et al. (2008) Alford, M. G., Schmitt, A., Rajagopal, K., & Schäfer, T. 2008, RvMP, 80, 1455, doi: 10.1103/RevModPhys.80.1455
  * Alvarez-Castillo et al. (2019) Alvarez-Castillo, D. E., Blaschke, D. B., Grunfeld, A. G., & Pagura, V. P. 2019, PhRvD, 99, 063010, doi: 10.1103/PhysRevD.99.063010
  * Anglani et al. (2014) Anglani, R., Casalbuoni, R., Ciminale, M., et al. 2014, RvMP, 86, 509, doi: 10.1103/RevModPhys.86.509
  * Annala et al. (2022) Annala, E., Gorda, T., Katerini, E., et al. 2022, PhRvX, 12, 011058, doi: 10.1103/PhysRevX.12.011058
  * Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., et al. 2018, RPPh, 81, 056902, doi: 10.1088/1361-6633/aaae14
  * Benic et al. (2015) Benic, S., Blaschke, D., Alvarez-Castillo, D. E., Fischer, T., & Typel, S. 2015, A&A, 577, A40, doi: 10.1051/0004-6361/201425318
  * Blaschke et al. (2020) Blaschke, D., Ayriyan, A., Alvarez-Castillo, D. E., & Grigorian, H. 2020, Univ, 6, 81, doi: 10.3390/universe6060081
  * Bombaci et al. (2021) Bombaci, I., Drago, A., Logoteta, D., Pagliara, G., & Vidaña, I. 2021, PhRvL, 126, 162702, doi: 10.1103/PhysRevLett.126.162702
  * Bonanno & Sedrakian (2012) Bonanno, L., & Sedrakian, A. 2012, A&A, 539, A16, doi: 10.1051/0004-6361/201117832
  * Brodie & Haber (2023) Brodie, L., & Haber, A. 2023, PhRvC, 108, 025806, doi: 10.1103/PhysRevC.108.025806
  * Christian & Schaffner-Bielich (2021) Christian, J.-E., & Schaffner-Bielich, J. 2021, PhRvD, 103, 063042, doi: 10.1103/PhysRevD.103.063042
  * Christian & Schaffner-Bielich (2022) Christian, J.-E., & Schaffner-Bielich, J. 2022, ApJ, 935, 122, doi: 10.3847/1538-4357/ac75cf
  * Christian et al. (2019) Christian, J.-E., Zacchi, A., & Schaffner-Bielich, J. 2019, PhRvD, 99, 023009, doi: 10.1103/PhysRevD.99.023009
  * Doroshenko et al. (2022) Doroshenko, V., Suleimanov, V., Pühlhofer, G., & Santangelo, A. 2022, NatAs, 6, 1444
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, ApJL, 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Glendenning & Kettner (2000) Glendenning, N. K., & Kettner, C. 2000, A&A, 353, L9.  https://arxiv.org/abs/astro-ph/9807155
  * Klähn et al. (2013) Klähn, T., Łastowiecki, R., & Blaschke, D. B. 2013, PhRvD, 88, 085001, doi: 10.1103/PhysRevD.88.085001
  * Lalazissis et al. (2005) Lalazissis, G. A., Niksic, T., Vretenar, D., & Ring, P. 2005, PhRvC, 71, 024312, doi: 10.1103/PhysRevC.71.024312
  * Lattimer (2023) Lattimer, J. M. 2023, Particles, 6, 30, doi: 10.3390/particles6010003
  * Li & Sedrakian (2019) Li, J. J., & Sedrakian, A. 2019, ApJL, 874, L22, doi: 10.3847/2041-8213/ab1090
  * Li & Sedrakian (2023) —. 2023, ApJ, 957, 41, doi: 10.3847/1538-4357/acfa73
  * Li et al. (2020) Li, J. J., Sedrakian, A., & Alford, M. 2020, PhRvD, 101, 063022, doi: 10.1103/PhysRevD.101.063022
  * Li et al. (2021) —. 2021, PhRvD, 104, L121302, doi: 10.1103/PhysRevD.104.L121302
  * Li et al. (2023a) —. 2023a, PhRvD, 107, 023018, doi: 10.1103/PhysRevD.107.023018
  * Li et al. (2023b) —. 2023b, ApJ, 944, 206, doi: 10.3847/1538-4357/acb688
  * Malfatti et al. (2020) Malfatti, G., Orsaria, M. G., Ranea-Sandoval, I. F., Contrera, G. A., & Weber, F. 2020, PhRvD, 102, 063008, doi: 10.1103/PhysRevD.102.063008
  * Margueron et al. (2018) Margueron, J., Hoffmann Casali, R., & Gulminelli, F. 2018, PhRvC, 97, 025806, doi: 10.1103/PhysRevC.97.025806
  * Masuda et al. (2013) Masuda, K., Hatsuda, T., & Takatsuka, T. 2013, ApJ, 764, 12, doi: 10.1088/0004-637X/764/1/12
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJL, 887, L24, doi: 10.3847/2041-8213/ab50c5
  * Miller et al. (2021) —. 2021, ApJL, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Montana et al. (2019) Montana, G., Tolos, L., Hanauske, M., & Rezzolla, L. 2019, PhRvD, 99, 103009, doi: 10.1103/PhysRevD.99.103009
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, RvMP, 89, 015007, doi: 10.1103/RevModPhys.89.015007
  * Paschalidis et al. (2018) Paschalidis, V., Yagi, K., Alvarez-Castillo, D., Blaschke, D. B., & Sedrakian, A. 2018, PhRvD, 97, 084038, doi: 10.1103/PhysRevD.97.084038
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21, doi: 10.3847/2041-8213/ab481c
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJL, 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Romani et al. (2022) Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2022, ApJL, 934, L17, doi: 10.3847/2041-8213/ac8007
  * Sagun et al. (2023) Sagun, V., Giangrandi, E., Dietrich, T., et al. 2023, ApJ, 958, 49, doi: 10.3847/1538-4357/acfc9e
  * Sedrakian (2023) Sedrakian, A. 2023, Particles, 6, 713, doi: 10.3390/particles6030044
  * Sedrakian et al. (2023) Sedrakian, A., Li, J. J., & Weber, F. 2023, PrPNP., 131, 104041, doi: 10.1016/j.ppnp.2023.104041
  * Tan et al. (2022) Tan, H., Dore, T., Dexheimer, V., Noronha-Hostler, J., & Yunes, N. 2022, PhRvD, 105, 023018, doi: 10.1103/PhysRevD.105.023018
  * Tan et al. (2020) Tan, H., Noronha-Hostler, J., & Yunes, N. 2020, PhRvL, 125, 261104, doi: 10.1103/PhysRevLett.125.261104
  * Tsaloukidis et al. (2023) Tsaloukidis, L., Koliogiannis, P. S., Kanakis-Pegios, A., & Moustakidis, C. C. 2023, PhRvD, 107, 023012, doi: 10.1103/PhysRevD.107.023012
  * Typel & Wolter (1999) Typel, S., & Wolter, H. H. 1999, NuPhA, 656, 331, doi: 10.1016/S0375-9474(99)00310-3
  * Zdunik & Haensel (2013) Zdunik, J. L., & Haensel, P. 2013, A&A, 551, A61, doi: 10.1051/0004-6361/201220697