Cooling of dark-matter admixed neutron stars with density-dependent equation of state

We propose a dark-matter (DM) admixed density-dependent equation of state where the fermionic DM interacts with the nucleons via Higgs portal. Presence of DM can hardly influence the particle distribution inside neutron star (NS) but can significantly affect the structure as well as equation of state (EOS) of NS. Introduction of DM inside NS softens the equation of state. We explored the effect of variation of DM mass and DM Fermi momentum on the NS EOS. Moreover, DM-Higgs coupling is constrained using dark matter direct detection experiments. Then, we studied cooling of normal NSs using APR and DD2 EOSs and DM admixed NSs using dark-matter modified DD2 with varying DM mass and Fermi momentum. We have done our analysis by considering different NS masses. Also DM mass and DM Fermi momentum are varied for fixed NS mass and DM-Higgs coupling. We calculated the variations of luminosity and temperature of NS with time for all EOSs considered in our work and then compared our calculations with the observed astronomical cooling data of pulsars namely Cas A, RX J0822-43, 1E 1207-52, RX J0002+62, XMMU J17328, PSR B1706-44, Vela, PSR B2334+61, PSR B0656+14, Geminga, PSR B1055-52 and RX J0720.4-3125. It is found that APR EOS agrees well with the pulsar data for lighter and medium mass NSs but cooling is very fast for heavier NS. For DM admixed DD2 EOS, it is found that for all considered NS masses, all chosen DM masses and Fermi momenta agree well with the observational data of PSR B0656+14, Geminga, Vela, PSR B1706-44 and PSR B2334+61. Cooling becomes faster as compared to normal NSs in case of increasing DM mass and Fermi momenta. It is infered from the calculations that if low mass super cold NSs are observed in future that may support the fact that heavier WIMP can be present inside neutron stars.


