# Sexaquark dilemma in neutron stars and its solution by quark deconfinement

M. Shahrbaf<EMAIL_ADDRESS>Institute of Theoretical Physics,
University of Wroclaw, 50-204 Wroclaw, Poland  D. Blaschke
<EMAIL_ADDRESS>Institute of Theoretical Physics, University of
Wroclaw, 50-204 Wroclaw, Poland  Bogoliubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, 141980 Dubna, Russia  National Research
Nuclear University (MEPhI), 115409 Moscow, Russia  S. Typel<EMAIL_ADDRESS>darmstadt.de Technische Universität Darmstadt, Fachbereich Physik, Institut
für Kernphysik, Schlossgartenstraße 9, D-64289 Darmstadt, Germany  GSI
Helmholtzzentrum für Schwerionenforschung GmbH, Theorie, Planckstraße 1,
D-64291 Darmstadt, Germany  G. R. Farrar<EMAIL_ADDRESS>Center for Cosmology
and Particle Physics, Department of Physics, New York University, NY, NY
10003, USA  D. E. Alvarez-Castillo<EMAIL_ADDRESS>Bogoliubov
Laboratory of Theoretical Physics, Joint Institute for Nuclear Research,
141980 Dubna, Russia  Henryk Niewodniczański Institute of Nuclear Physics,
Polish Academy of Sciences, 31-342 Cracow, Poland

###### Abstract

Following the idea that a stable sexaquark state with quark content (uuddss)
would have gone unnoticed by experiment so far and that such a particle would
be a good dark matter candidate, we investigate the possible role of a stable
sexaquark in the physics of compact stars given the stringent constraints on
the equation of state that stem from observations of high mass pulsars and
GW170817 bounds on the compactness of intermediate mass stars. We find that
there is a “sexaquark dilemma” (analogous to the hyperon dilemma) for which
the dissociation of the sexaquark in quark matter is a viable solution
fulfilling all present constraints from multi-messenger astronomy. The
parameters needed to model the hybrid star including sexaquarks are in line
with parameters of pre-existing quark- and hadronic-matter models. We find
that current constraints – tidal deformability in accordance with GW170817 and
maximum mass above the lower limit from PSR J0740+6620 – can be satisfied two
ways: with early quark deconfinement such that neither sexaquarks nor hyperons
are present in any NS interiors, or with later deconfinement such that a
neutron-sexaquark shell surrounds the inner quark matter core.

###### pacs:

05.30.-d, 12.39.-x, 25.75.Nq, 21.60.Gx, 24.85.+p

##  I Introduction

The present work is devoted to consideration of the consequences for neutron
star (NS) phenomenology that would follow from the existence of a possible
stable sexaquark (S) state with the quark content $uuddss$.

The S is an electrically neutral spin-less boson with baryon number $B_{S}=2$
and strangeness $S_{S}=-2$ in a flavor-singlet state. If it is light enough to
be stable against weak decay ($m_{S}<m_{p}+m_{\Lambda}+m_{e}=2054$ MeV), the S
is a good dark matter candidate Farrar (2003, 2018) and experiments to date
would not have been sensitive to it Farrar (2022). The sexaquark has to be
distinguished from the H-dibaryon (H) with the same quark content, which was
introduced by Jaffe Jaffe (1977) and estimated using the MIT Bag Model to have
a mass of 2150 MeV – greater than $m_{\Lambda}+m_{p}+m_{e}$ and hence weak-
interaction unstable with a lifetime $\mathcal{O}(10^{-10}s$). As illustrated
in Fig. 1, both a molecule of two $\Lambda$ hyperons as well as a more tightly
bound state of three diquarks could exist.

Our specific goal is to investigate whether the existence of neutron stars
with masses above $2M_{\odot}$ and the evidence that NS radii vary slowly with
mass above $\sim 1.4M_{\odot}$, excludes the existence of a scalar strongly-
interacting dibaryon with mass in the range where it is a potential dark
matter candidate, $\approx 2m_{N}$ to 2054 MeV. To this end, we use a pre-
existing formalism for the hadronic equation of state and extend it to include
an S along with a 2-flavor constant-speed-of-sound equation of state (EoS) for
deconfined quark matter, and we consider two methods for interpolating between
them. With a non-exhaustive exploration of model parameters (without adjusting
beyond the normal range), we find that a stable S with mass above $\approx
1885$ MeV is compatible with current neutron star observations. We find
solutions of two types: i) with early deconfinement and neither S nor hyperons
present in neutron stars of any mass and ii) with a layer containing both S
and nucleons outside a quark matter core. When our understanding of the
transition between hadronic and quark matter is improved, knowledge of neutron
stars may provide constraints on the allowed mass of a sexaquark, or even
point to its existence. Indeed, we find that an S naturally explains the soft
EoS at low densities implied by observations of GW170817.

Figure 1: Comparison of two realisations of a color singlet hadron with the
quark content (uuddss). Left: hadronic molecule of two $\Lambda$(uds)
hyperons, corresponding to the H-dibaryon. Right: compact bound state of three
diquarks, bound by color forces, which corresponds to a possible structure of
the sexaquark.

##  II Overview

One of the puzzling questions in the physics of hadrons concerns the structure
of recently discovered multiquark clusters such as pentaquarks and
tetraquarks. Are these states molecules of color-neutral (hadronic)
subclusters analogous to the hydrogen molecule of atomic physics? Or are they
multiquark states bound by confining forces between color charged quarks and
diquarks? Or some of each? As an example, the recently discovered all-charm
tetraquark X(6900) Aaij  _et al._ (2020) can a priori be viewed as a molecule
of two J/$\psi$ mesons or a diquark-antidiquark bound state. Similarly, a
pentaquark could be a baryon-meson molecule or a bound state of two diquarks
and an antiquark. For more details, see the textbook Ali  _et al._ (2019).

While the H state has been investigated in lattice QCD calculations Green  _et
al._ (2021); Francis  _et al._ (2019); Beane  _et al._ (2013); Inoue  _et al._
(2011); Beane  _et al._ (2011) and a binding energy of a few to 80 MeV has
been found at still unphysically large quark mass, the much lighter and
possibly more compact S state is not yet accessible to lattice QCD
calculations Detmold (2021).

In a recent calculation within a constituent quark model, Buccella obtained
the mass of the S as a bound state of three diquarks to be $m_{S}=1883$ MeV
Buccella (2020) while Azizi et al. have obtained it to be about $m_{S}=1180$
MeV using QCD sum rules Azizi  _et al._ (2020). Buccella showed also that this
estimate does not depend on the choice of the parameters for the constituent
quark masses by replacing the latter with an expression using the heavier
nonet mesons, $m_{S}=2M_{K(1430)}-(1/2)M_{f^{0}(1370)}+40~{}{\rm
MeV}+(1/2)[M_{N}-M_{\Delta}]+M_{\Sigma}-M_{Y^{*}}=1876.5$ MeV, where the first
three terms correspond to the constituent mass and the other terms to the
contribution of the chromomagnetic interaction Buccella (2020). An S in this
mass range is sufficiently heavy to not induce instability in the deuteron (D)
and other nuclei Farrar and Wang (2022); Farrar (2022)

For $m_{S}\leq m_{D}=1876.122$ MeV, the S is absolutely stable due to baryon
number conservation, while as long as $m_{S}\leq
m_{p}+m_{e}+m_{\Lambda}=2054.466$ MeV its decay requires $\Delta S=2$ and is
hence doubly-weak, and its lifetime exceeds the age of the universe Farrar and
Zaharijas (2004); Farrar and Wang (2022). In this study we suppose that the S
may be a deeply bound state with low enough mass to be absolutely or
effectively stable; in that case it is an attractive dark matter (DM)
candidate Farrar (2018, 2022). The observed dark matter to baryon ratio is
$\Omega_{DM}/\Omega_{B}=5.3\pm 0.1$ Ade  _et al._ (2016); Tanabashi  _et al._
(2018) and a successful model for dark matter should account for this value.
In fact, an abundance of S dark matter (SDM) in agreement with this
observation follows Farrar (2018) from statistical mechanics using quark
masses as determined from lattice QCD, given the effective freezeout
temperature $T_{\rm eff}\approx 150$ MeV for $m_{S}$ in the expected
range.111The predicted $\Omega_{DM}/\Omega_{B}$ is very insensitive to the
assumed freezeout temperature; the relevant range is motivated by the recent
result of $T_{c}=156.5\pm 1.5$ MeV for the pseudocritical temperature obtained
in lattice QCD simulations Bazavov  _et al._ (2019). Note a plotting error in
Farrar (2018) is corrected in v2 of Farrar  _et al._ (2020); the formulae
given in Farrar (2018) are correct.

Another environment for the production of the S, apart from the hadronization
transition in the early Universe, is the “Little Bang” in ultrarelativistic
heavy-ion collisions at the LHC, where not only traditional hadrons but also
light nuclei and antinuclei are abundantly produced. Recently, the ratio of
light sexaquarks relative to the deuteron to be expected under these
conditions has been estimated within a thermal statistical model Blaschke  _et
al._ (2021) and found to be sizeable, of the order one. Abundant production of
sexaquarks is not sufficient for their discovery – there remains the problem
of the detection of the S and discriminating it from the far more abundant
neutrons; see Farrar (2022) for an analysis of detection strategies.

In the present work, we consider S in the mass range of
$1885~{}\mbox{MeV}<m_{S}<2054$ MeV.222The lower end of the range is sufficient
that condensation of S below saturation density can be avoided and above the
upper end of the range the S decays in $\mathcal{O}(10^{-10}s)$ and is not a
Dark Matter candidate. We study the possible relevance of the S for the
properties of neutron stars in light of recent multi-messenger observations
that constrain the mass-radius diagram and tidal deformability. When the S is
present in cold dense baryonic matter in neutron stars, it forms a Bose
condensate as soon as the baryo-chemical potential in the center of the star
fulfills $\mu_{b}=m_{S}/2$, unless a transition to non-hadronic degrees of
freedom occurs first. For a density-independent mass, for instance
$m_{S}=1941\,(2054)$ MeV, this occurs for a star with
$M=0.21~{}M_{\odot}\,(0.7~{}M_{\odot})$ respectively. The mass $m_{S}=1941$
MeV is chosen as an illustrative case, because for a density-independent mass
and $m_{S}=1941$ MeV, the Bose condensation occurs exactly at saturation
density. Saturation density is roughly the lowest density at which
condensation can be tolerated due to the existence of nuclei. Therefore, due
to the saturation of the pressure once Bose condensation occurs,
$0.21~{}M_{\odot}$ is roughly the maximum NS mass that can be reached for a
non-interacting S with mass less than 1941 MeV, unless a transition to new
degrees of freedom replaces hadrons. Such a low maximum mass is in clear
contradiction with the observation of pulsars as massive as $2~{}M_{\odot}$
like PSR J0740+6620 Fonseca  _et al._ (2021) or PSR J0348+0432 Antoniadis  _et
al._ (2013) and many, many lower-mass neutron stars.

The problem with S Bose condensation is alleviated by the plausible assumption
that the effective mass of the S is medium-dependent and increases with
density, similar to the behavior of the other baryon masses at supersaturation
densities. However, as we will show in this work, simply a medium-dependent
mass for the S is insufficient to allow for neutron star masses consistent
with the observational constraint $M_{\rm max}\gtrsim 2~{}M_{\odot}$. This
constitutes the “sexaquark dilemma” to which this work is devoted. The same is
true for hyperons – the so-called hyperon puzzle. It is known that the hyperon
puzzle can be solved by quark deconfinement; see, e.g., Shahrbaf  _et al._
(2020a, b) and references therein. Therefore, we explore in the present paper
various scenarios for the sexaquark dilemma and the viability of solution by
quark deconfinement.

In a recent work, McDermott et al. McDermott  _et al._ (2019) argued that a
deeply bound S – in the mass range considered here – would be incompatible
with the delayed neutrino signal from supernova 1987A. However the inter-
conversion amplitude between S and two $\Lambda$’s is naturally small, and the
value needed to not impact SN1987a cooling is comfortably in the expected
range Farrar and Wang (2022); Farrar (2022), so SN1987A cooling is not the
show-stopper for sexaquark dark matter envisaged in McDermott  _et al._
(2019).333Similarly the conclusions of of Kolb and Turner (2019) do not apply
to sexaquark dark matter because the suppressed interconversion between
sexaquark and two baryons inhibits the destruction of sexaquarks in the high-
temperature hadronic phase. While the conditions relevant to the analysis of
McDermott  _et al._ (2019) are unlike the steady state conditions we consider,
we note that their assumed hyperon-rich environment and medium independence of
the masses, and neglect of the possibility of quark deconfinement would not be
applicable for our solution of the S dilemma.

