The impact of asymmetric dark matter on the thermal evolution of nucleonic and hyperonic compact stars

: We investigate the impact of asymmetric fermionic dark matter (DM) on the thermal evolution of neutron stars (NSs), considering a scenario where DM interacts with baryonic matter (BM) through gravity. Employing the two-fluid formalism, our analysis reveals that DM accrued within the NS core exerts an inward gravitational pull on the outer layers composed of BM. This gravitational interaction results in a noticeable increase in baryonic density within the core of the NS. Consequently, it strongly affects the star’s thermal evolution by triggering the early onsets of the direct Urca (DU) processes, causing enhanced neutrino emission and rapid star cooling. Moreover, the photon emission from the star’s surface is modified due to a reduction of radius. We demonstrate the effect of DM gravitational pull on nucleonic and hyperonic DU processes that become kinematically allowed even for NSs of low mass. We then discuss the significance of observing NSs at various distances from the Galactic center. Given that the DM distribution peaks toward the Galactic center, NSs within this central region are expected to harbor higher fractions of DM, potentially leading to distinct cooling behaviors.


Introduction
The remarkable compactness of neutron stars (NSs) makes them the densest objects accessible by direct observations.Due to this, NSs serve as ideal laboratories to study the properties of strongly interacting matter, as well as to probe General Relativity and explore physics beyond the Standard Model [1,2].Throughout the entire lifetime of an NS, a sizeable amount of dark matter (DM) could be accrued within the star's interior, affecting the matter distribution, mass, radius, etc. [3][4][5][6][7][8][9][10][11][12].An NS typically originates from a main sequence star of mass ranging from 8 to 20 M ⊙ that at the end of its lifetime undergoes a supernova explosion [13].On the other hand, this stellar journey starts from the gravitational collapse of molecular cloud regions surpassing the Jeans limit.It is worth highlighting that these proto-clouds may already harbor traces of DM, fastening up the collapse and ultimately creating newly formed stars admixed with a significant amount of DM [14].As the star evolves, it undergoes a core-collapse supernova, during which DM may potentially be created and/or accrued in the resulting remnant, namely, an NS [15].Once an NS is born, it can accumulate DM particles from the surrounding galactic medium, increasing the DM fraction within the object [16,17].The latter case is valid for asymmetric DM, which refers to models characterized by an asymmetry in the number densities of particles and antiparticles.
An estimation of the amount of DM being accreted through the Bondi accretion in the Galactic center gives up to 0.01% of the star's total mass [4].In addition, based on the predictions of many cosmological models, a rapid DM accumulation leading to a higher DM fraction inside the star could occur while passing through an extremely dense subhalo with primordial DM clumps [18].Thus, primordial density perturbations could result in a large fraction of DM forming gravitationally collapsed objects or clumps residing in subhalos [19,20] or even the formation of dark compact objects [21][22][23].
As DM is gravitationally trapped by an NS, it can form different configurations depending on its properties, e.g.either sinking to the center and forming a dense core, or creating a dilute halo embedding the whole NS.In the former scenario, heavy DM particles tend to create a dense inner region inside a baryonic star.The gravitational pull of the DM core leads to a denser and more compact NS configuration characterized by smaller maximum gravitational mass and radius compared to a purely baryonic NS.Consequently, these configurations can be perceived as the ones characterized by an apparent softening of the baryonic matter (BM) equation of state (EoS), or, equivalently, lower pressure at a certain baryonic density [8].
NSs with a dense DM core are more compact and harder to be deformed by an external gravitational field, a phenomenon that can be probed with forthcoming gravitational wave (GW) detections of binary NS mergers and black hole-NS mergers by assessing the tidal polarizability parameter Λ during the inspiral phase of the merger [5,6,8,11].Moreover, the first studies show that DM might manifest itself with an additional peak or pronounced oscillation mode in the post-merger GW spectrum [24], generation of an exotic waveform [25], modification in the kilonova ejecta [26], formation of a one-arm spiral instability [27], etc. [28][29][30].These smoking-gun signals are expected to become feasible with the next generation of GW detectors, i.e. the Einstein Telescope [31], Cosmic Explorer [32], and NEMO [33].On the other hand, when the radius of the DM component surpasses the one of the BM, a halo structure ensues, enveloping the entire NS [7].This hints at an increase of the NS gravitational mass, mimicking a stiffening of the BM EoS [4,34].Diluted DM distributed around the NS is more susceptible to deformation, consequently impacting the tidal polarizability parameter Λ of the star itself [5,11].
Another approach to probing the presence of DM within NSs involves the study of their thermal evolution.The standard cooling of NSs includes a combination of thermal radiation emitted from the surface and neutrino emission originating in the NS interior.This offers an opportunity to relate the particle composition, the effect of nucleonic and quark superfluidity/superconductivity, as well as the other properties of matter to the NS surface temperature through X-ray observations [35,36].Notably, thermal evolution depends on several different phenomena, occurring throughout the NS thermal history.The latter can be split into three primary stages: a newly formed NS with a thermally decoupled core and crust (age ≲ 100 years) [37], the subsequent phase dominated by neutrino emission (age 100-10 5 years), and a final phase dominated by photon emission (age ≳ 10 5 years) [38,39].The initial phase corresponds to the duration required for the core and crust of an NS to reach a thermal equilibrium, referred to as the thermal relaxation time.During this stage, the surface temperature of the NS remains constant as neutrinos gradually diffuse through the optically thick medium, supplying enough energy to counterbalance the cooling.As neutrinos reach the stellar surface, the star's temperature drops as a consequence of the thermal connection between the core and crust.
In the following stage, NS cooling is primarily determined by the particle composition of the NS core, which is defined by the BM EoS.Thus, the high proton fraction and the existence of hyperons in the NS interior are especially interesting as they allow the occurrence of β-decay, and its inverse process, known as the direct Urca (DU) process [40][41][42].The DU process, once kinematically allowed by the Fermi momenta triangle inequality condition, leads to strong neutrino emission, and, consequently, to a rapid stellar cooling [43][44][45].In the case where the threshold condition is not met, the less efficient modified Urca (MU) process becomes the dominant one.This process does not have a threshold, as the presence of a spectator nucleon facilitates the reaction.The MU together with the nucleonnucleon bremsstrahlung and the Cooper pair breaking and formation (PBF) mechanism contribute to the slow/intermediate cooling [44].
At density n = 2n 0 − 2.5n 0 , with n 0 being the nuclear saturation density, hyperons become new degrees of freedom [46][47][48].The presence of hyperons fundamentally alters the EoS, resulting in its significant softening and the corresponding reduction in the maximum mass achievable for NSs [49,50].As a consequence, hyperonic EoS barely provide the NS masses of 2 M ⊙ [51][52][53][54].While providing the existence of heavy NSs with masses above ∼ 2 M ⊙ , stiffer hadronic EoSs (see e.g.Ref. [55]) are discriminated by the observational constraints on the tidal deformability [56].These difficulties are naturally overcome within the scenario of early deconfinement of quark matter before the hyperon onset [57][58][59].In this case, all the observational constraints can be fulfilled simultaneously [60].
The presence of hyperons affects the NS thermal evolution as well.As a matter of fact, hyperonic DU takes place at a very low fraction of hyperons, typically of the order of several percent.In other words, the onset of the process is at densities slightly above the point where hyperons become energetically favorable [61,62].As it will be discussed in Sec. 4, in comparison to it, the nucleonic DU is highly dependent on the proton fraction and, therefore, might be allowed only at much higher baryonic densities, depending on the density dependence of the symmetry energy.
After ∼ 10 5 years since the NS formation, the neutrino emission processes become less relevant, and photon emission from the surface starts to take the dominant role.Further on, an NS keeps on cooling down until it becomes undetectable by X-ray telescopes, and the peak of its black body radiation shifts to longer wavelengths.They cool down so drastically that they enter the sensitivity band of the James Webb Space Telescope (JWST).Thus, the JWST becomes an important tool in searching for old NSs [63][64][65].However, deviations from this cooling scenario may occur due to additional heating or cooling channels.Accretion of matter from a companion star [66], magnetic field decay [67], rotochemical heating [68] as well as decay or annihilation of DM particles to Standard model states inside the star increase its surface temperature [69][70][71][72][73].The latter could be probed in both middle-and late-time heating scenarios [68,74].Old NSs may exhibit a plateau in their cooling curves, corresponding to a balance between the energy loss and heating caused by DM self-annihilation [75,76].The late-time heating is another sensitive probe of DM as the influence of nuclear superfluidity/superconductivity, magnetic field, and other factors that may wash out the effect, become irrelevant at this age.Accumulated asymmetric DM also contributes to a star's heating by depositing the kinetic energy gained during the infall, so-called "dark kinetic heating" [63,64].The primary contribution to the cooling of old NSs is photon emission from the surface, as neutrino cooling is suppressed, making it feasible to probe the heating effects of DM.The ongoing, i.e.XMM-Newton and Chandra [77], NICER [78,79], and future X-ray surveys, i.e.ATHENA [80], eXTP [81], and STROBE-X [82], are expected to increase the number of mass, radius, and surface temperature determinations, including the old NS detections, shedding light on the possible DM impact on their thermal evolution.
On the other hand, light DM particles emitted from an NS carry away energy, further cooling down the star [83].In comparison to the heavy DM particles with mass ≥MeV that could be evaporated from the NS's surface [84], light particles like axions could be emitted through the star.Thus, various studies have explored the effects of axion emission on NS and proto-NS thermal evolution [85], where axions produced within the NS core in nucleon bremsstrahlung or PBF mechanism, will contribute to the star's temperature drop [86][87][88].
This study focuses on asymmetric DM consistent with the ΛCDM model predicting cold (non-relativistic) collisionless DM particles [89].Accumulated asymmetric DM, which interacts with BM only through gravity, is considered to have no impact on the thermal evolution of NSs as it does not annihilate, is heavy, and does not participate in any processes involving BM.However, this study shows that accrued DM by modifying the BM distribution and compactness of the star, affects the DU processes responsible for rapid cooling.As a result, it is triggered in stars of lower mass.In this work, we extend our previous analysis [90] by considering hyperons.We investigate an interplay between hyperons and DM in the inner NS region and the modification of the stars' thermal evolution.This analysis assumes that a DM fraction remains unchanged or accretion is negligibly small after the NS formation, and, therefore, dark kinetic heating does not affect the cooling.
The paper is organized as follows.In Section 2, we present the considered EoSs for BM and DM.In Section 3 we discuss the two-fluid formalism, while Section 4 reviews the NS thermal evolution and all the nucleonic processes involved in it.In Section 5 we provide the results of the cooling simulations for different DM fractions.Conclusions are presented in Section 6.Throughout the article, we use the unit system in which h = c = G = 1.

