# The effects of self-interacting bosonic dark matter on neutron star
properties

Edoardo Giangrandi CFisUC, Department of Physics, University of Coimbra, Rua
Larga P-3004-516, Coimbra, Portugal  Violetta Sagun CFisUC, Department of
Physics, University of Coimbra, Rua Larga P-3004-516, Coimbra, Portugal
Oleksii Ivanytskyi Incubator of Scientific Excellence—Centre for Simulations
of Superdense Fluids, University of Wrocław, 50-204, Wroclaw, Poland
Constança Providência CFisUC, Department of Physics, University of Coimbra,
Rua Larga P-3004-516, Coimbra, Portugal  Tim Dietrich Institut für Physik und
Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, Potsdam,
Germany  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Am Mühlenberg 1, Potsdam D-14476, Germany

###### Abstract

We propose a model of asymmetric bosonic dark matter (DM) with self-repulsion.
By adopting the two-fluid formalism, we study different DM distribution
regimes, either, fully condensed inside the core of a star, or, otherwise,
distributed in a dilute halo around a neutron star (NS). We show that for a
given total gravitational mass, DM condensed in a core leads to a smaller
radius and tidal deformability compared to a pure baryonic star. This effect
may be interpreted as an effective softening of the equation of state. On the
other hand, the presence of a DM halo increases the tidal deformability and
total gravitational mass. As a result, an accumulated DM inside compact stars
could mimic an apparent softening/stiffening of strongly interacting matter
EoS and constraints we impose on it at high densities.

We limit the model parameter space by confronting the cross-section of the DM
self-interaction to the constraint extracted from the analysis of the Bullet
Cluster. Furthermore, from the performed analysis of the effect of DM
particles, interaction strength, and relative DM fractions inside NSs we
obtained a rigorous constraint on model parameters. To identify its impact on
NSs we consider the DM fraction may reach up to 5%, which could be considered
too high in several scenarios. Finally, we discuss several pieces of smoking
gun evidence of the presence of DM that is free from the abovementioned
degeneracy between the effect of DM and properties of strongly interacting
matter. These signals could be probed with future and ongoing astrophysical
and gravitational wave (GW) surveys.

Neutron Stars(1108) — Dark Matter(353) — Gravitational Waves(678)

##  1 Introduction

Since the first detection of the binary neutron star (NS) merger, GW170817
(Abbott et al., 2017a), which was accompanied by the observation of
electromagnetic signals originating from the same source, GRB170817A and
AT2017gfo (Abbott et al., 2017b), we have been witnessing exciting
breakthroughs in our understanding of compact stars and their merger dynamics.
In fact, gravitational wave (GW) astronomy and multi-messenger astrophysics
became new tools to extract information about the internal structure of NSs
from GW and electromagnetic observations (Bauswein et al., 2017; Annala et
al., 2018; Hinderer et al., 2019). Thus, from the combination of the analyses
of the GW170817 signal measured by the advanced LIGO and advanced Virgo
detectors, the constraint on the tidal deformability parameter of NS matter
$\Lambda_{1.4}\leq 800$ was extracted (Abbott et al., 2018). The second binary
NS merger event, GW190425 (Abbott et al., 2020a), provided constraints
consistent with GW170817, but due to its lower signal-to-noise ratio did not
deepen our knowledge about the NS Equation of State (EoS). In addition to GW
observations, also X-ray observations by NICER (Miller et al., 2019, 2021;
Riley et al., 2019, 2021; Raaijmakers et al., 2020), radio measurements of the
heaviest pulsars, e.g. PSR J0348+0432 of mass $2.01\pm 0.04$ $M_{\odot}$
(Antoniadis et al., 2013), PSR J0740+6620 of $2.08^{+0.07}_{-0.07}$
$M_{\odot}$ (Fonseca et al., 2021) and optical observations of the black widow
pulsars, e.g., PSR J1810+1744 of $2.13\pm 0.04$ $M_{\odot}$ (Romani et al.,
2021), and PSR J0952-0607 of $2.35\pm 0.17$ $M_{\odot}$ (Romani et al., 2022)
constrain properties of NSs.

While all the mentioned analyses and models assume that NSs are embedded in a
pure vacuum and do not contain dark matter (DM), they, indeed, could
accumulate a sizable amount of DM in their interior and surroundings. Due to
high compactness, NSs can effectively trap DM particles, which will rapidly
thermalize and become accrued inside the stars, altering their properties. The
presence of DM affects the internal structure and compactness of compact
stars. Thus, as it was shown, e.g., Ciarcelluti & Sandin (2011); Ellis et al.
(2018b); Nelson et al. (2019); Ivanytskyi et al. (2020); Das et al. (2020) and
Sagun et al. (2022), DM may either form an extended halo or a dense core
inside a NS. Depending on the mass of DM particles, their self-interaction
strength, and its relative abundance inside the star, one of the
abovementioned scenarios takes place. Since DM halos are invisible for typical
astrophysical observations, we would see only the baryonic matter (BM) radius,
independent of the fact that the outermost radius can extend further than its
BM component (Karkevandi et al., 2022). On the contrary, a DM core formation
will lead to a reduction in the NS radius.

Moreover, DM will affect tidal deformability parameters and the merger
dynamics (Ellis et al., 2018a; Bezares et al., 2019; Bauswein et al., 2023;
Leung et al., 2022). Nowadays, while there are studies investigating possible
alternative scenarios beyond standard compact binary mergers described by
general relativity in pure vacuum (Abbott et al., 2019a), the models used to
analyze GW signal do not account directly for DM.

Thus, to understand the effect of DM on the coalescence of NSs, numerical-
relativity simulations for different DM fractions, particle mass, and
interaction strength are required. As a step in this direction, there have
been the first two-fluid 3D simulations of coalescing binary NS admixed with
DM with the following studies of GW emission of the merger remnant, e.g.
Bauswein et al. (2023) and Emma et al. (2022). By considering different binary
masses and EoSs, Bauswein et al. (2023) showed that the GW frequency of the
orbiting DM components scales with the compactness of NSs. Moreover, the
relations between the DM GW frequency and the dominant post-merger GW
frequency of the stellar fluid or the tidal deformability were found, which
opens a possibility to probe the EoS effects during the binary inspiral. Emma
et al. (2022) studied the effect of mirror DM concentrated inside the core on
the deceleration of the inspiral phase, as well as on a modification of the
ejecta and debris disk formation.

Depending on whether DM has particle-antiparticle asymmetry, we will refer to
it as asymmetric or symmetric matter. Symmetric DM particles can self-
annihilate leaving a possibility of its detection via X-ray, $\gamma$-ray, or
neutrino telescopes (Kouvaris, 2008). Moreover, as studied in Pérez-García &
Silk (2012), self-annihilating DM in the inner regions of NSs may have a
significant impact on the kinematic properties, namely, velocity kicks and
rotation patterns.

Another possible effect of DM particle annihilation inside the NS core is
related to the late-time heating, which could be detected from observations of
the surface temperature of the old part of the NS population (de Lavallaz &
Fairbairn, 2010; Hamaguchi et al., 2019). Unfortunately, nowadays, our
database of old NSs is still quite limited.

Contrary to the annihilating DM, asymmetric DM will become accumulated inside
a star. Models that consider such scenario should allow old NSs to exist.
Especially, it is important for bosonic DM particles, which at zero
temperature could form Bose-Einstein condensate (BEC) leading to the
gravitational collapse of the bosonic DM to a black hole (Kouvaris, 2013a).

Light DM particles, such as axions, could contribute as an additional cooling
channel in compact stars. Thus, in the NS core axions can be produced either
in nucleon bremsstrahlung or in Cooper pair breaking and formation processes
(Sedrakian, 2016, 2019; Buschmann et al., 2022), causing an alteration of the
surface temperature and thermal evolution of a star. In addition, most of the
existing models are constrained by the results of neutrino emission coming
from the supernova observation SN 1987A (Chang et al., 2018) and existing NS
cooling data. The results of NS merger simulations (Dietrich & Clough, 2019)
show that axions produced in nucleon-nucleon bremsstrahlung do not lead to a
measurable change in the emitted GW signal, ejecta mass, as well as the
temperature profile of the merger remnant.

The fraction of DM in the compact star interior depends on different factors.
Thus, DM can be captured by NSs from the surrounding medium. Following
scattering processes, the kinetic energy of DM particles is transferred to the
star (Kouvaris, 2013b; Bell et al., 2020, 2021), and becomes gravitationally
bound with a star. The amount of DM accrued by an ordinary accretion
throughout a stellar evolution will depend on the position of the considered
NS in the Galaxy (Kouvaris & Tinyakov, 2010). As the DM density in the
Galactic Center is many orders of magnitude greater than in its arms, we may
expect a higher DM fraction in compact stars toward the Milky Way center (Del
Popolo et al., 2020). Furthermore, NSs in globular clusters may contain a
significant amount of DM (Bertone & Fairbairn, 2008). Moreover, we should not
forget that NS is the final stage of star’s evolution preceded by the
progenitor, main-sequence star, and supernova explosion with the formation of
a proto-NS. These and all other mechanisms are discussed in detail in
Karkevandi et al. (2022). As was estimated by Ivanytskyi et al. (2020), the
amount of accumulated DM in the most central part of the Galaxy accrued by a
spherically symmetric accretion scenario during the main-sequence star and
equilibrated NS stages is 0.01%. However, additional scenarios could lead to
high DM factions inside compact stars, e.g., DM production during a supernova
explosion, accretion of DM clumps formed at the early stage of the universe,
or initial star formation on a preexisting dark core. As we are interested in
identifying possible signatures on NSs, we have concentrated our study on
fractions up to 5%, which cannot be accumulated by a spherically symmetric DM
accretion followed by thermalization via the interaction with BM, and requires
the alternative mechanisms mentioned above.

Moreover, some local non-homogeneity of DM distribution may contribute to an
increase in DM fraction, leading even to dark compact objects (Dengler et al.,
2022) and dark stars (Kouvaris & Nielsen, 2015; Maselli et al., 2017).

Since DM properties are still unknown, different models have been employed,
considering its fermionic (Goldman et al., 2013; Gresham et al., 2017;
Ivanytskyi et al., 2020) and bosonic (Colpi et al., 1986; Petraki & Volkas,
2013; Karkevandi et al., 2022) nature. As it was discussed by Bramante et al.
(2013) to be consistent with the observations of old NSs, bosonic DM has to be
either self-interacting, decaying, or self-annihilating. Considering
asymmetric DM a repulsive self-interaction is required due to zero degeneracy
pressure. At the moment when accumulated bosonic asymmetric DM exceeds the
Chandrasekhar mass, nothing can prevent its gravitational collapse and the
formation of a black hole inside the NS, which could potentially disrupt the
star (Kouvaris, 2013a; Zurek, 2014).

Using an analog of visible matter and the Standard Model particles, we see
that all interactions have an exchange character, an interaction between
particles occurs due to an exchange of a mediator, e.g., the interaction
between nucleons is mediated by pions. In the present article, we extend this
approach for a dark sector by formulating a model of self-interacting
asymmetric bosonic DM, which includes vector interaction mediated by a real
$\omega$-field coupled to the scalar one. We model DM-admixed compact stars by
considering the mixed system of two fluids with different relative fractions.

