Constraining self-interacting fermionic dark matter in admixed neutron stars using multimessenger astronomy

We investigate the structure of admixed neutron stars with a regular hadronic component and a fraction of fermionic self-interacting dark matter. Using two limiting equations of state for the dense baryonic interior, constructed from piecewise generalised polytropes, and an asymmetric self-interacting fermionic dark component, we analyse different scenarios of admixed neutron stars depending on the mass of dark fermions $m_\chi$, interaction mediators $m_\phi$, and self-interacting strengths $g$. We find that the contribution of dark matter to the masses and radii of neutron stars leads to tension with mass estimates of the pulsar J0453+1559, the least massive neutron star, and with the constraints coming from the GW170817 event. We discuss the possibilities of constraining dark matter model parameters $g$ and $y \equiv m_\chi/m_\phi$, using current existing knowledge on neutron star estimations of mass, radius, and tidal deformability, along with the accepted cosmological dark matter freeze-out values and self-interaction cross-section to mass ratio, $\sigma_\mathrm{SI}/m_\chi$, fitted to explain Bullet, Abell, and dwarf galaxy cluster dynamics. By assuming the most restrictive upper limit, $\sigma_\mathrm{SI}/m_\chi<0.1$ cm$^2$/g, along with dark matter freeze-out range values, the allowed $g$-$y$ region is $0.01 \lesssim g \lesssim 0.1$, with $0.5 \lesssim y \lesssim 200$. For the first time, the combination of updated complementary restrictions is used to set constraints on self-interacting dark matter.