The structure of the present paper is as follows. In Sec. III we present the
formalism of calculation the equation of state (EoS) of hadronic matter based
on the density dependent (DD) relativistic mean field (RMF) model in neutron
star matter. In Sec. IV the EoS of quark matter is obtained and the theory of
two different approaches for constructing the hybrid stars are given. In Sec.
V we present our results for the properties of purely hadronic stars including
S particles and discuss their dilemmas. The deconfinement solution for this
dilemma is discussed in Sec. VI. Finally, the summary and conclusion are given
in Sec. VII.

##  III A relativistic density functional approach to hypernuclear matter with
sexaquark

The study of nuclear matter based on relativistic approaches with a Lagrangian
density including baryons and mesons as degrees of freedom has a long history
Walecka (1974); Serot and Walecka (1986); Reinhard (1989); Ring (1996). The
original treatment as full-fledged field theories was subsequently replaced by
the interpretation as effective field theories Rusnak and Furnstahl (1997) and
more recently as density functional theories Dreizler and Gross (2012). In
this work the EoS of hadronic matter is obtained from a generalized
relativistic density functional (GRDF) with baryon-meson couplings that depend
on the total baryon density of the system. The original density functional for
nucleonic matter considers the isoscalar $\sigma$ and $\omega$ mesons and the
isovector $\rho$ meson as exchange particles that describe the effective in-
medium interaction. The density dependence of the couplings is adjusted to
describe properties of atomic nuclei Typel and Wolter (1999); Typel (2005). It
has been confirmed that such GRDFs are successful in reproducing the
properties of nuclear matter around nuclear saturation Klahn  _et al._ (2006).

The explicit introduction of mesons in relativistic density functionals is not
necessary for the description of nuclear matter since the theoretical
formulation can be based on the baryonic degrees of freedom and their
densities only. This is realized, e.g., in relativistic point-coupling models,
see Sun  _et al._ (2019) and references therein. However, the use of mesons is
convenient from a practical point of view as they represent a simple means to
describe the relevant components of the effective interaction.

GRDFs with different functional forms of the density dependent couplings have
been studied, e.g., in Typel (2018), and various parameterizations are
available in the literature Dutra  _et al._ (2014). In the present work the
parameterization DD2 Typel  _et al._ (2010) is used for the $\sigma$,
$\omega$, and $\rho$ couplings. It was obtained by fitting properties of
finite nuclei: binding energies, charge and diffraction radii, surface
thicknesses and spin-orbit splittings. It predicts characteristic nuclear
matter parameters that are consistent with recent constraints Oertel  _et al._
(2017). In particular, a saturation density of $0.149065$ fm-3, a binding
energy per nucleon of $16.02$ MeV, and an incompressibility of $K=242.7$ MeV
in symmetric nuclear matter are found. The isospin dependence can be
characterized by a symmetry energy of $J=31.67$ MeV at saturation with a slope
parameter of $L=55.04$ MeV. Because no constraints were imposed at supra-
saturation densities in determining the DD2 parameters, the EoS constitutes a
pure extrapolation in this density range. The DD2 model leads to a rather
stiff EoS at high baryon densities with a maximum neutron star mass of $2.4$
solar masses in a pure nucleonic scenario of the strongly interacting system.
This approach with nucleons has been extended to include light and heavy
clusters as quasiparticles modified by medium effects Typel  _et al._ (2010);
Pais and Typel (2017) which are important in the finite-temperature EoS for
astrophysical applications, e.g., the simulation of core-collapse supernovae
or neutron-star mergers Oertel  _et al._ (2017). The internal structure of
these clusters is not taken in to account and they are treated as point-like
particles.

In the description of neutron-star matter one has to consider that new
baryonic degrees of freedom can become active with increasing density as the
chemical potentials rise. Neglecting any interaction, a new species appears
when the corresponding chemical potential crosses the particle mass. This
usually leads to a softening of the EoS as it is known for hyperons - the so-
called ’hyperon puzzle’ which has been discussed in Shahrbaf and Moshfegh
(2019); Shahrbaf  _et al._ (2020a, b) and the literature cited therein.
However, if the interaction with an additional $\phi$ meson, which couples
only to strangeness-carrying hyperons, is included in the GRDF, a sufficient
stiffening of the EoS can be achieved. In our present study, all particle in
the octet of spin $1/2$ baryons are included in the model as degrees of
freedom. Their interaction is described by the exchange of $\sigma$, $\omega$,
$\rho$, and $\phi$ mesons with appropriately adjusted couplings to the
baryons. The S particle is a further baryonic degree of freedom in the GRDF
with a density dependent mass shift that models the effect of the interaction
with the medium. For the description of neutron star matter electrons and
muons are added in the GRDF to achieve charge neutrality. The electron and
muon densities, assuming identical lepton chemical potentials, for given
baryon density are determined by the condition of $\beta$ equilibrium.

In the present work with the application to neutron stars we can restrict
ourselves to the case of matter at zero temperature. At densities below
saturation there are no hyperons or sexaquarks and the unified crust EoS of
the original GRDF-DD2 model with clusters is used. It contains the well-known
sequence of nuclei in a body-centered cubic lattice with a uniform background
of electrons and a neutron gas above the neutron drip line. The transition to
homogeneous matter just below the nuclear saturation density is described
consistently within the same approach. The main modification is to include the
new degrees of freedom at supersaturation densities in the GRDF model. This
will be described in the following subsections. All equations follow the
traditional convention of $\hbar=c=1$ of nuclear physics.

###  III.1 NS matter EoS with hyperons and sexaquark

The hadronic part of the core in a neutron star is assumed to consist of
homogeneous matter composed of strongly interacting baryons and charged, non-
interacting leptons in full thermodynamic equilibrium. Explicitly, these
degrees of freedom are protons ($p$), neutrons ($n$), hyperons ($\Lambda$,
$\Sigma^{+}$, $\Sigma^{0}$, $\Sigma^{-}$, $\Xi^{0}$, $\Xi^{-}$), the sexaquark
(S), electrons ($e$) and muons ($\mu$). All information on the thermodynamic
properties of the system can be obtained from a grand canonical thermodynamic
potential density $\Omega(\\{\mu_{i}\\})$ that depends only on the chemical
potentials $\mu_{i}$ of the individual particles at zero temperature. These
can be expressed as

$\mu_{i}=B_{i}\mu_{b}+Q_{i}\mu_{q}+S_{i}\mu_{s}+L_{i}\mu_{l}$ (1)

with the individual baryon ($B_{i}$), charge ($Q_{i}$), strangeness ($S_{i}$)
and lepton ($L_{i}$) numbers. Independent quantities are the baryon
($\mu_{b}$), charge ($\mu_{q}$), strangeness ($\mu_{s}$) and lepton
($\mu_{l}$) chemical potentials (assuming identical electron lepton and muon
lepton chemical potentials). The condition of full equilibrium with respect to
strangeness-changing reactions corresponds to $\mu_{s}=0$ and in case of
$\beta$ equilibrium we have $\mu_{l}=0$. Thus only two independent chemical
potentials, $\mu_{b}$ and $\mu_{q}$, remain. These have to be determined for a
given baryon density

$\displaystyle n_{b}$ $\displaystyle=$ $\displaystyle\sum_{i}B_{i}n_{i}^{(v)}$
(2) $\displaystyle=$ $\displaystyle
n_{p}^{(v)}+n_{n}^{(v)}+n_{\Lambda}^{(v)}+n_{\Sigma^{+}}^{(v)}+n_{\Sigma^{0}}^{(v)}+n_{\Sigma^{-}}^{(v)}$
$\displaystyle+n_{\Xi^{0}}^{(v)}+n_{\Xi^{-}}^{(v)}+2n_{S}^{(v)}$

and total charge density

$\displaystyle n_{Q}$ $\displaystyle=$ $\displaystyle\sum_{i}Q_{i}n_{i}^{(v)}$
$\displaystyle=$ $\displaystyle
n_{p}^{(v)}+n_{\Sigma^{+}}^{(v)}-n_{\Sigma^{-}}^{(v)}-n_{\Xi^{-}}^{(v)}-n_{e}^{(v)}-n_{\mu}^{(v)}$

with the particle number (or vector) densities $n_{i}^{(v)}$, see below.
Requiring local charge neutrality corresponds to $n{q}=0$ and thus the charge
chemical potential $\mu_{q}$ will be fixed too. Finally, the EoS of neutron
star matter depends only on the baryon chemical potential $\mu_{b}$.

All constituent particles with vacuum rest masses $m_{i}$ are considered as
quasiparticles in the medium with effective masses $m_{i}^{\ast}=m_{i}-S_{i}$
and effective chemical potentials $\mu_{i}^{\ast}=\mu_{i}-V_{i}$. These
contain scalar and vector potentials, $S_{i}$ and $V_{i}$, that describe the
interaction of the particles in the medium. They arise from the coupling of
the meson fields $\sigma$, $\omega$, $\rho$, and $\phi$ with coupling
strengths $\Gamma_{i\sigma}$, $\Gamma_{i\omega}$, $\Gamma_{i\rho}$, and
$\Gamma_{i\phi}$ to nucleons and hyperons or they are modeled effectively by a
mass shift $\Delta m_{S}$ for the sexaquark. The couplings are assumed to
depend on the total density of nucleons and hyperons

$n_{cpl}=n_{b}-2n_{S}$ (4)

and the mass shift depends on the baryon density (2). Using the meson field
names also for the field strengths, the potentials for nucleons and hyperons
have the form

$S_{i}=\Gamma_{i\sigma}\sigma$ (5)

and

$V_{i}=\Gamma_{i\sigma}\sigma+\Gamma_{i\omega}\omega+\Gamma_{i\rho}\rho+\Gamma_{i\phi}\phi+B_{i}V^{(r)}+W_{i}^{(r)}$
(6)

whereas they are given by

$S_{S}=-\Delta m_{S}$ (7)

and

$V_{S}=W_{S}^{(r)}$ (8)

for the sexaquark. The rearrangement contributions $V^{(r)}$, $W_{i}^{(r)}$,
see below, are required for thermodynamic consistency. The coupling of a
particle $i$ to a meson $m$ can be written as

$\Gamma_{im}=g_{im}\Gamma_{m}(n_{cpl})$ (9)

with prefactors $g_{im}$ and density dependent coupling functions
$\Gamma_{m}(n_{cpl})$. The dependence of these functions on the baryon density
$n_{b}$ is given by the functional forms as introduced in Typel and Wolter
(1999) and used for the DD2 parameterization, i.e.,

$\Gamma_{m}(n_{cpl})=\Gamma_{m}(n_{\rm sat}^{(v)})f_{m}(x)$ (10)

with $x=n_{cpl}/n_{\rm sat}^{(v)}$, the couplings at saturation
$\Gamma_{m}(n_{\rm sat}^{(v)})$ and functions

$f_{m}(x)=a_{m}\frac{1+b_{m}(x+d_{m})^{2}}{1+c_{m}(x+d_{m})^{2}}$ (11)

for $m=\sigma,\omega$ and

$f_{\rho}(x)=\exp[-a_{\rho}(x-1)]$ (12)

for the $\rho$ meson. See Typel  _et al._ (2010) for the actual parameters and
masses of the mesons. The coupling $\Gamma_{\phi}$ of the $\phi$ meson is
assumed to be constant in the present model with
$\Gamma_{\phi}=\Gamma_{\omega}(n_{\rm sat}^{(v)})$. In the following, also
derivatives $\Gamma_{m}^{\prime}=d\Gamma_{m}/dn_{cpl}$ of the couplings will
appear.

The coupling factors $g_{im}$ in eq. (9) of nucleons and hyperons to the
vector mesons can be expressed with help of the quantum numbers $B_{i}$,
$Q_{i}$, and $S_{i}$ of a particle $i$ as

$\displaystyle g_{i\omega}$ $\displaystyle=$ $\displaystyle
B_{i}+\frac{S_{i}}{3}$ (13) $\displaystyle g_{i\rho}$ $\displaystyle=$
$\displaystyle 2Q_{i}-B_{i}-S_{i}$ (14) $\displaystyle g_{i\phi}$
$\displaystyle=$ $\displaystyle\frac{\sqrt{2}}{3}S_{i}$ (15)