Baryonic matter
To account for the BM uncertainties at high densities, we adopt three different models characterized by different particle compositions, stiffness, and nuclear matter properties at the saturation density.Specifically, we have chosen these models based on the different values of their symmetry energy slope L and the incompressibility factor K 0 = 9( ∂p ∂n B ) n 0 , at the normal nuclear density.The first of these models is the Induced Surface Tension (IST) EoS, which is formulated considering the hard-core repulsion among particles.Notably, Sagun et al. [91] demonstrated that in a dense medium, a short-range repulsion between particles gives rise to an additional contribution to the single-particle energy, the IST term [92].At high densities, the IST contribution vanishes leading to a correct transition from the excluded volume approach to the proper volume regime.In the limit of a dilute gas, the IST EoS accurately recovers the first four virial coefficients of hard spheres.The values of the hard-core radius of nucleons were obtained from the fit of the particle ratios created in heavy-ion collision experiments in a wide range of the center of mass collision energies [92,93].Furthermore, the IST EoS reproduces the nuclear liquid-gas phase transition, including its critical endpoint [94].The generalization of the model to account for the long-term attraction and asymmetry energy was formulated in agreement with the NSs observations and tidal deformability from the GW170817 binary NS merger [95].In this work, we utilize Set B of the IST EoS developed in [96].
The second considered model is the nucleonic relativistic mean-field FSU2R EoS [46], which is the generalization of the FSU2 EoS [97].The FSU2R model is characterized by a smaller symmetry energy and neutron matter pressure, allowing it to describe NSs with smaller radii compared to the original FSU2 EoS.The model parameters were determined by fitting the binding energies, charge radii, and monopole response of atomic nuclei across the periodic table.Importantly, the FSU2R model reproduces equally well the properties of nuclear matter and finite nuclei.
To study the effects of hyperons on the NS thermal evolution, we adopt the FSU2H EoS in which the parameters of the isospin-symmetric matter were adjusted to solve the hyperon puzzle, providing a stiffer EoS and fulfilling the heavy NS constraints [46].The adjustments were made only above two saturation densities, i.e. after hyperonic onset, to preserve the properties of the model at saturation density and to fulfill finite nuclei and stellar radii.The model is parameterized with the values of the hyperons couplings and potentials (Λ, Σ, and Ξ) as described by Fortin et al. [62].However, following Ref.[98], the Σ potential U (N) Σ was set at 30 MeV and not spanned in a wide range of values.As a result, this parametrization does not allow Ξ 0 hyperons to appear.Table .1 summarizes the main model parameters.
The upper part of the left panel of Fig. 1 shows the relative fraction of proton for the IST and FSU2R EoS that have only n, p, e − as degrees of freedom.For comparison, the FSU2H EoS includes also µ, Λ, Ξ − , and Σ − hyperons.All particle fractions of the FSU2H EoS are depicted on the upper part of the right panel of Fig. 1.These particles appear as a result of the attractive nature of the ΛN and ΞN interactions [49,99].
Another important effect of hyperons on the matter properties is related to modifications of the particle fractions upon their onset.As shown on the upper right panel of Fig. 1, as Λ hyperons start to appear, their decay into protons leads to an increase of the latter ones.1. Parameters of the IST, FSU2R and FSU2H models.The table includes the saturation density n 0 , energy per baryon E/A, incompressibility factor K 0 , symmetry energy E sym , and its slope L at saturation density, as well as the maximum gravitational mass M max , and radius of the 1.4 M ⊙ star.All the BM EoSs were supplemented by the Haensel-Zdunik (HZ) EoS for the outer crust and the Negele-Vautherin (NV) EoS for the inner crust [100,101].

