Rapid neutron star cooling triggered by dark matter

We study the effect of asymmetric fermionic dark matter (DM) on the thermal evolution of neutron stars (NSs). No interaction between DM and baryonic matter is assumed, except the gravitational one. Using the two-fluid formalism, we show that DM accumulated in the core of a star pulls inwards the outer baryonic layers of the star, increasing the baryonic density in the NS core. As a result, it significantly affects the star's thermal evolution by triggering an early onset of the direct Urca process and modifying the photon emission from the surface caused by the decrease of the radius. Thus, due to the gravitational pull of DM, the direct Urca process becomes kinematically allowed for stars with lower masses. Based on these results, we discuss the importance of NS observations at different distances from the Galactic center. Since the DM distribution peaks towards the Galactic center, NSs in this region are expected to contain higher DM fractions that could lead to a different cooling behavior.


INTRODUCTION
Extremely high gravitational field and compactness inside neutron stars (NSs) make them a perfect laboratory to study the strongly interacting matter, test General Relativity and physics beyond the Standard Model (Baym et al. 2018a;Kramer et al. 2021).Throughout their entire lifetime, stars could accumulate a sizeable amount of dark matter (DM) in their interior, which will impact the matter distribution, masses, radii, evolution, etc. (Scott et al. 2008;Lopes & Lopes 2019).At the end of its evolution, a main sequence star of 8-20 M ⊙ undergoes a supernova explosion, creating an NS (Heger et al. 2003).The former is from the gravitational collapse of molecular cloud regions, which exceed the Jeans limit.The proto-cloud may already present traces of DM, facilitating the collapse and giving rise to newly born stars with a sizeable amount of DM (Yang et al. 2020).Once the star is born, DM particles could be further accreted from a surrounding medium, leading to an even higher DM fraction inside the object (Brito et al. 2015;Kouvaris & Tinyakov 2011).A steady DM accretion requires a long timescale and a high DM fraction in the surrounding medium.The maximum amount of the total accreted mass of DM is 0.01% of the star's total mass in the most central part of the Galaxy (Ivanytskyi et al. 2020).On the other hand, rapid accumulation and a significant increase of the DM fraction inside the star could occur while passing through an extremely dense region of primordial DM clumps in a subhalo (Bramante et al. 2022).As predicted by many cosmological models, primordial density perturbations could result in a large fraction of DM forming gravitationally ★ afonso.avila@student.uc.pt † edoardo.giangrandi@uni-potsdam.de ‡ violetta.sagun@uc.pt§ oleksii.ivanytskyi@uwr.edu.pl¶ cp@uc.ptcollapsed objects residing in subhalos (Erickcek & Sigurdson 2011;Buckley & DiFranzo 2018).
At the end of the stellar evolution, the star eventually reaches the iron-core stage, undergoing a core-collapse supernova explosion.During this incredibly energetic event, DM might be created and further accrued inside the remnant, i.e. an NS (Meyer & Petrushevska 2020).Once DM is trapped in the gravitational field of an NS, it may lead to different configurations depending on the DM properties: a core or halo configuration.In the former scenario, DM forms a compact core in the inner regions of an NS.A stronger gravitational pull by the inner core leads to more compact and denser configurations, characterized by smaller maximum gravitational masses and radii compared to the purely baryonic star.Thus, these configurations may be seen as an apparent softening of the baryonic equation of state (EoS) (Giangrandi et al. 2023).NSs with a compact DM core are also harder to deform, an effect that can be tested by future gravitational wave (GW) detections via the tidal polarizability Λ (Giangrandi et al. 2023;Sagun et al. 2022b;Karkevandi et al. 2022;Dengler et al. 2021;Diedrichs et al. 2023;Thakur et al. 2024).It could manifest itself by a supplementary peak or strong oscillation mode in the post-merger GW spectrum (Ellis et al. 2018;Bezares et al. 2019), production of an exotic waveform (Giudice et al. 2016), or modification of the kilonova ejecta (Emma et al. 2022).The latter probes are expected to be possible with the next generation of GW detectors, i.e., the Einstein Telescope (ET) (Punturo et al. 2010), Cosmic Explorer (CE) (Mills et al. 2018), and NEMO (Ackley et al. 2020).On the other hand, when the radius of the DM component exceeds the baryonic one, a halo structure is formed, fully embedding the star.This leads to an increase in the star's gravitational mass, mimicking a stiffening of the baryonic EoS (Sagun et al. 2022a).The halo is easier to be deformed due to the diluted DM distribution surrounding the star (Ivanytskyi et al. 2020;Shakeri & Karkevandi 2022), affecting the tidal polarizability Λ of the star (Sagun et al. 2022b;Diedrichs et al. 2023).
Another way to probe the presence of DM in NSs is by studying their thermal evolution.Standard NS cooling1 occurs via a combination of thermal radiation from the surface and neutrino emission from the interior, making it possible to test the particle composition and properties of matter through x-ray observations.In fact, the thermal evolution of compact stars contains very rich and complicated phenomena at different stages.Thus, it can be divided into three stages: newly born NS with the thermally decoupled core and crust (age ≲ 100 yr) (Sales et al. 2020), neutrino emission dominant stage (age 100 − 10 6 yr) and the photon emission dominant stage (age ≳ 10 6 yr) (Page et al. 2004).
The first stage corresponds to the time required for the core and crust of the newly born NS to become thermally equilibrated, the thermal relaxation time.A different composition of the NS's core and crust results in a substantial variance in their thermal conductivity and specific heat (Lattimer et al. 1994).Neutrinos, emitted from the core, create a cold front that advances toward the surface.When it reaches the surface, its temperature suddenly drops, marking the start of a core-crust thermal connection.Prior to this, the surface temperature of the star remains unchanged as neutrinos slowly diffuse from the core towards the surface, supplying energy that counterbalances cooling.
Further on, NS cooling is mainly defined by the particle composition of the NS core.Particularly, the amount of protons (or, equivalently, electrons and muons together) defines whether the star undergoes a rapid cooling governed by the direct Urca (DU) process, or cools down slowly/intermediately via the nucleon-nucleon bremsstrahlung, modified Urca (MU), and Cooper pair breaking and formation (PBF) processes (Page et al. 2006).In stars where the DU process takes place, a significant temperature drop is observed, heavily modifying their thermal evolution.The underlying EoS determines whether the DU process of neutron -decay and its inverse process can occur in the NS interior (Lattimer et al. 1991;Page et al. 2006;Potekhin et al. 2015).For these reactions to take place, the Fermi momenta of the involved particles have to satisfy the kinematic restriction of the triangle inequality,    +   ≥   given in terms of the Fermi momenta of protons, electrons and neutrons (for n,p,e matter).This condition ensures that for strongly degenerate fermions the reaction is constrained by the Pauli blocking principle, meaning that it can only take place when the energies of the particles involved are close to their Fermi energies.By considering charge neutrality and the relation between the Fermi momenta and the number density of each particle, the proton fraction,   , should be above ∼ 11% (Lattimer et al. 1991).When the threshold condition is not satisfied, the less effective MU process is the dominant neutrino emission process due to the presence of a spectator nucleon that mediates the reaction.
When the star is cool enough so that neutrino emission becomes less relevant, the photon emission from its surface becomes the dominant cooling process.Typically, the star continues to cool down until it becomes invisible to X-ray telescopes and the peak of the black body radiation shifts towards longer wavelengths.However, additional contributions that go beyond the minimal cooling paradigm could alter the above-mentioned scenario by contributing with additional heating or cooling channel.While the accretion of matter from a companion star (Wĳnands et al. 2017) and magnetic field decay (Aguilera et al. 2008) only deposit energy into the star, the presence of DM, depending on the considered candidate, could contribute to either cooling or heating.Through emission-evaporation of light DM from the NS core and/or surface, DM could carry away energy further cooling the star (Kumar et al. 2022).While the emission of DM from the star is mostly considered for light particles that can freely escape, e.g.axions, heavy DM (with mass above the MeV scale) could also evaporate from the star's surface (Garani & Palomares-Ruiz 2022).Many studies have been conducted to model the effect of axion emission on NS and proto-NS (Dietrich & Clough 2019) thermal evolution.Axions produced within the NS cores in, e.g., nucleon bremsstrahlung or PBF processes, leave the star, contributing to its cooling (Sedrakian 2016;Buschmann et al. 2021Buschmann et al. , 2022)).
On the other hand, DM may heat the star during the accretion or self-annihilating, depositing energy in the system (Baym et al. 2018b;Motta et al. 2018a,b;Berryman et al. 2022).DM via scattering with the Standard Model particles may deposit the kinetic energy gained falling into a steep NS gravitational potential, the so-called "dark kinetic heating" (Baryakhtar et al. 2017;Raj et al. 2018).Moreover, trapped symmetric DM could annihilate in the NS interior (de Lavallaz & Fairbairn 2010), or decay, leading to a further heating channel.Depending on the model, this process can be observed in either middle-time (Ángeles Pérez-García et al. 2022) or late-time heating (Hamaguchi et al. 2019).As shown by Kouvaris (2008) the NS cooling curves at the late stage may have a plateau corresponding to DM self-annihilation.This scenario appears to be the easiest to test, as the effect of nuclear superfluidity/superconductivity and magnetic field are negligible at this age.The major contribution to the cooling of old NSs comes from the photon emission, while the neutrino cooling stage sensitive to a particle composition is suppressed.Therefore, a possible heating mechanism of NSs due to DM self-annihilation could be probed by the increasing statistics on observational data of old NSs, in the range from soft x-ray to infrared bands, including the operating James Webb Space Telescope (Chatterjee et al. 2023) and the forthcoming Thirty Meter Telescope (TMT) and Extremely Large Telescope (ELT) (Skidmore et al. 2015).
Alternatively, asymmetric DM considered in this work interacts with BM only through gravity, and, therefore, does not contribute to either neutrino, photon emission, or self-annihilation.It is assumed that DM is accrued during the previous stages of a star's evolution and remains unchanged during the NS lifetime.Accounting for the accretion of DM during the NS equilibrated stage with consequent kinetic heating requires a more extensive study that we leave for future work.
The DM fraction up to a few percent could be accrued during the previous stages of star evolution by several mechanisms, including the DM production during a supernova explosion, the rapid accumulation in the main-sequence star and proto-NS while passing the extremely dense subhalos (Bramante et al. 2022) or be present in the star progenitor due to the primordial clump of DM formed during the DM-dominated era (Yang et al. 2020).
The paper is organized as follows.In Section 2, we present models for the BM and DM components.In Section 3, we discuss the main processes ruling the NS thermal evolution and how we implemented the second fluid in our calculations.In Section 4 the main results are presented, including the DM effects on the DU threshold and the cooling curves.Section 5 includes conclusions and discussions of the smoking gun signal of the presence of DM that could be tested in the near future.Throughout the article, we utilize the unit system in which ℏ =  =  = 1.  1. Parameters of the IST and FSU2R models.The table includes the saturation density  0 , energy per baryon /, incompressibility factor  0 , symmetry energy   , and its slope at the saturation density , as well as the maximum gravitational mass   , radius of the 1.4 M ⊙ star, and the NS mass, baryonic density and proton fraction that characterize the onset of the DU process.