The assumption of cold self-interacting DM provides a good agreement with the
large-scale structure of the universe and cosmology. It reconciles the success
of the cold DM (CDM) model with the non-observation of cuspy density profiles
of dwarf galaxies predicted by CDM N-body simulations and known as the core-
cusp problem (Moore, 1994). In comparison to an alternative mechanism to
flatten the central density profile by supernova-driven episodes of gas
removal, self-interacting DM is a more favorable one (Burger et al., 2022). On
the other hand, from the observed mass profiles of galaxies (Ahn & Shapiro,
2005) and Bullet Cluster observations (Clowe et al., 2006; Randall et al.,
2008) the DM self-interaction cross section per unit mass has an upper limit
of $\sigma/m<$1.25 cm2 g-1 (68% confidence level). The DM model considered in
this paper is in line with the above assumption, and therefore, provides
consistency with the state of the art of modern cosmology. Moreover, below we
explicitly account for the above constraint on $\sigma/m$ in order to limit
the model parameters. An implication of the proposed EoS and tests against
astrophysical and GW observations are performed in this work.

The paper is organized as follows. In Section 2, we present models for the BM
and DM components, with a detailed derivation provided in the Appendix A.
Section 3 is dedicated to the equilibrium configurations of DM-admixed compact
stars. In Section 4, we discuss how the speed of sound and the tidal
deformability are affected by the presence of DM. In Section 5 the main
results are presented, including the constraints on mass and interaction scale
of DM particles. In Section 6, we discuss the smoking gun signals of the
presence of DM that could be tested in the nearest future before concluding in
Section 7. In Appendix A, we show the full derivation of the DM EoS, with a
focus on the effective speed of sound for a DM-admixed NS in Appendix B. In
Appendix C, we show the scan over the model parameters and the obtained
constraints. Throughout the article, we utilize a unit system in which
$\hbar=c=G=1$.

##  2 Models of dark and baryonic matter

###  2.1 DM EoS

We consider the model of massive spinless DM particles carrying a conserved
charge. Such particles are described by a complex scalar field, have mass
$m_{\chi}$ and chemical potential $\mu_{\chi}$. At sufficiently low
temperatures bosonic DM exists in the form of the BEC. In the absence of
interaction such BEC has zero pressure and is mechanically unstable against
gravitational compression. We stabilize the BEC of DM by introducing repulsive
interaction mediated by a real vector field coupled to the scalar one. The
minimal Lagrangian representing this model is given in the Appendix A. It is
equivalent to a massive U(1) gauge theory of scalar particles, i.e., scalar
electrodynamics with massive photons. This Lagrangian implies a Noether
current corresponding to the invariance of action with respect to global U(1)
transformations. If the vector field was not a Yukawa but a gauge one, local
U(1) symmetry would also be respected and another Noether current could be
introduced (Brading & Brown, 2000). Given a quantum treatment, expectation
values of these two currents produce the same conserved charge, which is not
the case within the used mean-field approximation corresponding to a classical
treatment of the vector field (see Appendix A for details). We use the Noether
current resulting from global U(1) transformations, which leave the action
invariant even at the mean-field level.

In this work, we assume the vanishing temperature of the DM, being totally
converted to the BEC. In the considered case thermal fluctuations are
suppressed and mean-field approximation can be applied in order to derive the
corresponding EoS. Chemical potentials of the BM and DM components of NS scale
proportionally (for more details see Section 3). This significantly simplifies
solving two coupled Tolman–Oppen- heimer–Volkoff (TOV)–like equations for BM
and DM components, as shown by Ivanytskyi et al. (2020). Therefore, it is
convenient to formulate the DM EoS in the grand canonical ensemble (GCE),
where $\mu_{\chi}$ is an independent variable. Appendix A includes details of
the corresponding derivation for the interval of physical values of
$\mu_{\chi}\in[0,\sqrt{2}m_{\chi}]$ performed in the locally flat spacetime,
provided by small gradients of metrics and absence of the anisotropy issues
(see Karkevandi et al. (2022) for details). The corresponding pressure and
energy density are

$\displaystyle p_{\chi}$ $\displaystyle=$
$\displaystyle\frac{m_{I}^{2}}{4}\left(m_{\chi}^{2}-\mu_{\chi}\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}\right),$
(1) $\displaystyle\varepsilon_{\chi}$ $\displaystyle=$
$\displaystyle\frac{m_{I}^{2}}{4}\left(\frac{\mu_{\chi}^{3}}{\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}}-m_{\chi}^{2}\right),$
(2)

for $\mu_{\chi}\in[m_{\chi},\sqrt{2}m_{\chi}]$ and
$p_{\chi}=\varepsilon_{\chi}=0$ for $\mu_{\chi}\in[0,m_{\chi}]$. The parameter
$m_{I}$ has the unit of mass and controls the interaction strength. It is
proportional to the vector meson mass $m_{\omega}$ and inversely proportional
to its coupling $g$. Thus, large $m_{I}$ corresponds to weak interaction and
vice versa. At a first glance, the present EoS in the weak coupling regime
paradoxically leads to an infinite pressure due to $m_{I}\rightarrow\infty$.
This, however, is not the case since in the considered regime chemical
potential of the DM BEC $\mu_{\chi}$ coincides with its mass $m_{\chi}$
leading to the vanishing of the brackets in Eqs. (1) and (2). In the case of
$p_{\chi}$ the bracket vanishes faster than $m_{I}^{2}$ yielding to a zero
pressure $\sim m_{I}^{-2}$, while for $\varepsilon_{\chi}$ the bracket behaves
as $\sim m_{I}^{-2}$ providing a finite energy density of the DM BEC
$m_{\chi}n_{\chi}$. In the strong coupling regime, $m_{I}\rightarrow 0$
chemical potential of DM converges to $\sqrt{2}m_{\chi}$. As a result, the
bracket in Eq. (1) becomes equal to $m_{\chi}^{2}$ and the pressure vanishes
as $m_{I}^{2}m_{\chi}^{2}/4$. The corresponding bracket in Eq. (2) diverges as
$\sim m_{I}^{-2}$ leading to finite energy density $\sqrt{2}m_{\chi}n_{\chi}$.
Remarkably, weak and strong coupling limits of the present EoS are similar,
since DM pressure vanishes in both these cases. At $m_{I}\rightarrow\infty$ it
is due to the absence of repulsion. The limit $m_{I}\rightarrow 0$ is
equivalent to the case of the massless vector field, which does not have a
nontrivial mean-field solution needed to stiffen the EoS. Detailed analysis of
the weak and strong coupling limits of the present EoS is performed in
Appendix A.

Figure 1: Left panel: Scaled pressure $p_{\chi}/p_{\infty}$ (black solid
curve), energy density $\varepsilon_{\chi}/p_{\infty}$ (black dashed curve)
and speed of sound squared $c_{s,\chi}^{2}$ (red dotted curve) of DM as
functions of its chemical potential $\mu_{\chi}$ given in units of $m_{\chi}$.
Right panel: Scaled pressure $p_{\chi}/p_{\infty}$ (black solid curve) and
speed of sound $c_{s,\chi}^{2}$ (red dotted curve) of DM as functions of
scaled energy density $\varepsilon_{\chi}$ given in units of $p_{\infty}$.

A remarkable feature of the present EoS is that at infinite density its
pressure is limited by the value $p_{\infty}=m_{I}^{2}m_{\chi}^{2}/4$. This
regime is reached at $\mu_{\chi}=\sqrt{2}m_{\chi}$. Thus, the compressibility
of DM vanishes at asymptotically high densities regardless of $m_{\chi}$ and
$m_{I}$. The same conclusion holds for the speed of sound
$c_{s,\chi}^{2}=dp_{\chi}/d\varepsilon_{\chi}$. In other words, high-density
configurations of bosonic DM are gravitationally unstable at any strength of
the repulsive interaction. The left panel of Fig. 1 shows the pressure, energy
density, and speed of sound of the considered DM EoS as functions of the
corresponding chemical potential. It is worth mentioning, that the square of
the speed of sound is limited from above by the value $1/9$, which is reached
at $\mu_{\chi}=\sqrt{3/2}~{}m_{\chi}$ and does not depend on $m_{\chi}$ and
$m_{I}$. Thus, $c_{s,\chi}^{2}$ is bounded by quite small values and
corresponds to the soft EoS of DM. The right panel of Fig. 1 shows this EoS as
a function of energy density.

Observational data on the colliding clusters of galaxies 1E 0657-56 (the
Bullet Cluster) enable probing dynamics of the DM fluid on the cosmological
scale. This dynamic is determined by the DM self-interaction, which is
controlled by the corresponding cross section $\sigma_{\chi}$. Evaluation of
this quantity requires the invariant matrix element of scattering of two on-
shell DM particles from the initial state with three momenta ${\bf k}_{1}$ and
${\bf k}_{2}$ to the final one with momenta ${\bf k}_{1}^{\prime}$ and ${\bf
k}_{2}^{\prime}$. At the tree level, this matrix element includes
contributions from the $t$\- and $u$-channels

$\displaystyle i\mathcal{M}$ $\displaystyle=$
$\displaystyle\begin{array}[]{c}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a); \vertex[left =
of a] (b) {${\bf k}_{1}$}; \vertex[right = of a] (c) {${\bf k}_{1}^{\prime}$};
\vertex[below = of a] (d); \vertex[left = of d] (e) {${\bf k}_{2}$};
\vertex[right = of d] (f) {${\bf k}_{2}^{\prime}$}; \diagram{(b) -- [fermion]
(a) -- [fermion] (c)}; \diagram{(a) -- [boson, edge label = $q$] (d)};
\diagram{(e) -- [fermion] (d) -- [fermion] (f)};
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{array}$ (5)
$\displaystyle+$ $\displaystyle\begin{array}[]{c}\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(a); \vertex[left =
of a] (b) {${\bf k}_{1}$}; \vertex[below right = of a] (c) {${\bf
k}_{2}^{\prime}$}; \vertex[below = of a] (d); \vertex[left = of d] (e) {${\bf
k}_{2}$}; \vertex[above right = of d] (f) {${\bf k}_{1}^{\prime}$};
\diagram{(b) -- [fermion] (a) -- [fermion] (c)}; \diagram{(a) -- [boson, edge
label = $q$] (d)}; \diagram{(e) -- [fermion] (d) -- [fermion] (f)};
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\end{array}\,.$

These diagrams are formalized by the Feynman rules of scalar quantum
electrodynamics with massive photons. Each vertex corresponds to a factor
$-ig$. The annihilation of a particle with four-momentum $k$ and the creation
of another one with four-momentum $k^{\prime}$ in a vertex produces the factor
$-k_{\mu}-k_{\mu}^{\prime}$, while in the case of antiparticles, the momenta
entering this factor have the opposite signs. Wavy lines stand for the vector
field propagator $-ig^{\mu\nu}(q^{2}-m_{\omega}^{2})^{-1}$ with $q$ being a
transferred momentum. The conservation of energy and momentum ensures that in
the t\- and u-channels, which are represented by the upper and lower graphs in
Eq. (5), squared transferred momentum coincides with the Mandelstam variables
$t=(k_{1}-k_{1}^{\prime})^{2}=(k_{2}-k_{2}^{\prime})^{2}$ and
$u=(k_{1}-k_{2}^{\prime})^{2}=(k_{2}-k_{1}^{\prime})^{2}$, respectively. The
Lorentz indexes $\mu$ and $\nu$ appearing in the vertexes and in the vector
field propagator are dummies. With this, we arrive at

$\displaystyle\mathcal{M}=g^{2}\left[\frac{(k_{1}+k_{1}^{\prime})(k_{2}+k_{2}^{\prime})}{t-m_{\omega}^{2}}+\frac{(k_{1}+k_{2}^{\prime})(k_{2}+k_{1}^{\prime})}{u-m_{\omega}^{2}}\right].\quad$
(6)