INTRODUCTION
In addition to ordinary matter, our Universe is believed to be populated by a dark component known as dark matter (DM).Within the standard Λ-CDM cosmological model, DM accounts for approximately 80% of the total matter content in the Universe.Over the past few decades, extensive searches have been conducted to find this elusive massive particle candidate with null results (Bertone & Hooper 2018).The analysis of multi-messenger events involving radiation (DeRocco et al. 2019), neutrinos (Rembiasz et al. 2018), and gravitational waves (GW) (Badurina et al. 2021) has been explored to uncover hints of a new particle sector beyond the Standard Model.
Numerous candidates have been proposed to explain the existence of DM, including weakly (strongly) interacting massive particles, WIMPs (SIMPs), and feebly interacting neutrinos (Datta et al. 2021).In recent years, significant attention has been given to exploring additional DM candidates, particularly those in the light or ultra-light mass sector, such as weakly interacting axions or axion-like particles, ALPs.Despite extensive searches, the absence of any detection has led to highly constrained parameter spaces (Schumann 2019).Axions, in particular, have gained prominence due to their potential ★ E-mail: mmariani@fcaglp.unlp.edu.ar as a solution to the charge-parity problem in Quantum Chromodynamics (QCD) (Peccei 2008).They have also been investigated for their potential impact on phenomena like cooling and the generation of broadband radio signals in pulsars (Prabhu 2021).Additionally, self-interacting DM (SIDM) has been proposed as a solution to the core-cusp problem, which refers to the discrepancy between the inferred density profiles of low-mass galaxies and the predictions of cosmological N-body simulations (Spergel & Steinhardt 2000).SIDM has been explored as a means to reconcile these inconsistencies.These examples represent only a fraction of the wide array of DM candidates that have emerged from extensions of the Standard Model.
The clumpy nature of DM in the Universe allows for the existence of regions such as globular clusters with enhanced DM mass densities,   ∈ [10, 10 5 ]   0 , being   0 ∼ 0.4 GeV/cm 3 the solar neighbourhood DM density (Read 2014;Cautun et al. 2020).Nevertheless, detecting DM particles directly is an extremely challenging task due to the fact that they rarely interact with ordinary matter.So far direct, indirect and collider searches have tried to set up different experimental strategies, in addition to gravitational effects, to discern a signal that reveals the fundamental nature of DM.However, an alternative approach lies in indirectly detecting DM through accumulation effects in astrophysical bodies and, in particular, in dense neutron stars (NSs).In this sense, astrophysical bodies characterised by mass, , and radius, , play a significant role in this endeavour through the compactness, defined as  ≡  /( 2 ), with NSs or quark stars being especially efficient at this task (Sandin & Ciarcelluti 2009;Mukhopadhyay & Schaffner-Bielich 2016).
Despite the extremely weak interactions between DM particles and ordinary matter, the high density of NSs enables them to effectively capture and retain DM particles that may pass through (Cermeño et al. 2016).Once trapped within the star, these particles undergo successive collisions with the nucleons in the star, gradually thermalizing and eventually accumulating in the core over a sufficiently long period of time (Singh et al. 2023).This accumulation would lead to several detectable signatures, such as the heating up of cold NSs (Bertone & Fairbairn 2008;Bell et al. 2018;Pérez-García et al. 2022b), pulsar scintillation (Pérez-García et al. 2013), or changes in the rotational patterns (Kouvaris & Pérez-García 2014), making them detectable in future observations.Combining the efforts in the electromagnetic bands, with infrared capabilities of James Webb Space Telescope, the proposed gamma-ray telescopes like e-ASTROGAM (Angelis et al. 2018) and AMEGO, the upcoming Square Kilometer Array in radio astronomy, ATHENA in X-ray astronomy, along with 3rd generation gravitational wave detectors such as advanced LIGO, Virgo, and KAGRA, Cosmic Explorer (Reitze et al. 2019) or Einstein Telescope (Branchesi et al. 2023), represents a promising way to detect potential signals from a plethora of DM induced phenomenology (Boddy et al. 2022).In extreme cases, the accumulation of DM could trigger the collapse of the star into a black hole (McDermott et al. 2012;Zurek 2014;Singh et al. 2023) or the conversion into a quark star (Herrero et al. 2019), leading to bursting radiation signals (Zenati et al. 2023).
Besides accumulation in the NS core due to the stellar gravitational potential well and interactions with ordinary matter, we foresee another scenario where DM particles can appear as final or intermediate products in decay or creation processes, such as massive sterile neutrinos or massive scalars (Albertus et al. 2015;Rembiasz et al. 2018;?;Fornal & Grinstein 2018).In either case, the evidence (or lack thereof) of DM accumulation in NSs can provide valuable clues on where to direct experimental efforts to detect this type of matter.
Several studies have investigated SIDM in both the weakly and the strongly interacting regimes within compact objects (Narain et al. 2006;Leung et al. 2011;Deliyergiyev et al. 2019;Das et al. 2019Das et al. , 2021;;Husain & Thomas 2021;Kain 2021;Das et al. 2022;Miao et al. 2022;Rutherford et al. 2023;Routaray et al. 2023;Diedrichs et al. 2023).In this work, we model NSs including a fermionic SIDM component using a generalised piecewise polytropic (GPP) equation of state (EoS) for hadronic matter adjusted to reproduce, in a reasonably accurate way, observables like mass, radius and moment of inertia that would be obtained with realistic hadronic EoSs (O'Boyle et al. 2020).SIDM is incorporated through the two-fluid formalism (Sandin & Ciarcelluti 2009).
The simplest model of SIDM particles includes a Yukawa-type potential and a force carrier  mediating between two DM particles , with mass   .The key parameter governing the likelihood of self-interactions among DM particles is the SI cross-section ( SI ) to mass ratio,  SI /  .Furthermore, the separation between the DM halo of the Abell galaxy cluster and its stars can be explained in terms of self-interaction cross-sections of DM (see, for example, Kahlhoefer et al. 2015, and references therein).Additionally, SIDM provides a consistent cross-section that matches the DM halo profile of dwarf galaxies.However, there are discrepancies between this fit and the results from galaxy merger studies.This is why the analysis of velocity-dependent self-interactions related to freeze-out, ⟨ ann  rel ⟩ -where  ann is the DM annihilation cross-section, and  rel is the relative velocity between the annihilating particles-, which avoids the constraints posed by galaxy clusters, becomes relevant in such cases (Hayashi et al. 2021).Furthermore, it is crucial that DM selfinteractions are not so strong as to disrupt the elliptical shape of spiral galaxies or displace sub-clusters, as seen, for instance, in the Bullet Cluster.Moreover, the latter can also be employed to establish restrictions on the SIDM cross-section (Robertson et al. 2016).All these constraints place limits on the strength of DM self-interactions; for more details, see Rocha et al. (2013).
On the other hand, astronomical observations of NS received a boost in the past two decades.Detection in double pulsar systems -PSR J1614-2230 (Demorest et al. 2010), PSR J0348+0432 (Antoniadis et al. 2013), and PSR J0740+6620 (Cromartie et al. 2020;Fonseca et al. 2021)-of the ∼ 2  ⊙ NSs requires an acceptable EoS able to support such high masses.These observations posed the first strong restrictions to the behaviour of matter inside such compact stars.Moreover, multimessenger astronomy with GWs detected from the binary NS merger event GW170817 and its electromagnetic counterpart allowed to restrict the value of the dimensionless tidal deformability, Λ, of NSs (Abbott et al. 2017;Annala et al. 2018;Most et al. 2018;Raithel et al. 2018;Abbott et al. 2018;Capano et al. 2020) and ejecta properties (Pérez-García et al. 2022a).Additionally, the Neutron Star Interior Composition Explorer (NICER) has measured the mass and radius of the isolated millisecond-pulsars PSR J0030+0451 (Riley et al. 2019;Miller et al. 2019) and (together with XMM-Newton data) PSR J0740+6620 (Riley et al. 2021;Miller et al. 2021) with great precision.The latter study showed that despite having a mass ∼ 40 % larger, the radii of PSR J0740+6620 and PSR J0030+0451 are of the same order (Riley et al. 2021;Miller et al. 2021).
Our goal is to analyse how the current astrophysical constraints associated with NSs and cluster dynamics can set limits on the mass of the DM particle,   , and on parameters of the SIDM model, such as the self-interaction (SI) coupling strength constant, , or the mass scale of the interaction or, equivalently, the associated generic mediator mass,   .Furthermore, for the first time, constraints coming from multi-messenger astronomy of NSs are used, combined with the restrictions of SI cross-section from galaxy clusters, dwarf galaxies and the thermally averaged annihilation cross-section related to the cosmological DM freeze-out value, to set constraints to SIDM.
The paper is organised in the following manner.Section 2 is devoted to describing the theoretical framework used in this paper.In particular, Subsection 2.1 provides the description of both the hadronic EoS and fermionic SIDM EoS used in this study.Subsection 2.2 offers a brief review of the two-fluid formalism employed to investigate NSs with a DM component using the two aforementioned EoS.Our main findings are presented in Section 3, while Section 4 contains a summary and a discussion of the astrophysical implications of our results.Unless stated otherwise, we use units where  =  = ℏ = 1 throughout the paper.

Generic hadronic matter and fermionic SIDM EoS
First, let us discuss the description of the hadronic content in NSs.To ensure that our conclusions are not dependent on any specific EoS, we constructed two different hadronic EoSs using the GPP formalism presented in O' Boyle et al. (2020).
These two EoSs, referred to as the soft and stiff, are designed so that the resulting mass-radius relationships encompass a family of EoSs consistent with chiral effective field theory (EFT) up to a mass density of  = 1.1  nuc , being  nuc = 0.16 fm −3 the nuclear saturation density (Hebeler et al. 2013;Annala et al. 2020).In both cases, the crust is described by a GPP fit to the SLy4 crust EoS, presented in O' Boyle et al. (2020).Within these parameterised EoSs, the pressure (), energy density  (), and speed of sound   () are continuous functions.This is a crucial aspect of the GPP formalism as the dimensionless tidal deformability explicitly depends on the value of the speed of sound (for more details, see, for example, Leung et al. 2022, and references therein).
For each interval in mass density range [ −1 ,   ] the EoS adopts a power-law form where the parameters   , Γ  ,   , Λ  characterise the fit and ensure continuity and differentiability of both () and  () at the dividing densities  0 ,  1 and  2 , where  0 <  1 <  2 .This leads to continuous speed of sound,  s .By imposing subsequent mathematical relations among these quantities, we obtain seven parameters used to construct each hadronic EoS, as presented in Table 1.
The soft and the stiff EoSs act as an envelope for numerous microscopic hadronic EoSs in the literature that satisfy current astrophysical constraints on NSs.This approach has recently been used to study both macroscopic and oscillating properties of compact objects independently of any particular theoretical model (see, for example, Ranea-Sandoval et al. 2022;Gonçalves & Lazzari 2022;Saes & Mendes 2022;Lugones et al. 2023;Lenzi et al. 2023;Ranea-Sandoval et al. 2023a;Ranea-Sandoval et al. 2023b).
In addition to considering hadronic matter, we also incorporate an additional fraction of massive fermionic DM particles that can inhabit the NS.This arises in our scenario as a result of the dark component being gravitationally bound to the compact star.
The capability of a NS to capture DM from an external distribution is mainly determined by its gravitational pull, opacity, and the kinematics of incoming particles from the surrounding dark environment (Press & Spergel 1985;Gould 1987).It is important to note that the DM capture rates in NSs have been calculated based on scattering off nucleons () (Kouvaris 2008).Typically, the capture rate   can be approximated as follows -see Eq. (1) in ?-, where   denotes the probability of at least one scattering event between a DM particle () and a  taking place within the NS.At this point, it is worth emphasising that in our treatment the - interaction is considered secondary with respect to gravitational effects.The previous expression considers NS mass and radius values, along with General Relativity corrections, for typical benchmark compact stars with masses around  ∼ 1.4 ⊙ and radii of approximately  ∼ 10 km.This is accurate to within factors of order unity, and any associated uncertainties are inherently included in the fraction of DM populating the star.
Setting a minimum upper limit on the DM-nucleon cross-section,   ≳ 10 −46 cm 2 , results in   ∼ 1, indicating that the NS saturates its capability to capture and bind DM.Therefore, the amount of DM inside the NS is a tiny fraction compared to the baryonic number, which is of the order of 10 58 .Given the current rates and the fact that the oldest NS lifetimes are  ∼ 10 9 yr, it is likely that an additional mechanism is necessary to incorporate DM to the few-per-cent level, as assumed in most works, including ours.Now, using a simplified scattering model, we will consider the contributing terms for the   →   process.We will assume that the SI terms are described by generic Lorentz scalar  and vector   mediator fields.Explicitly, the Lagrangian density for DM interaction involving , DM anti-particles χ and ,   fields is written as where   and   are the scalar and vector coupling strengths, respectively, and   are the Dirac matrices.Generically, we assume a bosonic mediator, , which we can take to be either a Lorentz scalar or a vector to illustrate, with mass scale   ∼   , where   is the interaction mass scale.In its simplest form, this scenario assumes a unified fine-structure constant   ∼  , (where  , are, respectively, the scalar and vector fine structure constants) and   ≪   .
Since we consider asymmetric DM coupled to light , , force mediators, if they are sufficiently light, then the interaction between DM particles becomes long-range.More specifically, for an interaction described in the non-relativistic regime by a sort of Yukawa potential, V , = − ,  − ,  /, where  , =  2 , /4, long-range effects occur if the mediator mass is smaller than the Bohr momentum,  , ≲  ,   /2.
Following the work of Mukhopadhyay & Schaffner-Bielich ( 2016), we model SIDM as a Fermi gas of SI particles with mass   .The explicit expressions for the fermionic SIDM EoS are given by Mukhopadhyay & Schaffner-Bielich (2016), In order to handle dimensionless quantities, we have defined   =  , /  as the ratio of DM Fermi momentum to the DM particle mass and a strength parameter  =   /  , the quotient between the DM particle mass and the interaction mass scale.The terms that are proportional to  2 account for the SI effects, which are included analogously to the vector interaction terms in Nambu-Jona-Lasinotype models for quark matter (Orsaria, M. and Rodrigues, H. and Weber, F. and Contrera, G. A. 2014).In our fermionic SIDM model, we can approximate the self-scattering cross-section  SI using dimensional arguments (Girmohanta & Shrock 2022).This allows us to express the self-interaction terms in Eq. ( 6) in the following form Since we do not expect large values of accreted or produced DM in the interior of the NS and assuming a two-body self-interaction for simplicity, we note that the contribution to the energy density or pressure is proportional to the square of the scalar density  , or the vector density   .Note that   =   /  is the  number density, and in this context, For   =  ,  /  → 0, we find that  ,  can be approximated by Thus, in our scenario,  ,  ≃   .These scalar and vector particle number densities contribute with the same strength to   and   , affecting the DM EoS.
In the non-relativistic limit, the scalar SI cross-section for   →   scattering can be expressed as Therefore, the ratio of the SI cross-section to the mass of the  particle becomes where we have introduced the notation   =   /  .It is worth noting that in the non-relativistic limit, a vector interaction yields the same expression for the SI cross-section to mass ratio, i.e.,  ,→ ∼   , → .Henceforth, we will consider a generic coupling constant   =   =  and impose a single mass scale   =   =   =   to constrain the product  or equivalently /  , based on the combined properties of NS mass-radius, dimensionless tidal deformability, Λ, thermal DM freeze-out value and the SI cross-section obtained from galactic dynamics.
In the current DM paradigm, to achieve a finite DM relic density, it is necessary to consider the self-annihilation processes involving χ fields.This introduces a correction to the scattering SI cross-section, which can be expressed in a generic form as Therefore, it is important to assess the extent to which the DM in our SIDM model is consistent with the observed DM relic abundance, thereby determining the value of thermally averaged cross-section ⟨ ann  rel ⟩.The commonly used canonical value for a generic weakly interacting DM candidate is typically stated as ⟨ ann  rel ⟩ ∼ 3 × 10 −26 cm 3 s −1 , with unspecified uncertainty, and is assumed to be independent of the  mass.Recent studies on the search for annihilation products of DM suggest that 2.2 × 10 −26 cm 3 /s ≲ ⟨ ann  rel ⟩ ≲ 5.2 × 10 −26 cm 3 /s, with a weak dependence on   > 10 GeV (Steigman et al. 2012).The irreducible annihilation channel   →  and tree-level cross-section reads It should be noted that this expression can be modified by O (1) pre-factors, which we will neglect at this point, depending on factors such as whether  is a Majorana or a Dirac fermion and whether  is a scalar or a gauge boson.Additionally, even in simple models, ⟨ ann  rel ⟩ typically receives contributions from other annihilation the -halo branch corresponds to the configurations with a DM halo, where   >  mat , while the -core branch includes those with a DM core, where   <  mat .The limiting configuration between branches is indicated with a circle, where  eq =   =  mat .In the cases in which  eq exists,  eq =  ( eq ).The coloured bars and contoured clouds correspond to astrophysical constraints for NSs (see the text for details).channels.In our scenario, this is described by the process χ  →  .At the lowest order, it can be approximated as where  ann =  χ→  explicitly.Furthermore, in the context of the process  → χ , the decay width is given by where  ≡ .To ensure perturbative consistency, it is necessary to verify that Γ → χ /  ∼  ≪ 1.By considering  ∼ 0.1, we find that  ≲ 1.8, allowing for a positive and loosely constrained coupling strength.Note that for  > 1/2 this bound is still valid as derived in Peskin & Schroeder (1995).

Stellar configurations in the two fluid formalism
Considering that DM interacts only weakly with ordinary matter, it is possible to obtain equilibrium configurations of compact objects formed by an admixture of ordinary and DM in the so-called twofluid formalism in which hadronic (mat) and dark () components interact only gravitationally (see Mukhopadhyay & Schaffner-Bielich (2016) and references therein).Thus, in this model, the pressure and energy density for each type of matter are assumed to be essentially decoupled, so that pressure can be written as follows and energy density can be expressed as Since in this formalism, the DM-ordinary matter interaction is assumed to be negligible, and both components interact only gravitationally, the classical Tolman-Oppenheimer-Volkoff (TOV) equations are replaced by the following four coupled differential equations for pressure and gravitational mass ) where  = mat, .Note that, in this description, the total gravitational mass is defined as  () =  mat () +   ().The admixed TOV equations are simultaneously solved numerically using the prescribed EoSs of Table 1 for ordinary matter and the DM EoS described in Section 2. Along with the specified EoSs, it is necessary to make explicit the boundary conditions at the centre of the star,  = 0, setting  mat (0) =   (0) = 0. We also need to specify the central pressures of the two fluids.For this purpose, we define the central pressure of the matter,  mat (0), and the DM fraction,   , which is given by from which the central value of the dark fluid pressure,   (0), can be determined.The radius of the stellar configuration is defined as  = max( mat ,   ) where the radius of each component is determined when the condition   (  ) = 0 is satisfied.Due to the monotonicity of the radial coordinate, a stellar configuration could exist on the stable curve of the admixed TOV solutions in which dark and ordinary matter have the same size, i.e.,  eq =   =  mat .For con-figurations with  >  eq a DM halo forms, and for  <  eq there is a DM core.In such case, if there is a stable configuration with  =  eq , its gravitational mass can be defined as  eq =  ( eq ).