following the usual SU(6) coupling scheme, see, e.g., Weissenborn  _et al._
(2012). The coupling factors for the $\sigma$ meson to the nucleons are given
by $g_{p\sigma}=g_{n\sigma}=1$. For hyperons, they are determined by fixing
their in-medium potential $U_{Y}$ with $Y=\Lambda,\Sigma,\Xi$ in symmetric
nuclear matter at saturation density $n_{\rm sat}^{(v)}$ of the DD2
parameterization. Explicitly they are given by

$\displaystyle g_{Y\sigma}$ $\displaystyle=$
$\displaystyle\left[\left(g_{Y\omega}\Gamma_{\omega}(n_{\rm
sat}^{(v)})+\Gamma_{\omega}^{\prime}(n_{\rm sat}^{(v)})n_{\rm
sat}^{(v)}\right)\omega_{\rm sat}\right.$ (16)
$\displaystyle\left.-\Gamma_{\sigma}^{\prime}(n_{\rm sat}^{(v)})n_{\rm
sat}^{(s)}-U_{Y}\right]\left(\Gamma_{\sigma}(n_{\rm
sat}^{(v)})\sigma\right)^{-1}$

with the vector density $n_{\rm sat}^{(v)}$ and the scalar density $n_{\rm
sat}^{(s)}$ at saturation of the DD2 parameterization.

In principle, the interaction of the S particle with the medium can be
described by a coupling to the mesons as for the other baryonic degrees of
freedom. However, this requires to introduce several unknown meson-sexaquark
couplings. Since the sexaquark will exist as a condensate at zero temperature,
only the difference $V_{S}-S_{S}$ of the vector and scalar potentials appears
in the equations. This effective potential rises, to lowest order,
proportional to the baryon density. Thus it is sufficient to describe the
interaction with the medium simply by a mass shift

$\Delta m_{S}=m_{S}x_{S}\frac{n_{b}}{n_{0}},$ (17)

giving the effective mass

$m^{*}_{S}=m_{S}-S_{S}=m_{S}+\Delta m_{S}.$ (18)

For simplicity in a first exploration, we describe this with a single,
adjustable parameter $x_{S}$ and fixed constants $m_{S}$ and $n_{0}$. The mass
shift (17) rises linearly with the baryon density and represents an effective
repulsive interaction with the medium for positive $x_{S}$. Due to this
density dependence, there appears a rearrangement contribution in the vector
potential (8) of the S particle. It depends linearly on the scalar density of
the sexaquark and is thus considerably smaller than the mass shift as long as
this density is small compared to the baryon density.

Obviously, electrons and muons are treated as free particles since they don’t
participate in the strong interaction and there is no electric potential in
uniform matter.

The meson fields themselves are found from the meson field equations

$\displaystyle m_{\sigma}^{2}\sigma$ $\displaystyle=$
$\displaystyle\Gamma_{\sigma}n_{\sigma}$ (19) $\displaystyle
m_{\omega}^{2}\omega$ $\displaystyle=$
$\displaystyle\Gamma_{\omega}n_{\omega}$ (20) $\displaystyle m_{\rho}^{2}\rho$
$\displaystyle=$ $\displaystyle\Gamma_{\rho}n_{\rho}$ (21) $\displaystyle
m_{\phi}^{2}\phi$ $\displaystyle=$ $\displaystyle\Gamma_{\phi}n_{\phi}$ (22)

with the source densities

$\displaystyle n_{\sigma}$ $\displaystyle=$
$\displaystyle\sum_{i}g_{i\sigma}n_{i}^{(s)}$ (23) $\displaystyle n_{\omega}$
$\displaystyle=$ $\displaystyle\sum_{i}g_{i\omega}n_{i}^{(v)}$ (24)
$\displaystyle n_{\rho}$ $\displaystyle=$
$\displaystyle\sum_{i}g_{i\rho}n_{i}^{(v)}$ (25) $\displaystyle n_{\phi}$
$\displaystyle=$ $\displaystyle\sum_{i}g_{i\phi}n_{i}^{(v)}~{}.$ (26)

The rearrangement contributions

$V^{(r)}=\Gamma_{\omega}^{\prime}n_{\omega}\omega+\Gamma_{\rho}^{\prime}n_{\rho}\rho-\Gamma_{\sigma}^{\prime}n_{\sigma}\sigma$
(27)

and

$W_{i}^{(r)}=n_{S}^{(s)}\frac{\partial\Delta m_{S}}{\partial n_{i}^{(v)}}$
(28)

in eqs. (6) and (8) arise due to the dependence of the couplings and mass
shifts, respectively, on the densities.

###  III.2 Particle densities and relativistic density functional

The densities of the particles depend on their effective chemical potentials
$\mu_{i}^{\ast}$ and effective masses $m_{i}^{\ast}$. In a relativistic model,
vector and scalar densities have to be distinguished. The vector density of a
fermion is given by

$\displaystyle n_{i}^{(v)}$ $\displaystyle=$ $\displaystyle
g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}\>\theta(p_{i}-p)=\frac{g_{i}}{6\pi^{2}}p_{i}^{3}$
(29)

with degeneracy factor $g_{i}$ and an integration up to the Fermi momentum