Introduction
Neutron stars are excellent celestial laboratories for investigating the supradense nuclear matter which is otherwise inaccessible to terrestrial laboratories. The density inside neutron stars is several times the saturated nuclear density, hence exotic particles like hyperons [1,2], pion or kaon condensate [3,4] and quarks [5] are believed to be present inside the core. Dark matter (DM) particles may also be captured and accumulated inside neutron stars [6][7][8]. Exotic particles soften the equation of state (EOS) and reduce the tidal deformability of the neutron star [9]. Exact nature of the matter is a challenging task and yet to be known. Any model suggested should not only describe the superdense matter but also reproduce the properties of matter observed at saturation densities [10,11]. Recently, the unprecedented joint detection of neutron star merger GW170817 by Advanced LIGO and Virgo observatories has put stronger constraints on the equation of state by constraining tidal deformability of NSs [12,13]. Using Shapiro delay measurements, a very massive neutron star has been found in the form of PSRJ0740+6620 with mass 2.17 + 0. 11 − 0.10 [14]. This can put stringent constraint on the equation of state.
Nowadays there are various cosmological and astrophysical indications for the existence of dark matter in the Universe like rotation curves of spiral galaxies, large-scale structures of the Universe, anisotropies of cosmic microwave background radiation (CMBR), gravitational lensing etc. The detection of dark matter is attempted following three different ways i.e. direct detection, indirect detection and collider searches (LHC). However, till now no experimental signature of dark matter has been discovered. Direct detection experiments put upper bounds on the dark matter-nucleon elastic scattering cross-sections for different DM masses. In the literature, many theoretical particle dark matter models are proposed to indirectly detect the dark matter and to explain the existence of few unsolved phenomenological evi-dences such as gamma ray excesses observed by Fermi-LAT gamma ray telescope [15,16], positron excesses measured by PAMELA [17], AMS-02 [18], DAMPE [19] experiments etc. Till now many particle candidates of dark matter are proposed like Weakly Interacting Massive Particles (WIMPs) [20][21][22][23][24], Axions [25,26], Febbly Interating Massive Particles (FIMPs) [27,28], Fuzzy dark matter [29,30], neutralino [20], Kaluza Klein dark matter [31] etc. In this work our proposed particle candidate of dark matter is WIMP. In the early Universe, WIMPs are produced thermally and initially they are at thermal equilibrium but when the temperature drops below the WIMPs mass they are decoupled at a particular temperature ∼ M χ 20 called freeze-out temperature. After decoupling, WIMP would possibly be a relic particle and may constitute a particle candidate of cold dark matter (CDM). WIMPs can cluster with stars gravitationally and also form a background density in the universe. Several studies have indicated that neutron stars being highly compact objets can capture more dark matter particles during the formation stage in the supernova explosion as compared to the non-compact objects [32]. Recently, it is shown that the admixture of DM inside NSs softens the equation of state and hence tidal deformability is reduced [7]. It has been proven in [33] that the DM capture could be highly improved if it happens in binary pulsars. Since the DM present inside DM admixed NSs can possibly change the global properties of neutron stars, this open another indirect window to study DM apart from the other numerous ways. The structures of DM admixed NSs have been studied recently. It is shown that mass and radius of NSs can be remarkably affected by mirror DM [34]. It has been shown that fermionic DM could soften the equation of state and hence reduce the maximum mass supported by the NS [7]. This effect is sensitive to the mass of DM particle and the selfinteraction within the dark matter. Since the normal matter and DM are believed to interact gravitationally, presence of DM can hardly influence the the particle distribution inside NS but can significantly affect the structure as well as EOS of NS.
Cooling of neutron stars have been well studied by several authors [35][36][37][38][39]. Some study has been done on the effect of DM on cooling of NSs [40][41][42]. It has been found that the heating due to dark matter annihilation can affect the temperature of the stars older than 10 7 years and consequently flattening out the temperature at 10 4 K for the neutron stars [41]. Moreover, recently it has been found that slowdown in the pulsar rotation can drive the NS matter out of beta equilibrium and the resultant imbalance in chemical potentials can induce late-time heating, named as rotochemical heating which can heat a NS up to 10 6 K for t = 10 6 −10 7 years [42]. In Ref. [40], the authors have studied the cooling of DM admixed NS with dark matter mass ranging from 0.1 GeV to 1.3 GeV. In the present work, we have considered low as well as high dark matter masses (upto 500 GeV) and also varied the dark matter Fermi momenta for the cooling calculations. For these calculations, we have considered dark-matter modified density-dependent (DD2) EOS [43,44] and the results are compared with the observational data. Also, our work is different from the above mentioned works [41,42] because we don't consider heating due to WIMP annihilation owing to the very small annihilation cross-section. We consider indirect effect on cooling of NS stars due to change in the neutron star structure in presence of dark matter. With the introduction of dark matter, cooling properties can change significantly as compared to the normal NSs mainly because of changes in neutrino emissivity, neutrino luminosity and heat capacity. For given mass, neutron emissivity will be different due to significant change in stellar structures and consequently, neutrino luminosity will also be different. Heat capacity related to EOS will be different for normal NSs and DM admixed NSs because of softening of the EOS in case of the latter. Thus, normal NSs can be distinguished from DM admixed NSs using astronomical observation data related to surface temperature and age of pulsars. We have considered the complete set of cooling data of both young and cool and old and warm neutron stars namely Cas A, RX J0822-43, 1E 1207-52, RX J0002+62, XMMU J17328, PSR B1706-44, Vela, PSR B2334+61, PSR B0656+14, Geminga, PSR B1055-52 and RX J0720.4-3125. We adopted the temperature data sets of the above mentioned pulsars from the Refs. [7,45,48] and the luminosity data from the Ref. [49] (Table 1 and 3). As representatives for late time-cooling, a group of above mentioned three NSs (PSR B0656+14, Geminga and PSR B1055-52) is chosen forming a class of nearby objects that allows spectral fits to their X-ray emission ( [50][51][52][53] and references therein). We studied NS cooling of both normal NSs using DD2 EOS [43,44] and Akmal-Pandharipande-Ravenhall (APR) EOS [54] and DM admixed NSs using DD2 EOS modified with DM sector. It is important to mention here that although DD2 is marginally allowed by the tidal deformability constraint obtained from the analysis of GW170817 with Phenom PNRT model [55], DM admixed DD2 will be softened and might be considerably allowed by the GW170817 constraints. Earlier DM admixed NSs are studied by some groups [6][7][8] where they adopted σ -ω-ρ model but our approach differs from theirs in the sense that meson-nucleon couplings are density-dependent in our model which gives rise to an extra term called rearrangement term [44,56] in the nucleon chemical potential. This paper is organised as follows. In Sect. 2, we describe baryonic EOS model and DM admixed baryonic EOS model. We constrain dark matter-Higgs coupling parameter from the direct detection experiments as discussed in Sect. 3. In Sect. 4, we discuss the cooling mechanism of neutron star. Furthermore, the results and calculations are presented in Sect. 5. Finally, we conclude our work in Sect. 6.