Dark matter
Hereby, DM is described as a fermionic gas of relativistic particles with spin 1/2 [4,5,102].This model is the simplest realization of the asymmetric DM that due to the Pauli blocking principle can resist a gravitational collapse and form stable configurations inside NSs [103].Thus, DM particles are characterized by a small number of parameters, i.e. mass m DM and degeneracy factor g DM = 2, making it a widely used model [4].
The pressure p DM and energy density ε DM of the system can be written down as follows where µ DM , n DM = ∂p DM /∂µ DM and k DM = µ 2 DM − m 2 DM θ(µ DM − m DM ) the DM chemical potential, number density, and Fermi momentum, respectively.
In this work, we consider the DM particle mass m DM =1 GeV.This choice of the DM particle's mass is motivated by the intriguing similarity with nucleons.A study of the broad particle mass range of the fermionic DM in the MeV-GeV scale can be found in Ref. [4].
The Einstein's field equations can be written down as The stress-energy tensor has a form where ε, p, and u µ are the proper energy density, pressure, and four-velocity of the fluid component.Due to a negligibly weak nongravitational interaction between BM and DM the only nonvanishing components in the perfect fluid approximation in the stress-energy tensor are Consequently, each fluid individually satisfies the equations of motion of a single fluid, hence the energy-momentum conservation, ∇ µ T i µν = 0. Thus, the Tolmann-Oppenheimer-Volkoff (TOV) equation [108,109] that describes a non-rotating spherically symmetric star in hydrostatic equilibrium can be split into two coupled first-order ODEs, one for each component [4].
Here p tot ≡ p BM + p DM and M tot = M BM (R BM ) + M DM (R DM ) are the total pressure and total gravitational mass enclosed in a sphere of radius r which can be written as The radii R i are evaluated using the zero-pressure condition at the surface p i (R i ) = 0.For a given EoS and the boundary conditions m i (0) = 0, p i (0) = p c i with p c i being the pressure of each of the components in the center of a star we can obtain the total gravitational mass of the system, radius profile, etc.At the same time, the pressure of each of the components is fixed by the corresponding chemical potential.Therefore, in this work, we used the boundary conditions for the chemical potentials µ i (0) = µ c i with µ c i being the corresponding central values.In this case, the system of Eqs. ( 5) can be converted to obtained in Ref. [4].It follows from these relations that the chemical potentials of BM and DM scale proportionally according to their central values.
It is convenient to define the DM mass fraction as The mass-radius curves for the aforementioned hadronic models are represented in Fig. 2, including the DM-admixed configurations.They were obtained by solving the TOV Eqs. ( 7) at varying the central value of the BM chemical potential and adjusting the central value of the DM one so that its fraction is kept equal to the values indicated in the figure.Total gravitational mass of DM-admixed and baryonic NSs as a function of their baryonic radius R B calculated for the DM particle mass m DM =1 GeV.Solid curves correspond to pure BM stars described by the IST (blue curve), FSU2R (red curve), and FSU2H (green curve) EoSs.Dashed and dotted curves depict the M-R relations obtained for the relative DM fractions equal to 2% and 4%, respectively.Orange and dark yellow bands represent 1σ constraints on the mass of PSR J0348+0432 [52] and PSR J1810+1744 [110].While olive green and light orange contours show the NICER measurements of PSR J0030+0451 [111,112], purple and magenta contours correspond to the PSR J0740+6620 measurement [79,113].Observations of GW170817 [114] and GW190425 [115] binary NS mergers by LIGO-Virgo collaboration are shown in blue and green.The 2σ contour of HESS J1731-347 [116] is plotted in light red.The shaded gray region is excluded by the rotation of the fastest spinning pulsar PSR J1748-2446ad [117].
Fig. 2 shows the solid M-R curves for baryonic stars described by the IST, FSU2R, and FSU2H EoSs in blue, red, and green, respectively.The dashed and dotted curves represent the DM-admixed configurations with DM fractions of 2% and 4% for the corresponding BM EoSs.