$p_{i}=\left\\{\begin{array}[]{lll}\sqrt{[\mu_{i}^{\ast}]^{2}-[m_{i}^{\ast}]^{2}}&\mbox{if}&\mu_{i}^{\ast}>m_{i}^{\ast}\\\
0&\mbox{if}&\mu_{i}^{\ast}\leq m_{i}^{\ast}\end{array}\right.$ (30)

at zero temperature. The expression for the scalar density of a fermion

$\displaystyle n_{i}^{(s)}$ $\displaystyle=$ $\displaystyle
g_{i}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{m_{i}^{\ast}}{E(p,m_{i}^{\ast})}\>\theta(p_{i}-p)$
(31) $\displaystyle=$
$\displaystyle\frac{g_{i}}{4\pi^{2}}\left[\mu_{i}^{\ast}p_{i}-(m_{i}^{\ast})^{2}\ln\frac{\mu_{i}^{\ast}+p_{i}}{m_{i}^{\ast}}\right]$

contains an additional factor $m_{i}^{\ast}/E$ in the integral as compared to
the vector density with the energy

$E(p,m_{i}^{\ast})=\sqrt{p^{2}+[m_{i}^{\ast}]^{2}}$ (32)

depending on the momentum $p$, thus $\mu_{i}^{\ast}=E(p_{i},m_{i}^{\ast})$
with the Fermi momentum $p_{i}$. For the bosonic sexaquark (S), there is only
a finite density at zero temperature when the condition
$\mu_{S}^{\ast}=m_{S}^{\ast}$ is satisfied and we have Bose-Einstein
condensation. In this case the vector and scalar densities are identical and
can be written as

$n_{S}^{(v)}=n_{S}^{(s)}=\left\\{\begin{array}[]{lll}g_{S}\xi_{S}&\mbox{if}&\mu_{S}^{\ast}=m_{S}^{\ast}\\\
0&\mbox{else}&\end{array}\right.$ (33)

with degeneracy $g_{S}=1$ and a factor $\xi_{S}$. The latter can be determined
from the given baryon density $n_{b}$ as

$\displaystyle\xi_{S}$ $\displaystyle=$
$\displaystyle\left(n_{b}-n_{n}^{(v)}-n_{p}^{(v)}-n_{\Lambda}^{(v)}-n_{\Sigma^{+}}^{(v)}-n_{\Sigma^{0}}^{(v)}-n_{\Sigma^{-}}^{(v)}\right.$
(34) $\displaystyle\left.-n_{\Xi^{0}}^{(v)}-n_{\Xi^{-}}^{(v)}\right)/g_{S}$

at the condensation point.

The relativistic density functional assumes the form of a grand canonical
thermodynamic potential density in the present approach. It is given
explicitly by the expression

$\Omega(\\{\mu_{i}\\})=\sum_{i\in\mathcal{F}}\Omega_{i}+\Omega_{S}+\Omega_{\rm
meson}-\Omega_{\rm meson}^{(r)}-\Omega_{\rm mass}^{(r)}$ (35)

with a contribution

$\Omega_{i}=-\frac{1}{4}\left[\mu_{i}^{\ast}n_{i}^{(v)}-m_{i}^{\ast}n_{i}^{(s)}\right]\>,$
(36)

of the fermionic quasi particles
$i\in\mathcal{F}=\left\\{p,n,\Lambda,\Sigma^{+},\Sigma^{0},\Sigma^{-},\Xi^{0},\Xi^{-}\right\\}$.
The condensate contribution of the bosonic sexaquarks is formally given by

$\Omega_{S}=g_{S}\xi_{S}\left(m_{S}^{\ast}-\mu_{S}^{\ast}\right)$ (37)

and the meson contribution has the form

$\displaystyle\Omega_{\rm meson}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(\Gamma_{\omega}n_{\omega}\omega+\Gamma_{\rho}n_{\rho}\rho+\Gamma_{\rho}n_{\rho}\rho-\Gamma_{\sigma}n_{\sigma}\sigma\right)\>.$

Finally, the meson rearrangement contribution

$\Omega_{\rm meson}^{(r)}=V^{(r)}n_{cpl}$ (39)

and the mass shift rearrangement contribution

$\Omega_{\rm mass}^{(r)}=\sum_{i=n,p,S}n_{i}^{(v)}W_{i}^{(r)}$ (40)

appear in eq. (35). With the above definition of $\Omega$ the standard
relations

$n_{i}^{(v)}=-\left.\frac{\partial\Omega}{\partial\mu_{i}}\right|_{\\{\mu_{j}\\}_{j\neq
i}}\qquad n_{i}^{(s)}=\left.\frac{\partial\Omega}{\partial
m_{i}}\right|_{\\{\mu_{j}\\}}$ (41)

are valid as required for thermodynamic consistency. It is easily noted that
$\Omega_{S}=0$ because either $\xi_{S}=0$ for $m_{S}^{\ast}\neq\mu_{S}^{\ast}$
or $m_{S}^{\ast}=\mu_{S}^{\ast}$ with $\xi_{S}>0$ if condensation occurs.
Generally, the condition $\mu_{S}^{\ast}\leq m_{S}^{\ast}$ applies to the
model with sexaquarks. Since S is uncharged and does not couple to the mesons,
this corresponds to

$B_{S}\mu_{b}\leq m_{S}-S_{S}+W_{S}^{(r)}=m_{S}+\Delta
m_{S}+n_{S}^{(s)}\frac{\partial\Delta m_{S}}{\partial n_{S}^{(v)}}$ (42)

with baryon number $B_{S}=2$. Thus there is a condition on the baryon chemical
potential $\mu_{b}$ that is limited from above. If this maximum is reached,
the baryon density can only rise by an increase of the condensate density of
the sexaquark. With an increase of its effective mass due to the mass shift,
the condensation point moves to higher baryon chemical potentials.

###  III.3 Parameters and thermodynamic properties

The rest masses $m_{i}$ of nucleons, hyperons and mesons used in the GRDF are
given in tables 1 and 2.

Table 1: Rest masses of nucleons and hyperons in the generalized relativistic density functional. The rest mass of the sexaquark depends on the scenario, see text. particle $i$ |  $m_{i}$ [MeV] | particle $i$ |  $m_{i}$ [MeV]
---|---|---|---
p | 938.272081 | $\Sigma^{0}$ | 1192.642
n | 939.565413 | $\Sigma^{-}$ | 1197.449
$\Lambda$ | 1115.683 | $\Xi^{0}$ | 1314.86
$\Sigma^{+}$ | 1189.37 | $\Xi^{-}$ | 1321.71
Table 2: Rest masses of mesons in the generalized relativistic density functional. meson $i$ | $\sigma$ | $\omega$ | $\rho$ | $\phi$
---|---|---|---|---
$m_{i}$ [MeV] | 546.212459 | 783 | 763 | 1019.461

The baryon masses follow the recommendations of the Particle Data Group Zyla
_et al._ (2020) except for the sexaquark that depends on the considered
scenario, see section V. Here we assume the lowest mass of S that can be
considered is determined by the value that avoids a condensation at baryon
densities below the nuclear saturation density in symmetric nuclear matter in
the DD2 model. The masses of the $\omega$ and $\rho$ mesons are standard
values used in relativistic density functionals and $m_{\sigma}$ is a fit
parameter of the DD2 interaction. The mass of the $\phi$ meson as well as
those of the electron and muon are taken from Zyla  _et al._ (2020). Values of
the baryon-meson couplings $\Gamma_{m}$ at the reference density $n_{\rm
sat}^{(v)}$ and the coefficients in the functions $f_{m}$, cf. eqs. (11) and
(12), can be found in Typel  _et al._ (2010) and we recall that
$\Gamma_{\phi}=\Gamma_{\omega}(n_{\rm sat}^{(v)})$. The reference density
$n_{\rm sat}^{(v)}=0.149065$ fm-3 is the nuclear saturation density of the DD2
parameterization corresponding to a scalar density of $n_{\rm
sat}^{(s)}=0.139650$ fm-3. The values $U_{\Lambda}=-28$ MeV, $U_{\Sigma}=30$
MeV, and $U_{\Xi}=-14$ MeV from Raduta  _et al._ (2020) are adopted for the
hyperon potentials at saturation in symmetric nuclear matter to calculate the
coupling factors (16). The parameter $x_{S}$ in the mass shift (17) is varied
in the model to explore changes of the EoS caused by different onsets for the
condensation of sexaquarks.

All thermodynamic properties of the system can be derived from the grand
canonical thermodynamic potential density (35), e.g., the pressure $P$ is
simply given by $P=-\Omega$. The free energy density $\varepsilon$ is equal to
the internal energy density $f$ and can be calculated easily as

$f=\varepsilon=\Omega+\sum_{i}\mu_{i}n_{i}^{(v)}\>.$ (43)

The speed of sound in unit of the speed of light is found from

$c_{s}^{2}=\frac{dP}{d\varepsilon}$ (44)

and should not exceed one to be physically permitted.

##  IV High density equation of state

###  IV.1 Quark matter EoS

At low temperatures and very high densities, it has been proven that the
favorite state of three flavor color-superconducting quark matter would be the
color-flavor-locking (CFL) phase Alford  _et al._ (1999). Then one could ask
what will happen to the quark matter state at lower densities, in particular
in the region of the quark-hadron phase transition. In order to answer this
question, one has to solve self-consistently the coupled gap equations for
quark masses and diquark pairing gaps in the competing channels. Such a task
has been attacked within the Nambu–Jona-Lasinio (NJL) model for dense quark
matter Buballa (2005). The solution which for the case of three flavors and
three colors was presented for the first time in Refs. Ruester  _et al._
(2005); Blaschke  _et al._ (2005); Abuki and Kunihiro (2006) shows that there
is a corridor for the 2-flavor color superconducting (2SC) phase between
hadronic matter and the CFL phase.

We consider here the scenario that at high enough density S particles are
dissociated into a 2SC phase of quark matter which is microscopically
described by non-local NJL (nlNJL) model while its formulation has been fitted
Antić  _et al._ (2021) to the simple parameterization of constant speed of
sound (CSS) Zdunik and Haensel (2013) at higher densities and zero
temperature. The microscopic nlNJL model is a covariant extension of the NJL
model in which the quark fields interact via momentum dependent vertices
Schmidt  _et al._ (1994); Gomez Dumm  _et al._ (2006); Blaschke  _et al._
(2007). It has been shown in Antić  _et al._ (2021) that the EoS obtained from
this method is fitted very well (regarding the $\chi^{2}$ value) to CSS
parameterization for quark matter EoS. Indeed, a mapping from the
2-dimensional space of nlNJL model to the 3-dimensional space of CSS model has
been performed in this paper using a simple functional. With this functional
that is introduced in Antić  _et al._ (2021), the EoS of color superconducting
quark matter can be obtained for each value of vector meson coupling
($\eta_{V}$) and diquark coupling ($\eta_{D}$). The resulting EoS which has a
microscopic justification, is fitted to the following parameterization for CSS
quark matter

$P_{QM}(\mu)=A(\mu/\mu_{0})^{(1+1/c_{s}^{2})}-B,$ (45)

where $\mu_{0}=1000$ MeV defines a scale for the chemical potential. In (45),
there are three free parameters which have to be defined. The squared speed of
sound $c_{s}^{2}=dP/d\varepsilon$ has to be large enough so that the maximum
mass $M_{\rm max}$ of the corresponding neutron star sequence exceeds the
observational lower bound on it, which presently amounts to $2.01~{}M_{\odot}$
at 68.3% credibility Fonseca  _et al._ (2021). The prefactor $A$ could change
the slope of the $P-\mu$ line and the effective bag pressure $B$ imposes the
confinement effect at low densities in quark matter EoS and makes the pressure
to have negative value in this region. Once the pressure as a function of
chemical potential is obtained, one could calculate the baryon density
($n_{b}=dP/d\mu$) as well as the energy density ($\varepsilon=n_{b}\mu-P$) for
the quark matter. For the relation between the nlNJL and CSS models of quark
matter, see also the recent work Contrera  _et al._ (2022).

In Antić  _et al._ (2021), for nlNJL parameters $0.70<\eta_{D}<0.80$ and
$0.11<\eta_{V}<0.18$, the parameter mapping to the space of the CSS model
resulted in the following range of parameters: $0.449<c_{s}^{2}<0.541$,
$91.484<A[{\rm MeV/fm}^{3}]<101.116$ and $82.437<B[{\rm MeV/fm}^{3}]<92.290$.
In the current work, we use several different sets of these parameters in
order to find out which of them is modeling the hybrid star with S in the
outer core in such a way to fulfill all observational constraints.

We consider four different scenarios for solving the S dilemma in compact
stars denoted as follows, where the subscript $i$ labels two mass choices for
the S ($i=1885,1941$), in MeV:

  * •

Scenario DSiY: An early deconfinement (D) before sexaquark (Si) and hyperon
(Y) onset,

  * •

Scenario SiDY: deconfinement after sexaquark but before hyperon onset,

  * •

Scenario SiYD: deconfinement after both sexaquark and hyperon onset when S is
prior,

  * •

Scenario YSiD: like previous case but hyperon is prior.

In this work, the value of $\eta_{D}$ has been taken to be $0.75$ in
accordance to the Fierz transformation of a one-gluon exchange interaction
Buballa (2005).

For the present hadronic matter EoS and the quark matter EoS, the transition
from hadronic matter to deconfined quark matter occurs at low chemical
potential. Moreover, because of the softening of the hadronic matter EoS at
higher chemical potential after the appearance of hyperons, we encounter the
reconfinement problem Zdunik and Haensel (2013) which will be ignored in
applying the Maxwell construction following the argument Shahrbaf  _et al._
(2020a) that once the deconfinement occurs at a certain critical density, the
hadronic matter EoS is not valid anymore beyond that point.

For moving the transition point to higher chemical potentials, we may apply an
extra bag pressure to the quark matter EoS which could be either a constant
one ($B_{0}$) or a $\mu-$dependent one. The $\mu-$dependent bag pressure
$B_{\rm eff}(\mu)$ is efficient for making a stronger phase transition with a
big jump in density (a bigger difference in the slope of hadronic matter and
quark matter EoS at the transition point) at the transition point.

The larger the density jump, the stronger the effect of softening the EoS
which leads to more compact hybrid stars at the transition. This may help to
fulfill the demanding constraint on the tidal deformability,
$70<\Lambda_{1.4}<580$ at $1.4~{}M_{\odot}$, deduced from the gravitational
wave signal of the binary neutron star merger GW170817 Abbott  _et al._
(2018a). For a sufficiently large jump in energy density
($\Delta\varepsilon\gtrsim\varepsilon_{c}$, see Alford  _et al._ (2013) for
details), a gravitational instability may occur which, after stability is
recovered for sufficiently stiff quark matter at higher densities, leads to
the formation of an alternate branch of hybrid stars in the $M-R$ diagram
Gerlach (1968). This potentially observable feature of a strong phase
transition is accompanied with the phenomenon of mass-twin stars Glendenning
and Kettner (2000). Due to the constraint of causality on the speed of sound,
however, a parameterization of $B_{\rm eff}(\mu)$ with a sufficiently large
value of $\Delta\varepsilon$ may not be admissible.

The $\mu$-dependence of the bag pressure has been introduced in Alvarez-
Castillo  _et al._ (2019) in the following form

$B_{\rm eff}(\mu)=B_{1}f_{<}(\mu)~{},$ (46)

with the switching functions

$f_{<}(\mu)=\frac{1}{2}\left[1-\tanh\left(\frac{\mu-\mu_{<}}{\Gamma_{<}}\right)\right]~{}.\\\
$ (47)

While in Alvarez-Castillo  _et al._ (2019) the parameterization of $B_{\rm
eff}(\mu)$ resulted in mass-twins, in the present work (as also in Ref.
Contrera  _et al._ (2022)), it will suffice to accommodate the tidal
deformability constraint but not lead to mass-twins. It is worth mentioning
that in (45), the value of $B$ reads

$B=B_{0}+B_{\rm eff}(\mu)$ (48)

in which $B_{0}$ is a constant value and $B_{\rm eff}$ is defined according to
(46). We have changed the parameters in above definition as well as the
parameters of quark matter EoS several times to find the best sets (to the
best of our knowledge) of parameters which not only fulfills all observational
constraints but also results in an appropriate transition point to be a proper
solution for one of the considered scenarios in this work.

###  IV.2 Hybrid EoS

Within the two-phase approaches to hybrid neutron star matter, the most common
phase transition constructions are the Maxwell construction (MC) and the
recently developed crossover interpolations, see Baym  _et al._ (2018) for a
recent review. While the MC has been in use already since the early days of
discussing quark deconfinement in neutron stars Baym and Chin (1976) with
varying but physically equivalent formulations, for the crossover
interpolations a word of caution may be in order.

The idea of the interpolation construction pioneered in Masuda  _et al._
(2013a, b) is to facilitate a thermodynamically consistent transition from a
relatively soft hadronic EoS that could be trusted up to and slightly beyond
saturation density, to a stiff quark matter EoS with a region of validity
above 3-5 times saturation density. Strictly speaking, in the crossover
transition region neither of the input EoS are trustworthy. Therefore, no
Maxwell construction would be applicable and alternatives have been developed
to interpolate between soft hadronic and stiff quark matter EoS in the
crossover region. The first construction Masuda  _et al._ (2013b) was not
formulated in natural variables of the thermodynamic potential, which was
corrected in the second version Masuda  _et al._ (2013a), but both
interpolations were defined by a “mixing” of the input EoS, even in the
crossover region where they are strictly speaking not applicable. Furthermore,
the weight functions for the mixing were chosen as Fermi functions with the
unphysical implication of a nonzero probability for quark matter at low and
hadronic matter at high densities. An alternative to facilitate treatment of
the hadron-to-quark matter crossover is the replacement interpolation
construction (RIC) developed in Ayriyan  _et al._ (2018); Abgaryan  _et al._
(2018) to describe the situation of a mixed phase due to pasta structures in
the hadron-to-quark matter transition. The RIC treatment was applied to the
stiff-soft transition case in the context of solving the hyperon puzzle
Shahrbaf  _et al._ (2020a) without introducing an unphysical high-density
hadronic component. A more general two-zone interpolation method that followed
the intentions discussed in Baym  _et al._ (2018) has recently been developed
in Ayriyan  _et al._ (2021) and used within a Bayesian analysis of neutron
star constraints on the EoS.

Here we employ both the MC as well as RIC for the first order phase transition
from DD2Y-T+S EoS to CSS quark matter EoS, to investigate the appearance of S
particles in hybrid stars. To do this, we try to find a stable hybrid star
with the 2SC phase of quark matter in the core surrounded by a layer of
nuclear matter or hypernuclear matter which includes S.

####  IV.2.1 Maxwell construction

In the MC, the Gibbs conditions for phase equilibrium have to be fulfilled
globally, i.e.

$\displaystyle\mu_{H}$ $\displaystyle=$ $\displaystyle\mu_{Q}=\mu_{c}~{},$
(49) $\displaystyle T_{H}$ $\displaystyle=$ $\displaystyle T_{Q}=T_{c}~{},$
(50) $\displaystyle P_{H}(\mu_{b},\mu_{e})$ $\displaystyle=$ $\displaystyle
P_{Q}(\mu_{b},\mu_{e})=P_{c}~{}.$ (51)

The above conditions guarantee the chemical, thermal and mechanical phase
equilibrium at the transition point between hadronic phase and quark matter
phase for which the critical values of the thermodynamic variables are denoted
with the subscript $c$. We apply the conditions for the zero temperature case
in the present work.

Within the MC scheme, the reconfinement phenomenon Zdunik and Haensel (2013)
can occur which consists in a second (unphysical) crossing of the hadronic and
quark matter EoS. It can be dealt with in the MC by ignoring the unphysical
crossing. The reconfinement situation can also be removed by applying a RIC
with negative value of $\Delta_{p}$, as it has been discussed in Shahrbaf  _et
al._ (2020b).

####  IV.2.2 Replacement interpolation construction

Within the RIC, we assume that neither the hadronic matter nor the quark
matter EoS are reliable in a corridor around their unphysical crossing, e.g.,
in the reconfinement region. Therefore, two boundary points, i. e., $\mu_{H}$
before the reconfinement and $\mu_{Q}$ after the reconfinement, are defined so
that the interpolated pressure could be described by a parabolic form between
them,

$P_{M}(\mu)~{}=~{}\alpha_{2}(\mu-\mu_{c})^{2}~{}+~{}\alpha_{1}(\mu-\mu_{c})~{}+~{}(1+\Delta_{P})P_{c},$
(52)

where the $\mu_{c}$ and $p_{c}$ corresponds to the critical point at which the
wrong Maxwell transition from quark matter to hadronic matter has occurred.
The $\alpha_{1}$, $\alpha_{2}$ as well as $\mu_{H}$ and $\mu_{Q}$ could be
obtained from the continuity conditions at the borders of the mixed-phase

$\displaystyle P_{H}(\mu_{H})$ $\displaystyle=$ $\displaystyle
P_{M}(\mu_{H})~{},$ (53) $\displaystyle P_{Q}(\mu_{Q})$ $\displaystyle=$
$\displaystyle P_{M}(\mu_{Q})~{},$ (54)
$\displaystyle\frac{\partial}{\partial\mu}P_{H}(\mu_{H})$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\mu}P_{M}(\mu_{H})~{},$ (55)
$\displaystyle\frac{\partial}{\partial\mu}P_{Q}(\mu_{Q})$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\mu}P_{M}(\mu_{Q})~{}.$ (56)