RESULTS
To obtain the main results of this work, some of the SIDM EoS space parameters are considered within the following ranges.The DM particle mass is explored in the region 100 MeV ≤   ≤ 10 4 MeV ; while the mediator mass, whose value for our fermionic SIDM model is not determined, lies in the interval Within these parameter ranges, we construct a variety of admixed NSs using the previously mentioned soft and stiff hadronic EoSs.The results for the mass-radius plane, dimensionless tidal deformability, second Love number  2 , and DM parameter spaces   -  and - are presented in the following subsections.

Mass-radius curve for admixed NSs with fermionic SIDM
First, in Fig. 1, we present a schematic mass-radius relationship, illustrating an example of an admixed NS with a DM component.In the diagram, we define several relevant quantities: the -halo branch, corresponding to configurations with a DM halo,   >  mat ; the -core branch, corresponding to configurations with a DM core,   <  mat , and the limiting stellar configuration between these branches, denoted by  eq =   =  mat and  eq =  ( eq ).While the -core branch generally exhibits the typical morphology of hadronic NSs, the halo branch deviates from this shape, showing a sudden change of / beyond  eq toward larger radii.In this figure, we also present in detail the current astrophysical constraints  2021).The cyan bar is associated with the least massive neutron star in a double pulsar system known to date, J0453+1559, with a mass of  = 1.174 ± 0.004 ⊙ (Martinez et al. 2015).The coloured contoured clouds correspond to  and  determinations by GW or X-ray observations.The pink and purple clouds correspond to constraints derived from the observations of the GW170817 event, at 90% confidence (Abbott et al. 2017(Abbott et al. , 2018)); the dark and light blue clouds correspond to constraints from the GW190425 event, also at 90% confidence (Abbott et al. 2020).Both GW events are displayed for the low-spin prior scenario.The brown (Miller et al. 2019) and yellow (Riley et al. 2019) clouds correspond to the restrictions derived from NICER observations of the pulsar J0030+0451.The dark green (Miller et al. 2021) and light green (Riley et al. 2021) clouds represent the constraints obtained from the joint observation of PSR J0740+6620 by NICER and XMM-Newton.For both PSR J0030+0451 and PSR J0740+6620, two contoured clouds are presented (labelled in the figure as 1 and 2) since both observations were analysed and reduced by two independent research groups, yielding different results.In these last two constraints, the outer (inner) contours correspond to 95% (68%) confidence level.All current astrophysical constraints detailed above should be satisfied by admixed NSs.
The results for the mass-radius of the admixed NSs are shown in Fig. 2 -for the soft hadronic EoS-and in Fig. 3 -for the stiff hadronic EoS-.In both figures, panels () displays all the studied EoS -considering the entire mentioned ranges of the   ,   and   parameters-, while panels () show only the results satisfying all the astrophysical constraints that we will now detail.The impact of including DM in the stellar configurations has been analysed by imposing the three most relevant astrophysical constraints as filters to obtain our results: the lower limit of the most massive double pulsar measured to date  max ≳ 2.01  ⊙ (light red horizontal bar) (Fonseca et al. 2021), the upper limit to the mass of PSR J0453+1559, the less massive NS in a binary system measured to date,  min < 1.178  ⊙ (cyan horizontal bar) (Martinez et al. 2015), and the data coming from the binary NS merger associated to the GW170817 event (purple and pink clouds) (Abbott et al. 2017(Abbott et al. , 2018)).Our results reveal that these three constraints are the most restrictive ones.Once these three are met, all other current astrophysical constraints for NSs are satisfied; in this sense, the other constraints are encompassed by these three filters.After applying the constraint filters, we discarded all the EoS either having too massive -halo branch, producing  min > 1.178  ⊙ or laying over the GW170817 clouds, or yielded an excessive amount of DM, resulting in  max < 2.01  ⊙ .The colour of the curves in panel (a) represents the   / ratio for each stellar configuration, while panel (b) displays the central DM particle number density,    .
Although the wide array of curves displays tangled behaviour, each panel presents a variety of curves that show different features in our model results.In each panel we observe the existence of distinct -core branches, reminiscent of typical purely hadronic mass-radius curves, but with varying arc length extensions.Thus, as the -mass component increases, the maximum mass and radii decrease.On the other hand, the most notable effect of the DM population is the emergence of the -halo branches that detach from the -core branches.These -halo branches exhibit a more horizontal trend towards larger radii with minimal mass variation.The value of  eq where these branches join also depends strongly on the model parameters and we will discuss these dependencies in Subsection 3.3.Nevertheless, it worth be noting here that the values of  eq do not appear to correlate with the   / ratio, as there is a significant variation in the  eq values unrelated to the increase or decrease of   /.In particu-   2. For completeness, the black curves correspond to purely hadronic EoSs.While the similarities or differences of the pure -core configurations with the hadronic EoS depend on the particular selection of the parameters, the curves with significant -halo branches present an important deviation from the hadronic curve.lar, for the soft model cases, an increase in the   / leads to a decrease in the maximum mass value, resulting in the exclusion of most of the curves with larger values of   / when the filter associated with the most restrictive mass constraint  max ≥ 2.01  ⊙ is applied (Fonseca et al. 2021), panel () of Fig. 2.Although there are some stellar configurations for some particular EoS that contain   / ≳ 0.1, the majority of EoSs, as well as all stellar configurations with astrophysical interest, that is, with  min ≳ 1  ⊙ , have   / < 0.06.As mentioned earlier, the increasing of DM in the -core objects produces a smaller radius.This effect implies that larger fractions of DM could lead to extremely stiff hadronic EoSs that satisfy the constraint related to GW170817 event.Considering the opposite case, soft models results suggest that for high fractions of DM, some baryonic EoSs could be excluded since they would fail to satisfy the GW190425 constraint (dark-blue and light-blue clouds in the figures).
Regarding the behaviour of    in the panel () of both figures, it is important to emphasise that we integrate the admixed TOV equations while considering the relationship between the central values of both the hadronic and DM EoS through the parameter   .Consequently, it is expected that configurations with higher mass will be obtained as the central density    -along with  mat  -increases, similar to what occurs in the pure hadronic scenario.In both the soft and stiff models, there are a few particular EoSs that suggest that the GW170817 event could potentially be produced by -halo objects.