DM-ADMIXED STARS
In this section, we first summarize the properties of the two BM models chosen, IST and FSU2R, and the DM model.We next review the TOV equations for two-fluid configuration and show the massradius curves for the BM models chosen with different fractions of DM.

BM models
To address the uncertainties of the BM EoS we chose two models with different nuclear matter properties at the saturation density.
Particularly, with the different values of the symmetry energy slope L and incompressibility factor  0 = 9(

𝜕 𝑝 𝜕𝑛 𝐵𝑀
)  0 at saturation.The first model is based on the induced surface tension (IST) approach, formulated with an explicit account of the hard-core repulsion among the particles.It was shown by Sagun et al. (2014) that in a dense medium, a short-range repulsive interaction among the particles induces an additional contribution to the single-particle energy, the IST term.At high densities, the IST contribution is negligibly small compared to other terms in the single-particle energy, and the excluded volume treatment of hard-core repulsion is switched to the proper volume regime.In the dilute gas limit, this approach recovers the first four virial coefficients of hard spheres.In the IST EoS the hard-core radius of nucleons was obtained from the fit of the heavy-ion collision data (Sagun et al. 2018;Bugaev et al. 2021).Application of the IST EoS to the nuclear liquid-gas phase transition with its critical endpoint allowed constraining the parameters of the eigensurface tension of nucleons (Sagun et al. 2017), while the attraction and symmetry energy terms were fitted to the NS observables (Sagun et al. 2019).Hereby, we use the Set B of the IST EoS developed in Sagun et al. (2020).
The second considered model is the nucleonic relativistic meanfield FSU2R EoS (Tolos et al. 2017), which is a further development of the FSU2 approach (Chen & Piekarewicz 2014) with a softer symmetry energy and neutron matter pressure.This allows the FSU2R EoS to describe NS with radii smaller than in the case of FSU2 EoS.The parameters of this model were fitted to the binding energies, charge radii, and monopole response of atomic nuclei across the periodic table.It equally well reproduces the properties of nuclear matter and finite nuclei.
For the realistic description of the outer layers, the IST and FSU2R EoSs are supplemented by the Haensel-Zdunik (HZ) EoS for the outer crust and the Negele-Vautherin (NV) EoS for the inner crust (Haensel & Zdunik 1990;Negele & Vautherin 1973).