In the BEC case, three momenta of the incoming and outgoing particles vanish,
which produces $4m_{\chi}^{2}$ in the numerator of each fraction in the
previous expression and corresponds to $t=u=0$. This yields
$\mathcal{M}=-8m_{\chi}^{2}/m_{I}^{2}$ and differential cross section of the
DM self-interaction $d\sigma_{\chi}/d\Omega=|\mathcal{M}|^{2}/64\pi^{2}E_{\rm
cm}$, where $E_{\rm cm}=2m_{\chi}$ is the center-of-mass energy of the
incoming particles with ${\bf k}_{1}={\bf k}_{2}=0$. This differential cross
section is independent of the angle variables. Therefore, the total one is
obtained by multiplying $d\sigma_{\chi}/d\Omega$ by $4\pi$. Finally, we obtain

$\displaystyle\frac{\sigma_{\chi}}{m_{\chi}}=\frac{2m_{\chi}}{\pi m_{I}^{4}}.$
(7)

Numerical simulations of the Bullet Cluster combined with the results from
X-ray, strong and weak lensing, optical observation set an upper limit on this
ratio Randall et al. (2008). Within the 68 % confidence interval
$\sigma_{\chi}/m_{\chi}<$1.25 cm2 g-1, while assuming equal mass-to-light
ratios in the subcluster and the main cluster prior to the merger yields even
more stringent constraint $\sigma_{\chi}/m_{\chi}<$0.7 cm2 g-1. In order to
keep the parameter space of the model as wide as possible, we use a more
relaxed version of this constraint. Thus, we require

$\displaystyle m_{I}~{}[{\rm MeV}]>18.24~{}\sqrt[4]{m_{\chi}~{}[{\rm MeV}]}.$
(8)

This requirement obviously discredits the strong coupling limit
$m_{I}\rightarrow 0$ with respect to the above cosmological constraint on the
DM self-interaction cross section. Further analysis is limited to the region
of the model parameter space, which respects Eq. (8).

###  2.2 Baryon matter EoS

In order to thoroughly study the impact of DM on compact stars made of mostly
BM, we consider two EoSs of different stiffness. One of them is the induced
surface tension (IST) EoS, formulated on the basis of the hard-core approach.
Thus, nucleons are characterized by an effective hard-core radius that
provides a short-range repulsion between the particles of different species.
This part of the model was fixed from the fit of heavy-ion collision data
(Sagun et al., 2018), while the IST contribution was implemented by accounting
for an interparticle interaction at high density. The corresponding parameters
were fitted to reproduce the nuclear matter ground-state properties, correct
behavior of the nuclear liquid-gas phase transition and its critical point
(Sagun et al., 2017) and proton flow constraint (Ivanytskyi et al., 2018).
Furthermore, in Sagun et al. (2019b) the model was generalized to describe NSs
showing a big application range of the unified IST approach Sagun et al.
(2019a). In the present work, we consider the Set B described in detail in
Sagun et al. (2020), while the crust is modeled in a simplified way by the
polytropic EoS with adiabatic index $\gamma=4/3$. This model parameterization
reproduces GW170817 and GW190425 tidal deformability limit (Abbott et al.,
2018, 2020a), NICER mass-radius measurements (Miller et al., 2019, 2021; Riley
et al., 2019, 2021; Raaijmakers et al., 2020), as well as the maximum mass
constraint.

In addition, we consider the DD2 EoS (Typel & Wolter, 1999; Typel et al.,
2010) with and without $\Lambda$ hyperons. The DD2 is a mean-field
relativistic nuclear model with density-dependent couplings, whose parameters
were fitted to the ground-state properties of nuclei. Hyperons have been
included in several works. In the present study, the density dependence of the
hyperon couplings to the $\sigma$, $\omega$ and $\rho$ mesons are considered
to be the same as one of the nucleons. For the $\phi$ coupling, the density
dependence of the $\omega$ meson is considered. The couplings of the $\sigma$
meson to the $\Lambda$ and $\Xi$ have been taken from Fortin et al. (2017) and
Fortin et al. (2020), respectively, and have been fitted to the binding energy
of $\Lambda$ and $\Xi$ hypernuclei. The coupling to the $\Sigma$ hyperon was
chosen so that the $\Sigma$ potential in symmetric nuclear matter is $+30$
MeV, see Gal et al. (2016) for a discussion. For the vector mesons, the quark
model predictions are used,

$\displaystyle g_{\omega\Lambda}=g_{\omega\Sigma}=\frac{2}{3}g_{\omega
N},\quad g_{\omega\Xi}=\frac{1}{3}g_{\omega N},$ $\displaystyle
g_{\phi\Lambda}=g_{\phi\Sigma}=-\frac{\sqrt{2}}{3}g_{\omega N},\quad
g_{\phi\Xi}=-\frac{2\sqrt{2}}{3}g_{\omega N}.$

Finally, the effective $\rho$-meson coupling is determined by the product of
the hyperon isospin with the $\rho$ meson-nucleon coupling. The DD2 EoS
without $\Lambda$ hyperons reproduces the maximum mass constraint
$(M_{max}=2.4M_{\odot})$ and NICER data, while $\Lambda_{1.4}$ is 700. At the
same time, the DD2 EoS with $\Lambda$ hyperons gives the maximum mass of
$2~{}M_{\odot}$. Further on the DD2 EoS with $\Lambda$ hyperons will be
referred as the DD2$\Lambda$.

The complete NS EoS contains, besides the core EoS, the BPS EoS Baym et al.
(1971) for the outer crust, and the inner crust was calculated within a
Thomas-Fermi calculation taking DD2 as the underlying model and allowing for
the appearance of several geometries as discussed in Grill et al. (2014). The
inner crust EoS has been published in Fortin et al. (2016).

##  3 Mixed system of two components

We assume no interaction between DM and BM, except through gravity. This
assumption is fully justified by the latest constraints coming from the DM
direct detection experiments and Bullet Cluster (Clowe et al., 2006; Randall
et al., 2008), showing that the DM-BM cross section is many orders of
magnitude lower than the typical nuclear one, $\sigma_{\chi}\sim 10^{-45}\
\mathrm{cm}^{2}\ll\sigma_{N}\sim 10^{-24}\ \mathrm{cm}^{2}$.

Therefore, the stress-energy tensors of both components are conserved
separately, leading to the system of the TOV equations with split components
(Oppenheimer & Volkoff, 1939; Tolman, 1939)

$\frac{dp_{i}}{dr}=-\frac{(\epsilon_{i}+p_{i})(M_{\mathrm{tot}}+4\pi
r^{3}p_{\mathrm{tot}})}{r^{2}\left(1-{2M_{\mathrm{tot}}}/{r}\right)},$ (9)

which describes the relativistic hydrostatic equilibrium of a DM-admixed NS.
In Eq. (9), the subscript index refers both to the BM and DM, i.e., $i=B,D$,
while $p_{\mathrm{tot}}\equiv p_{B}+p_{\chi}$ and $M(r)$ are the total
pressure and gravitational mass enclosed inside a sphere of radius $r$,
respectively,

$M_{i}(r)=4\pi\int^{r}_{0}\varepsilon_{i}(r^{\prime})r^{\prime 2}dr^{\prime}.$
(10)

Using Eq. (10), we define the total gravitational mass as the sum of the two
components, $M_{\mathrm{tot}}=M_{B}(R_{B})+M_{D}(R_{D})$, where the radii
$R_{i}$ are evaluated using the zero-pressure condition at the surface

$p_{i}(R_{i})=0.$ (11)

After having the total mass of the system, it is possible and convenient to
write the fraction of the accumulated DM as

$f_{\chi}=\frac{M_{D}}{M_{\mathrm{tot}}}.$ (12)

It is worth noting, that we refer to the microscopic/thermodynamic DM
parameters as $\chi$, while the macroscopic ones have an index $D$.

It is easy to obtain directly from Eq. (9) the relation between the chemical
potentials of the BM and DM. In fact, Ivanytskyi et al. (2020) showed that

$\frac{d\ln\mu_{B}}{dr}=\frac{d\ln\mu_{\chi}}{dr}=-\frac{M_{\mathrm{tot}}+4\pi
r^{3}p_{\mathrm{tot}}}{r^{2}(1-2M_{\mathrm{tot}}/r)},$ (13)

which yields the conclusion that the two chemical potentials are proportional
to each other. The value their ratio attains in the center of the star is the
proportionality constant, which can be used to simplify the model

$\mu_{\chi}=\left(\frac{\mu_{\chi}}{\mu_{B}}\right)_{r=0}\mu_{B}.$ (14)

|  |
---|---|---

|  |
---|---|---

Figure 2:  Left column: total gravitational mass of the DM-admixed NS as a
function of its visible radius $R$ obtained for $m_{\chi}$=100 MeV,
$m_{I}$=250 MeV (upper panel) and $m_{\chi}$=1000 MeV, $m_{I}$=1000 MeV (lower
panel). Black solid, dashed-dotted, and dotted curves correspond to pure BM
stars described by the IST EoS, DD2 EoS, and DD2 EoS with hyperons. Red, blue,
and green colours depict relative DM fractions equal to 1%, 3%, and 5%,
correspondingly. Green, gray, and teal bands represent 1$\sigma$ constraints
on mass of PSR J0348+0432 (Antoniadis et al., 2013), PSR J1810+1744 (Romani et
al., 2021), and PSR J0952-0607 (Romani et al., 2022). Pink and beige contours
show the NICER measurements of PSR J0030+0451 (Miller et al., 2019; Riley et
al., 2019), while orange and blue contours depict the PSR J0740+6620
measurements (Miller et al., 2021; Riley et al., 2021). LIGO-Virgo
observations of GW170817 (Abbott et al., 2018) and GW190425 (Abbott et al.,
2020a) binary NS mergers are shown in blue and magenta. Middle column: energy
density profiles for the BM (dotted curves) and DM (dashed curves) components
are shown for the DD2 EoS. The solid black curve represents the profile for
pure BM $1.4~{}M_{\odot}$ NS, while the other profiles were sampled to have
the same total gravitational mass. The upper panel is obtained for
$m_{\chi}$=100 MeV, $m_{I}$=250 MeV and the lower one for $m_{\chi}$=1000 MeV,
$m_{I}$=1000 MeV. Right column: the same as on the middle column, but for the
IST EoS.

By solving the TOV Eq. (9) with the boundary conditions and accounting for the
relation between both components from Eq. (14), we calculate the M-R relations
for DM-admixed NSs for different values of DM fractions $f_{\chi}$, particle’s
mass $m_{\chi}$, and the interaction scale $m_{I}$. To better understand the
impact of each parameter we consider light and heavy DM particles with
$m_{\chi}$=100 MeV and $m_{\chi}$=1000 MeV (see the left column of Fig. 2).
Moreover, to address our ignorance of the EoS for baryonic component we
studied the effect of DM on the soft IST EoS, depicted as a solid black curve
on the left panels of Fig. 2, as well as on the stiff DD2$\Lambda$ EoS (dotted
black curve) and DD2 EoS (dashed-dotted black curve). The chosen EoSs
represent different sides of mass and radius region allowed by the recent
astrophysical, GW, and nuclear physics constraints, and therefore, provide
good coverage of BM parameters. As it can be seen, the DD2$\Lambda$ EoS
(dotted black curve) and DD2 EoS coincide until $\sim 1.4~{}M_{\odot}$, a
point where the onset of hyperons happens. Further, hyperon production softens
the EoS leading to a smaller total maximum mass and star’s radius.

The left panels of Fig. 2 show the effect of DM with different relative
fractions inside a star on its mass and radius. Thus, we see a reduction of
$M_{\rm max}$ and radius of stars for larger DM fractions caused by a DM core
formation. In fact, the formation of more compact objects for an outside
observer would look like a softening of the BM EoS. This degeneracy between
the effect of DM and possible change of the strongly interacting matter
properties at high density will be discussed in Section 6.