Effect of fermionic SIDM on NS tidal polarisability
It is known that tidal effects in the late inspiralling phase of a binary NS merger could be detected by GW detectors.To describe this phenomenology the individual NS tidal deformability is defined as Λ = 2 3  2  −5 with  2 the second Love number (Hinderer 2008).In Fig. 4, we present the Love number  2 as a function of the compactness, , for soft, panel (), and stiff, panel (b), hadronic EoSs, using the selection of Table 2.We also include the results for purely hadronic EoS for comparison.We selected four particular admixed EoSs to consider and compare combinations of the soft and the stiff scenarios with two of them (EoS 1 and 3) having no -halo branch, and the other two (EoS 2 and 4) featuring a substantial presence of a -halo branch.
Table 2. EoSs chosen to analyse the behaviour of the second Love number,  2 , as a function of the compactness of admixed NSs.EoS 1 and 3 have an absence of -halo branch, while EoS 2 and 4 have a significant presence of a -halo branch.
As observed in Fig. 4, the cases with -halo branches exhibit significant differences from the purely hadronic curve, consistently being smaller.However, the extent of the differences between the latter and the cases where the -halo branch is absent depends on the specific admixed EoS chosen.This suggests that the study of DM through the determination of  2 is influenced by the underlying hadronic models, although certain trends can still be identified.This becomes evident when comparing the results of the soft EoS 1 and the purely hadronic EoS, as they are indistinguishable from each other, as is shown in panel ().Furthermore, this characteristic appears to be independent of the DM fraction   .In our setting, the SIDM model is more sensitive to smaller values of DM particles in the region where  ≲ 0.15.
On the other hand, the constraints from GW170817 event are related to the total gravitational mass of the binary system,  ≈ 2.74 ⊙ , and the individual NS gravitational masses  1 ≈ (1.36 − 1.6)  ⊙ and  2 ≈ (1.17 − 1.36)  ⊙ .An analysis was performed for this event, which was marginalised over the selection methods (Abbott et al. 2018), leading to the discovery of an upper bound for the effective dimensionless tidal deformability of the binary, Λ ≤ 800 at a 90% confidence level.This analysis assumed a low-spin prior, which disfavours EoSs predicting larger radii for stars.These findings arise due to the individual Λ 1 , Λ 2 that are inbuilt in the actual definition of Λ( 1 ,  2 , Λ 1 , Λ 2 ) (Abbott et al. 2018).Tighter constraints were found in a follow-up reanalysis (Abbott et al. 2019) obtaining Λ = 300 +420 −230 (using the 90% highest posterior density interval), under minimal assumptions about the nature of the compact star binary system.In our work, we use a more refined value of the tidal deformability of a 1.4  ⊙ NS obtained from  and 3.The colour of the curves represents the ratio between the DM mass and the total gravitational mass,   /.The nearly vertical branches in the curves with a higher value of Λ correspond to configurations where a DM halo exists,   >  mat .The branches displaying a traditional hadronic morphology correspond to configurations with a DM core,   <  mat .The green vertical segment corresponds to the constraint imposed for a 1.4  ⊙ NS based on the analysis of the GW170817 event (Abbott et al. 2018).Similar to the effect of   / observed on Figs. 2 and 3, as shown in panel (), the inclusion of DM component may contribute to satisfying the dimensionless tidal deformability constraint for some of the otherwise discarded EoSs.Abbott et al. (2018) and estimated to be Λ 1.4 = 190 +390 −120 at a 90% credible level when a common EoS is imposed.Note that this value is actually more constraining than the one obtained when assuming a binary NS merger involving similar merging NSs, in which case Λ = Λ.
The results for the Λ-mass plane are presented in Fig. 5, for the soft, panel () and stiff, panel (b), hadronic EoS.In these figures, we only show the filtered family of stars satisfying the astrophysical constraints, corresponding to panels () of Figs. 2 and 3.The colour of the curves represents the ratio between the DM mass and the total gravitational mass of the respective star, denoted by   /.It can be observed that soft models satisfy the GW170817 constraint (indicated by the green segment) regardless of the DM parameters.For stiff models, configurations with a small proportion of DM, resembling the purely hadronic configuration, barely satisfy this constraint.Moreover, as the DM ratio increases, the GW170817 constraint is easily satisfied, consistent with the behaviour observed in the massradius plane, where the radius decreases with increasing DM ratio.