NS thermal evolution
The thermal evolution of isolated NSs is determined by the energy balance equation which can be written down as Here C v is the total specific heat of the stellar matter, T ∞ s , L ∞ ν , L ∞ γ are the red-shifted surface temperature, neutrino luminosity, and photon luminosity, respectively [44].The last term on the right-hand side of Eq. ( 9) refers to extra mechanisms contributing to heating (+) or cooling (−) the compact star.In this work, no additional cooling/heating mechanisms were considered, hence H ∞ ≡ 0. As photons are emitted from the stellar surface, the photon luminosity L γ = 4πR 2 B σT s 4 strongly depends on the stellar baryonic radius R B .Here, σ is the Stefan-Boltzmann constant.Once the photon luminosity in the local frame is obtained, it can be red-shifted to obtain the corresponding value in the frame of an observer placed at infinity L ∞ γ .This operation can be performed by multiplying the quantity by the factor e Φ , i.e. by the g tt component of the metric, derived by integrating the following ODE Here ε tot = ε BM + ε DM is the total energy density.
On the other hand, neutrinos are emitted throughout the star as the NS matter is optically thin for these particles.The neutrino luminosity is highly sensitive to various factors, such as Cooper pairing between particles, particle composition, their fractions, etc.The latter ones highly depend on the underlying BM EoS [43][44][45].In fact, the enhanced cooling due to the neutron β-decay and its inverse process, known as the DU process, occur when the Fermi momenta of the involved particles satisfy the triangle inequality condition, p F,i + p F,j ≥ p F,k , where i, j and k are the involved particles and their cyclic permutations.In Eq. ( 11) ℓ stands for a lepton, e − or µ.The triangle inequality condition ensures that, in the presence of strongly degenerate fermions, it can only take place when the energies of the involved particles are close to their Fermi energies.In the case of proton, electron, and neutron composition the triangle condition, reads as p Fp + p Fe ≥ p Fn , with p Fp , p Fe , and p Fn being the Fermi momenta of particles.Taking into account charge neutrality and the connection between Fermi momenta and the number density of each particle, the proton fraction (Y p ) must be above approximately 11% [43].The latter value is obtained for n,p,e − matter.While muons are present, the minimal proton fraction attains a larger value where x e = n e /(n e + n µ ) and n e , n µ are densities of electrons and muons.The triangle condition might not be met for some models, meaning that up to the TOV mass, the nucleonic DU process is not triggered in the stellar core.In our study, all three considered models, the IST, FSU2R, and FSU2H, have the onset of the DU at the central densities and corresponding masses shown in Tables 2, 3. Further on, we will refer to n DU and mass M DU as the baryonic density at which the DU process onsets and the corresponding mass of the star with the central baryonic density equal to n DU .On the other hand, after hyperonization of matter in the inner core, the neutrino luminosity is amplified, triggered by the hyperonic DU [118].Each of the hyperonic DU channels is open after the appearance of all the species involved in the process and the fulfillment of the momentum conservation condition [41,119].Thus, among all hyperonic DU processes, only the following are allowed by the FSU2H EoS