DM model
DM is modeled as a relativistic Fermi gas of non-interacting particles with spin one-half, which has been extensively studied in the literature, e.g.Nelson et al. (2019); Ivanytskyi et al. (2020); Sagun et al. (2022b).The expressions for the pressure and energy density in the grand canonical ensemble can be written as where ) are the DM degeneracy factor, chemical potential, number density, and Fermi momentum, respectively.We consider DM to be at zero temperature.It is well motivated considering that DM is accrued long enough to thermalize with the NS matter and reach the same vanishing temperature.

Two-fluid configurations
Since DM interacts with BM only gravitationally, the stress-energy tensors of the two fluids are conserved separately, which allows us to write down the two coupled Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman 1939;Oppenheimer & Volkoff 1939) where the subscript index  = BM, DM labels the components,  tot =    +   with   / = 4 2   is the total gravitational mass enclosed to the sphere of radius  and  tot =    +   is the total pressure.
By fixing the values of the central chemical potentials of each of the components we are able to obtain NSs of different mass admixed with a different amount of DM.It is convenient to work in the grand canonical ensemble, as the chemical potentials of two components are related to each other, as The full derivation of this relation is presented in Ivanytskyi et al. (2020).It shows, that for a given stellar configuration the chemical potentials of two components scale proportionally, i.e.   ∝    , which simplifies solving the system of two coupled TOV equations.
After integration of the TOV equations, the gravitational masses of each of the components are recovered, and it is possible to define the DM mass fraction as    =     tot .The mass-radius relations computed for the IST EoS (solid curves) and FSU2R EoS (dashed curves) are shown in Figure 1.The effect of increasing the DM faction from 2% to 4% is shown in different colors.As is seen, the presence of the DM core reduces the gravitational maximum mass of NSs and their radii, being similar to an apparent softening of the BM EoS.Also included in the figure are the constraints obtained from different NS observations as identified in the caption.
For our study, we chose a single value of the DM particle mass,    =1 GeV, that leads to the formation of the distinctive DM core.However, for completeness, the scans over the relative DM fraction and particle's mass are performed in Section 4.4.In Fig. 2 the depicted number densities for BM and DM as a function of radius for the IST EoS (upper panel) and FSU2R EoS (bottom panel) show the relation between the compactness of the DM core and BM redistribution.Thus, with an increase of the DM fraction, a more compact DM core pulls BM inward, leading to a smaller star radius.This effect will be further discussed in the context of the photon luminosity in Section 4.2.The radii of the 1.9 M ⊙ DM-admixed stars with the respective profiles shown in Fig. 2 are given in Table 2.

