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 fermioinc 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 three pulsars namely PSR B0656+14, Geminga and PSR B1055-52. 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, we found that in case of medium and heavier mass NSs, all chosen DM masses and fermi momenta agree well with the observational data but for lower mass NSs, all DM fermi momenta and high DM masses barely agree with the observations. Furthermore, only lower DM mass agrees with observations in case of lighter NSs. 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], quarks [5] and possisbly the dark matter (DM) particles [6]- [8] are believed to be present inside the core. 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 evidences 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]- [23], Axions [24,25], Febbly Interating Massive Particles (FIMPs) [26,27], Fuzzy dark matter [28,29], neutralino [30], 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 [33]. It has been proven in [34] 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 [35]. It has been shown that fermionic DM could soften the equation of state and hence reduce the maximum mass supported by the NS [33]. This effect is sensitive to the mass of DM particle and the self-interaction 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 [36]- [38]. However, very little study has been done on the effect of DM on cooling of NSs [39]. 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 cooling data of PSR B0656+14, Geminga and PSR B1055-52 [40] in our work . We studied NS cooling of both normal NSs using density-dependent (DD2) EOS [41,42] and Akmal-Pandharipande-Ravenhall (APR) EOS [43] and DM admixed NSs using DD2 EOS modified with DM sector. It is important to mention here although DD2 is marginally allowed by the tidal deformability constraint obtained from the the analysis of GW170817 with Phenom PNRT model [44], 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]- [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 [42,45] in the nucleon chemical potential. This paper is organised as follows. In section 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 section 3. In section 4, we discuss the cooling mechanism of neutron star. Further, the results and calculations are presented in section 5. Finally, we conclude our work in section 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 [41,42] is given by where ψ B denotes the nucleon fields, τ B is the isospin operator and gs represent the densitydependent meson-baryon couplings. These couplings are determined by following the prescription adopted by Typel et al. [42,46]. The functional dependence of the couplings on density was introduced for the first time in [47] 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 [47] 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(n 0 ) and the coefficients a α , b α , c α and d α [42,46]. The fit provides the saturation density n 0 = 0.149065f 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 Section 2.1.  Table 1: Meson masses and parameters of meson-nucleon couplings in DD2 EOS.

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 Section 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.35 [48] 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 assumptions that the average density of nuclear matter inside the neutron star is 10 3 times greater than the average dark matter density (ρ DM ) and the fermi momentum of dark matter is constant [49] 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 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 ψ = ψ p ψ n , scalar meson (σ), vector meson (ω µ ) and isovector meson (ρ µ ), DM particle (χ) and Higgs boson h can be derived from Eq. (1) and Eq. (4) as where masses of DM paticle and Higgs particle are denoted by M χ and M h = 125.09 GeV respectively. Σ B = Σ 0 B + Σ r B 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 term i.e. Σ r which appears because of density-dependence of meson-nucleon couplings [45]. 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 [42,46]. 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 µ n = µ p + µ e , respectively. Here, the chemical potentials µ e and µ µ are given as whereas the nucleon chemical potentials contain the rearrangement term also because of densitydependence 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) [50] equations of hydrosatic equilibrium to generate the mass-radius and pressure-radius profiles as shown in Figures 2-3. In Figure 1, we present EOSs for different DM masses and fermi momenta along with APR and DD2 EOSs. 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 Figure  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 χ = 500GeV, the EOS corresponding to k DM F = 0.06 becomes softest among all the cases of DM masses. This sudden softening of EOS for k DM F = 0.06 at M χ = 500GeV might be due to dominance of dark matter over baryonic matter at such extreme parameters of DM and consequently, dominance of gravitational attraction over internucleon repulsion leading to a reduced pressure. In Figure 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 Figure 1. Here also mass-radius profile for k DM F = 0.06 at M χ = 500GeV follows a trend contrary to other combinations of k DM F and M χ . In this case the majority of mass contribution is from DM and hence leading to smaller radius than other cases due to enhanced gravitational contraction. Figure 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 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( mq v ) . 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 = MnMχ Mn+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 [51], PandaX-II [52], LUX [53] and DarkSide-50 [54] experimental bounds for constraining the parameter "y". We accordingly fix "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 [55] 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 zero. The photon luminosity is calculated using the Stefan-Boltzmann law [56] 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 [57] is utilised in the present work for calculating the neutrino and photon luminosities. Several neutrino emitting processes contribute in the cooling of neutron stars [55,58,59]. 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 . Modified Urca process will become more dominant provided the proton fraction is below the threshold. The modified 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 . Cooper pairing of nucleons are the other set of neutrino emiitting processes as n + n → [nn] + ν + ν, p + p → [pp] + ν + ν. These are medium processes where luminosity varies with temperature as L f ast ν ∝ T 7 9 . There are several other neutrino emitting processes involved in the cooling as follows e − + e + → ν + ν (electron-positron pair annihilation), e − → e − + ν + ν (electron synchrotron), γ + e − → e − + ν + ν (photoneutrino emission), e − + Z → e − + Z + ν + ν (electron-nucleus bremsstrahlung), n + n → n + n + ν + ν (neutron-neutron bremsstrahlung) and n + p → n + p + ν + ν (neutron-nucleus bremsstrahlung).

Calculations and Results
We utilised NSCool Numerical code for studying cooling of NSs adopting different EOSs like APR, DD2 and DM admixed DD2. We considered different neutron star masses namely 1.0 M , 1.4 M and 2.0 M for the calculations. In case of DM admixed DD2, we explored the effect of variation of DM mass (50 GeV,200 GeV and 500 GeV) and DM fermi momentum k DM F (0.0 GeV, 0.01 GeV, 0.02 GeV, 0.03 GeV and 0.06 GeV) on the cooling of NSs. It is important to mention here that k DM F = 0 GeV means dark matter density is zero but effective mass of nucleons will be effected due to non-zero Higgs-nucleon Yukawa coupling (Eq. (4)). For demonstrating DM-effect on neutron star cooling we plot the variations of luminosity with time (Figures 4-6) and temperature with time (Figures 7-9). As seen in all the plots (Figures 4-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 Figures 4-6, luminosity versus time are plotted for different NS masses and for every NS mass diffferent EOSs are considered. For lower NS masses, cooling with DD2 is fastest and with APR it is slowest. Moreover, for fixed DM mass, cooling is faster for higher values of DM fermi momentum. But in heavier NSs, cooling with APR becomes fastest which might be due to appearance of extra neutino emitting channels and variation due to k DM F is the same as previously. For both the medium mass as well as heavier DM admixed NSs all chosen DM fermi momenta are considerably consistent with cooling data of pulsars namely PSR B0656+14 and Geminga but for lower mass NS all DM fermi momenta barely agree with the observations (PSR B0656+14, Geminga). The effect of k DM F on cooling of heavier NSs becomes more and more prominent at higher values of DM mass as is evident from the right most panels of Figures 5 and 6. Figures 7-9

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 three pulsars namely PSR B0656+14, Geminga and PSR B1055-52. We found that APR EOS agrees well with the pulsar data for lighter and medium mass NSs whereas DD2 agrees well for medium and heavier mass NSs and marginally for low mass NSs. For DM admixed DD2 EOS, we found that in case of medium and heavier NSs, all chosen DM fermi momenta agree well with the observational data but for lower mass NS all DM fermi momenta barely agree with the observations. Furthermore, it is found that all chosen DM masses agree well with the observations in case of medium and heavier NSs but only lower DM mass agrees with observations in case of lighter NSs. 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.