FSU2H EoS
Further on in the text, we will refer to these processes as np, Λp, and Σ − Λ.The remaining processes and their inverse, i.e.
are not kinematically allowed.After the onset of hyperons, the minimum proton fraction for nucleonic DU showed in Eq. ( 12) changes to where and n j (j = e, µ, n, p, Σ − , Σ + , Ξ − ) are the densities of the corresponding sort of particles [120].As the isospin of Σ − hyperon is one, its occurrence is strongly affected by the density dependence of the symmetry energy [120].The DU threshold densities and the star's masses described by the FSU2H model are presented in Table 3.
The triggered nucleonic and hyperonic DU processes have the emissivity of ε ∼ T 6 leading to an enhanced star cooling.When the threshold conditions are not met, the less efficient MU, i.e.N + n → N + p + ℓ + ν ℓ , N + p + ℓ → N + n + ν ℓ (with N = n, p, ℓ = e − , µ), and nucleon bremsstrahlung processes become the predominant mechanisms for neutrino emission with ε ∼ T 8 [121,122].
Finally, at lower temperatures, the Cooper pairing of nucleons contributes to the neutrino emissivity [45].The existence of an attractive force among particles results in the formation of pairs, where the particle excitations are gapped.Thus, when the temperature falls below the critical value for neutron superfluidity and proton superconductivity, denoted as T ≪ T c , the emission of neutrinos is drastically suppressed, following the Boltzmann factor of e −∆/T , where ∆ corresponds to the energy gap.Once activated at the critical temperature, the effect of the PBF mechanism leads to the neutrino emissivity of the order of ε ∼ T 7 [44].In the inner crust, neutrons are expected to create pairs in the singlet 1 S 0 state and the triplet 3 P 2 state in the core.Cooper pairs of protons occur in the singlet 1 S 0 state in the NS core.We examine the above-mentioned nucleon pairing in this article.For simplicity, in this work, no hyperon superfluidity was considered.Based on the previous studies, e.g. in Ref. [50], we expect that incorporation of the hyperonic superfluidity would alter the onset of the hyperonic DU processes towards higher densities.Thus, the process would not be triggered right after the onset of hyperons.This could have a significant impact on the star's cooling whereas nucleonic DU is only permitted at densities exceeding the onset of hyperons, as observed in the three models studied here.However, another interplay between the nucleonic and hyperonic DU is discussed in [50].
To model the thermal evolution of NSs we adopted the public code NSCool [123], using the implementation method described by Ávila et al. [90].The specific heat is calculated as the sum of the contributions from its constituent particles: neutrons, protons, electrons, and hyperons (if present) using the standard NSCool implementation.The impact of nucleon superfluidity is accounted for by introducing control functions R c which depend on the corresponding critical temperature [124,125].
To account for the DM effects, all the NS profiles were obtained using the two-fluid formalism.Despite accumulated DM does not interact with BM, it leads to BM redistribution.Consequently, the DM affects not only the profiles of such quantities as the total pressure, energy density, and gravitational mass, which explicitly include a DM term, but also the profiles of the metric functions, baryon density, and particle fractions.At the same time, the quantities, that do not have a DM contribution, can be directly used as an input for the single-fluid cooling formalism of the NSCool.This allows us to treat a two-fluid cooling problem within the single-fluid framework.

Hyperonic and nucleonic DU onsets
We start our analysis by studying three different mass configurations representing the low-, middle-, and heavy-mass stars (see Table 4).Purely baryonic stars of 1.9 M ⊙ modeled by the IST and FSU2R EoSs do not exhibit a rapid cooling due to the DU process as it is triggered in heavier stars where the central density is higher (see Table 2 for details).However, heavy DM particles of ≳ GeV scale tend to form a dense core inside a star, pulling BM inwards from the outer layers, and leading to BM redistribution.As a consequence, the baryonic density in the inner core increases, triggering the nucleonic and hyperonic DU processes as shown in Figs.3-4.Thus, the onset of the enhanced neutrino emission occurs at the same particle fractions and central baryonic density for stars with and without DM.However, as it can be seen in Fig. 1 and Table 2, with the increase of the DM fraction the total gravitational mass at which the DU processes are kinematically allowed is shifted towards lower masses.Thus, for the FSU2H EoS the mass of the star at the nucleonic DU onset, M DU , drops from 1.85 M ⊙ to 1.75 M ⊙ , 1.71 M ⊙ , and 1.66 M ⊙ for 2%, 3%, and 4% of DM, respectively (see Table 3).Moreover, due to hyperonic DU processes, such as Λp and Σ − Λ, triggered before the nucleonic DU process, np, an accumulation of DM allows for an NS to exhibit a fast cooling behavior at 1.24 M ⊙ , for 4% DM (as shown in Fig. 1).
The values for the central baryonic densities and total gravitational masses at which the hyperonic and nucleonic DU processes become active are presented in Table 3.
As was discussed in Ref. [90], another consequence of the presence of the DM core is related to the reduction of the stellar radius due to the stronger gravitational pull compared to a pure BM configuration.This effect is presented in Table 4.As the photon emission is related to the stellar radius, this effect alters the photon luminosity that becomes important for stars of ≳ 10 5 years [39].