Due to the fact that in the considered model at
$\mu_{\chi}\rightarrow\sqrt{2}m_{\chi}$ energy density diverges at finite
pressure, DM falls under the Schwarzschild radius forming a black hole. It
takes place for the high-mass stars for which the DM chemical potential in the
center reaches the limit (see the upper left panel of Fig. 2).

The panels on the middle and right columns of Fig. 2 demonstrate the split
energy density profiles of DM (dashed curves) and BM (dotted curves). The
solid black curve depicts the energy density profile for the $1.4~{}M_{\odot}$
star. The profiles for DM-admixed NSs are shown for stars with the same total
gravitational mass as the pure BM NS. As the onset of hyperons occurs after
$1.4M_{\odot}$, two formulations of the DD2 EoS give the same prediction for
the matter distribution inside the stars. Therefore, in Fig. 2 we show
profiles only for the DD2 EoS.

For heavy bosons a compact DM core is formed, which is seen from the high
values of the $\epsilon_{D}$, being an order of magnitude above $\epsilon_{B}$
(see the middle and right panels of the low row isn Fig. 2). Furthermore, the
$\epsilon_{D}$ drops to zero at radius $\sim$2 km corresponding to the size of
a DM core.

For the DM fraction 5% and $m_{\chi}=100$ MeV, $m_{I}=$ 250 MeV (see the
middle upper panel of Fig. 2) a DM halo is formed with the radius of 13.0 km.

##  4 Tidal deformability of DM-admixed NSs

The tidal deformability parameter $\lambda$ quantifies the response of an
object to a static external quadrupolar tidal field $\mathcal{E}_{ij}$ by
relating it to a quadrupolar moment
$\mathcal{Q}_{ij}=-\lambda\mathcal{E}_{ij}$. For a given stellar configuration
of the total mass $M_{\mathrm{tot}}$ and radius $R$ this tidal deformability
can be expressed through the Love number $k_{2}$ as $\lambda=2k_{2}R^{5}/3$
and is commonly mapped to the dimensionless
$\Lambda=\lambda/M_{\mathrm{tot}}^{5}$ (Hinderer, 2008). In the two-component
case, $R$ should be understood as the outermost radius, i.e., $R=R_{B}$ in the
DM core scenario and $R=R_{D}$ in the DM halo one. The Love number is defined
through the solution of an ordinary differential equation (ODE) appearing as a
leading order expansion of the Einstein equations with a metric perturbed by
the external gravitational field (Regge & Wheeler, 1957). The microscopic
properties of matter are encoded into this ODE through the change of total
pressure $p_{\mathrm{tot}}\equiv p_{B}+p_{\chi}$ caused by perturbation of the
total energy density
$\varepsilon_{\mathrm{tot}}\equiv\varepsilon_{B}+\varepsilon_{\chi}$. This
change is quantified by the derivative
$dp_{\mathrm{tot}}/d\varepsilon_{\mathrm{tot}}$. In the barotropic one-fluid
case, this derivative represents the corresponding speed of sound. In the two-
fluid case, the speed of sound derivation as
$dp_{\mathrm{tot}}/d\varepsilon_{\mathrm{tot}}$ is mathematically identical to
the expression obtained by Das et al. (2022). Therefore, in what follows, we
refer to it as the effective speed of sound of the two-fluid system. It can be
expressed through the speed of sound of baryonic $c_{s,B}^{2}$ and dark
$c_{s,\chi}^{2}$ components as

$\displaystyle c_{s,\mathrm{eff}}^{2}=\eta
c_{\mathrm{s,B}}^{2}+(1-\eta)c_{s,\chi}^{2}$ (15)

with $\eta\in[0,1]$. The lower and upper edges of this interval correspond to
the cases of pure DM and BM, respectively. Appendix B provides the derivation
of Eq. (15) and parameter $\eta$. This expression demonstrates that the
effective speed of sound lies between the ones of pure components.

In Fig. 3 we show the effective speed of sound for different
$\xi=\frac{\mu_{\chi}}{\mu_{B}}$ values, as well as the speed of sound for
pure BM and DM components. A relation between the parameters $\xi$ and $\eta$
is given in Eq. (B3) in Appendix B. The upper panel of Fig. 3 indicates how
the effective speed of sound behaves with DM accumulated in a core of a
compact star. Note, that it is in between the speed of sound values for pure
components. On the lower panel of Fig. 3 we see that the effective speed of
sound follows the BM, and only in the outer crust the DM component stars
dominate, which is related to a halo configuration.

Figure 3: The effective speed of sound for a mixture of BM and DM as a
function of total energy density. Upper panel: the curves were obtained for
$m_{\chi}$=750 MeV and $m_{I}$=250 MeV, which represents a DM core
configuration. Lower panel: the same as on the upper panel, but for
$m_{\chi}$=100 MeV and $m_{I}$=250 MeV illustrating a DM halo configuration.
The horizontal line at low densities corresponds to the polytropic EoS for the
crust. Figure 4: Tidal deformability as a function of total gravitational mass
calculated for pure BM stars (black curves) and DM-admixed NSs with relative
DM fractions 1%, 3%, and 5%, in red, blue, and green, correspondingly. Solid,
dash-dotted and dotted curves represent the IST EoS, DD2 EoS, and DD2$\Lambda$
EoS. The colors and symbols coincide with the ones used in Fig. 2 for a better
comparison. The figure is obtained for $m_{\chi}$=1000 MeV, $m_{I}$=1000 MeV.
Green, gray, and teal bands represent 1$\sigma$ constraints on mass of PSR
J0348+0432 (Antoniadis et al., 2013), PSR J1810+1744 (Romani et al., 2021),
and PSR J0952-0607 (Romani et al., 2022). The magenta area visualizes the
constraints obtained from GW170817 (Abbott et al., 2018).

As can be seen in Fig. 4, for the given total gravitational mass DM condensed
in a core leads to a smaller tidal deformability parameter compared to a pure
baryonic star. A similar effect has been shown in Fig. 2 for the radius. For a
distant observer, these effects will be perceived as an effective softening of
the EoS. On the other hand, the presence of a DM halo leads to a significant
increase in the outermost radius that goes beyond the BM component, an
increase of the tidal deformability parameter, and consequent effective
stiffening of the EoS. The considered IST, DD2, and DD2$\Lambda$ EoSs make us
conclude that the soft EoS, being on the lower limit of the GW170817 90% CL
region (see the magenta area in Fig. 4), provides a stringent constraint on a
DM core scenario, while the stiff EoS, being on the upper border of it, allows
much higher DM fractions, and disfavors an extended halo configuration. This
degeneracy between the effect of DM and strongly interacting matter properties
at high densities possesses limitations on DM detection, except for several DM
smoking guns that are going to be discussed in Section 6. Despite it, we have
to be aware of the fact that observational data on compact stars could be
affected by accumulated DM and, consequently, constraints we put on strongly
interacting matter at high densities.

##  5 Results

To study an interplay between boson mass and the interaction scale, as well as
to put constraints on the DM fraction, we perform a scan over those parameters
for the IST EoS (upper row), DD2 EoS (middle row), and DD2$\Lambda$ EoS
(bottom row) for fixed DM fractions of 1%, 3%, and 5% (see Fig. 5 in Appendix
C). The color maps represent the total maximum gravitational mass of DM-
admixed NSs. The white curve on each panel corresponds to
$M_{max}=1.4~{}M_{\odot}$, whereas the red curve represents
$M_{max}=2.0~{}M_{\odot}$. In the case $2.0~{}M_{\odot}$ configurations are
not reachable, we indicate $1.9~{}M_{\odot}$ stars with a green curve. As one
can see from the upper row on Fig. 5, the increase of the DM fraction narrows
the range of the values of the interaction scale $m_{I}$ consistent with the
masses of the heaviest known pulsars. On the other hand, the existence of the
high-mass stars with a significant amount of heavy DM requires low values of
the interaction scale.

We see the same dependence between $m_{\chi}$ and $m_{I}$ values. In fact,
lower $m_{I}\equiv m_{\omega}/g$ values correspond to the higher coupling
constant $g$ or, equivalently, stronger repulsion between the DM particles.
The IST EoS for any DM fraction is always in agreement with the tidal
deformability constraint, independently of $m_{\chi}$ and $m_{I}$ (see the
upper row in Fig. 5). At the same time, only for 1% and 3% of DM the total
maximum mass of DM-admixed NSs can reach $2.0~{}M_{\odot}$. Thus, to
simultaneously reproduce $2.0~{}M_{\odot}$ and GW170817 tidal deformability
constraints the boson mass and interaction scale are restricted to the values
shown in yellow. The shaded areas correspond to the non-allowed regions of
parameters that cannot simultaneously provide the heaviest pulsars and GW
constraints.

For 3% and 5% of DM the DD2 EoS reproduces both constraints in a wide range of
parameters disfavouring megaelectronvolt mass range of bosonic DM with low
values of the interaction strength. The black curve in the middle and bottom
rows of Fig. 5 depicts the GW170817 tidal deformability constraint
$\tilde{\Lambda}_{1.36}=720$ (Abbott et al., 2019b) above which the model is
consistent with the GW170817 merger. The dashed area corresponds to a non-
allowed range of parameters, including 1% of DM for the DD2 EoS. For the
DD2$\Lambda$ EoSs there are no $m_{I}$ and $m_{\chi}$ values that
simultaneously reproduce the heaviest pulsars and GW constraints. In fact,
only one of these criteria was reproduced for considered values of DM
fractions. This is directly related to the fact that at the onset of $\Lambda$
hyperons the EoS becomes softer in addition to the DM softening effect in a
core configuration.

From this analysis, we can conclude that, contrary to the stiff BM EoS (the
DD2 EoS, as an example), the soft BM EoS (the IST EoS, as an example) provides
a weaker limit on DM particle mass and interaction strength. This is related
to the fact that the pure baryonic DD2 or DD2$\Lambda$ EoSs are on the upper
border of the $\Lambda_{1.4}$ constraint from GW170817. Any decrease in
$\Lambda_{1.4}$ due to a DM core will not violate this condition, whereas a
small DM halo configuration will do it. As can be seen in Fig. 4, the IST EoS
is located on the lower limit of the magenta area favoring a halo formation.

It is worth noting that this result is obtained under the assumption of a
similar DM fraction in all galaxies. As a matter of fact, an application of
the GW170817 tidal deformability result and multi-messenger data as a
universal constraint on the amount of DM is questionable. Each galaxy could be
characterized by a different DM profile, as well as have local DM
inhomogeneities. Strictly speaking, GW170817 probes an amount of DM only in a
part of the NGC 4993, the host galaxy for this particular merger. Therefore, a
larger sample size of NS-NS and NS-BH mergers is required to constrain the DM
properties.

Due to current uncertainties of the BM EoS at high density, we cannot
discriminate between the effect of DM and the properties of BM. As it will be
discussed in the following Section 6, we expect a higher DM fraction inside
compact stars toward the Galactic Center. If so, the compact star population
would follow the scenarios presented from the left to right panels on Fig. 5,
i.e., from low to high DM fraction.

##  6 Discussions

As described above, there are various effects of DM on compact stars. A
natural question arises: how we can narrow down the proposed DM models and
constrain the DM properties using NSs? Can compact stars provide a smoking gun
evidence for the presence of DM? There are several different approaches:

(i) By measuring the mass, radius, and moment of inertia of NSs with few-
percent accuracy. Nowadays, NICER (Miller et al., 2019; Raaijmakers et al.,
2020; Miller et al., 2021; Raaijmakers et al., 2021) and in the near future
ATHENA (Cassano et al., 2018), eXTP (in ’t Zand et al., 2019), and STROBE-X
(Ray et al., 2019) are expected to measure $M$ and $R$ of NSs with a high
accuracy. Using the synthetic data for the STROBE-X telescope, and assuming
two NSs of the same mass and BM EoS, Rutherford et al. (2023) concluded that a
measurement of radii with a 2% accuracy would be enough to draw a conclusion
about the presence of DM in star’s interior. However, the existence of the
deconfinement phase transition in a core would exhibit in the same way,
leading to a degeneracy between the effect of DM and the phase transition. The
main drawback of this approach is that the effect of DM could mimic the
softening/stiffening of BM at high density and vice versa. Current
uncertainties of the baryonic EoS do not allow discrimination of two effects.
In addition, radio telescopes, e.g., MeerKAT (Bailes et al., 2018), SKA (Watts
et al., 2015), and ngVLA (Bower et al., 2018) plan to increase radio pulsar
timing and discover Galactic Center pulsars. A mass reduction of NSs toward
the Galactic Center or variation of mass, radius, and moment of inertia in
different parts of the Galaxy could shed light on the amount of accumulated DM
in compact stars. In fact, we could see a paucity of old millisecond pulsars
in the Galactic Center either due to light extinction on dust, or the collapse
of DM-admixed NSs into black holes after exceeding the Schwarzschild limit
(Bramante & Linden, 2014).

(ii) By performing binary numerical-relativity simulations and kilonova ejecta
for DM-admixed compact stars for different DM candidates, mass of particles,
interaction strength, and fractions with the further comparison to GW and
electromagnetic signals. The smoking gun of the presence of DM could be a
supplementary peak in the characteristic GW spectrum of NS mergers (Ellis et
al., 2018a), exotic waveforms (Giudice et al., 2016), modification of the
kilonova ejecta, or the presence of a strong oscillation mode in the waveforms
during the post-merger stage (Bezares et al., 2019). The next generation of GW
detectors, i.e., the Cosmic Explorer (CE) (Mills et al., 2018) and Einstein
Telescope (ET) (Punturo et al., 2010) will open another perspective of
detection of post-merger regimes and probing an internal composition of
compact stars.

(iii) By detecting a new feature in the binary Love relation (Yagi & Yunes,
2016). Thus, as it was shown in Fig. 3, DM could produce a bump, or any other
irregular behavior, in the effective speed of sound that would affect the
binary Love relation. Similar, as it was demonstrated for the strongly
interacting matter by Tan et al. (2022). This mark may be revealed by the next
generation of GW detectors that are planned to have the measurement precision
of $\delta\Lambda\sim 5$ for a GW170817-like event.

(iv) By detecting objects that go in contradiction with our understanding. A
potential candidate for DM-admixed NS could be the secondary component of
GW190814 (Abbott et al., 2020b). While likely being a black hole (Essick &
Landry, 2020; Tews et al., 2021), this compact object with the mass of $\sim
2.6~{}M_{\odot}$ raised debates about its nature (Tsokaros et al., 2020) as a
pure baryon matter EoS would not be able to explain a compact star of $\sim
2.6~{}M_{\odot}$. Hence, if not being a black hole, the compact object would
have to be supplemented either with exotic degrees of freedom, such as
hyperons and/or quarks (Tan et al., 2020; Dexheimer et al., 2021), an early
deconfinement phase transition (Ivanytskyi & Blaschke, 2022), very fast
rotation (Zhang & Li, 2020), or extra stiffening of the EoS at high densities
(Fattoyev et al., 2020). An alternative explanation of this puzzle would be a
DM-admixed NS (Di Giovanni et al., 2022), which could also explain a formation
of a black hole of so low mass as a collapsed DM-admixed NS (Bramante &
Linden, 2014).

The recently announced measurement of the central compact object within the
supernova remnant HESS J1731-347 (Doroshenko et al., 2022) is another object
that puzzles our understanding. This lightest and smallest compact star ever
observed could be explained as a NS admixed with DM (Sagun et al., 2023).

(v) Modification of the pulsar pulse profile due to the extra light-bending
(Miao et al., 2022) and/or gravitational microlensing in the case of the
existence of a dark halo.

(vi) Modification of the cooling rate of compact stars (de Lavallaz &
Fairbairn, 2010; Hamaguchi et al., 2019; Ángeles Pérez-García et al., 2022;
Buschmann et al., 2022). We want to note, that this effect is the most
inaccurate among the abovementioned ones. Thus, NSs need to have a well-
measured surface luminosity and age. In addition to it, uncertainties related
to particle composition, EoS, magnetic field, superfluidity/superconductivity,
NS masses, the chemical composition of an atmosphere, etc., could wash out an
effect of DM. Old NSs are less affected by the mentioned effects, as a photon
cooling stage starts to dominate over a neutrino cooling stage that is very
sensitive to a particle composition and superfluidity/superconductivity (Page
et al., 2004). The magnetic field is also expected to be unimportant for old
isolated NSs. Therefore, a possible heating mechanism of NSs due to DM
annihilation could be probed by increasing statistics on observational data of
old NSs.

##  7 Conclusions

We proposed a model of bosonic DM represented by a complex scalar field
coupled to the vector one through the covariant derivative, which is
equivalent to scalar electrodynamics with massive photons. The model describes
DM existing in the form of BEC with repulsive interaction. Pressure of the
present EoS saturates at asymptotically high densities leading to the
vanishing speed of sound and compressibility at this regime. From the
thermodynamic requirements, the chemical potential of DM existing as such BEC
is limited to the interval $\mu_{\chi}\in[m_{\chi},\sqrt{2}m_{\chi}]$, with
$m_{\chi}$ being the DM particle mass. In the weak and strong coupling limits,
this interval shrinks to its lower and upper bounds, respectively, while
pressure vanishes even at any density. This spectacular feature of the present
model makes its weak and strong coupling limits qualitatively similar and
requires further clarification.

At the same time, the strong coupling limit is shown to be inconsistent with
the Bullet Cluster cosmological data. Confronting the model prediction on the
cross section of the DM self-interaction to the results of numerical
simulations and observations allowed us to constrain the interaction scale
$m_{I}$ from below depending on the DM particle mass $m_{\chi}$.

DM-admixed compact stars were modeled by considering the mixed system of two
fluids with different relative fractions. The performed derivation of the
effective speed of sound for a two-fluid system allowed us to calculate the
tidal deformability parameter for compact stars admixed with different amounts
of DM. We argue that the one-fluid approach cannot be applied to a mixed
system of several components with the different proper speed of sound values.

To account for a discrepancy related to the baryonic component the soft IST
EoS and stiffer DD2 EoS with and without hyperons were considered. For
different DM particle’s mass, its relative fraction, and interaction scale we
found the conditions of DM core formation. We argue that in the framework of
the considered model only a small DM halo is possible, with the outermost
radius around twice the baryonic one. However, the total maximum gravitational
mass of this configuration is below 2 $M_{\odot}$.

We performed a thorough analysis of the effect of DM particle mass in the
100-1000 MeV mass range and self-interacting scale on maximum total
gravitational mass and tidal deformabilities of NSs for several fixed DM
fractions. We found that for 1%, 3% of DM for the IST EoS and 3%, 5% of DM for
the DD2 EoS the model can simultaneously reproduce the heaviest pulsars and
GW170817 tidal deformability constraint. The obtained allowed region of boson
mass $m_{\chi}$ and interaction scale $m_{I}$ for a fixed DM fraction shows an
anticorrelated dependence between these parameters, i.e., a high $m_{\chi}$
value favors a low $m_{I}$ value. For the DD2$\Lambda$ EoS no allowed region
of parameters was found due to the inability to simultaneously reproduce both
constraints.

The considered DM fractions up to 5% were chosen to demonstrate the effects on
compact star properties, and as it was discussed in the Introduction, are
higher than could be accumulated by Bondi accretion. While up-to-date
calculations are based on the interaction between DM and BM, the self-
interaction of DM could lead to enhanced DM accretion, hence the DM fraction
here proposed. However, we leave this for future studies.

In Section 6, we discussed the possible smoking gun signatures of DM in
compact stars that could be probed in the near future, e.g., alteration of
maximum total gravitational mass and radius of compact stars as a function of
a distance from the Galactic Center; modification of the surface temperature
(an additional heating or cooling mechanism) of NSs towards the Galactic
Center; lack of old millisecond pulsars in the Galactic center; the presence
of supplementary peak(s) in the GW signal from NS-NS and/or NS-BH mergers,
exotic waveform, or modification of the kilonova ejecta; gravitational-lensing
effect or alteration of the pulsar pulse profile due to the extra light
bending in a dark halo. Moreover, such objects as a secondary component of the
GW190814 event and central compact object within the supernova remnant HESS
J1731-347 challenge the existing models of compact stars and black holes,
offering the possibility of these objects being a DM-admixed NS.

We argue that compact stars and their mergers provide a novel sensitive
indirect method of detection and constraining the DM properties. Based on the
performed analysis it is clear that the present data analysis of X-ray, radio,
and GW observations without accounting for an accumulated DM could miss a
valuable piece of information as well as give an incorrect prediction about
the strongly interacting matter properties at high density.

The work of E.G., C.P., and V.S. was supported by national funds from FCT –
Fundação para a Ciência e a Tecnologia, I.P., within Projects No.
UIDB/04564/2020, UIDP/04564/2020, EXPL/FIS-AST/0735/2021. E.G. also
acknowledges the support from Project No. PRT/BD/152267/2021. C.P. is
supported by the Project No. PTDC/FIS-AST/28920/2017. V.S. also acknowledges
the PHAROS COST Action CA16214. The work of O.I. was supported by the Polish
National Science Center under the grant No. 2019/33/BST/03059.

##  Appendix A Lagrangian Model of the DM EOS

Here, we derive the DM EoS, presented with Eqs. (1) and (2), and we analyze it
in the weak and strong coupling regimes. The minimal Lagrangian describing the
chosen model of DM should include mass and kinetic terms of the complex scalar
$\chi$ and real vector $\omega^{\mu}$ fields, which are coupled through the
covariant derivative $D^{\mu}=\partial^{\mu}-ig\omega^{\mu}$ with $g$ being
the corresponding Yukawa coupling constant. Thus, we start with

$\mathcal{L}=(D_{\mu}\chi)^{*}D^{\mu}\chi-
m_{\chi}^{2}\chi^{*}\chi-\frac{\Omega_{\mu\nu}\Omega^{\mu\nu}}{4}+\frac{m_{\omega}^{2}\omega_{\mu}\omega^{\mu}}{2},$
(A1)

where $m_{\chi}$ and $m_{\omega}$ are masses of the scalar and vector fields,
respectively, and
$\Omega_{\mu\nu}=\partial_{\mu}\omega_{\nu}-\partial_{\nu}\omega_{\mu}$.
Before going further we would like to discuss the Neother current resulting
from the Lagrangian (A1)

$j^{\mu}=i\left[\chi^{*}\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\chi^{*})}-\chi\frac{\partial\mathcal{L}}{\partial(\partial_{\mu}\chi)}\right]=i\left[\chi^{*}D^{\mu}\chi-\chi(D^{\mu}\chi)^{*}\right]=i(\chi^{*}\partial^{\mu}\chi-\chi\partial^{\mu}\chi^{*})+2g\omega^{\mu}\chi^{*}\chi.$
(A2)

This current is equivalent to the one corresponding to local U(1) symmetry of
$\mathcal{L}$ (Brading & Brown, 2000). The density of conserved charge
associated to the DM particles is obtained by averaging the zeroth component
of this current. The second term in the expression for $j^{\mu}$ vanishes
under such averaging, since it includes an odd number of creation and
annihilation operators of the vector field. Thus, $n_{\chi}=\langle
j^{0}\rangle=i\langle\chi^{*}\partial^{0}\chi-\chi\partial^{0}\chi^{*}\rangle$.
This relation is artificially violated in the case when $\omega^{\mu}$ is
treated classically, e.g., within the mean-field approximation. This
discrepancy can be removed by using the Noether current resulting from the
invariance with respect to global U(1) transformation (Brading & Brown, 2000)