The RIC allows us to have a solution for which the quark core of NS is
surrounded by a mixture of usual nuclear matter, sexaquark and even hyperons
if the upper boundary of the mixed phase goes beyond the hyperon onset.
Therefore, we can have solutions for all scenarios in this way. Moreover,
changing the parameters of the RIC enable us to have a larger part of the
sequence in the sexaquark-nuclear phase than in the MC case.

##  V Purely hadronic stars and their problems

Since numerous abbreviations have been used in this work, we collect them all
together for the convenience of the readers in Appendix A.

With the generalized relativistic density functional (GRDF) for hadronic
matter at hand, it is possible to explore the EoS and corresponding properties
of neutron stars. There are various scenarios to be distinguished in the
following. The simplest case with only nucleons and leptons corresponds to the
original DD2 model as presented in Typel  _et al._ (2010). Adding hyperons the
model is called DD2Y-T to distinguish it from the similar model DD2Y
introduced in Marques  _et al._ (2017). The DD2Y-T predictions were compared
already to other EoS models with hyperons in Stone  _et al._ (2021). Finally,
after including also the sexaquark, there is the full model which will be
denoted DD2Y-T+S in the following.

To start, the mass of the sexaquark is taken to be constant. When the mass of
the sexaquark is equal to $2054$ MeV, the onset of the S (sexaquark) occurs at
about 0.25 fm-3 (more than $1.5$ times the saturation density) with an
immediate appearance of the BEC and the neutron star will collapse with
increasing central density because the pressure remains constant.

The obvious suggestion is including the repulsive interaction of S by
considering an increasing density-dependent mass for S particles, as done for
hyperons. Considering the S particle as a composite system of diquarks
embedded in the baryonic medium suggests that a dissociation of the bound
state into a many-body correlation in the continuum could be studied in a
microscopic model. But this requires intricate calculation that are not in the
scope of the present work and we leave it for the future. In this work, we
only consider the S particle as a point-like boson with a medium dependent
mass. As it was mentioned in subsection A, we suggest that all possible
substructures can lead to a medium dependent mass of S particles.

Hereafter, we consider a linear density dependence of the mass of the S in
which the mass shift is positive so that the instability problem of NS is
prevented since the corresponding increase of the chemical potential following
the condensation criterion allows a rising pressure. At zero density, the
effective mass of S particles, i.e., $m_{S}^{\ast}$ is given by the vacuum
mass of S particle, $m_{S}$, and the rate of increase that is given by a
factor $x_{S}$ times the baryon density $n_{b}$ in the unit of the reference
density $n_{0}=0.15$ fm-3,

$m^{*}_{S}(n_{b},x_{S})=m_{S}\left(1+x_{S}\frac{n_{b}}{n_{0}}\right)\>.$ (57)

This assumption results in an increase of the onset density of S formation and
subsequent condensation, so that there is still an increase of the pressure at
higher densities.

An estimate of reasonable values for $x$ can be derived from the change of the
effective potential $U_{i}=V_{i}-S_{i}$ of the octet baryons with the baryon
density $n_{b}$ as

$x_{i}=\frac{n_{0}}{m_{i}}\frac{dU_{i}}{dn_{b}}\>.$ (58)

Figure 2: The effective slope of the mass shift for octet baryons calculated
based on the effective potentials the in DD2Y-T model.

These effective slopes calculated in the DD2Y-T model are depicted in Fig. 2.
The slopes for the various hyperons are typically large and positive near
nuclear saturation density and transition to being slowly increasing beyond
$n_{b}\approx 0.4\,{\rm fm}^{-3}$. As we shall see, without a transition to a
deconfined phase the observed range of neutron star masses cannot be
reproduced, so only a portion of the range shown in Fig. 2 will prove to be
physically relevant. It should also be noted that the constant-slope
approximation we use for the S in the present work may be an over-
simplification as it is only applicable for hyperons for relatively small
ranges of $n_{b}$.

Figure 3: Panel (a) shows double baryonic chemical potential as a function of
baryon density while panel (b) shows the pressure as a function of baryon
density for DD2, DD2Y-T and DD2Y-T+S for three different slopes of the mass
shift of S particles.

The onset density for the condensation of S depends on its rest mass in vacuum
and the slope parameter $x$. One can see this onset by plotting $2\mu_{b}$ as
a function of density for different parameter choices, as shown in the panel
(a) of Fig. 3. The onset density occurs as soon as the baryo-chemical
potential of the hadronic matter fulfills $2\mu_{b}=m_{S}$. The full and
dashed black lines in figure 3 show the value of $2\mu_{b}$ as a functions of
$n_{b}$ for the DD2 and DD2Y-T models. The effective mass of S is also plotted
in the same figure, for three values of the vacuum rest mass and several slope
parameters. When the vacuum mass of S is $\leq 1885$ MeV, Bose condensation
occurs for any density in the absence of repulsive interactions between S and
nucleons and more careful modeling would be required. Thus we do not consider
lower masses for the S here.

For $m_{S}=1941$ MeV, condensation occurs exactly at nuclear saturation
density if $x=0$, while for $S_{2054}$ the intersection point with the black
lines occurs at a higher baryon density of about $0.25$ fm-3. Obviously, the
onset of condensation occurs at higher densities as $x$ increases. Therefore,
varying the value of $x$ results in different scenarios for NS regarding the
priority of Y onset and S onset if the density-dependence of hyperon masses is
taken as given.

Different behaviors of the effective mass of S are reflected in how the
pressure changes with baryon density, as shown in panel (b) of Fig. 3 for
illustrative cases. The pressure is constant beyond the onset density for Bose
condensation causing an instability of the neutron star calculated with the
corresponding equations of state. For positive $x$, the increase of the
effective mass causes a rise of the pressure; this could also be attributed to
the effect of a repulsive potential at short distances.

Figure 4: Particle fraction of hadronic matter including an increasing S mass
in $\beta$-equilibrium for NS matter. The solid lines correspond to S1885
while the dashed lines corresponds to S1941. DD2Y-T, DD2Y-T+S with $x=0.05$
and DD2Y-T+S with $x=0.08$ are shown in panel (a), (b) and (c) respectively.

The particle fractions as a function of baryon density in $\beta$ equilibrated
matter are shown in Fig. 4 where S is taken to have an effective mass that
increases with the density. The dashed lines show the results for the heavier
S, i.e., S1941 ($m_{S}=1941$ MeV) and the solid lines correspond to the
lighter S, i.e., S1885 ($m_{S}=1885$ MeV). As the figure shows, for $x=0.05$
the S particles appear before hyperons but for $x=0.08$, S onset occurs after
the hyperon onset. However, for the lower mass of S, the hyperon onset is very
close to the S onset and the only hyperon which appears before S, is the
lightest one, i.e., $\Lambda$. It is concluded from Fig. 4 that the appearance
of S affects the hyperon onset. When there is a delay in S1941 onset because
of the larger mass, hyperons appear at lower densities compared to the S1885
case.

The speed of sound $c_{s}$ is a quantity that can be derived directly from the
pressure and energy density, see equation (44). It has to be smaller than the
speed of light to have a causal EoS.

Figure 5: The square of speed of sound as a function of energy density for the
purely hadronic EoS cases.

The square of $c_{s}$ for the hadronic equations of state is depicted in Fig.
5 showing that the causality condition is fulfilled in the whole range of
energy densities for all combinations of S mass and slope parameter $x$.

###  V.1 Sexaquark and hyperon dilemma in neutron stars

Each combination of $m_{S}$ and $x$ leads to an equation of state of neutron
star matter and corresponding mass-radius relation that can be compared to the
predictions of the original DD2 and DD2Y-T models without sexaquark. In Fig.
6, the mass-radius (panel (a)) and mass-density (panel (b)) relations for
neutron stars obtained from DD2 and DD2Y-T as well as DD2Y-T+S are shown when
the S particles are supposed to have a mass that increases with density.

Figure 6: Mass-radius relation (panel (a)) and mass versus central density
(panel (b)) for compact stars in which the S particle is supposed to have an
effective mass $m_{S}^{\ast}(n_{b},x)$ given by equation (57) with
$m_{S}=1941$, and $1885$ MeV and $x=0.03,0.05,0.08$. The solid lines
correspond to S1885 while the dashed lines correspond to S1941. For a
comparison the new $1-\sigma$ mass-radius constraints from the NICER analysis
of observations of the massive pulsar PSR J0740+6620 Fonseca  _et al._ (2021)
are indicated in red Riley  _et al._ (2021) and blue Miller  _et al._ (2021)
regions. Additionally, the green bar marks the radius of a $1.4$ solar mass
neutron star from a joint analysis of the gravitational-wave signal GW170817
with its electromagnetic counterparts at $90\%$ confidence Dietrich  _et al._
(2020). The cyan region is from the NICER mass-radius measurement on PSR
J0030+0451 Miller  _et al._ (2019). The gray and light orange regions
corresponds to the estimates of the components of the binary system labeled as
$M_{1}$ and $M_{2}$ of the GW170817 merger Abbott  _et al._ (2018b).

The S dilemma for the DD2Y-T+S models, i.e., the decrease of the maximum
neutron star masses due to the softening of the EoS with the appearance of S,
is also seen. As the figure shows, for $x=0.03,0.05$ the maximum mass of the
neutron stars does not reach 2M⊙ and more repulsive effects are needed in
hadronic matter to reach the lower bound of the maximum mass of neutron stars,
as exemplified by $x=0.08$.

Another important constraint for an EoS is the tidal deformability of a
neutron star, deduced by the LIGO and Virgo collaborations from GW170817
Abbott  _et al._ (2018a) and translated to a multi-messenger constraint on the
radius at $1.4~{}M_{\odot}$ in Dietrich  _et al._ (2020), shown as green bar
in panel (a) of Fig. 6. The only line which crosses it is the blue solid line
corresponding to DD2Y-T+S1885 with $x=0.05$.

Figure 7: (a): Tidal deformability parameters $\Lambda_{1}$ and $\Lambda_{2}$
of the high- and low-mass components of the binary merger. The results have
been shown for $m_{S}=1941,1885$ MeV and $x=0.050,0.08$. The solid lines
correspond to S1885 while the dashed lines correspond to S1941. The green
region is the placed constraints on the tidal effects by LIGO and Virgo
collaboration from GW170817 Abbott  _et al._ (2018a). (b): Dimensionless tidal
deformability $\Lambda$ as a function of the star mass. The green line shows
the $\Lambda_{1.4}$ constraint from the low-spin prior analysis of GW170817
Abbott  _et al._ (2018a).

The tidal deformability parameters $\Lambda_{1}$ and $\Lambda_{2}$ of the
high- and low-mass components of the binary merger event have been constructed
from the $\Lambda(M)$ relation Abbott  _et al._ (2018a). We show in Fig. 7 the
observational constraints for the tidal deformability parameters along with
predictions for the various equations of state.