NS COOLING
The evolution of NSs is governed by the thermal balance equation (Page et al. 2006) where  v is the total specific heat of the stellar matter,  ∞  ,  ∞  ,  ∞  are the redshifted surface temperature, neutrino luminosity, and photon luminosity, respectively.The last term in Eq. ( 4) accounts for any additional source of heat (sign "+") or carrying energy away (sign "−") discussed in Introduction.In this study, no additional source of heating or cooling is considered, therefore,  ∞ ≡ 0.
The photon luminosity depends on the star's radius as   = 4R 2   4 .The redshifted functions are obtained by multiplication by the  Φ factor, being the   component of the Schwarzschild metric.It includes the metric function Φ, which dependence on the radial coordinate , and can be obtained by solving the following ODE where  tot =   +    ,  tot =   +    are the total pressure, and energy density, respectively.In comparison to the photon luminosity that follows the typical black body radiation law, the neutrino emissivity in each process depends on factors such as density, temperature, and the Cooper pairing between nucleons.When the DU process is not allowed, the MU process makes a dominant contribution in removing heat in the form of  emission, unless  and  are in a paired state.The latter substantially suppresses the rates of neutrino emission.Thus, neutrons/protons in the NS interior exist in superfluid/superconducting states by forming Cooper pairs.Despite the suppression of the rates of neutrino emission, the breaking of Cooper pairs is a new source of  balancing the former one.As a result, Cooper pairs are constantly breaking and forming (the PBF process) providing the medium cooling rate (Potekhin et al. 2015).This process gets activated for  and  when the temperature reaches their respective critical values.The pairing strength is defined by the gap parameter, its form, and the peak of the PBF emissivity, varying for the adopted gap model.Free neutrons in the inner crust and protons in the core pair in a singletstate ( 1  0 ), while neutrons in the core are expected to undergo a triplet-state pairing ( 3  2 ) (Bardeen et al. 1957;Page et al. 2006).Consequently, the application of various models of nucleon pairing to explore the cooling of NSs also offers the opportunity to gain deeper insights into the properties of NS matter.
Implementation.We use the publicly available thermal evolution code NScool (Page 2016).To account for the gravitational impact of DM, profiles were generated for each target mass using the twofluid formalism.As asymmetric DM does not directly contribute to a star's cooling, all particle species remain unchanged.Hence, there are two types of variables used by the one-fluid framework of NSCool: the ones related to the BM EoS that do not include the DM contribution, i.e. baryon density and particle fractions, and the ones that include such a term, i.e., total pressure, total energy density, total gravitational mass and the metric functions as the two fluids exist in the same spacetime.