$j^{\mu}=i(\chi^{*}\partial^{\mu}\chi-\chi\partial^{\mu}\chi^{*}).$ (A3)

Below we use this expression for the DM four current.

Within the mean-field approximation operator of the vector field in this
Lagrangian is replaced by its constant expectation value. This value can be
obtained from the corresponding Euler-Lagrange equation

$(\partial^{2}+m_{\omega}^{2}+2g^{2}\chi^{*}\chi)\omega^{\mu}-\partial^{\mu}\partial^{\nu}\omega_{\nu}+gj^{\mu}=0.$
(A4)

In the hydrostatic case only the zeroth component of the DM four-current
attains nonzero mean value being nothing else as the DM particle number
density $n_{\chi}$, i.e., $\langle j^{\mu}\rangle=\eta^{\mu 0}n_{\chi}$ with
$\eta^{\mu\nu}$ standing for the Minkowski metric tensor. Replacing
$\omega^{\mu}$ by its constant mean value we eliminate derivatives of the
vector field in Eq. (A4). Furthermore, we replace $\chi^{*}\chi$ by the scalar
field condensate $\langle\chi^{*}\chi\rangle$ and DM current by its mean
value. Thus, the Euler-Lagrange equation of the vector field under the mean-
field approximation can be given a form of

$\langle\omega^{\mu}\rangle=\eta^{\mu 0}\omega\quad{\rm
with}\quad\omega=-\frac{gn_{\chi}}{m_{\omega}^{2}+2g^{2}\langle\chi^{*}\chi\rangle}.$
(A5)

At finite particle number densities, DM has finite chemical potentials
$\mu_{\chi}$, which serves as a Lagrange multiplier in
$\mathcal{L}+\mu_{\chi}j^{0}$ ensuring that mean value of $j^{0}$ coincides
with $n_{\chi}$. Inserting Eq. (A5) into Eq. (A1), one gets

$\mathcal{L}+\mu_{\chi}j^{0}=\partial_{\mu}\chi^{*}\partial^{\mu}\chi-
M_{\chi}^{2}\chi^{*}\chi+\nu_{\chi}j^{0}+\frac{m_{\omega}^{2}\omega^{2}}{2}.$
(A6)

The first three terms on the right-hand side of this equation represent free
quasi-particles with the effective mass and chemical potential defined as

$\displaystyle M_{\chi}^{2}$ $\displaystyle\equiv$ $\displaystyle
m_{\chi}^{2}-g^{2}\omega^{2},$ (A7) $\displaystyle\nu_{\chi}$
$\displaystyle\equiv$ $\displaystyle\mu_{\chi}+g\omega,$ (A8)

respectively. At zero temperature these bosonic quasi-particles condense to
zero mode with $\langle\chi^{*}\chi\rangle=\zeta^{2}$ and $\zeta$ being an
amplitude of this mode. This BEC contributes to the DM pressure
$\zeta^{2}(\nu_{\chi}^{2}-M_{\chi}^{2})$ (see, e.g., chapter 2.4 of Kapusta &
Gale (2006)). The last term in Eq. (A6) does not include any dynamical
variables, and therefore, simply renormalizes the pressure, which becomes

$p_{\chi}=\zeta^{2}(\nu_{\chi}^{2}-M_{\chi}^{2})+\frac{m_{\omega}^{2}\omega^{2}}{2}.$
(A9)

Mean value of the vector field defined by Eq. (A4) maximizes the pressure. The
same is the case for $\zeta$, i.e., $\frac{\partial
p_{\chi}}{\partial\omega}=\frac{\partial p_{\chi}}{\partial\zeta}=0$. These
conditions should be supplemented with a definition of the DM particle number
density given by the thermodynamic identity $n_{\chi}=\frac{\partial
p_{\chi}}{\partial\mu_{\chi}}$ yielding the following system of equations,
which should be solved in order to find $\omega$ and $\zeta$ and construct the
DM EoS:

$\displaystyle 2g\zeta^{2}(\nu_{\chi}+g\omega)+m_{\omega}^{2}\omega=0,$ (A10)
$\displaystyle 2\zeta(\nu_{\chi}^{2}-M_{\chi}^{2})=0,$ (A11) $\displaystyle
n_{\chi}=2\zeta^{2}\nu_{\chi}.$ (A12)

Before solving this system, we want to demonstrate that it is consistent with
Eq. (A5). For this we formally express $\nu_{\chi}$ from Eq. (A12) and insert
the result to Eq. (A10) yielding
$gn_{\chi}+(2g^{2}\zeta^{2}+m_{\omega}^{2})\omega=0$. Then, replacing
$\zeta^{2}$ by $\langle\chi^{*}\chi\rangle$ this relation can be written in
the desired form.

In the BEC, amplitude of the bosonic zero mode is $\zeta\neq 0$. Therefore,
Eq. (A11) requires $\nu_{\chi}=M_{\chi}$, which can be solved with respect to
the vector field as

$\omega=\frac{-\mu_{\chi}\pm\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}}{2g}.$ (A13)

From this condition we immediately conclude that physical values of the DM
chemical potential are limited to the range
$-\sqrt{2}m_{\chi}\leq\mu_{\chi}\leq\sqrt{2}m_{\chi}$. At positive
$\mu_{\chi}$ BEC is constituted by the DM particles, while negative
$\mu_{\chi}$ corresponds to antiparticles. For definiteness in this work, we
consider the case of the DM particles, which is equivalent to requiring
$\mu_{\chi}\geq 0$. Eq. (A13) ensures that the first term in Eq. (A9) vanishes
and DM pressure $p_{\chi}\propto\omega^{2}$. At zero $n_{\chi}$ this pressure
should vanish provided by $\omega=0$. At non-negative chemical potential of DM
this condition can be fulfilled only at $\mu_{\chi}=m_{\chi}$ if the sign
“$+$” is chosen in Eq. (A13). In order to obtain the DM particle number
density we first express $\zeta^{2}$ from Eq. (A10) and then insert the result
into Eq. (A12). Thus, at $\mu_{\chi}\in[m_{\chi},\sqrt{2}m_{\chi}]$ the DM EoS
becomes

$\displaystyle p_{\chi}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(\frac{m_{\omega}}{g}\right)^{2}\left(m_{\chi}^{2}-\mu_{\chi}\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}\right),$
(A14) $\displaystyle n_{\chi}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{m_{\omega}}{g}\right)^{2}\frac{\mu_{\chi}^{2}-m_{\chi}^{2}}{\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}},$
(A15)

while at $\mu_{\chi}\in[0,m_{\chi}]$ one gets $p_{\chi}=0$ and $n_{\chi}=0$.
It is seen from these expressions that vector field mass $m_{\omega}$ and
coupling $g$ do not enter the DM EoS independently but appear as the ratio
$m_{I}\equiv m_{\omega}/g$, which is a relevant parameter. With this notation
and thermodynamic identity $\varepsilon_{\chi}=\mu_{\chi}n_{\chi}-p_{\chi}$ we
arrive at Equations (1) and (2).

The weak and strong coupling limits of the present EoS are obtained at
$g\rightarrow 0$ and $g\rightarrow\infty$, respectively. This corresponds to
$m_{I}\rightarrow\infty$ and $m_{I}\rightarrow 0$. In order to consider these
limits we treat the DM particle density as an independent quantity. For this
we first write Eq. (A15) as a quadratic equation for
$\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}$ and solve it as

$\displaystyle\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}=\sqrt{m_{\chi}^{2}+\frac{n_{\chi}^{2}}{m_{I}^{4}}}-\frac{n_{\chi}}{m_{I}^{2}}.$
(A16)

This solution allows us to express $\mu_{\chi}$ and expand it up to the
leading order in $m_{I}^{-2}$ or $m_{I}^{2}$

$\displaystyle\mu_{\chi}=\left\\{\begin{array}[]{l}m_{\chi}+\frac{n_{\chi}}{m_{I}^{2}}+\mathcal{O}\left(m_{I}^{-4}\right),\quad\quad~{}~{}~{}m_{I}\rightarrow\infty\\\
\sqrt{2}m_{\chi}-\frac{m_{\chi}^{3}m_{I}^{4}}{8\sqrt{2}n_{\chi}^{2}}+\mathcal{O}\left(m_{I}^{8}\right),\quad
m_{I}\rightarrow 0.\end{array}\right.$ (A19)

From this we conclude that in the weak coupling limit chemical potential of DM
converges to the smallest value of the BEC interval
$\mu_{\chi}\in[m_{\chi},\sqrt{2}m_{\chi}]$, while in the strong coupling limit
it converges to the largest value. This conclusion holds for any $n_{\chi}$.
The DM pressure behaves as

$\displaystyle
p_{\chi}=\frac{m_{I}^{2}}{8}\left(\mu_{\chi}-\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}\right)^{2}=\left\\{\begin{array}[]{l}\frac{n_{\chi}^{2}}{2m_{I}^{2}}+\mathcal{O}\left(m_{I}^{-4}\right),\quad
m_{I}\rightarrow\infty\\\
\frac{m_{\chi}^{2}m_{I}^{2}}{4}+\mathcal{O}\left(m_{I}^{4}\right),\quad
m_{I}\rightarrow 0.\end{array}\right.$ (A22)

where on the second step $\sqrt{2m_{\chi}^{2}-\mu_{\chi}^{2}}$ was
approximated by $m_{\chi}-n_{\chi}/m_{I}^{2}$ at $m_{I}\rightarrow\infty$,
while at $m_{I}\rightarrow 0$ it was neglected compared to $\mu_{\chi}$. In
both of the considered limits the pressure vanishes leading to the energy
density mentioned in Subsection 2.1.

##  Appendix B Effective Speed of Sound of Two-Fluid System

In order to calculate $c_{s,\chi}^{2}$ we notice that in the GCE pressure and
energy density of each component are functions of the corresponding chemical
potential. With this we can write

$\displaystyle
c_{s,tot}^{2}=\frac{dp_{tot}}{d\varepsilon_{tot}}=\frac{\frac{\partial
p_{B}}{\partial\mu_{B}}+\frac{\partial
p_{\chi}}{\partial\mu_{\chi}}\frac{d\mu_{\chi}}{d\mu_{B}}}{\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}+\frac{\partial\varepsilon_{\chi}}{\partial\mu_{\chi}}\frac{d\mu_{\chi}}{d\mu_{B}}}=\frac{\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}c_{s,B}^{2}+\frac{\partial\varepsilon_{\chi}}{\partial\mu_{\chi}}\frac{d\mu_{\chi}}{d\mu_{B}}c_{s,\chi}^{2}}{\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}+\frac{\partial\varepsilon_{\chi}}{\partial\mu_{\chi}}\frac{d\mu_{\chi}}{d\mu_{B}}},$
(B1)

where on the last step we used identities $c_{s,B}^{2}=\frac{\partial
p_{B}}{\partial\mu_{B}}\bigl{/}\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}$
and $c_{s,\chi}^{2}=\frac{\partial
p_{\chi}}{\partial\mu_{\chi}}\bigl{/}\frac{\partial\varepsilon_{\chi}}{\partial\mu_{\chi}}$
in order to express derivatives of the pressures of two components with
respect to the corresponding chemical potentials. This expression can be given
the form of Eq. (15) with $\eta$ defined as