Analysis and constraints on the SIDM-model parameter space
After integrating the two-fluid TOV equations for the admixed EoSs, we study how applying the astrophysical filters mentioned above constrains the parameter space of SIDM.In Figs. 6 and 7 for the soft and stiff hadronic EoSs, respectively, we present the plane   -  for different values of   = 0.02, 0.05, 0.1, within the allowed mass ranges for   and   .The circles and diamonds indicate the sampling sets studied.The circle (diamond) indicate the sets that do (do not) fulfil the three restrictive astrophysical constraints considered.Each diamond is composed of three triangles representing the minimum mass 1.178  ⊙ limit (left triangle), the maximum mass 2.01  ⊙ (right triangle) limit, and the GW170817 restriction (bottom triangle).The colour green (red) indicates if the corresponding con-straint is (is not) satisfied.Although related to specific astronomical observations, these filters can be revealing since they are qualitative indicators of the DM effects on the maximum and minimum mass and on the radius of admixed NSs.The circles in Figs. 6 and 7 indicate parameter sets fulfilling the three constraints.The colour bar indicates the value of the corresponding  eq for each EoS set.To simplify the classification and visualisation of the results, we present with dark violet, assigned to  eq = 0, the cases in which no halo branch appears in the mass-radius curve of stable solutions, meaning that  eq does not exist.Panels (a), () and () of Figs. 6 and 7 show the upper left corner excluded according to astrophysical constraints on NSs.This zone corresponds to weakly SIDM or low  =   /  values.Depending on the hadronic EoS and the fraction   chosen, the excluded region corresponds to 0.2 <  < 0.6 (stiff EoS and   = 0.02) or 0.2 <  < 10 (soft EoS and   = 0.1).In the colour map area, it can be observed that most of the regions (dark purple zone) correspond to EoS with predominantly -core branches, as  eq becomes remarkably small or absent.Some sets in the intermediate   -low   region (corresponding to the colours greenish and yellow in the colour map) produce EoS with a non-negligible -halo branch contribution.As we already noted in the mass-radius relationships (Figs. 2 and 3), in the most extreme cases, the  eq reach high enough values that one or both objects involved in the merger that produced the GW170817 event could have been -halo objects.Although the sampled values of   differ by approximately a factor 5, it can be seen in the soft case, that a larger value for this parameter prevents compliance of the more restrictive constraint of 2.01  ⊙ , producing less massive stellar configurations.On the other hand, due to the higher masses generally achieved by stiff EoSs, our stiff models fulfil the maximum mass constraint for almost all studied sets and therefore exclude fewer sets than the soft models in this high   -low   region.
Finally, within the framework of the DM paradigm, in Fig. 8 (soft hadronic EoS) and Fig. 9 (stiff hadronic EoS) we study the  SI generic coupling constant,  -interaction parameter, , plane, considering the constraints on  SI /  and ⟨ ann  rel ⟩, and the upper bound of g, given by Eq. ( 13).For  SI /  , we apply  SI /  ≲ 0.1, 1.0, 10 cm 2 /g corresponding to order-of-magnitude values coming from Bullet-type massive clusters (Robertson et al. 2016), the Abell cluster galaxies (Kahlhoefer et al. 2015) and dwarf galaxies (Hayashi et al. 2021), respectively; the thermally averaged self-annihilation cross-section from the cosmological DM freeze-out is taken to be 2.2 × 10 −26 cm 3 /s ≲ ⟨ ann  rel ⟩ ≲ 5.2 × 10 −26 cm 3 /s (Steigman et al. 2012).In both figures, panels (), () and () correspond to   = 0.02, 0.05, 0.1, respectively.The dark blue contoured region in the - plane is constructed considering the  SI /  in the prescribed uncertainty range [0.1, 10] cm 2 /g, and using the constraints shown in Fig. 6 and 7, for the mediator mass, the dark fermion mass, considering the restrictions coming from the multimessenger astronomy of NSs.We colour with light blue the region where  SI /  < 0.1 cm 3 /s in accordance with the Bullet cluster constraint.In addition, we present in light pink the constraint for ⟨ ann  rel ⟩, also considering the restrictions in   and   coming from NS observations.In green colour we present the region  ≲ 1.8, arising from Eq. ( 13).
In order to establish these constraints over  and  magnitudes, we determine the intersection between the most restrictive limit of   /  < 0.1 and the pink area given by the ⟨ ann  rel ⟩ restriction.Although the coloured areas are slightly larger in Fig. 9 than in Fig. 8, it is possible to restrict the strength parameter  in the range 2 ≲  ≲ 200, while the generic coupling constant is allowed to vary between 0.01 ≲  ≲ 0.1.When considering results for the stiff hadronic EoS, Fig. 9, allowed  values are shifted to the lower side with the lowest  ∼ 0.5 for the   = 0.02 case.we have used two extreme cases of GPP EoS i.e. a soft and a stiff EoS compatible with both low-density chiral EFT calculations and modern NS astronomical observations.It is important to emphasise that the results presented here are intended to be representative of a wide variety of EoS that fall between these two extreme cases.A different choice of soft and stiff EoS would not qualitatively change the conclusions drawn from our work.The DM EoS was modelled considering a Fermi gas of self-interacting massive particles.Different families of admixed NSs considering a fixed fraction of DM were built through the integration of the TOV equations within the two-fluid formalism.
We have built the mass-radius and tidal deformability-mass relationships for the different families of stellar configurations with a dark component obtained.The DM parameters associated with the SIDM model considered were varied in a wide range of values to obtain the most general and representative results possible for our model.From the different sets of admixed EoS, we have set a filter to select only those admixed NSs satisfying the current NS constraints.This filtering process allowed us not only to study the different possible NS configurations but also to establish updated constraints over the SIDM parameters.
The results on the mass-radius plane suggest that the effect of the SIDM on admixed NSs could be quite relevant, producing long branches of halo-type objects and/or considerably decreasing the mass and radius of these compact objects.However, when the current astrophysical restrictions of NSs are applied as filters, much of the studied EoSs should be discarded, so these constraints could effectively help improve the DM models.In the cases where the available constraints from NS observations are fulfilled, the effect of fermionic SIDM is more evident for compact low-mass objects producing branches of configurations with approximately constant gravitational mass and a variable radius which is larger than the purely baryonic NS would have, the so-called -halo branch.In particular, our results do not discard the possibility that the GW170817 event could be produced by at least one -halo object.
Interestingly, in some cases, we found that the inclusion of a DM component displaces the radius of the obtained configurations towards lower values and contributes, in the case of stiff EoS, to satisfy the constraint of the GW170817 event.On the contrary, in the case of the soft hadronic EoS, this change in radius could cause the curves to fall outside of the GW190425 restriction.This behaviour is also evident in the tidal deformability-mass plane, where the DM contribution leads to a decrease in Λ and, thus, helps satisfy the GW170817 constraint in this plane in the stiff scenario.
In addition, we have explored the impact of the presence of SIDM on the NS structure, looking in particular at the change in the Love number  2 which can be estimated indirectly from GW data emitted during binary NS mergers.We find that low compactness admixed NSs present the best scenario for testing the presence of fermionic SIDM in NSs.
The characterisation and classification of the -halo and -core branches based on the astrophysical constraints allow us to set the most restrictive constraints when studying the DM effects on NSs.In addition to the thoroughly studied maximum mass constraint, the impact of DM on the branches of -halo and on the value of  eq could potentially be in tension with the minimum mass constraint and the mass and radius constraints coming from GW170817.In this regard, future studies, in the context of potential future observations, should not lose sight of this possible behaviour of admixed NSs.
Considering all the admixed EoSs constructed along with the current astrophysical constraints, we have studied and constrained the parameter space of our fermionic SIDM model.We found that weakly interacting fermionic SIDM -corresponding to low values for   and high values for   -i.e., low values for -is excluded for admixed NSs.This effect occurs due to a excessively massive  eq -thereby failing to satisfy  min < 1.178  ⊙ or GW170817 -or having a maximum mass  max that is too small, less than 2.01  ⊙ .In more detail, we have seen that the SIDM EoS constraints cover a wider range in the   −   plane if a soft hadronic EoS is considered.In this scenario,   ≳ 10 MeV is excluded for   ≲ 10 3 MeV.However, when considering the stiff hadronic EoS, this situation undergoes a quantitative change.In this scenario, the constraints coming from the observations of NSs are less stringent, again ruling out DM candidates with masses   of a few hundred MeV if   ≳ 10 MeV.Notably, in this case, the results are almost independent of the DM fraction,   .
The constraints on the values of DM SI cross-sections, obtained from the dynamics of Bullet and Abell clusters, dwarf galaxies, and cosmological DM freeze-out intervals, have been combined with multi-messenger NS astronomy to constrain fermionic SIDM.This wealth of information helps us constrain the - parameter space effectively.Specifically, we find that, when considering the most restrictive upper limit condition,  SI /  < 0.1 cm 3 /s, along with the range of DM freeze values for ⟨ ann  rel ⟩, the allowed - region falls within 0.01 ≲  ≲ 0.1 and 0.5 ≲  ≲ 200.
In a recent study by Shirke et al. (2023), similar SIDM constraints were derived, focusing solely on the DM vector interaction strength and based on only one of the conditions we have imposed, namely, the maximum mass values of NSs.Additionally, they worked within the traditional single-fluid TOV formalism.The updated constraints we present here, considering NS observations and SIDM cross-section restrictions simultaneously, can provide complementary and tighter constraints.
Finally, our results indicate that measuring the mass-radius of a low mass NS potentially in the -halo branch or detecting a GW NS merger event with determinations of Λ and  2 could offer insights into the presence of DM in NSs.As already mentioned, there are other observable quantities that could also indicate DM presence in compacts objects, such as gravitational waves from NS oscillations (Shirke et al. 2023), explosive kilonovae events (Bramante et al. 2014) or X-ray pulse profiles (Miao et al. 2022).We anticipate that this diverse range of observable predictions, coupled with the wealth of data expected from future experimental and observational projects, will contribute significantly to advancing our understanding of both the fundamental nature of DM and the internal structure and composition of NSs.