The DU onset
We start with a pure BM 1.9 M ⊙ NS, where DU is still not active in its center, and add DM particles of mass    =1 GeV.As shown schematically in Fig. 3, the accrued DM inside the core triggers the DU process.An increase of the DM fraction causes an extension of the region where the neutron -decay and its inverse process are allowed (these two regions are depicted in Fig. 3 in dark grey and light red, respectively).The rest of the star matter outside of the light red region cools down via the slow/medium processes.For a better comparison, in Fig. 3 all the radii are normalized to the outermost baryonic radius of the considered configuration.The physical values of the radii and total gravitational masses of the considered configurations in this figure and Figs.5-6 are presented in Table 2.
The asymmetric DM interacting only via gravity with the Standard Model particles does not directly affect the neutrino or photon emission.Therefore, the DU process is kinematically allowed at the same proton fraction and central baryonic density for stars with and without DM.As it can be seen in Fig. 4 the presence of DM alters the value of the total gravitational mass at which the DU process is kinematically allowed.With an increase of the DM fraction, we see a drastic reduction in the total gravitational mass for both models.The total gravitational mass of the star with the triggered DU gets dramatically reduced.Thus, for the FSU2R EoS the mass of the star at the DU onset,   , drops from 1.92 M ⊙ to 1.83 M ⊙ , 1.79 M ⊙ , and 1.75 M ⊙ for 2%, 3%, and 4% of DM, respectively.For pure BM stars modeled by the IST EoS, the DU processes are allowed at 1.91M ⊙ .The vertical solid (for the IST EoS) and dashed (for the FSU2R EoS) grey lines in Fig. 4 correspond to the central BM densities of the stars at which the DU process is activated.The points at which the vertical grey lines cross the sequence of curves in Fig. 4 indicate the total gravitational mass and proton fraction at the DU onset.The blue, red, and green curves indicate 2%, 3%, and 4% of DM, while the grey curves show the BM stars (for details see Fig. 4).For the IST EoS the enhanced cooling starts at 1.83M ⊙ , 1.80M ⊙ , and 1.76M ⊙ for 2%, 3%, and 4% of DM, respectively.The small differences observed between both models are due to the denser BM core in IST compared to FSU2R.As can be seen, this phenomenon has a drastic effect on NS cooling.These results could be compared to the cooling of pure baryonic stars described by the FSU2R EoS  2. All configurations correspond to NSs with a total gravitational mass of 1.9 M ⊙ .performed in Negreiros et al. (2018) and by the IST EoS in Tsiopelas & Sagun (2020, 2021).Fig. 4 shows that the proton fractions of the DU process at the onset are different for the IST and FSU2R EoSs.This is due to the fact that the former takes into account the effects of excluded volume in nuclear matter, which at a given value of baryon density leads to larger Fermi momenta of nucleons compared to the point-like case of the FSU2R EoS.