$\displaystyle\eta=\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}\left[\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}+\frac{\partial\varepsilon_{\chi}}{\partial\mu_{\chi}}\frac{d\mu_{\chi}}{d\mu_{B}}\right]^{-1},$
(B2)

The thermodynamic identities $n_{B}=\frac{\partial p_{B}}{\partial\mu_{B}}$
and $\varepsilon_{B}=\mu_{B}n_{B}-p_{B}$ can be used in order to obtain
$\frac{\partial\varepsilon_{B}}{\partial\mu_{B}}=\mu_{B}\frac{\partial
n_{B}}{\partial\mu_{B}}$. We similarly obtain
$\frac{\partial\varepsilon_{\chi}}{\partial\mu_{\chi}}=\mu_{\chi}\frac{\partial
n_{\chi}}{\partial\mu_{\chi}}$. Within stellar interiors chemical potential of
two components scale proportionally to each other (Ivanytskyi et al., 2020).
This allows us to conclude that
$\frac{d\mu_{\chi}}{d\mu_{B}}=\frac{\mu_{\chi}}{\mu_{B}}=\xi$. Thus, Eq. (B2)
becomes

$\displaystyle\eta=\frac{\partial n_{B}}{\partial\mu_{B}}\left[\frac{\partial
n_{B}}{\partial\mu_{B}}+\xi^{2}\frac{\partial
n_{\chi}}{\partial\mu_{\chi}}\right]^{-1}.$ (B3)

Mechanical stability of BM and DM requires $\frac{\partial
n_{B}}{\partial\mu_{B}}>0$ and $\frac{\partial n_{\chi}}{\partial\mu_{\chi}}$,
respectively. In this case $\eta$ by construction is limited to the interval
$\eta\in[0,1]$. The lower edge of this interval $\eta=0$ is obtained at
$\mu_{B}=0$, which corresponds to the absence of BM. On the other hand, the
case of pure BM is obtained at $\mu_{\chi}=0$ yielding $\eta=1$.

##  Appendix C Scan Over Boson Mass and the Interaction Scale

|  |
---|---|---

|  |
---|---|---

|  |
---|---|---

Figure 5:  Upper row: Parameter space in the $m_{I}-m_{\chi}$ plane calculated
for the IST EoS and different values of DM fraction, 1% (left panel), 3%
(middle panel), 5% (right panel). The color represents the total maximum
gravitational mass of DM-admixed NSs. The white, and red curves correspond to
stars with the total maximum gravitational mass equal to 1.4$~{}M_{\odot}$,
and 2$~{}M_{\odot}$, respectively. For a DM fraction 5% (see the right panel)
the maximum mass do not reach 2$~{}M_{\odot}$, thus, we show 1.9$~{}M_{\odot}$
configurations in green. Shaded areas correspond to the non-allowed regions of
parameters, whereas simultaneously the astrophysical and GW constraints are
not fulfilled. Middle row: The same as on the upper row, but calculated for
DD2 EoS. The tidal deformability constraint $\tilde{\Lambda}_{1.36}=720$ (90%
credible interval) is shown as a black curve. All the area above this curve is
consistent with the GW170817 tidal deformability constraint (Abbott et al.,
2019b). The 1% DM case is not consistent with the tidal deformability
constraint for all values of $m_{I}$ and $m_{\chi}$. Lower row: The same as on
the middle and upper rows, but calculated for DD2$\Lambda$ EoS.