Equation of state model
In this section, we utilize density-dependent relativistic hadron field theory for describing strongly interacting superdense nuclear matter inside neutron stars. Nucleon-nucleon strong interaction is mediated by the exchanges of scalar σ meson, responsible for strong attractive force, vector ω, responsible for strong repulsive force and ρ meson, responsible for symmetry energy. The Lagrangian density [43,44] is given by where ψ B denotes the nucleon fields, τ B is the isospin operator and gs represent the density-dependent meson-baryon couplings. These couplings are determined by following the prescription adopted by Typel et al. [44,57]. The functional dependence of the couplings on density was introduced for the first time in [58] as described here where ρ b is the total baryon density, 1+c α (x+d α ) for α = ω, σ . In order to reduce parameters, the functions are constrained as f σ (1) [58] is considered for the isovector meson ρ μ because g ρ B decreases at higher densities. Finite nuclear properties are fitted to determine the saturation density, the mass of σ meson, the couplings g α B (ρ 0 ) and the coefficients a α , b α , c α and d α [44,57]. The fit provides the saturation density ρ 0 = 0.149065 f m −3 , binding energy per nucleon as −16.02 MeV and compressibility factor K = 242.7 MeV. The nucleon mass M n is considered to be 939 MeV through out our work.
Leptons are treated as non-interacting particles and described by the free Lagrangian density where l = e − , μ − and m l = m e , m μ . The energy and pressure due to leptons will be explicitly mentioned in Sect. 2.1.