Thermal Evolution
Fig. 5 shows the redshifted surface temperature  ∞  as a function of stellar age for 1.2, 1.6, and 1.9 M ⊙ stars.The presented cooling curves show a fit of the observational data obtained considering 1  0 neutron and proton pairings, described by the SFB (Schwenk et al. 2003) and CCDK (Chen et al. 1993) models.The color grade in Fig. 5 depicts the different DM fractions.With an increase of the DM fraction, the curves are lighter.Thus, pure BM stars for the IST EoS (upper panel) and FSU2R EoS (bottom panel) are shown in red, blue, and green for 1.2, 1.6, and 1.9 ⊙ stars, respectively.As it can be seen, in both models the DU process does not operate at 1.9 ⊙ .However, an accumulation of DM particles with    = 1 GeV of    ≃ 0.161% (IST EoS) and    = 0.378% (FSU2R EoS) triggers the previously forbidden process to operate.
We also address the possibility of different envelope composition that affects the relation between the surface and core temperatures, and, consequently, the observed surface luminosity.Thus, the fraction of light elements is accounted for by the factor  = Δ/, whereas Δ is the mass of light elements in the upper envelope.The lightelement envelope, i.e. hydrogen, helium, with  = Δ/ = 10 −7 is depicted with the dashed curves in Figs.5-6, while the heavyelements envelope, mostly carbon, with  = Δ/ = 10 −16 , is shown with solid curves.Figs.5-6 depict the recently updated measurement of the HESS J1731-347 star (object number 8) which is reported to be the lightest and smallest compact object ever detected (Doroshenko et al. 2022).Our analysis shows that the surface temperature of the HESS J1731-347 compact star can be described by the light-element envelope, which is in agreement with the results of Sagun et al. (2023).

Cas A
The compact object in the center of Cas A, being about 356 years old, is an interesting and subject-to-debate object that shows an unusually rapid cooling (Ho & Heinke 2009;Heinke & Ho 2010;Shternin et al. 2011;Elshamouty et al. 2013;Wĳngaarden et al. 2019).The recent combined analysis of the x-ray spectra of Cas A suggests the surface temperature drop is 1.4-2.5% over two decades of observations, the mass of the star of  = 1.55±0.25 M ⊙ (Shternin et al. 2022).Many models have been suggested to explain this behaviour via a strong 3  2 pairing between neutrons in the core (Page et al. 2011;Ho et al. 2015), rapid cooling via the DU process (Taranto et al. 2016), medium-modified one-pion exchange in dense matter (Blaschke et al. 2012), beyond the Standard Model Physics (Hamaguchi et al. 2018), etc.As shown in Tsiopelas & Sagun (2020, 2021) the IST EoS does not necessarily require the inclusion of neutron superfluidity and/or proton conductivity to explain the Cas A temperature drop.However, the obtained mass of the star fitted to the data was found to be 1.96 M ⊙ (Tsiopelas & Sagun 2021).In fact, in many models, the DU process is allowed only at high masses due to the existing threshold.The inclusion of DM solves this problem and allows the DU process to be activated in medium/low mass stars.
The insets on Fig. 6 demonstrate the best fit of the Cas A observational data points.These curves were obtained by supplementing 1  0 neutron and proton pairings with 3  2 neutron pairing described by the T72 (Takatsuka 1972) model.Thus, the upper panel of Fig. 6 ) envelopes, correspondingly.The effect of neutron superfluidity in the 1  0 channel via the SFB model (Schwenk et al. 2003) and proton superconductivity in the 1  0 channel with the CCDK model (Chen et al. 1993) were taken into account.The target masses refer to total gravitational masses.
shows the cooling curves for stars of  = 1.2, 1.6, 1.9 M ⊙ calculated for the IST EoS and supplemented with n 1  0 (SFB model), p 1  0 (CCDK model), n 3  2 pairing (T72 model) with the maximum critical temperature   = 6.596 • 10 8 K.The best agreement with the observational data is obtained for the  = 1.9 M ⊙ star with    = 2 − 4%.For the FSU2R EoS, the best fit is achieved for the same combination of the pairing gaps as for the IST EoS with a slightly higher value of the maximum critical temperature for the n 3  2 pairing,   = 7.105 • 10 8 K.The light blue curve on the inset on the lower panel of Fig. 6 shows the  = 1.6 M ⊙ star with    = 4% and light-elements envelope that provides the best fit to the Cas A data and agrees with its estimated mass.