## References

  * Abbott et al. (2017a) Abbott, B. P., et al. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, Astrophys. J. Lett., 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2018) —. 2018, Phys. Rev. Lett., 121, 1611012013Sci…340..448A, doi: 10.1103/PhysRevLett.121.161101
  * Abbott et al. (2019a) —. 2019a, Phys. Rev. Lett., 123, 011102, doi: 10.1103/PhysRevLett.123.011102
  * Abbott et al. (2019b) —. 2019b, Phys. Rev. X, 9, 011001, doi: 10.1103/PhysRevX.9.011001
  * Abbott et al. (2020a) —. 2020a, Astrophys. J. Lett., 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2020b) Abbott, R., et al. 2020b, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Ahn & Shapiro (2005) Ahn, K.-J., & Shapiro, P. R. 2005, Mon. Not. Roy. Astron. Soc., 363, 1092, doi: 10.1111/j.1365-2966.2005.09492.x
  * Ángeles Pérez-García et al. (2022) Ángeles Pérez-García, M., Grigorian, H., Albertus, C., Barba, D., & Silk, J. 2022, Phys. Lett. B, 827, 136937, doi: 10.1016/j.physletb.2022.136937
  * Annala et al. (2018) Annala, E., Gorda, T., Kurkela, A., & Vuorinen, A. 2018, Phys. Rev. Lett., 120, 172703, doi: 10.1103/PhysRevLett.120.172703
  * Antoniadis et al. (2013) Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448, doi: 10.1126/science.1233232
  * Bailes et al. (2018) Bailes, M., et al. 2018, PoS, MeerKAT2016, 011, doi: 10.22323/1.277.0011
  * Bauswein et al. (2023) Bauswein, A., Guo, G., Lien, J.-H., Lin, Y.-H., & Wu, M.-R. 2023, Phys. Rev. D, 107, 083002, doi: 10.1103/PhysRevD.107.083002
  * Bauswein et al. (2017) Bauswein, A., Just, O., Janka, H.-T., & Stergioulas, N. 2017, Astrophys. J. Lett., 850, L34, doi: 10.3847/2041-8213/aa9994
  * Baym et al. (1971) Baym, G., Pethick, C., & Sutherland, P. 1971, Astrophys. J., 170, 299, doi: 10.1086/151216
  * Bell et al. (2021) Bell, N. F., Busoni, G., Motta, T. F., et al. 2021, Phys. Rev. Lett., 127, 111803, doi: 10.1103/PhysRevLett.127.111803
  * Bell et al. (2020) Bell, N. F., Busoni, G., Robles, S., & Virgato, M. 2020, JCAP, 09, 028, doi: 10.1088/1475-7516/2020/09/028
  * Bertone & Fairbairn (2008) Bertone, G., & Fairbairn, M. 2008, Phys. Rev. D, 77, 043515, doi: 10.1103/PhysRevD.77.043515
  * Bezares et al. (2019) Bezares, M., Viganò, D., & Palenzuela, C. 2019, Phys. Rev. D, 100, 044049, doi: 10.1103/PhysRevD.100.044049
  * Bower et al. (2018) Bower, G. C., et al. 2018, ASP Conf. Ser., 517, 793.  https://arxiv.org/abs/1810.06623
  * Brading & Brown (2000) Brading, K., & Brown, H. R. 2000, arXiV.  https://arxiv.org/abs/hep-th/0009058
  * Bramante et al. (2013) Bramante, J., Fukushima, K., & Kumar, J. 2013, Phys. Rev. D, 87, 055012, doi: 10.1103/PhysRevD.87.055012
  * Bramante & Linden (2014) Bramante, J., & Linden, T. 2014, Phys. Rev. Lett., 113, 191301, doi: 10.1103/PhysRevLett.113.191301
  * Burger et al. (2022) Burger, J. D., Zavala, J., Sales, L. V., et al. 2022, Phys. Rev. Lett., 129, 191103, doi: 10.1103/PhysRevLett.129.191103
  * Buschmann et al. (2022) Buschmann, M., Dessert, C., Foster, J. W., Long, A. J., & Safdi, B. R. 2022, Phys. Rev. Lett., 128, 091102, doi: 10.1103/PhysRevLett.128.091102
  * Cassano et al. (2018) Cassano, R., et al. 2018, arXiv.  https://arxiv.org/abs/1807.09080
  * Chang et al. (2018) Chang, J. H., Essig, R., & McDermott, S. D. 2018, JHEP, 09, 051, doi: 10.1007/JHEP09(2018)051
  * Ciarcelluti & Sandin (2011) Ciarcelluti, P., & Sandin, F. 2011, Phys. Lett. B, 695, 19, doi: 10.1016/j.physletb.2010.11.021
  * Clowe et al. (2006) Clowe, D., Bradač, M., Gonzalez, A. H., et al. 2006, The Astrophysical Journal, 648, L109, doi: 10.1086/508162
  * Colpi et al. (1986) Colpi, M., Shapiro, S. L., & Wasserman, I. 1986, Phys. Rev. Lett., 57, 2485, doi: 10.1103/PhysRevLett.57.2485
  * Das et al. (2022) Das, A., Malik, T., & Nayak, A. C. 2022, Phys. Rev. D, 105, 123034, doi: 10.1103/PhysRevD.105.123034
  * Das et al. (2020) Das, H. C., Kumar, A., Kumar, B., et al. 2020, Mon. Not. R. Astron. Soc., 495, 4893, doi: 10.1093/mnras/staa1435
  * de Lavallaz & Fairbairn (2010) de Lavallaz, A., & Fairbairn, M. 2010, Phys. Rev. D, 81, 123521, doi: 10.1103/PhysRevD.81.123521
  * Del Popolo et al. (2020) Del Popolo, A., Le Delliou, M., & Deliyergiyev, M. 2020, Universe, 6, 222, doi: 10.3390/universe6120222
  * Dengler et al. (2022) Dengler, Y., Schaffner-Bielich, J., & Tolos, L. 2022, Phys. Rev. D, 105, 043013, doi: 10.1103/PhysRevD.105.043013
  * Dexheimer et al. (2021) Dexheimer, V., Gomes, R. O., Klähn, T., Han, S., & Salinas, M. 2021, Phys. Rev. C, 103, 025808, doi: 10.1103/PhysRevC.103.025808
  * Di Giovanni et al. (2022) Di Giovanni, F., Sanchis-Gual, N., Cerdá-Durán, P., & Font, J. A. 2022, Phys. Rev. D, 105, 063005, doi: 10.1103/PhysRevD.105.063005
  * Dietrich & Clough (2019) Dietrich, T., & Clough, K. 2019, Phys. Rev. D, 100, 083005, doi: 10.1103/PhysRevD.100.083005
  * Doroshenko et al. (2022) Doroshenko, V., Suleimanov, V., Pühlhofer, G., & Santangelo, A. 2022, Nature Astronomy, doi: 10.1038/s41550-022-01800-1
  * Ellis et al. (2018a) Ellis, J., Hektor, A., Hütsi, G., et al. 2018a, Phys. Lett. B, 781, 607, doi: 10.1016/j.physletb.2018.04.048
  * Ellis et al. (2018b) Ellis, J., Hütsi, G., Kannike, K., et al. 2018b, Phys. Rev. D, 97, 123007, doi: 10.1103/PhysRevD.97.123007
  * Emma et al. (2022) Emma, M., Schianchi, F., Pannarale, F., Sagun, V., & Dietrich, T. 2022, Particles, 5, 273, doi: 10.3390/particles5030024
  * Essick & Landry (2020) Essick, R., & Landry, P. 2020, Astrophys. J., 904, 80, doi: 10.3847/1538-4357/abbd3b
  * Fattoyev et al. (2020) Fattoyev, F. J., Horowitz, C. J., Piekarewicz, J., & Reed, B. 2020, Phys. Rev. C, 102, 065805, doi: 10.1103/PhysRevC.102.065805
  * Fonseca et al. (2021) Fonseca, E., et al. 2021, Astrophys. J. Lett., 915, L12, doi: 10.3847/2041-8213/ac03b8
  * Fortin et al. (2017) Fortin, M., Avancini, S. S., Providência, C., & Vidaña, I. 2017, Phys. Rev. C, 95, 065803, doi: 10.1103/PhysRevC.95.065803
  * Fortin et al. (2016) Fortin, M., Providencia, C., Raduta, A. R., et al. 2016, Phys. Rev. C, 94, 035804, doi: 10.1103/PhysRevC.94.035804
  * Fortin et al. (2020) Fortin, M., Raduta, A. R., Avancini, S., & Providência, C. 2020, Phys. Rev. D, 101, 034017, doi: 10.1103/PhysRevD.101.034017
  * Gal et al. (2016) Gal, A., Hungerford, E. V., & Millener, D. J. 2016, Rev. Mod. Phys., 88, 035004, doi: 10.1103/RevModPhys.88.035004
  * Giudice et al. (2016) Giudice, G. F., McCullough, M., & Urbano, A. 2016, JCAP, 10, 001, doi: 10.1088/1475-7516/2016/10/001
  * Goldman et al. (2013) Goldman, I., Mohapatra, R. N., Nussinov, S., Rosenbaum, D., & Teplitz, V. 2013, Phys. Lett. B, 725, 200, doi: 10.1016/j.physletb.2013.07.017
  * Gresham et al. (2017) Gresham, M. I., Lou, H. K., & Zurek, K. M. 2017, Phys. Rev. D, 96, 096012, doi: 10.1103/PhysRevD.96.096012
  * Grill et al. (2014) Grill, F., Pais, H., Providência, C., Vidaña, I., & Avancini, S. S. 2014, Phys. Rev. C, 90, 045803, doi: 10.1103/PhysRevC.90.045803
  * Hamaguchi et al. (2019) Hamaguchi, K., Nagata, N., & Yanagi, K. 2019, Phys. Lett. B, 795, 484, doi: 10.1016/j.physletb.2019.06.060
  * Hinderer (2008) Hinderer, T. 2008, The Astrophysical Journal, 677, 1216, doi: 10.1086/533487
  * Hinderer et al. (2019) Hinderer, T., et al. 2019, Phys. Rev. D, 100, 06321, doi: 10.1103/PhysRevD.100.063021
  * in ’t Zand et al. (2019) in ’t Zand, J. J. M., et al. 2019, Sci. China Phys. Mech. Astron., 62, 029506, doi: 10.1007/s11433-017-9186-1
  * Ivanytskyi et al. (2018) Ivanytskyi, A. I., Bugaev, K. A., Sagun, V. V., Bravina, L. V., & Zabrodin, E. E. 2018, Phys. Rev. C, 97, 064905, doi: 10.1103/PhysRevC.97.064905
  * Ivanytskyi & Blaschke (2022) Ivanytskyi, O., & Blaschke, D. 2022, Phys. Rev. D, 105, 114042, doi: 10.1103/PhysRevD.105.114042
  * Ivanytskyi et al. (2020) Ivanytskyi, O., Sagun, V., & Lopes, I. 2020, Phys. Rev. D, 102, 063028, doi: 10.1103/PhysRevD.102.063028
  * Kapusta & Gale (2006) Kapusta, J. I., & Gale, C. 2006, Finite-Temperature Field Theory: Principles and Appli cations, 2nd edn., Cambridge Monographs on Mathematical Physics (Cambridge University Press), doi: 10.1017/CBO9780511535130
  * Karkevandi et al. (2022) Karkevandi, D. R., Shakeri, S., Sagun, V., & Ivanytskyi, O. 2022, Physical Review D, 105, doi: 10.1103/physrevd.105.023001
  * Kouvaris (2008) Kouvaris, C. 2008, Phys. Rev. D, 77, 023006, doi: 10.1103/PhysRevD.77.023006
  * Kouvaris (2013a) —. 2013a, The Dark Side of Neutron Stars, arXiv, doi: 10.48550/ARXIV.1308.3222
  * Kouvaris (2013b) —. 2013b, Adv. High Energy Phys., 2013, 856196, doi: 10.1155/2013/856196
  * Kouvaris & Nielsen (2015) Kouvaris, C., & Nielsen, N. G. 2015, Phys. Rev. D, 92, 063526, doi: 10.1103/PhysRevD.92.063526
  * Kouvaris & Tinyakov (2010) Kouvaris, C., & Tinyakov, P. 2010, Phys. Rev. D, 82, 063531, doi: 10.1103/PhysRevD.82.063531
  * Leung et al. (2022) Leung, K.-L., Chu, M.-c., & Lin, L.-M. 2022, Phys. Rev. D, 105, 123010, doi: 10.1103/PhysRevD.105.123010
  * Maselli et al. (2017) Maselli, A., Pnigouras, P., Nielsen, N. G., Kouvaris, C., & Kokkotas, K. D. 2017, Phys. Rev. D, 96, 023005, doi: 10.1103/PhysRevD.96.023005
  * Miao et al. (2022) Miao, Z., Zhu, Y., Li, A., & Huang, F. 2022, Astrophys. J., 936, 69, doi: 10.3847/1538-4357/ac8544
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, The Astrophysical Journal, 887, L24, doi: 10.3847/2041-8213/ab50c5
  * Miller et al. (2021) Miller, M. C., et al. 2021, Astrophys. J. Lett., 918, L28, doi: 10.3847/2041-8213/ac089b
  * Mills et al. (2018) Mills, C., Tiwari, V., & Fairhurst, S. 2018, Phys. Rev. D, 97, 104064, doi: 10.1103/PhysRevD.97.104064
  * Moore (1994) Moore, B. 1994, Nature, 370, 629, doi: 10.1038/370629a0
  * Nelson et al. (2019) Nelson, A. E., Reddy, S., & Zhou, D. 2019, J. Cosm. Astrop. Phys., 2019, 012, doi: 10.1088/1475-7516/2019/07/012
  * Oppenheimer & Volkoff (1939) Oppenheimer, J. R., & Volkoff, G. M. 1939, Phys. Rev., 55, 374, doi: 10.1103/PhysRev.55.374
  * Page et al. (2004) Page, D., Lattimer, J. M., Prakash, M., & Steiner, A. W. 2004, Astrophys. J. Suppl., 155, 623, doi: 10.1086/424844
  * Pérez-García & Silk (2012) Pérez-García, M. Á., & Silk, J. 2012, Phys. Lett. B, 711, 6, doi: 10.1016/j.physletb.2012.03.065
  * Petraki & Volkas (2013) Petraki, K., & Volkas, R. R. 2013, Int. J. Mod. Phys. A, 28, 1330028, doi: 10.1142/S0217751X13300287
  * Punturo et al. (2010) Punturo, M., et al. 2010, Class. Quant. Grav., 27, 194002, doi: 10.1088/0264-9381/27/19/194002
  * Raaijmakers et al. (2020) Raaijmakers, G., et al. 2020, Astrophys. J. Lett., 893, L21, doi: 10.3847/2041-8213/ab822f
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, Astrophys. J. Lett., 918, L29, doi: 10.3847/2041-8213/ac089a
  * Randall et al. (2008) Randall, S. W., Markevitch, M., Clowe, D., Gonzalez, A. H., & Bradac, M. 2008, Astrophys. J., 679, 1173, doi: 10.1086/587859
  * Ray et al. (2019) Ray, P. S., et al. 2019, arxiv.  https://arxiv.org/abs/1903.03035
  * Regge & Wheeler (1957) Regge, T., & Wheeler, J. A. 1957, Physical Review, 108, 1063, doi: 10.1103/PhysRev.108.1063
  * Riley et al. (2019) Riley, T. E., et al. 2019, Astrophys. J. Lett., 887, L21, doi: 10.3847/2041-8213/ab481c
  * Riley et al. (2021) —. 2021, Astrophys. J. Lett., 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Romani et al. (2021) Romani, R. W., Kandel, D., Filippenko, A. V., Brink, T. G., & Zheng, W. 2021, Astrophys. J. Lett., 908, L46, doi: 10.3847/2041-8213/abe2b4
  * Romani et al. (2022) —. 2022, Astrophys. J. Lett., 934, L18, doi: 10.3847/2041-8213/ac8007
  * Rutherford et al. (2023) Rutherford, N., Raaijmakers, G., Prescod-Weinstein, C., & Watts, A. 2023, Phys. Rev. D, 107, 103051, doi: 10.1103/PhysRevD.107.103051
  * Sagun et al. (2023) Sagun, V., Giangrandi, E., Dietrich, T., et al. 2023, In preparation
  * Sagun et al. (2022) Sagun, V., Giangrandi, E., Ivanytskyi, O., Lopes, I., & Bugaev, K. 2022, Proc. Sc., PANIC2021, 313, doi: 10.22323/1.380.0313
  * Sagun et al. (2019a) Sagun, V., Lopes, I., & Ivanytskyi, A. 2019a, Nucl. Phys. A, 982, 883, doi: 10.1016/j.nuclphysa.2018.10.024
  * Sagun et al. (2020) Sagun, V., Panotopoulos, G., & Lopes, I. 2020, Phys. Rev. D, 101, 063025, doi: 10.1103/PhysRevD.101.063025
  * Sagun et al. (2017) Sagun, V. V., Bugaiev, K. A., Ivanytskyi, A. I., Oliinychenko, D. R., & Mishustin, I. N. 2017, EPJ Web Conf., 137, 09007, doi: 10.1051/epjconf/201713709007
  * Sagun et al. (2019b) Sagun, V. V., Lopes, I., & Ivanytskyi, A. I. 2019b, The Astrophysical Journal, 871, 157, doi: 10.3847/1538-4357/aaf805
  * Sagun et al. (2018) Sagun, V. V., Bugaev, K. A., Ivanytskyi, A. I., et al. 2018, Eur. Phys. J. A, 54, 100, doi: 10.1140/epja/i2018-12535-1
  * Sedrakian (2016) Sedrakian, A. 2016, Phys. Rev. D, 93, 065044, doi: 10.1103/PhysRevD.93.065044
  * Sedrakian (2019) —. 2019, Phys. Rev. D, 99, 043011, doi: 10.1103/PhysRevD.99.043011
  * Tan et al. (2022) Tan, H., Dexheimer, V., Noronha-Hostler, J., & Yunes, N. 2022, Phys. Rev. Lett., 128, 161101, doi: 10.1103/PhysRevLett.128.161101
  * Tan et al. (2020) Tan, H., Noronha-Hostler, J., & Yunes, N. 2020, Phys. Rev. Lett., 125, 261104, doi: 10.1103/PhysRevLett.125.261104
  * Tews et al. (2021) Tews, I., Pang, P. T. H., Dietrich, T., et al. 2021, Astrophys. J. Lett., 908, L1, doi: 10.3847/2041-8213/abdaae
  * Tolman (1939) Tolman, R. C. 1939, Phys. Rev., 55, 364, doi: 10.1103/PhysRev.55.364
  * Tsokaros et al. (2020) Tsokaros, A., Ruiz, M., & Shapiro, S. L. 2020, Astrophys. J., 905, 48, doi: 10.3847/1538-4357/abc421
  * Typel et al. (2010) Typel, S., Ropke, G., Klahn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev. C, 81, 015803, doi: 10.1103/PhysRevC.81.015803
  * Typel & Wolter (1999) Typel, S., & Wolter, H. H. 1999, Nucl. Phys. A, 656, 331, doi: 10.1016/S0375-9474(99)00310-3
  * Watts et al. (2015) Watts, A., et al. 2015, PoS, AASKA14, 043, doi: 10.22323/1.215.0043
  * Yagi & Yunes (2016) Yagi, K., & Yunes, N. 2016, Class. Quant. Grav., 33, 13LT01, doi: 10.1088/0264-9381/33/13/13LT01
  * Zhang & Li (2020) Zhang, N.-B., & Li, B.-A. 2020, Astrophys. J., 902, 38, doi: 10.3847/1538-4357/abb470
  * Zurek (2014) Zurek, K. M. 2014, Phys. Rept., 537, 91, doi: 10.1016/j.physrep.2013.12.001