Mapping the nucleonic and hyperonic DU regions inside the star
To investigate an interplay between the size of the DM core and the star regions in which np, Λp and Σ − Λ DU processes are operating we show the pie charts, where each slice represents a different BM model.Thus, in Fig. 3 the top right, top left and bottom slices illustrate stars of 1.75 M ⊙ described by the IST, FSU2R, and FSU2H EoSs, respectively.The DM fraction continuously changes from 0% to 4.5%.All radii of the DM-admixed configurations are normalized to the baryonic radius of each configuration.The physical values of radii and total gravitational masses are listed in Table 4.
As it is shown in Tables 2, 3 none of the pure BM stars of 1.75 M ⊙ have the operating nucleonic DU in their interior.However, by increasing the DM fraction the np DU process is triggered at 3.87%, 4.17%, 1.96% for the IST, FSU2R, and FSU2H EoSs respectively.A much lower value for the FSU2H EoS is related to the fact that the nucleonic DU is triggered at a lower stellar mass, i.e. 1.85 M ⊙ , than for the IST and FSU2R EoSs, i.e. 1.908 M ⊙ and 1.921 M ⊙ (see Tables 2, 3).The nucleonic DU region is depicted with the light red color in Fig. 3.Moreover, parts of the star in which the Λp and Σ − Λ DU processes take place are shown in yellow and light green colors.4. Baryonic radii for different gravitational mass configurations for the three models used in this work and their modifications due to the addition of DM.

IST EoS
For a comparison, Fig. 4 shows stars of 1.90 M ⊙ described with the same models as in Fig. 3.In the case of the FSU2H EoS, not only the hyperonic DU processes but also the nucleonic DU are activated.Moreover, as it can be seen in Fig. 4, no DM-admixed configurations with a DM fraction higher than 2% exist in this model.This is related to the FSU2H EoS that has an upper limit on the baryonic density, n B = 1 fm −3 , at which the effective nucleonic mass vanishes [46].As higher DM fractions require higher central baryonic densities, which can not be obtained within the FSU2H EoS, such configurations do not exist.

Cooling curves
Fig. 5 depicts the red-shifted surface temperature as a function of age for stars of 1.2 M ⊙ (red curves), 1.6 M ⊙ (blue curves), and 1.9 M ⊙ (green curves) modeled within the IST EoS (left panel) and the FSU2R EoS (right panel).The cooling curves for the stars of 1.2 M ⊙ (red curves), 1.5 M ⊙ (blue curves), and 1.7 M ⊙ (green curves) modeled within the FSU2H EoS are shown in Fig. 6.This choice of lower mass values in comparison to Fig. 5 is related to the upper limit of the baryonic density and absence of heavy stars admixed with DM within the FSU2H model.
The cooling curves were obtained considering the neutron and proton 1 S 0 pairing, described by the SFB [126] and CCDK [127] models, as the ones that provide the best description of the observational data.Following [90], in Figs. 5, 6 the color grade represents the different DM fractions, whereas the higher DM fraction corresponds to a lighter shade.Pure BM stars, i.e. with 0% of DM, are shown with the darkest shade of each color.
The envelope composition plays an important role in modeling the surface luminosity as it affects the relation between the surface and core temperatures.In particular, the relation between the internal and surface temperature is very sensitive to the outermost layer composition, which plays the role of an insulator until very low temperatures.To address this question, we consider the fraction of light elements via the factor η = ∆M/M with ∆M being the mass of light elements in the envelope.In Figs. 5, 6 the light-element envelope, i.e. hydrogen, helium, with η = ∆M/M = 10 −7 is depicted with the dashed curves, while the heavy-elements envelope, mostly carbon, with η = ∆M/M = 10 −16 , is shown with solid curves.
On both panels in Fig. 5 the cooling curves for the 1.9M ⊙ stars with 0% of DM show a slow cooling, corresponding to the absence of the DU process.However, an accrued DM of m DM = 1 GeV of f DM ≃ 0.161% (IST EoS) and f DM = 0.378% (FSU2R EoS) triggers the previously forbidden process to operate.Thus, the green curves for f DM = 2%, 3% and 4% exhibit a rapid drop of the surface temperature.
As seen in Fig. 6 the FSU2H model provides a good agreement with the observational data that fully agrees with the results of Ref. [128].For the pure BM star of 1.7M ⊙ only the Λp DU process operates leading to a medium cooling.At f DM ≃ 0.61% the Σ − Λ DU process onsets causing a rapid temperature drop, while the np DU process takes place for much higher DM fractions (see light green curves in Fig. 6).