As can be seen, the only EoS which lies in the green credibility region, is
DD2Y-T+S when $m_{S}=1885$ MeV and $x=0.05$. However, this model does not
fulfill the constraint on the maximum mass, $M_{\rm max}\geq 2~{}M_{\odot}$
Fonseca  _et al._ (2021), see Figs. 6 and 7 panel (b). It should also be
stressed that the standard hadronic matter model, DD2Y-T, also fails to
satisfy the tidal deformability constraint.

Therefore we conclude that within our setting, where the stiff nucleonic RMF
parameterization ”DD2” (that violates the tidal deformability constraint) sets
the hadronic matter baseline, a purely hadronic EoS is ruled out – independent
of the existence of a sexaquark.

In Fig. 6, some special values for the slope factor $x$ and $m_{S}$ have been
selected based on the range which fulfills the constraints of the combination
of $m_{S}$ and $x$. These constraints for $x$ as a function of $m_{S}$ are
shown in Fig. 8. The boundaries are defined as follows:

  1. 1.

S appears below the maximum mass of NS in matter with only nucleons: region
above blue solid line is excluded.

  2. 2.

S appears below the maximum mass of NS in matter with nucleons and hyperons:
region above red solid line is excluded.

  3. 3.

Condensation of S occurs only above the saturation density: region below
violet solid line is excluded.

  4. 4.

Condensation of S doesn’t occur below the saturation density: region left of
yellow solid line is excluded.

  5. 5.

No decreasing mass of S (i.e., constant pressure above density of condensation
and below endpoint of neutron star $M$-$R$ relation): region below $x=0$ is
excluded.

  6. 6.

Avoiding an unstable S that decays in $10^{-10}$s: region above $2054$ MeV is
excluded ($m_{S}>m_{p}+m_{e}+m_{\Lambda}=2054$ MeV)

Therefore, only the white region remains. However, the border of the white
region which lies on the line $x=0$ is also strictly excluded within the
framework of our analysis because it leads to Bose-Einstein condensation (BEC)
of a constant-mass S which limits the maximum mass of purely hadronic neutron
stars to $M_{\rm max}\lesssim 0.7~{}M_{\odot}$. A scenario according to which
such a BEC of S triggers a catastrophic rearrangement of the neutron star
structure (star quake) which entails the conversion to a hybrid neutron star
with a color superconducting (CFL) quark matter core has recently been
considered in Blaschke  _et al._ (2022).

Figure 8: Allowed (white) and excluded (light and dark grey) regions in the
plane of $m_{S}$ and $x$ values as defined by constraints on the sexaquark
condensation. See text for details and explanation of the additional lines.
The region of masses $m_{S}>m_{p}+m_{e}+m_{\Lambda}=2054$ MeV is also excluded
by the requirement of sexaquark stability on cosmological timescales, if S is
the dark matter. For further details, see text.

In addition to these boundaries, three dotted vertical lines for specific
values of the sexaquark mass are shown. Furthermore the green line corresponds
to the condition of S condensation at the onset of hyperons. The sexaquark
appears before (after) hyperons for points below (above) this line. Dashed
blue and red lines, which are almost on top of each other, show S condensation
at the center of a neutron star, which includes nucleons and nucleons-
hyperons, respectively, with a mass of $M=1.4$ M⊙.

In the next sections we discuss hybrid EoS which include a deconfined phase of
matter to understand if they favor, exclude or constrain the existence of a
sexaquark.

##  VI Deconfinement solution for the hyperon and sexaquark dilemmas

It was mentioned in section IV that two different types of construction are
used for dealing with hybrid stars. We compare their results regarding the
maximum mass, radius, tidal deformability and the sequence of particles onset
for constructed hybrid stars in this section.

The resultant parameters for the EoS which have been obtained as a hybrid
solution for different scenarios using the MC and RIC are listed in Table 3
and Table 4, respectively. The quark matter parameters have been selected and
tuned in such a way to describe a desirable critical point for the transition
from hadronic phase to quark phase. (The corresponding hadronic matter and
quark matter EoS for each scenario which includes S1941 are given in Appendix
B in Fig. 14 in the P-$\mu$ plane while the ones which include S1885 are shown
in Fig. 15.) As we already mentioned, the “reconfinement” phenomenon at higher
chemical potentials has been ignored when a Maxwell construction is employed.

The effect of applying a RIC is to shift the deconfinement onset to higher
densities than in a Maxwell construction which, therefore, allows S particles
to exist in the core of higher mass neutron stars – around the mass of the
lightest observed NS. Around the reconfinement point, neither the hadronic EoS
nor the quark matter one is reliable. Indeed, nature prefers the phase with
the higher pressure (transition occurs from a stiff EoS to a softer one) but
at the same time a transition from deconfined quark matter to a confined
hadronic phase is nonphysical. In this situation, the RIC with a negative
value of $\Delta_{p}$ in (52) is a mathematical approach that allows to
describe a crossover transition from the softer hadronic phase to the stiffer
quark phase. It is worth mentioning that RIC with a negative value of
$\Delta_{p}$ results in an interpolated EoS which is stiffer than both the
hadronic and quark matter EoS’s. One should notice that for each set of
hadronic EoS and quark matter EoS, there is always a minimum value of
$\Delta_{p}$ for which causality is fulfilled, i.e.,
$c_{s}^{2}={(n/\mu)}/{(d^{2}P/d\mu^{2})}<1$. Therefore not only
$d^{2}P/d\mu^{2}$ must be positive but also it has to be greater than $n/\mu$.
The first value of $\Delta_{p}$ in Table 4 for each set is the minimum allowed
value.

The mass-radius relations of the hybrid stars corresponding to the EoS plotted
in Fig. 14 are shown in Fig. 9.

It can be seen that all obtained solutions fulfill all observational
constraints regarding the mass-radius of neutron stars except for the recent
radius constraint for the $1.4$ solar mass neutron star obtained by Dietrich
et al. Dietrich  _et al._ (2020) which is not fulfilled by any of MC lines in
upper panels. Only MC-S1941-5 crosses it marginally while all hybrid solutions
with S1885 in lower panels cross this green bar. We expect this issue will
result in remarkable difference in tidal deformability plots for upper and
lower panels.

Figure 9: The upper panels show the mass-radius relation for the hybrid star
in each scenario which includes S1941 based on the Table 3. Mass-radius
relation for the constructed hybrid stars based on the Table 3 when S1885 is
included in hadronic EoS and also based on Table 4 are shown in the lower
panels. The Blue diamond and maroon circle show the S onset and Y onset
respectively in all panels. Moreover, DD2Y-T and DD2Y-T+S lines have been
shown with blue and black solid lines in all panels. Colored regions are
described in the caption of Fig. 6.

Figure 10: Tidal deformability parameters $\Lambda_{1}$ and $\Lambda_{2}$ for
all MC hybrid stars based on Table 3 are shown in panel (a). In panel (b) the
tidal deformability parameters $\Lambda_{1}$ and $\Lambda_{2}$ correspond to
all RIC hybrid stars based on Table 4 are shown. In panel (b), the pure
hadronic lines are also depicted for comparison. The only hadronic line that
fulfills the $90\%$ credibility constraint is DD2Y-T+S1885 with $x=0.05$ which
can not reach the necessary maximum mass for fulfilling the mass constraint of
neutron stars, See Fig. 6.

The maximum mass is obtained via MC for MC-S1941-5 and MC-S1941-3 because the
stiffest quark matter EoS with $\eta_{v}=0.17$ which results in the maximum
value of the speed of sound, i.e., $c_{s}^{2}/c^{2}=0.52$, has been used for
constructing these two solutions. As the lower panels of Fig. 9 show, the sets
that include S1885 are mainly a solution for the early deconfinement scenario
and the second crossing corresponding to reconfinement is also visible in the
$M-R$ curves. If one neglects this early deconfinement and apply RIC to the
second crossing, the appearance of S is also possible at interpolated region.
Evidently, the onset of deconfinement has been shifted to higher mass for RIC
compared to MC in the lower panels of Fig. 9 and therefore, the constructed
hybrid stars are a solution for S1885DY scenario instead of DS1885Y. Comparing
the three lower panels of this figure, one sees that for fixed mass of S,
increasing the slope of the mass shift, $x$, increases the reconfinement
density and therefore allows the interpolated region to cover higher masses of
hybrid stars, around $1.4M_{\odot}$.

The softer EoS has the smaller radius and therefore, around $1.4~{}M_{\odot}$,
the softest EoS should have the best tidal deformability. As Fig. 9 shows, the
$M-R$ curves with S1885 which are shown in the lower panels, have the smaller
radius. The tidal deformability has also been calculated for all obtained
hybrid stars and the results are plotted in Fig. 10. In this figure, panel (a)
shows the results for MC while the results for RIC are plotted in panel (b).

It can be seen in Fig. 10 that the tidal deformability results for MC
solutions with S1941 are either totally out of the green credibility region or
marginally cross it, especially the YS1941D scenario which is strongly
disfavored by the results of tidal deformability.

However for the hybrid star solutions with S1885, not only the M-R constraints
are fulfilled but also the tidal deformability puzzle has been solved and the
observational constraints from GW170817 with the $90\%$ credibility is
fulfilled. Both MC and RIC with S1885 result in hybrid stars compatible with
both mass-radius constraints and tidal deformability constraints. Since
(neglecting repulsion between S and nucleons) condensation of S occurs at
saturation density or below if $m_{S}\leq 1885$ MeV, and $m_{S}=1885$ MeV is
the mass for which the condensation occurs exactly at saturation density,
$m_{S}$ near 1885 MeV may be most favorable for compatability with
observational constraints of neutron stars Dietrich  _et al._ (2020).

Figure 11: Panel (a): Pressure as a function of energy density for the
obtained EoS. The horizontal lines show the Maxwell construction. The gray
region corresponds to the EoS constraint from Hebeler  _et al._ (2013) while
the narrow region between black dotted lines corresponds to the analysis of
Miller  _et al._ (2021). Panel (b): Squared sound speed for these EoS.

Panel (a) of Fig. 11 shows that all obtained EoS are almost inside the gray
region which has been introduced as the accepted region for $P-\varepsilon$
lines in Hebeler  _et al._ (2013). From panel (b) of this figure one concludes
that the causality constraint $c_{s}\leq 1$ (in units of the speed of light)
is fulfilled for the all solutions which include S1885.

As can be seen in the lower panels of Figs. 9 and panels (a) and (b) of Fig.
10, we have found two different types of stable hybrid stars including S,
which fulfill all modern observational constraints for pulsars. Taking
$m_{S}=1885$ MeV with medium-dependent mass, and using both the MC as well as
the RIC, we have a transition from DD2Y-T+S on the exterior to a 2SC phase
quark matter on the interior. We display the radial structure of a number of
illustrative cases, with the composition in panels $(a_{1})$, $(b_{1})$ and
$(c_{1})$ and energy density and pressure in panels $(a_{2})$, $(b_{2})$ and
$(c_{2})$. The stars in Fig. 12 do not satisfy both the M-R and tidal
deformability constraints, while those in Fig. 13 satisfy all constraints.

Figure 12: The profile of the stars for DD2Y-T, DD2Y-T+S1885 and MC-S1941-4.
The particle fraction as a function of the distance from the center of star is
shown in panels $(a_{1})$ DD2Y-T, $(b_{1})$ DD2Y-T+S1885 and $(c_{1})$
MC-S1941-4. The pressure and the energy density of each star are shown in
panels $(a_{2})$, $(b_{2})$ and $(c_{2})$. The profiles correspond to
$M=2M_{\odot}$ and $M=1.4M_{\odot}$ for each star. The EoS for DD2Y-T+S1885
when $x=0.05$ is a soft one which cannot reach to $2M_{\odot}$ None of these
models satisfies both maximum mass and tidal deformability constraints.

Figure 13: Same as Fig. 12, but for MC-S1885-7 (panels $(a_{1})$ and
$(a_{2})$) and RIC-C-1 (panels $(b_{1})$ and $(b_{2})$), both of which satisfy
both the maximum mass and tidal deformability constraints. In the first model,
MC-S1885-7, outside the quark matter core there are only nucleons and leptons,
which in the RIC-C-1 model there is a layer outside the quark matter core
which has a significant sexaquark component.

In Fig. 12, panels $(a_{1})$ and $(a_{2})$ show the pure hadronic star from
DD2, panels $(b_{1})$ and $(b_{2})$ show the hyperonic star including S when
$m_{S}=1885$ MeV and $x=0.05$, and panels $(c_{1})$ and $(c_{2})$ show the
MC-S1941-4. Fig. 13 shows MC-S1885-7 and RIC-C-1 which are representative of
the different types of successful hybrid solutions with S. As expected there
is a jump in energy density for the MC solutions but we don’t see this jump in
RIC which is a crossover transition. Among the solutions for which the profile
of the star is plotted, the only solution which has a quark matter core for
$M=1.4M_{\odot}$ is MC-S1885-7 because of the early deconfinement.