Effect of dark matter on equation of state
A uniformly distributed fermionic dark matter (WIMP) is considered inside neutron star. Dark matter interacts with Higgs field h with coupling strength y. DM-Higgs coupling y is explicitly discussed in Sect. 3. Three different WIMP masses (M χ = 50 GeV, 200 GeV, 500 GeV) are considered in our calculations. Higgs field h interacts with the nucleons via effective Yukawa coupling f M n /v, where f denotes the nucleon-Higgs form factor and is estimated to be approximately 0.3 [59] and v = 246.22 GeV denotes Higgs vacuum expectation value (VEV). In the Higgs potential, terms higher than quadratic are dropped because they are negligible in the mean field approximation (MFA). Hence the dark sector and its interaction with nucleons and Higgs field is described by the Lagrangian density Here we consider the assumption that the average dark matter number density inside neutron star is 10 3 times smaller than saturated nuclear matter number density [6,60] and the Fermi momentum of dark matter is constant [6] through out the neutron star. With these assumptions, the fractional mass of dark matter inside neutron star for M χ = 200 GeV can be expressed as .
which gives k DM F ∼ 0.033 GeV. We vary k DM F in our calculations from 0.01 GeV to 0.06 GeV and dark matter densties ρ DM will also vary accordingly. Equations of motion for nucleon doublet scalar meson (σ ), vector meson (ω μ ) and isovector meson (ρ μ ), DM particle (χ ) and Higgs boson h can be derived from Eqs. (1) and (4) as where masses of DM paticle and Higgs particle are denoted by M χ and M h = 125.09 GeV respectively.
is the vector self energy in which the first term consists of the usual non-vanishing components of vector mesons i.e. 0 B = g ωB ω 0 −g ρ B τ 3B ρ 03 and the second term is the rearrangement ρ B which appears because of density-dependence of mesonnucleon couplings [56]. Here g α B = ∂g α B ∂ρ B where α = σ , ω, ρ and τ 3B is the isospin projection of B = n, p. Due to density dependence of nucleon-meson couplings, chemical potential of the nucleons takes the form In the mean-field approximation (MFA), fields are replaced by their expectation values and above equations are simplified as The effective masses of nucleons and dark matter are respectively given as The baryon density (ρ), scalar density (ρ s ) and dark matter density (ρ DM s ) are where k F and k DM F are the Fermi momenta for nucleonic matter and dark matter respectively and γ = 2 is the spin degeneracy factor of nucleons. The masses of mesons and meson-nucleon couplings at saturation density ρ 0 are given in Table 1 [44,57]. In order to get the density dependent profile for M * n and M * χ , Eqs. (6) and (8) should be solved self consistently. The energy and pressure i.e. EOS are provided by expectation values of energy-momentum tensor in the static case as = T 00 and P = 1 3 T ii . The total energy density and pressure for the combined Lagrangian L B + L DM are obtained as where ρ n and ρ p are the neutron and proton number densities and k n F and k p F are the corresponding Fermi momenta of neutron and proton, respectively. The nuclear matter inside the neutron star will be charge neutral and β-equilibrated. The conditions of charge neutrality and β-equilibrium are given as and whereas the nucleon chemical potentials contain the rearrangement term also because of density-dependence of couplings as mentioned earlier. The particle fractions of neutron, proton, electron and muon will be determined by the self consistent solution of Eqs. (11) and (12) for a given baryon density. The energy density and pressure due to the non-interacting leptons are given as So the total energy density and pressure of the charge neutral β-equilibrated neutron star matter are For all the EOSs considered in our work, we solve numerically Tolman-Oppenheimer-Volkoff (TOV) [61] equations of hydrosatic equilibrium to generate the mass-radius and pressure-radius profiles as shown in Figs. 2, 3. In Fig. 1, we present EOSs for different DM masses and Fermi momenta along with APR and DD2 EOSs. It is important to mention that all of these EOSs satisfy the causality condition i.e. c 2 s < 1. It is evident that for a fixed DM-Higgs coupling y and fixed DM mass, EOS becomes softer for higher values of DM Fermi momentum and is softest for k DM F = 0.06 GeV for lower to moderate values of density and APR is stiffest for higher values of density. Moreover, it is inferred from the comparison of three panels of Fig. 1 that for fixed values of y and k DM F , higher values of DM masses leads to softer EOS. It is important to mention here that for the higher DM mass M χ = 500 GeV, the EOS corresponding to k DM F = 0.06 GeV becomes softest among all the cases of DM masses. This sudden softening of EOS for k DM F = 0.06 GeV at M χ = 500GeV might be due to dominance of dark matter over baryonic matter at such extreme parameters of DM. Nevertheless all the neutron star configurations for the dark matter densities considered in this work are stable and will not undergo black hole formation [60]. In Fig. 2, mass-radius profile is plotted for NS masses 1M , 1.4M and 2M . These plots can be explained the same way as in case of  contribution is from DM and hence leading to smaller radius than other cases due to enhanced gravitational contraction. Fig. 3 shows the pressure-radius profile for different NS masses where it is evident that for fixed DM mass and NS mass, higher values of k DM F leads to lower pressure except for the case of M χ = 500 GeV and M N S = 2M where k DM F = 0.06 GeV leads to higher pressure in the inner region of the star. This is because the star becomes more centrally condensed at very high DM mass.