Thermal evolution of Cassiopeia A as a DM-admixed NS
The central compact object (CCO) in Cassiopeia A (Cas A) has been raising various discussions due to its unusually rapid cooling [129][130][131][132][133]. Recent combined analysis of the Xray spectra of Cas A suggests that the surface temperature drop is 2.2 ± 0.3 % over 10 years when the absorbing hydrogen column density N H is used as a free model parameter, and 1.6 ± 0.2 % for a constant N H value [134].The thermal evolution of such a young NS, with an estimated age around 356 years, opened up a pandora box of models, between which the 3 P 2 pairing of neutrons in the core [35,135], rapid cooling via the DU process [136], the impact of medium-modified one-pion exchange in dense matter [137] and beyond the Standard Model physics [138].As the estimated mass of the CCO in Cas A is around 1.55 M ⊙ [134], a realistic description should take place for a middle mass star.
As it was shown in Ref. [139], the IST EoS does not necessarily require the incorporation of the neutron superfluidity and/or proton conductivity to explain the Cas A temperature drop.The model reproduced the Cas A data with 1.66 M ⊙ and 1.91 M ⊙ stars with the inclusion of neutron and proton singlet pairings, as well as with the 1.96 M ⊙ star for unpaired matter [140].As it was shown in [90] the FSU2R EoS describes the Cas A cooling with a combination of n 1 S 0 (SFB model), p 1 S 0 (CCDK model), n 3 P 2 pairing (T72 model) with the maximum critical temperature T c = 7.105 • 10 8 K.In Ref. [128], the best fit of Cas A within the FSU2H EoS is obtained for the 1.88 M ⊙ star considering the proton and neutron single pairing together with the triplet neutron pairing characterized by the maximum critical temperature T c ∼ 1.41 • 10 9 K.For the same EoS, we also see a better agreement with Cas A data while varying the neutron 3 P 2 pairing.
As one of the results of our analysis, the best fit is obtained for the DM-admixed star of 1.2 M ⊙ with a DM fraction of 3%.The results of the fit are presented in Fig. 7 whereas the curves are obtained for the heavy-elements envelope η = ∆M/M = 10 −16 and a combination of n 1 S 0 (SFB model [126]), p 1 S 0 (CCDK model [127]), and n 3 P 2 pairing (T72 model [141]) with the maximum critical temperature T c ∼ 8.52 • 10 8 K.

Conclusions
Compact star cooling is an important probe of the particle composition and physics governing the internal dynamics.In this work, we extend a study of Ávila et al. [90] of the impact of asymmetric fermionic DM on the thermal evolution of NSs.Within this scenario, DM does not self-annihilate due to the particle-antiparticle asymmetry and, therefore, is accrued in the star's interior.We demonstrate that considering its interaction with BM through gravity the accumulation of DM leads not only to the well-known modification of the star's gravitational mass, radius, and tidal deformability but also affects its thermal evolution.Thus, accrued DM indirectly alters the cooling behavior of NSs by pulling inward BM from the outer layers, increasing the central density and consequently affecting the BM distribution.The latter triggers the onset of nucleonic and hyperonic DU processes in lower NS masses and rapid star cooling due to the enhanced neutrino emission.At the same time, the proton fraction corresponding to the nucleonic DU onset remains equivalent to the one of a pure BM star.The same occurs with the onset density of hyperons.These onset densities correspond to the NSs of smaller masses due to their compactification.) and light elements (η = ∆M/M = 10 −7 ), respectively.The impact of neutron superfluidity in the 1 S 0 channel, employing the SFB model [126], and proton superconductivity in the 1 S 0 channel, employing the CCDK model [127], is considered.The figure is adopted from [90].The utilized observational data are listed in the Appendix.
Another consequence of the gravitational pull of BM towards the center results in a reduction of the baryonic radius of the star.Therefore, the stellar surface decreases, leading to a lower photon luminosity.This effect becomes visible during the photon-dominated stage when neutrino emission plays a secondary role.
To thoroughly study this effect, we considered the nucleonic IST and FSU2R models that have the DU onset at 1.908 M ⊙ and 1.921 M ⊙ , respectively, and the hyperonic FSU2H model that kinematically allows both nucleonic and hyperonic DU processes.Accumulation of the DM particles with a mass of m χ = 1 GeV of f χ ≃ 0.161% (IST EoS) and f χ = 0.378% (FSU2R EoS) triggers the previously energetically forbidden process.
At the same time, the FSU2H EoS gives the possibility to study the sequence of the DU onsets with an increase of the DM fraction, both due to nucleonic and hyperonic DU processes.Thus, for the 1.7M ⊙ star with 0% of DM, only the Λp DU process operates leading to a medium cooling.At f DM ≃ 0.61% the Σ − Λ DU process onsets causing a rapid temperature drop, while the final contribution of np DU process occurs for f DM > 2%.
We show that the rapid cooling of the CCO in Cas A could be explained by the DMadmixed star of 1.2 M ⊙ described by the FSU2H EoS.In [90] the best fit of the Cas A cooling rate was obtained for the M = 1.6 M ⊙ star modeled within the FSU2R EoS with f DM = 4% and light-elements envelope.Thus, the accrued DM helps to reconcile the star mass at which the DU process is kinematically allowed with the observational data on the mass of the CCO in Cas A by lowering the mass compared to a pure BM star.We demonstrate that an increase of the DM fraction causes a shift of the DU onset towards lower gravitational masses of the star.This effect could serve as a distinctive signature of the presence of DM in compact stars.In Ref. [90] it was shown that a similar result can be obtained by considering heavier DM particles, leaving the DM fraction to be low.
Thus, low/middle mass NSs that are not expected to have the operated DU process, might show a rapid cooling due to the presence of DM.Consequently, stars of similar masses will display different cooling patterns depending on the accrued DM fraction in their interiors.As DM fraction increases towards the Galactic center, the thermal evolution of NSs might exhibit distinct features compared to stars further away, e.g. in the Solar vicinity and beyond.However, the measurements of the distant NSs are not so precise.Therefore, several radio telescopes such as FAST [142], SKA [143], and CHIME [144] will search for the old, isolated compact stars near the Earth.The recently launched JWST together with the forthcoming Thirty Meter Telescope (TMT) and Extremely Large Telescope (ELT) [145] Figure 6.The same as in Fig. 5, but for stars with the total gravitational masses of M = 1.2, 1.5, 1.7 M ⊙ modeled within the FSU2H EoS and supplemented with n 3 P 2 pairing described by the T72 model.
will provide a sensitive probe of old NSs [65].These searches together with studies of the NS matter EoS are important for understanding the interior composition of compact stars and possible degeneracy between the effects of strongly interacting matter and DM [8].