Scan over DM parameters
Fig. 4 shows how an increase in DM fraction for    =1 GeV particles leads to a decrease in the star's mass.As it was discussed by Giangrandi et al. (2023), the DM particle's mass and relative fraction inside the star have a comparable impact.Thus, similar DMadmixed configurations could be obtained by increasing the particle mass for the fixed value of the fraction, or vice versa.To see this effect on the DU onset and total gravitation mass for the IST EoS (upper panel) and FSU2R EoS (bottom panel) we perform scans shown in Fig. 7.The color maps denote the total gravitational mass of the stars at which the DU process is kinematically allowed.The black dash-dotted and solid curves depict the attained maximum total gravitational mass for these stars.From Fig. 7 we can see that, while the onset of the DU process for the pure BM stars occurs at 1.91 M ⊙ (IST EoS) and 1.92 M ⊙ (FSU2R EoS), it could decrease below 1.6 M ⊙ for    ≥ 3 GeV and/or DM fraction    ≥ 2 %.Therefore, to see a significant effect on NS cooling the fraction of heavy DM does not need to be high.If we consider an average estimated mass of Cas A central object as 1.6 M ⊙ (Shternin et al. 2022) an accumulation of DM in the range of fractions and masses is shown in the top right corners on both panels of Fig. 7.

CONCLUSIONS
In this study, we focused on the effects of asymmetric fermionic DM on the NS thermal evolution.Despite asymmetric DM that interacts with BM only gravitationally contributes neither to neutrino and photon emission directly nor deposits energy to the system, it alters the thermal evolution of NSs indirectly.The calculations were performed under the assumption of a negligibly low DM accretion rate during the NS evolution, or, equivalently, a fixed DM fraction after the NS formation.In this scenario, there is no effect of dark kinetic heating on the NS's thermal evolution.We demonstrate that an accumulated DM pulls inwards BM from the outer layers, significantly increasing the central density, hence modifying the BM distribution.Consequently, the onset of the DU process is triggered at lower NS masses, leading to a highly efficient and rapid cooling, which is substantially different from the case when it is forbidden.At the same time, the proton fraction corresponding to the DU onset remains the same, as for the pure BM star with the same central BM density.We show that despite the DU process is kinematically allowed only at 1.91 M ⊙ for the IST EoS and 1.92 M ⊙ for the FSU2R EoS, an accumulation of DM particles with    = 1 GeV of    ≃ 0.161% (IST EoS) and    = 0.378% (FSU2R EoS) triggers the previously forbidden process.An increase of the DM particle's mass    ≥ 3 GeV and/or DM fraction    ≥ 2 % shifts the DU onset even below 1.6 M ⊙ .
The effect of DM is also illustrated for the thermal evolution of the compact object in the center of Cas A. We find the best fit of the surface temperature drop of Cas A with the FSU2R EoS supplemented with neutron and proton singlet pairing, and triplet neutron pairing for the  = 1.6 M ⊙ star with a DM fraction of 4%.This result agrees very well with the estimated mass of the object of  = 1.55±0.25 M ⊙ (Shternin et al. 2022).
An additional effect of DM is related to the pull of BM inward, creating a more compact core and reduction of the baryonic radius.Thus, the total surface of the star is reduced leading to a lower photon luminosity.This effect is clearly visible at the photon-dominated stage when the neutrino emission takes a subdominant role.
We showed that either an increase of the DM fraction or the particle's mass causes a shift of the DU onset toward lower star gravitational masses.This effect could be a smoking gun signature of the presence of DM in compact stars.Thus, low/middle mass NSs that are not expected to have an operating DU process in fact might have it due to the presence of DM.As a result, stars of the same mass will show a different cooling pattern depending on the DM fraction accrued in their interior.This effect could be also used to map the DM distribution.In this context, present, e.g.XMM-Newton, NICER (Riley et al. 2021;Miller et al. 2021), and future, e.g.ATHENA (Cassano et al. 2018), eXTP (in 't Zand et al. 2019), and STROBE-X (Ray et al. 2019), x-ray observational programs look very promising as they expect to increase the number of mass, radius, and surface temperature determinations.As we expect a higher DM fraction inside compact stars toward the Galactic Center, their thermal evolution could exhibit a distinct feature from the stars in the solar vicinity.