One can see the sequence of the onset of different particles at different
radius in the profiles of the composition of each figure. The only star in
which $\Lambda$ hyperon could appear using the Maxwell construction is
MC-S1941-4; it gives a solution for YS1941D scenario but it is not compatible
with tidal deformability constraint. Employing RIC, a stable hybrid star with
a typical radius of about 12 km can be obtained; it has an inner core of quark
matter while a layer with S and neutrons connects it to a normal nuclear
matter crust. Not only does the sexaquark appear in this model but also all
constraints from neutron stars are well fulfilled.

Thus we have shown that – contrary to the expectation of McDermott  _et al._
(2019) that a deeply bound dibaryon would be incompatible with neutron stars –
the observed NS properties are compatible with the hadron spectrum containing
a stable sexaquark, for parameters similar to those commonly adopted for
hyperons and deconfined quark matter. Whether the S appears in neutron stars
or not depends on how early matter deconfines to quarks.

##  VII Conclusions

Motivated by a stable sexaquark444S=$uuddss$, with $m_{S}<2054$ MeV to assure
a lifetime greater than the age of the Universe Farrar (2022). being a good
dark matter candidate, we have investigated whether the maximum mass and tidal
deformability measurements of neutron stars are compatible with the existence
of such a particle. To address this question we extended a relativistic mean
field model of hadrons to include sexaquarks and used a causal
phenomenological EoS to describe the quark matter phase, restricting
parameters to the range suggested by two flavor superconductivity. A first-
principles treatment of the cross-over region between hadronic and quark
matter phases is not available for all EoS parameters of interest, so we
investigate two different approaches. When the conditions for applicability of
the Maxwell construction are met we use that, and otherwise we use a
Replacement Interpolation construction.

We find that the existence of a stable sexaquark is well-compatible with both
the maximum mass and the highly constraining tidal deformability of GW170817,
given present knowledge of the properties of hadronic and quark matter.
Whether or not a sexaquark exists, within our framework the most massive stars
must have a quark matter core with stiff EoS to support their high mass.
555There are purely hadronic models, some of them including hyperons, with a
softer EoS at high densities than the DD2 parameterization which can satisfy
the maximum mass and tidal deformability constraints. However, this requires a
modification of the high-density couplings with more parameters, e.g., DD2-F
Alvarez-Castillo  _et al._ (2016), or specifically adjusted hyperon-meson
couplings, e.g., Thapa  _et al._ (2021), or the introduction of additional
coupling mechanisms, e.g., Grigorian  _et al._ (2018). Among the successful
hybrid star solutions which we develop, we find two general types: those in
which the quark matter core is surrounded by nucleons, and others in which the
core and interpolated region is surrounded by a layer with a substantial
sexaquark (but negligible hyperon) fraction. This potentially may lead to an
observable signature of sexaquarks in the cooling curve or kilonova
properties.

Satisfying the tidal deformability constraints on neutron stars of $\approx
1.4~{}M_{\odot}$ while accommodating neutron star masses above $2~{}M_{\odot}$
is very challenging. At low density the EoS must be quite soft to produce
compact $\approx 1.4~{}M_{\odot}$ neutron stars, while at high density it must
be stiff to support high mass stars. Intriguingly, the sexaquark neatly solves
this problem by naturally producing the needed softening. However the
softening at low density needed to fit current deformability constraints can
also be achieved by an early onset of the quark matter phase. Thus within our
treatment we cannot judge whether neutron star properties may actually call
for a sexaquark or not. We also stress that our approach in this paper has
been entirely phenomenological; theoretical work is needed to decide whether
the parameters leading to the successful models with sexaquarks can emerge
from fundamental QCD, and to accurately treat the cross-over between the quark
and hadron phases. (The same caveats also apply to models without sexaquarks.)

## Acknowledgements

We thank Sanjay Reddy for valuable discussions. We also acknowledge
discussions with Kazem Azizi. This work was supported by the Polish National
Science Centre (NCN) under grant number 2019/33/B/ST9/03059. D.B. received
support from the Russian Fund for Basic Research (RFBR) under grant No.
18-02-40137 and from the Russian Federal Program ”Priority-2030”. The research
of G.R.F. is supported by NSF-PHY-2013199 and the Simons Foundation. D. A-C.
acknowledges support from the Bogoliubov-Infeld program for the collaboration
between JINR Dubna and Polish Institutions. This work is part of a project
that has received funding from the European Union’s Horizon 2020 research and
innovation program under the grant agreement STRONG – 2020 - No 824093. We are
grateful to the COST Actions CA15213 ”THOR” and CA16214 ”PHAROS” for their
support of our networking activities.

##  Appendix A Guide to terminology

For the convenience of the readers, we collect in this Appendix the
explanations of the numerous abbreviations that have been used in the main
text of this work.

DD2

Parameterization of a GDRF for hadronic matter including only nucleons

DD2Y-T

DD2 with hyperons by S. Typel

DD2Y-T+S

DD2Y-T with sexaquark

GRDF

Generalized relativistic density functional

MC

Maxwell Construction

RIC

Replacement Interpolation Construction

DSiY

Deconfinement (D) before sexaquark (Si) and hyperon (Y) onset, where i=1885,
1941, 2054 denotes to the mass of S in [MeV]

SiDY

Deconfinement after sexaquark but before hyperon onset

SiYD

Deconfinement after both sexaquark and hyperon onset when S is prior

YSiD

Deconfinement after both sexaquark and hyperon onset when hyperon is prior

CSSj, j=1,2,…,9

Constant speed of sound parameterization for quark matter EoS, where $j$
denotes to the number of the set in Table 3

MC-Si-j, j=1,2,…,9

Hybrid solution which has been obtained from Maxwell construction where $j$
denotes the number of sets in Table 3

CSS-A

Constant speed of sound parameterization for quark matter EoS where A denotes
the case with x = 0.03 in Table 4

CSS-B

Constant speed of sound parameterization for quark matter EoS where B denotes
the case with x = 0.04 in Table 4

CSS-C

Constant speed of sound parameterization for quark matter EoS where C denotes
the case with x = 0.05 in Table 4

RIC-A-k, $k=1,2,3$

Hybrid solution obtained with RIC when S1885 is in hadronic matter EoS. A is
defined the same as in CSS-A and $k$ corresponds to the different sets for
replacement interpolation construction in Table 4

RIC-B-k, $k=1,2$

Same as RIC-A-k, for A $\leftrightarrow$ B, corresponding to CSS-B

RIC-C-k, $k=1$

Same as RIC-A-k, for A $\leftrightarrow$ C, corresponding to CSS-C

##  Appendix B Hadronic and quark matter EoS in the P-$\mu$ plane

In this appendix we illustrate the construction of the hybrid EoS via MC and
RIC given the hadronic and quark matter EoS as input.

In Fig. 14, the MC is performed for different scenarios of hadronic and quark
matter EoS combinations for the case of $m_{S}=1941$ MeV. As was discussed in
subsection IV.1, the reconfinement phenomenon is ignored at higher chemical
potentials when a Maxwell construction is employed. The CSS$\\_$4 is almost
masqueraded with the hadronic matter EoS around the transition point. However,
a $\mu$-dependent bag pressure has been used for constructing CSS$\\_$3 and
CSS$\\_$4 to produce a bigger change in the slope at the transition point and
therefore a bigger jump at the density, but a bigger slope for the quark
matter EoS is not compatible with fulfilling the causality condition.

In Fig. 15, MC and RIC are performed for the case $m_{S}=1885$ MeV. Panel (d)
shows a MC when an early deconfinement happens while the other panels show RIC
with the possible minimum value of $\Delta_{p}$ for which the mechanical
stability is fulfilled. The described situation of the reconfinement and
crossover transition by RIC is depicted in this figure.

Figure 14: MC for different scenarios of hybrid stars. The solutions for the
S1941DY scenario when the deconfinement occurs after S1941 onset and before Y
onset has been shown in panel (a). Panel (b) shows the solution for the
YS1941D scenario when the deconfinement occurs after both S onset and Y onset
and panel (c) shows the solution for the DS1941Y scenario when the
deconfinement occurs before S1941 onset and Y onset. The curves labelled
CSS$\\_$1 to CSS$\\_$5 correspond to the quark matter EoS which after MC give
rise to the parameter sets MC-S1941-1 to MC-S1941-5 which are given in Table
3.

Figure 15: MC and RIC solutions when $S_{1885}$ is included in the hadronic matter EoS. The curves labelled CSS-A, CSS-B and CSS-C correspond to the quark matter EoS which after RIC give rise to the parameter sets RIC-A, RIC-B and RIC-C given in Table 3 while CSS-7, CSS-8 and CSS-9 correspond to the quark matter EoS which after MC give rise to the sets MC-S1885-7, MC-S1885-8 and MC-S1885-9 which are given in Table 4.  Table 3:  Parameter sets characterizing hybrid EoS which have been investigated as a solution for MC when S1941 or S1885 are included in hadronic matter.  set | $m_{S}$ | $~{}~{}x~{}~{}$ | $~{}~{}\eta_{V}~{}~{}$ | c${}_{s}^{2}$/c2 | A | B | B0 | B1 | $\mu_{<}$ | $\Gamma_{<}$ | scenario
---|---|---|---|---|---|---|---|---|---|---|---
| [MeV] |  |  |  | [MeVfm-3] | [MeVfm-3] | [MeVfm-3] | [MeVfm-3] | [MeV] | [MeV] |
MC-S1941-1 | 1941 | 0.05 | 0.11 | 0.45 | 93.49 | 81.73 | 9.0 | 0 | 0 | 0 | S1941DY
MC-S1941-2 | 1941 | 0.05 | 0.14 | 0.48 | 94.99 | 85.50 | 5.7 | 0 | 0 | 0 | S1941DY
MC-S1941-3 | 1941 | 0.05 | 0.17 | 0.52 | 96.85 | 90.18 | 0 | 7.0 | 1020 | 110 | S1941DY
MC-S1941-4 | 1941 | 0.08 | 0.14 | 0.48 | 94.99 | 85.50 | 0 | 10.0 | 1250 | 100 | YS1941D
MC-S1941-5 | 1941 | 0.05 | 0.17 | 0.52 | 96.85 | 90.18 | 0 | 0 | 0 | 0 | DS1941Y
MC-S1885-6 | 1885 | 0.030 | 0.14 | 0.48 | 94.99 | 85.50 | 0 | 0 | 0 | 0 | S1885DY
MC-S1885-7 | 1885 | 0.035 | 0.13 | 0.47 | 94.45 | 84.14 | 0 | 0 | 0 | 0 | DS1885Y
MC-S1885-8 | 1885 | 0.040 | 0.13 | 0.47 | 94.45 | 84.14 | 1.8 | 0 | 0 | 0 | DS1885Y
MC-S1885-9 | 1885 | 0.050 | 0.13 | 0.47 | 94.45 | 84.14 | 3.5 | 0 | 0 | 0 | DS1885Y
Table 4:  Parameter sets characterizing hybrid EoS which have been investigated as a solution for RIC when S1885 is included in hadronic matter. In different sets, A corresponds to $x=0.03$, B corresponds to $x=0.04$ and C corresponds to $x=0.05$.  set | $m_{S}$ | $~{}~{}x~{}~{}$ | $~{}~{}\eta_{V}~{}~{}$ | c${}_{s}^{2}$/c2 | A | B | B0 | $\Delta_{p}$ | $\mu_{H}$ | $\mu_{Q}$ | scenario
---|---|---|---|---|---|---|---|---|---|---|---
| [MeV] |  |  |  | [MeVfm-3] | [MeVfm-3] | [MeVfm-3] |  | [MeV | MeV] |
RIC-A-1 | 1885 | 0.03 | 0.14 | 0.48 | 94.99 | 85.50 | 0 | -3$\%$ | 1027.76 | 1121.92 | S1885DY
RIC-A-2 | 1885 | 0.03 | 0.14 | 0.48 | 94.99 | 85.50 | 0 | -5$\%$ | 1019.32 | 1261.74 | S1885YD
RIC-A-3 | 1885 | 0.03 | 0.14 | 0.48 | 94.99 | 85.50 | 0 | -7$\%$ | 1012.81 | 1418.92 | S1885YD
RIC-B-1 | 1885 | 0.04 | 0.14 | 0.48 | 94.99 | 85.50 | -4 | -7$\%$ | 1064.10 | 1335.80 | S1885YD
RIC-C-1 | 1885 | 0.05 | 0.12 | 0.46 | 93.95 | 82.88 | 0 | -3$\%$ | 1122.39 | 1286.74 | S1885YD
RIC-C-2 | 1885 | 0.05 | 0.12 | 0.46 | 93.95 | 82.88 | 0 | -5$\%$ | 1097.33 | 1398.25 | S1885YD