Appendix: Observational data
Observational data of the CCO in the Cas A supernova remnant are depicted as source 0 in Figs.5,6,7.It shows the temperature measurements using Chandra ACIS in GRADED and FAINT modes with 1σ error bars [35].The surface temperature is subject to variation across different observations and the hydrogen column density N H may also vary [134], leading to different T ∞ s values.Thus, the plotted data are the outcome of examining both the GRADED and FAINT modes mode of the ACIS detector.One data set involves varying the hydrogen column density, while the other considers it fixed.

Figure 1 .
Figure 1.Left panel: The proton fraction, Y p , (top) and total gravitational mass of the DM-admixed and baryonic NSs (bottom) as a function of the central baryonic density, n B,c , for the IST and FSU2R EoSs.On the bottom panel, the royal blue, red, and green curves represent the relative DM fractions equal to 2%, 3%, and 4%, respectively.The vertical solid and dashed gray lines correspond to the central BM densities of the stars at which the DU process is activated for the IST and FSU2R EoSs, respectively.The intersection points depict the nucleonic DU threshold.Right panel: The particle fractions, Y j , (top) and total gravitational mass of the DM-admixed and baryonic NSs (bottom) as a function of the central baryonic density, n B,c , for the FSU2H EoS.The vertical solid, dashed, and dashed-dotted gray lines correspond to the central BM densities at which the DU processes are activated, i.e. nucleonic (np) and hyperonic (Λp and Σ − Λ), respectively.

Figure 2 .
Figure 2.Total gravitational mass of DM-admixed and baryonic NSs as a function of their baryonic radius R B calculated for the DM particle mass m DM =1 GeV.Solid curves correspond to pure BM stars described by the IST (blue curve), FSU2R (red curve), and FSU2H (green curve) EoSs.Dashed and dotted curves depict the M-R relations obtained for the relative DM fractions equal to 2% and 4%, respectively.Orange and dark yellow bands represent 1σ constraints on the mass of PSR J0348+0432[52] and PSR J1810+1744[110].While olive green and light orange contours show the NICER measurements of PSR J0030+0451[111,112], purple and magenta contours correspond to the PSR J0740+6620 measurement[79,113].Observations of GW170817[114] and GW190425[115] binary NS mergers by LIGO-Virgo collaboration are shown in blue and green.The 2σ contour of HESS J1731-347[116] is plotted in light red.The shaded gray region is excluded by the rotation of the fastest spinning pulsar PSR J1748-2446ad[117].

Figure 3 .
Figure 3. Stellar configurations with different DM fractions for the IST EoS (top right slice), FSU2R EoS (top left slice), and FSU2H EoS (bottom slice).The size of the np, Σ − Λ, Λp DU regions, and DM core are depicted in light red, light green, yellow, and dark gray, respectively.For a better comparison, the radii are normalized to the baryonic radius of each configuration and are given in Table4.All configurations correspond to NSs with a total gravitational mass of 1.75 M ⊙ .

Figure 4 .
Figure 4.The same as in Fig.3, but for the total gravitational mass of 1.9 M ⊙ .The empty part of the bottom slice corresponds to the non-existing configurations (see details in the text).

Figure 5 .
Figure 5. Cooling curves for stars with total gravitational masses of M = 1.2, 1.6, and 1.9 M ⊙ described by the IST EoS (left) and FSU2R EoS (right) are shown.The considered DM fraction is f DM = 0%, 2%, 3% and 4% for particle's mass m DM = 1 GeV.The solid and dashed curves correspond to envelopes composed of heavy elements (η = ∆M/M = 10 −16) and light elements (η = ∆M/M = 10 −7 ), respectively.The impact of neutron superfluidity in the 1 S 0 channel, employing the SFB model[126], and proton superconductivity in the 1 S 0 channel, employing the CCDK model[127], is considered.The figure is adopted from[90].The utilized observational data are listed in the Appendix.

Figure 7 .
Figure 7.The results of the fit of the surface temperature as a function of age for the CCO in Cas A measured by Chandra ACIS-S in GRADED and FAINT modes.The red and magenta data points correspond to variable and fixed absorbing hydrogen column density N H = 1.656 • 10 22 cm −2 in the FAINT mode, while the green and blue data points depict the same data, but in the GRADED mode.
km] Table

Table 2 .
Baryonic densities and total gravitational masses at the onset of the nucleonic DU process for the IST and FSU2R EoSs.

Table 3 .
The same as in Table2, but for the FSU2H model.The columns correspond to different hyperonic and nucleonic DU onsets.