Figure 1 .
Figure1.Total gravitational mass of the DM-admixed NS as a function of its baryonic radius   calculated for the DM particle mass    =1 GeV.Solid and dashed black curves correspond to pure BM stars described by the IST EoS and FSU2R EoS, respectively.Royal blue, red, and green curves characterize relative DM fractions equal to 2%, 3%, and 4%, respectively.Green, gray, and cyan bands represent 1 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).Red and black contours show the NICER measurements of PSR J0030+0451(Miller et al. 2019;Riley et al. 2019), while yellow and blue contours correspond to the PSR J0740+6620 measurement(Raaĳmakers et al. 2021;Miller et al. 2021).LIGO-Virgo collaboration observations of GW170817(Abbott et al. 2018) and GW190425(Abbott et al. 2020) binary NS mergers are shown in blue and pink.The 1 and 2 contours of HESS J1731-347(Doroshenko et al. 2022) are plotted in dark and light orange.

Figure 2 .
Figure 2. The split number density profiles for BM and DM for the 1.9  ⊙ star.The impact of DM of different fractions on the star's profile is shown for the IST EoS (upper panel) and FSU2R EoS (bottom panel).The calculations are performed for DM particle mass    = 1 GeV.The red arrow on the y-axis indicates the density of the DU threshold.

Figure 3 .
Figure 3. Stellar configurations with different DM fractions for the IST EoS (upper part) and FSU2R EoS (lower part).The size of the DU region and DM core are depicted in dark grey and light red, respectively.For comparison, the radii are normalized to the outermost baryonic radius of each configuration and are given in Table2.All configurations correspond to NSs with a total gravitational mass of 1.9 M ⊙ .

Figure 4 .
Figure 4.The proton fraction,   , and total gravitational mass of NSs as a function of the central baryonic density,  , .The curves are equivalent to those in Fig. 1.The vertical solid (for the IST EoS) and dashed (for the FSU2R EoS) grey lines correspond to the central BM densities of the stars at which the DU process is activated.The intersection points correspond to the minimal total gravitational mass at which the nucleonic DU process threshold is satisfied.

Figure 5 .
Figure 5. Cooling curves for stars of different mass  = 1.2, 1.6 and 1.9 M ⊙ .The calculations are performed for the IST EoS (upper panel) and FSU2R EoS (bottom panel) and the relative DM fractions of    = 0%, 2%, 3% and 4%.The considered mass of DM particles is    = 1 GeV.The axes correspond to the redshifted surface temperature vs. year.The solid and dashed lines refer to heavy ( = Δ/ = 10 −16 ) and light-element ( = Δ/ = 10 −7 ) envelopes, correspondingly.The effect of neutron superfluidity in the 1  0 channel via the SFB model(Schwenk et al. 2003) and proton superconductivity in the 1  0 channel with the CCDK model(Chen et al. 1993) were taken into account.The target masses refer to total gravitational masses.

Figure 6 .
Figure6.The same as Fig.5supplemented with the triplet neutron pairing described by the T72 model(Takatsuka 1972).The insets on both panels show the observational data and the best fit within the IST and FSU2R EoSs of the Cas A measured by Chandra ACIS-S in GRADED and FAINT modes.The green and blue data points correspond to the variable and fixed absorbing hydrogen column density   = 1.656 • 10 22  −2 in the FAINT mode, whereas the green and blue data points depict the same data for the GRADED mode taken fromShternin et al. (2022).

Figure 7 .
Figure 7. Scan over the particle's mass,    , and fraction,    , of DMadmixed stars, shown for the IST EoS (upper panel) and FSU2R EoS (bottom panel).The color shows the star's total gravitational mass at which the DU process starts to operate.Dash-dotted and solid black curves are the contour lines showing the maximum gravitational mass obtainable for the two models.

Table 2 .
Parameters of the considered DM-admixed stars in Figs. 5 and 6.