## References

  * Farrar (2003) G. R. Farrar, Int. J. Theor. Phys. 42, 1211 (2003).
  * Farrar (2018) G. R. Farrar, (2018), arXiv:arXiv:1805.03723 [hep-ph] .
  * Farrar (2022) G. R. Farrar, (2022), arXiv:2201.01334 [hep-ph] .
  * Jaffe (1977) R. L. Jaffe, Phys. Rev. Lett. 38, 195 (1977), [Erratum: Phys.Rev.Lett. 38, 617 (1977)].
  * Aaij  _et al._ (2020) R. Aaij _et al._ (LHCb), Sci. Bull. 65, 1983 (2020), arXiv:2006.16957 [hep-ex] .
  * Ali  _et al._ (2019) A. Ali, L. Maiani, and A. D. Polosa,  _Multiquark Hadrons_ (Cambridge University Press, 2019).
  * Green  _et al._ (2021) J. R. Green, A. D. Hanlon, P. M. Junnarkar, and H. Wittig, Phys. Rev. Lett. 127, 242003 (2021), arXiv:2103.01054 [hep-lat] .
  * Francis  _et al._ (2019) A. Francis, J. R. Green, P. M. Junnarkar, C. Miao, T. D. Rae, and H. Wittig, Phys. Rev. D 99, 074505 (2019).
  * Beane  _et al._ (2013) S. R. Beane, E. Chang, S. D. Cohen, W. Detmold, H. W. Lin, T. C. Luu, K. Orginos, A. Parreno, M. J. Savage, and A. Walker-Loud (NPLQCD), Phys. Rev. D 87, 034506 (2013).
  * Inoue  _et al._ (2011) T. Inoue, N. Ishii, S. Aoki, T. Doi, T. Hatsuda, Y. Ikeda, K. Murano, H. Nemura, and K. Sasaki (HAL QCD), Phys. Rev. Lett. 106, 162002 (2011).
  * Beane  _et al._ (2011) S. R. Beane _et al._ (NPLQCD), Phys. Rev. Lett. 106, 162001 (2011).
  * Detmold (2021) W. Detmold, (2021), private communication.
  * Buccella (2020) F. Buccella, PoS CORFU2019, 024 (2020).
  * Azizi  _et al._ (2020) K. Azizi, S. S. Agaev, and H. Sundu, J. Phys. G 47, 095001 (2020), arXiv:1904.09913 [hep-ph] .
  * Farrar and Wang (2022) G. R. Farrar and Z. Wang, (2022), in preparation.
  * Farrar and Zaharijas (2004) G. R. Farrar and G. Zaharijas, Phys. Rev. D 70, 014008 (2004).
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016).
  * Tanabashi  _et al._ (2018) M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Bazavov  _et al._ (2019) A. Bazavov _et al._ (HotQCD), Phys. Lett. B 795, 15 (2019), arXiv:1812.08235 [hep-lat] .
  * Farrar  _et al._ (2020) G. R. Farrar, Z. Wang, and X. Xu, (2020), arXiv:2007.10378 [hep-ph] .
  * Blaschke  _et al._ (2021) D. Blaschke _et al._ , Int. J. Mod. Phys. A 36, 2141005 (2021).
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021).
  * Antoniadis  _et al._ (2013) J. Antoniadis _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Shahrbaf  _et al._ (2020a) M. Shahrbaf, D. Blaschke, A. G. Grunfeld, and H. R. Moshfegh, Phys. Rev. C 101, 025807 (2020a).
  * Shahrbaf  _et al._ (2020b) M. Shahrbaf, D. Blaschke, and S. Khanmohamadi, J. Phys. G 47, 115201 (2020b).
  * McDermott  _et al._ (2019) S. D. McDermott, S. Reddy, and S. Sen, Phys. Rev. D 99, 035013 (2019).
  * Kolb and Turner (2019) E. W. Kolb and M. S. Turner, Phys. Rev. D 99, 063519 (2019), arXiv:1809.06003 [hep-ph] .
  * Walecka (1974) J. D. Walecka, Annals Phys. 83, 491 (1974).
  * Serot and Walecka (1986) B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16, 1 (1986).
  * Reinhard (1989) P. G. Reinhard, Rept. Prog. Phys. 52, 439 (1989).
  * Ring (1996) P. Ring, Prog. Part. Nucl. Phys. 37, 193 (1996).
  * Rusnak and Furnstahl (1997) J. J. Rusnak and R. J. Furnstahl, Nucl. Phys. A 627, 495 (1997).
  * Dreizler and Gross (2012) R. M. Dreizler and E. K. Gross, “Density functional theory: an approach to the quantum many-body problem,” (Springer Science & Business Media, 2012).
  * Typel and Wolter (1999) S. Typel and H. H. Wolter, Nucl. Phys. A 656, 331 (1999).
  * Typel (2005) S. Typel, Phys. Rev. C 71, 064301 (2005).
  * Klahn  _et al._ (2006) T. Klahn _et al._ , Phys. Rev. C 74, 035802 (2006).
  * Sun  _et al._ (2019) B. Y. Sun, Z. W. Liu, and R. Y. Xing, AIP Conf. Proc. 2127, 020020 (2019), arXiv:1904.10372 [nucl-th] .
  * Typel (2018) S. Typel, Particles 1, 3 (2018).
  * Dutra  _et al._ (2014) M. Dutra, O. Lourenço, S. S. Avancini, B. V. Carlson, A. Delfino, D. P. Menezes, C. Providência, S. Typel, and J. R. Stone, Phys. Rev. C 90, 055203 (2014), arXiv:1405.3633 [nucl-th] .
  * Typel  _et al._ (2010) S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010).
  * Oertel  _et al._ (2017) M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017).
  * Pais and Typel (2017) H. Pais and S. Typel, “Comparison of equation of state models with different cluster dissolution mechanisms,” in  _Nuclear Particle Correlations and Cluster Physics_ (World Scientific, 2017) Chap. 4, pp. 95–132, arXiv:arXiv:1612.07022 [nucl-th] .
  * Shahrbaf and Moshfegh (2019) M. Shahrbaf and H. R. Moshfegh, Annals Phys. 402, 66 (2019).
  * Weissenborn  _et al._ (2012) S. Weissenborn, D. Chatterjee, and J. Schaffner-Bielich, Phys. Rev. C 85, 065802 (2012), [Erratum: Phys.Rev.C 90, 019904 (2014)].
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Raduta  _et al._ (2020) A. R. Raduta, M. Oertel, and A. Sedrakian, Mon. Not. Roy. Astron. Soc. 499, 914 (2020).
  * Alford  _et al._ (1999) M. G. Alford, K. Rajagopal, and F. Wilczek, Nucl. Phys. B 537, 443 (1999).
  * Buballa (2005) M. Buballa, Phys. Rept. 407, 205 (2005), arXiv:hep-ph/0402234 .
  * Ruester  _et al._ (2005) S. B. Ruester, V. Werth, M. Buballa, I. A. Shovkovy, and D. H. Rischke, Phys. Rev. D 72, 034004 (2005), arXiv:hep-ph/0503184 .
  * Blaschke  _et al._ (2005) D. Blaschke, S. Fredriksson, H. Grigorian, A. M. Oztas, and F. Sandin, Phys. Rev. D 72, 065020 (2005), arXiv:hep-ph/0503194 .
  * Abuki and Kunihiro (2006) H. Abuki and T. Kunihiro, Nucl. Phys. A 768, 118 (2006), arXiv:hep-ph/0509172 .
  * Antić  _et al._ (2021) S. Antić, M. Shahrbaf, D. Blaschke, and A. G. Grunfeld, (2021), arXiv:arXiv:2105.00029 [nucl-th] .
  * Zdunik and Haensel (2013) J. L. Zdunik and P. Haensel, Astron. Astrophys. 551, A61 (2013).
  * Schmidt  _et al._ (1994) S. M. Schmidt, D. Blaschke, and Y. L. Kalinovsky, Phys. Rev. C 50, 435 (1994).
  * Gomez Dumm  _et al._ (2006) D. Gomez Dumm, D. B. Blaschke, A. G. Grunfeld, and N. N. Scoccola, Phys. Rev. D 73, 114019 (2006).
  * Blaschke  _et al._ (2007) D. B. Blaschke, D. Gomez Dumm, A. G. Grunfeld, T. Klahn, and N. N. Scoccola, Phys. Rev. C 75, 065804 (2007).
  * Contrera  _et al._ (2022) G. A. Contrera, D. Blaschke, J. P. Carlomagno, A. G. Grunfeld, and S. Liebing, (2022), arXiv:2201.00477 [nucl-th] .
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018a), arXiv:1805.11581 [gr-qc] .
  * Alford  _et al._ (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Gerlach (1968) U. H. Gerlach, Phys. Rev. 172, 1325 (1968).
  * Glendenning and Kettner (2000) N. K. Glendenning and C. Kettner, Astron. Astrophys. 353, L9 (2000), arXiv:astro-ph/9807155 .
  * Alvarez-Castillo  _et al._ (2019) D. E. Alvarez-Castillo, D. B. Blaschke, A. G. Grunfeld, and V. P. Pagura, Phys. Rev. D 99, 063010 (2019).
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Baym and Chin (1976) G. Baym and S. A. Chin, Phys. Lett. B 62, 241 (1976).
  * Masuda  _et al._ (2013a) K. Masuda, T. Hatsuda, and T. Takatsuka, PTEP 2013, 073D01 (2013a), arXiv:1212.6803 [nucl-th] .
  * Masuda  _et al._ (2013b) K. Masuda, T. Hatsuda, and T. Takatsuka, Astrophys. J. 764, 12 (2013b), arXiv:1205.3621 [nucl-th] .
  * Ayriyan  _et al._ (2018) A. Ayriyan, N. U. Bastian, D. Blaschke, H. Grigorian, K. Maslov, and D. N. Voskresensky, Phys. Rev. C 97, 045802 (2018), arXiv:1711.03926 [nucl-th] .
  * Abgaryan  _et al._ (2018) V. Abgaryan, D. Alvarez-Castillo, A. Ayriyan, D. Blaschke, and H. Grigorian, Universe 4, 94 (2018), arXiv:1807.08034 [astro-ph.HE] .
  * Ayriyan  _et al._ (2021) A. Ayriyan, D. Blaschke, A. G. Grunfeld, D. Alvarez-Castillo, H. Grigorian, and V. Abgaryan, Eur. Phys. J. A 57, 318 (2021), arXiv:2102.13485 [astro-ph.HE] .
  * Marques  _et al._ (2017) M. Marques, M. Oertel, M. Hempel, and J. Novak, Phys. Rev. C 96, 045806 (2017).
  * Stone  _et al._ (2021) J. R. Stone, V. Dexheimer, P. A. M. Guichon, A. W. Thomas, and S. Typel, Mon. Not. Roy. Astron. Soc. 502, 3476 (2021).
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021).
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021).
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019).
  * Abbott  _et al._ (2018b) B. P. Abbott, R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, _et al._ , Physical review letters 121, 161101 (2018b).
  * Blaschke  _et al._ (2022) D. Blaschke, O. Ivanytskyi, and M. Shahrbaf (2022) arXiv:2202.05061 [nucl-th] .
  * Hebeler  _et al._ (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013).
  * Alvarez-Castillo  _et al._ (2016) D. Alvarez-Castillo, A. Ayriyan, S. Benic, D. Blaschke, H. Grigorian, and S. Typel, Eur. Phys. J. A 52, 69 (2016), arXiv:1603.03457 [nucl-th] .
  * Thapa  _et al._ (2021) V. B. Thapa, M. Sinha, J. J. Li, and A. Sedrakian, Phys. Rev. D 103, 063004 (2021), arXiv:2102.08787 [astro-ph.HE] .
  * Grigorian  _et al._ (2018) H. Grigorian, D. N. Voskresensky, and K. A. Maslov, Nucl. Phys. A 980, 105 (2018), arXiv:1808.01819 [astro-ph.HE] .