Direct detection
Dark matter direct detection experiments do not show any signatures of collision events as of now. These experiments give an upper bound on the elastic scattering cross-section as a function of dark matter mass. In the present scenario, fermionic dark matter (WIMP) can undergo an elastic collision with the detector nucleus (quark level) by the Higgs exchange. Therefore, the effective Lagrangian contains the scalar operatorχχqq and can be written as where q represent the valence quarks and α q = y f ( This scalar operator contributes to the spin independent (SI) scattering cross-section for the fermionic dark matter candidate and can be expression as In Eq. (19), m r = M n M χ M n +M χ is the reduced mass. We calculate the SI scattering cross-section using Eq. (19) and then constrain the parameter "y" using the direct detection experiments in such a way that calculated scattering cross-section for different dark matter masses are below the experimental bounds. In the present scenario, we use XENON-1T [62], PandaX-II [63], LUX [64] and DarkSide-50 [65] experimental bounds for constraining the parameter "y". We checked that by varying the value of the parameter "y" in the cooling calculations, no significant differences are found and hence, we accordingly fixed "y" to be 0.01 and calculate the corresponding scattering cross-section for three chosen DM masses as tabulated in Table 2.

Cooling mechanism of neutron stars
It is well known that the surface temperature of the neutron star decreases with time which is the direct indication of cooling. In order to calculate the thermal evolution of the neutron star, one needs to solve the the energy balance equation for the neutron star which can be expressed as [66] where E th represents the thermal energy content of the star, T and T e are the internal and effective temperatures of the star, respectively, C v is the heat capacity of the core and H is the source term which includes different "heating mechanisms" important in the later stage of neutron star evolution. In Eq. (20), L ν and L γ denote the neutrino and photon luminosities, respectively. H (T ) is considered here to be zero. The photon luminosity is calculated using the Stefan-Boltzmann law [67] as This relation is obtained using T e ∝ T 0.5+α (α 1), where R is the radius of the star and σ denotes the Stefan-Boltzmann constant. NSCool code [68] is utilised in the present work for calculating the neutrino and photon luminosities.
Several neutrino emitting processes contribute in the cooling of neutron stars [53,66,69]. Direct Urca processes and modified Urca processes are the two main neutrino emitting processes for the cooling. The direct Urca processes are n → p + e − + ν e , p + e − → n + ν e .
These are possible in neutron stars only if the proton fraction crosses a critical threshold. The two processes are fast and the luminosity varies with the temperature as L f ast ν ∝ T 6 9 . The baryon direct Urca processes are the strongest and open for some EOSs satisfying the threshold condition at sufficiently high densities. Cooling of neutron stars becomes rapid (enhanced) if these processes are allowed. If these processes are forbidden, then baryon modified Urca and baryon bremsstrahlung processes are the main neutrino emitting processes which lead to slow (standard) cooling. These reactions are abundant and the number of this type of open reactions grows quickly as density increases. However, one should not be bothered about this wealth of reactions because as the density increases direct Urca process will be determining the neutrino luminosity since the modified Urca and bremsstrahlung processes are negligible compared to it. Inside the neutron star core, direct Urca is the basic process bringing nucleons into the state of beta-equilibrium with chemical potentials satisfying the equality μ n = μ p + μ e . In equilibrium both reactions take place at the same rate, and the direct Urca process leaves the nucleon composition of matter unchanged. Threshold is the most important feature of direct Urca process. Since the process is several orders of magnitude more efficient than other neutrino processes, it is very important to know the threshold exactly. The reaction is allowed only if the Fermi momenta k n F , k p F and k e F satisfy the triangle condition (can be sides of one triangle). In other words, the value of each Fermi momentum should be smaller than the sum of two others. In neutron star matter k n F is larger than k p F and k e F , and the triangle condition reads: k n F < k p F + k e F . For ρ ∼ ρ 0 one typically has k n F ∼ 340MeV/c, k e F ∼ k p F ∼ (60−100) MeV/c and the condition is invalid, i.e., the direct Urca is forbidden [69]. However, k p F and k e F may grow with density faster than k n F , and the process can be open at higher densities, a few times ρ 0 . The direct Urca process in the npe matter is allowed for the EOSs having large symmetry energy at densities several times the nuclear saturation density. For APR EOS, the direct Urca process is allowed only for M 1.97M [36].
Modified Urca process will become more dominant provided the proton fraction is below the threshold for direct Urca. The modified Urca processes are n + n → n + p + e − + ν e , n + p + e − → n + n + ν e , p + n → p + p + e − + ν e and p + p + e − → p + n + ν e .
These are slow processes and luminosity varies with the temperature as L slow ν ∝ T 8 9 . The neutron and proton branches of modified Urca process have similar emissivity with main difference in the threshold for the proton branch which is allowed at k n F < 3k p F + k e F . The neutron branch doesn't have the threshold. In the npe matter, this inequality is equivalent to k n F < 4k e F , i.e., to the proton fraction Y p exceeding the critical value Y cp = 1/65 = 0.0154 [69]. The latter condition is satisfied almost anywhere in the neutron star core. It can be violated only for the equations of state with very low symmetry energy at ρ ρ 0 , forbidding the proton branch in the outermost part of the core [69]. Contrary to the case of the direct Urca process, the emissivity increases smoothly from zero while the density exceeds the threshold value. The proton process is especially efficient at higher densities, near the threshold of the direct Urca process. Both branches of the modified Urca process are the leading standard (slow) neutrino generating mechanisms in non-superfluid neutron star cores, provided the direct Urca processes are forbidden.
For understanding the cooling mechanism of neutron stars, neutron star in Cassiopi A (Cas A) supernova remnant is very important. It is the youngest known thermally emitting isolated neutron star in our galaxy and is also the first neutron star for which cooling has been observed directly. Moreover, it is among very few isolated neutron stars whose age and surface temperatures are very well determined and hence is important in understanding the thermal evolution and interior properties of neutron stars. Nearly 20 years of monitoring of this neutron star, since it was discovered by Chandra X-ray Observatory in 1999, has shown that there is a decrease in it's temperature by 2-3% per decade. This cooling rate is significantly faster than that can be explained by standard neutron star cooling theories. This provides a strong evidence for the existence of superfluidity in neutron star cores [70][71][72]. Once the temperature is below the critical temperature of superfluid transition, neutron superfluidity and proton superconductivity opens a new channel for neutrino emission by continuous breaking and formation of cooper pairs. Hence the cooling of the neutron star is enhanced. This rapid cooling is expected to continue for several more decades. In future, Chandra observations will be more reliable and temperature measurements of this neutron star will be more accurate.
Free neutrons in the crust and both neutrons and protons in the core of a neutron star are likely to be in the superfluid (SF) state. We assume the singlet-state pairing of the protons, and either singlet-state or triplet-state pairing of the neutrons. In    Fig. 4 but for M χ = 500 GeV a uniform, low-density matter (near the core-crust interface) the neutron pairing is known to be of the singlet type, but it switches to the triplet-state type at higher densities. Once the superfluidity is taken into account, the neutrino emissivity is highly suppressed at low temperature. The reason of this suppression is the gap in the energy spectrum which suppresses the excitation around the Fermi surface [73]. But at the same time another important neutrino emitting channel so called the Cooper pair breaking and formation (PBF) comes in play. It is associated with the Cooper pair breaking due to the thermal disturbance and its subsequent reformation. The triplet pairing in the core dominates the neutron PBF. The PBF neutrino emitting processes are given as These are medium processes where luminosity varies with temperature as L medium ν ∝ T 7 9 . It has been shown earlier that superfluidity in the crust affects the cooling curves at the initial thermal relaxation stage, while superfluidity in the core accelerates cooling at later stages. The duration of the thermal relaxation stage is greatly reduced by the effect of superfluidity on heat capacity of free neutrons in the stellar crust [69]. When thermal relaxation is over and the isothermal state is established throughout the star, the cooling is mainly   Fig. 7 but for M χ = 500 GeV regulated by the neutrino luminosity and heat capacity of the stellar core. The neutrino processes and heat capacity of the crust cease to play a significant role except for the very lowmass stars with large crusts. However, it has been found that for some superfluid models the neutrino luminosity of the crust may affect the cooling for a short period of time during the transition from the neutrino cooling era to the photon era [69].

Calculations and results
We utilised NSCool Numerical code [68] for studying cooling of NSs adopting different EOSs like APR, DD2 and DM  7,8,9). As seen in all the plots (Figs. 4,5,6,7,8,9,10,11,12,13,14), shortly after birth, NS cooling becomes dominated by neutrino emitting processes as mentioned earlier. When the internal temperature has sufficiently dropped in nearly about 10 4 −10 5 year then the cooling is dominated by photon emission from the NSs surface. In Figs. 4, 5, 6, luminosity vs time are plotted for dif-   10,11,12,13), effect on cooling due to varying DM masses becomes more evident for lower mass NSs except for the case of k DM F = 0.06 GeV where variation due to DM mass is prominent even for heavier mass NS (Fig. 14) which is due to very high DM Fermi momentum. For all masses of DM admixed NSs, all chosen DM Fermi momenta are considerably consistent with cooling data of pulsars namely PSR B0656+14, Geminga, Vela, PSR B1706-44 and PSR B2334+61 as evident from the Figs. 4,5,6,7,8,9,10,11,12,13,14. It is noted that Cas A barely agrees with the cooling curves corresponding to higher values of dark matter Fermi momenta and DM mass. These observed pulsars might contain dark matter (WIMP) with lower to moderate mass. Furthermore, as seen from left most panels of Figs. 10, 11, 12, 13, it is evident that if small mass and super cold NSs are found in future astronomical cooling observations, we can say that heavier WIMPS may actually exist inside NSs.
For demonstrating the effect of superfluidity on cooling, we plotted luminosity vs time and temperature vs time profiles in Fig. 15 considering different pairing gap models inside the NSCool code because the actual value of neutron 3 P 2 gap is unknown [74]. In Fig. 15, 'pairing 0' means no pairing is considered and 'pairing a, b, c' correspond to [74] and the references therein. It is evident that cooling is faster when pairing is considered and among the three pairing models, cooling with 'pairing a' is slightly faster. Hence, we have considered 'pairing a' model throughout the work.