Figure 1 .
Figure1.Schematic mass-radius relationship for a typical admixed NS with a SIDM component.We indicate over the curve the two possible branches: the -halo branch corresponds to the configurations with a DM halo, where   >  mat , while the -core branch includes those with a DM core, where   <  mat .The limiting configuration between branches is indicated with a circle, where  eq =   =  mat .In the cases in which  eq exists,  eq =  ( eq ).The coloured bars and contoured clouds correspond to astrophysical constraints for NSs (see the text for details).

Figure 2 .
Figure 2. Mass-radius relationship of different families of soft admixed NSs.Panel () shows the results with all the considered EoSs.Panel () displays the filtered EoS that fulfil the astrophysical constraints.The colour on the curves represents the value of the ratio   / (left panel) and of the central  number density,   (right panel) for each star.In both panels, it can be seen the diverse morphology of the -halo branches and the variation of the  eq .In panel (a), it is evident that increasing   / reduces the maximum mass and radii in the -core branches.In panel (b), similar to the effect of increasing the central baryon density, as central   increases, the mass also increases.The coloured bars and contoured clouds represent the current astrophysical constraints; see Fig.1and text for more details.The inclusion of the SIDM in the NS models may satisfy these constraints or result in the exclusion of certain baryonic EoSs (see discussion in the text).Note that solutions in this plot, by construction, do always include a finite fraction of DM.
. The bars correspond to mass constraints derived by different pulsar observations: solid and hatched red bars correspond, respectively, to pulsars J1614-2230 and J0348+0432, originally measured by Antoniadis et al. (2013) and Demorest et al. (2010); the constraint for pulsar J1614-2230 was subsequently improved in 2018 by Arzoumanian et al. (2018).The light red solid bar corresponds to the restriction in the mass of the pulsar J0740+6620 (Cromartie et al. 2020), later refined by Fonseca et al. (

Figure 4 .
Figure 4. Love number,  2 , as a function of the compactness for the soft, panel () and stiff, panel () hadronic EoSs of the Table2.For completeness, the black curves correspond to purely hadronic EoSs.While the similarities or differences of the pure -core configurations with the hadronic EoS depend on the particular selection of the parameters, the curves with significant -halo branches present an important deviation from the hadronic curve.

Figure 5 .
Figure 5. Λ-mass plane for the soft, panel () and stiff, panel (b), admixed NSs satisfying the astrophysical NSs constraints, panels () of Figs.2and 3. The colour of the curves represents the ratio between the DM mass and the total gravitational mass,   /.The nearly vertical branches in the curves with a higher value of Λ correspond to configurations where a DM halo exists,   >  mat .The branches displaying a traditional hadronic morphology correspond to configurations with a DM core,   <  mat .The green vertical segment corresponds to the constraint imposed for a 1.4  ⊙ NS based on the analysis of the GW170817 event(Abbott et al. 2018).Similar to the effect of   / observed on Figs.2 and 3, as shown in panel (), the inclusion of DM component may contribute to satisfying the dimensionless tidal deformability constraint for some of the otherwise discarded EoSs.

Figure 6 .
Figure 6.Plane of   -  for soft hadronic EoSs and for different values of   = 0.02, 0.05, 0.1 -panels (), () and (), respectively-.The circles and diamonds indicate different EoSs within the ranges of the SIDM parameters chosen.In each panel, circles (diamonds) indicate the stellar configurations that do (do not) fulfil the three main astrophysical constraints considered (see text for details).The colour map in the circle regions indicates the value of the mass  eq for the stellar configuration with   =  mat (see details in the text).The three-sided diamonds indicate which of the three constraints is (in green) or is not (in red) satisfied.More specifically, the right and left sides of each diamond correspond to the NS maximum and minimum mass restrictions, respectively; the lower side of each diamond corresponds to the limitations imposed by GW170817.

Table 1 .
Parameters of the selected hadronic EoSs constructed with the prescription described in O'Boyle et al. (2020).