Summary and conclusions
In this work, we first prepared dark-matter admixed DD2 equation of state and explored the effect of dark matter mass and Fermi momentum on the neutron star equation of state. Dark matter-Higgs coupling is constrained using dark matter direct detection experiments namely XENON-1T, PandaX-II, LUX and DarkSide-50. Then, we studied cooling of normal NSs using APR and DD2 equation of states and DM admixed neutron stars using dark-matter modified DD2 with varying dark matter mass and Fermi momentum for fixed DM-Higgs coupling. We have done our analysis by considering three neutron star masses one each from the lighter (1.0 M ), medium (1.4 M ) and heavier (2.0 M ) NSs. We demonstrate our results by choosing three different DM masses namely 50 GeV, 200 GeV and 500 GeV and different Fermi momenta k DM F namely 0.01 GeV, 0.02 GeV, 0.03 GeV and 0.06 GeV. We calculated the variations of luminosity and temperature of the above mentioned neutron star masses with time and compared our calculations with the observed astronomical cooling data of pulsars namely Cas A, RX J0822-43, 1E 1207-52, RX J0002+62, XMMU J17328, PSR B1706-44, Vela, PSR B2334+61, PSR B0656+14, Geminga, PSR B1055-52 and RX J0720.4-3125. We found that APR EOS agrees well with the pulsar data for lighter and medium mass NSs whereas for DM admixed DD2 EOS, it is found for all considered NS masses, all chosen DM Fermi momenta are consistent with the observational data of PSR B0656+14, Geminga, Vela, PSR B1706-44 and PSR B2334+61. Cooling becomes faster as compared to normal NSs in case of increasing DM masses and Fermi momenta. It is observed from the calculations that if low mass super cold NSs are observed in future that may support the fact that heavier WIMP can be present inside neutron stars.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .