Scattering searches for dark matter in subhalos: Neutron stars, cosmic rays, and old rocks

Scattering Searches for Dark Matter in Subhalos: Neutron Stars, Cosmic Rays, and Old Rocks Joseph Bramante,1,2,* Bradley J. Kavanagh ,3,†and Nirmal Raj 4,‡ 1The Arthur B. McDonald Canadian Astroparticle Physics Research Institute, Department of Physics, Engineering Physics, and Astronomy, Queen’s University, Kingston, Ontario K7L 2S8, Canada 2Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada 3Instituto de Física de Cantabria (IFCA, UC-CSIC), Avenida de Los Castros s/n, 39005 Santander, Spain 4TRIUMF, 4004 Wesbrook Mall, Vancouver, British Columbia V6T 2A3, Canada (Received 20 September 2021; revised 3 March 2022; accepted 6 May 2022; published 9 June 2022) In many cosmologies dark matter clusters on subkiloparsec scales and forms compact subhalos, in which the majority of Galactic dark matter could reside. Null results in direct detection experiments since their advent four decades ago could then be the result of extremely rare encounters between the Earth and these subhalos. We investigate alternative and promising means to identify subhalo dark matter interacting with standard model particles: (1) subhalo collisions with old neutron stars can transfer kinetic energy and brighten the latter to luminosities within the reach of imminent infrared, optical, and ultraviolet telescopes; we identify new detection strategies involving single-star measurements and Galactic disk surveys, and obtain the first bounds on self-interacting dark matter in subhalos from the coldest known pulsar, PSR J2144–3933; (2) subhalo dark matter scattering with cosmic rays results in detectable effects; (3) historic Earth-subhalo encounters can leave dark matter tracks in Paleolithic minerals deep underground. These searches could discover dark matter subhalos weighing between gigaton and solar masses, with corresponding dark matter cross sections and masses spanning tens of orders of magnitude. DOI: 10.1103/PhysRevLett.128.231801 Introduction.—The enigma of dark matter persists [1].In the hope that its microscopic properties will be revealed through its scattering and annihilation into standard model (SM) states, extensive experimental efforts are under way to detect it [2–4]. Both direct and indirect dark matter (DM) searches often assume that DM is smoothly distributed over the Galactic halo. However, cold DM is expected to form subhalo structures down to Earth-mass scales (∼10−6M⊙) [5–7]. Such substructure should persist on length scales below a kiloparsec if small DM halos remain undisrupted during hierarchical clustering, as suggested by N-body simulations [8–10]. Furthermore, in many cosmological scenarios the small-scale power is enhanced—via, e.g., an early matter-dominated era or DM self-interactions augmenting primordial density perturbations—resulting in large fractions of DM surviving in subhalos [11–18], where certain collisionless DM cosmologies predict subhalo masses ranging from 10−28 −107M⊙[11,19,20]. It is known that subhalos could boost Galactic indirect detection fluxes, since the DM annihilation rate scales as the density squared [21–25]. Due attention has thereby been paid to annihilation signals, but it may be that DM only reveals its SM interactions via scattering processes, such as when it is asymmetric [26] or its annihilation is p-wave dominated [27–29]. Should the interacting component of DM reside in subhalos, direct detection (DD) experiments [30–32] may have observed no conclusive signal simply because the Earth has yet to encounter a subhalo since their inception. Thus other, novel probes of DM scattering are needed; in this Letter we outline three strategies. The first is to observe heating of old, nearby neutron stars (NSs) by subhalo DM. These NSs will travel through subhalos and capture constituent DM particles, in a process that transfers DM kinetic energy of order the DM mass, because of the steep gravitational potential around NSs. This method was first proposed to detect smooth halo DM [33], which would result in infrared emission from NSs observable by upcoming telescopes, and is currently an active area of research [34–57]. Here, we show that subhalo DM can impart larger NS luminosities than diffuse DM, bringing imminent optical and ultraviolet observations into play. Moreover, for DM with strong self-interactions, we show that due to enhanced accretion there is parameter space already constrained by observations of the coldest known NS, PSR J2144–3933 [58]. All in all, NS luminosity searches can probe gigaton to solar mass subhalos, for DM-nucleon cross sections as low as 10−45 cm2. Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW LETTERS 128, 231801 (2022) 0031-9007=22=128(23)=231801(9) 231801-1 Published by the American Physical Society Two other strategies to detect distant DM subhalos involve terrestrial probes that do not rely on Earth’s spatial position or the usual years-long runtime of terrestrial experiments. Specifically, encounters of subhalo DM with Galactic cosmic rays and with geological minerals deep underground can be discerned in several ways discussed below. Under the assumption of a smooth halo, cosmic rays and minerals have already been utilized to constrain strongly interacting DM [59–65]; here we will use these to set limits on DM subhalo masses and radii. Current constraints on subhalos come from gravitational microlensing surveys, applicable for asteroid to solar mass subhalos with sizes ranging up to stellar radii [66,67], and additional gravitational probes include pulsar timing arrays [68] and accretion glows in molecular clouds [69]. Direct detection of clustered DM up to the Planck mass was studied in Ref. [17], of higher mass composite DM in Refs. [70–77], and similar sensitivities with ancient minerals in Refs. [78–80]. Asteroid-mass subhalos colliding with stars were studied in Ref. [81], assuming geometric cross sections. We do not make this assumption and treat more general microscopic DM cross sections. Finally, DM substructure could impact DD experiments [82–84] and the biosphere [85]. Neutron star heating and cooling.—In this study we consider subhalos of uniform density of mass Mand maximum radius Rsh; other density profiles are expected to impact these results at the Oð1Þlevel. For simplicity, we assume all NSs are of the same mass MNS and radius RNS, with benchmark parameters MNS ¼1.5M⊙;R NS ¼10 km:ð1Þ The rate of NS-subhalo encounters for a single NS with MNS ≫Mis given by (see Supplemental Material [86]): ΓmeetðrÞ¼4ffiffiffi π pfχ ρχðrÞ MσvðrÞ ×Rsh þRNS ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þðvesc hvreliÞ2 r2 ;ð2Þ where ρχðrÞand σvðrÞare the DM density and speed dispersion at Galactic position r, and the escape velocity at the NS surface is vesc ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2GMNS=RNS p. Here, we have averaged the rate over the relative velocity vrel of the encounters, which we assume to follow a MaxwellBoltzmann distribution with dispersion σvrel ðrÞ¼ ffiffiffi 2 pσv ([100], Problem 8.8). For the solar position r⊙¼8.3kpc, we have ρ⊙¼0.4GeV=cm3,σv¼156 km=s[101] and the average subhalo-NS relative speed hvreli≃350 km=s, and for Rsh ≪RNS: Γmeetðr⊙Þ¼ fχ 1.9Gyr 10−15 M⊙ M:ð3Þ For concreteness, in our numerical results we fix the fraction of DM in subhalos as fχ¼1. Assuming that all of the incident DM flux in the subhalo is captured by the NS, the DM kinetic energy deposited in the NS after subhalo passage is Emeet ¼zM min1;Rco Rsh2;ð4Þ where the blueshift of DM falling into the neutron star increases the DM kinetic energy according to 1þz¼1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1−2GMNS=RNS p¼1.34;ð5Þ and the effective gravitational collection radius Rco of the neutron star is Rco ¼vesc hvrelil RNSð1þzÞ;ð6Þ where l¼1if the DM in the subhalo is effectively collisionless [102], and l¼2for DM that behaves like a fluid, in which case Rco will be the Bondi radius [103–106]. The latter expression assumes that the sound speed in the subhalo does not exceed the NS-subhalo relative velocity. In Eq. (6), the 1þzfactor accounts for the apparent NS radius to distant observers. DM in subhalos will heat an NS from internal temperature ˜ Tcold to ˜ Thot given by Z˜ Thot ˜ Tcold d ˜ TcVð˜ TÞ¼Emeet:ð7Þ Here and throughout this Letter, the tilde denotes temperatures measured by a distant observer: ˜ T≡T=ð1þzÞ, where Tis the NS local core temperature. The NS heat capacity cVis given by [107] (see the Appendix for more on NS cooling) cVð˜ TÞ¼4.8×1026 J=K˜ T 104K ¼2.7×10−21 M⊙=K˜ T 104K:ð8Þ Putting Eqs. (7) and (8) together, we obtain ˜ Thot 104K2 ¼˜ Tcold 104K2 þEmeet 6.2×10−18 M⊙ :ð9Þ Between subhalo encounters the NS cools by emitting photons from its surface and neutrinos from its interior; the latter mode is subdominant for ˜ T≲108K, which corresponds to a surface temperature (generally different from ˜ T PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-2 due to an insulating “envelope”on the NS) of ˜ Ts≲106K. We take the NS cooling time tcool from Ref. [108]; see the Appendix for further notes. Frequent NS encounters.—If subhalos heat NSs often enough, the NSs will maintain a steady-state temperature. Specifically, if the NS cooling time tcoolð˜ TsÞexceeds the timescale between subhalo encounters tmeet ¼Γ−1 meet, and if the encounters are energetic enough to impart temperatures > ˜ Ts, the NS will maintain a temperature of at least ˜ Ts.To find this temperature in terms of subhalo parameters, we set ˜ T2 hot ¼2 ˜ T2 cold in Eq. (9) to estimate the Emeet required to attain a steady-state temperature in the vicinity of ˜ Tcold, which we call ET meet. Then for a given subhalo radius and NS accretion rate, ET meet corresponds to some Min Eq. (4), which also implies a unique tmeet in Eq. (2). In Fig. 1we plot as a function of ˜ Ts(corresponding to the ET meet in the top xaxis) both tcool and tmeetðr⊙Þfor various Rsh. The condition tmeet <t cool then defines the NS temperatures for which a steady steady state scenario is possible. It can be seen that the subhalo encounter timescale does not exceed the cooling time for ˜ Ts≲3000 K in the limits where the subhalo is much smaller or much larger than the NS, whereas for Rsh comparable to RNS the condition could produce hotter NSs. This follows from Eqs. (2) and (4).For Rsh ≪RNS both tmeet and Emeet ∝M, and for Rsh ≫RNS both tmeet and Emeet ∝MR−2 sh , thus the tmeet vs ET meet curves are coincident in these limits. For intermediate Rsh the Sommerfeld enhancement factors for Γmeet and Emeet differ in such a way that greater energies can be deposited in less frequent encounters. The coldest NS observed, PSR J2144–3933, has a temperature upper limit of 2.9×104K for our benchmark NS [58], whereas our largest achievable steady-state temperature is ˜ Ts<6×103K (for Rsh ¼104km), for DM that accretes onto NSs in the collisionless regime, i.e., l¼1in Eq. (6). In the top panel of Fig. 2we show, in the space of Mvs Rsh, the subhalo parameters that can be explored in the future by measuring ˜ Tsdown to 1000 K using upcoming infrared telescopes like JWST, ELT, and TMT, within achievable times [33,35]. This region is bounded on the left by a black-dashed contour of ˜ Thot ¼ 1000 K [obtained from Eq. (9) for ˜ Tcold ≪ ˜ Thot], and on the right by a black-dashed contour of tmeet ¼tcoolð1000 KÞ¼ 1.2×107yr, and would be excluded if an NS at r¼r⊙is measured to be below 1000 K. (For a similar region corresponding to 104K, see Supplemental Material [86]). On the other hand, if DM has substantial self-interactions then PSR J2144–3933 already places a strong bound on subhalo DM. In this case we consider fluid accretion of the DM onto NSs, especially relevant for theories of dissipative and strongly self-interacting DM [109–113], for which DM can have a sizable sound speed (much like the interstellar medium), if collected into dense subhalos. The effective DM collection radius of the NS will now be enhanced by another factor of vesc=hvreli[i.e., l¼2in Eq. (6)] [103–106], and in a large parametric region Bondi accretion imparts enough energy for all NSs near Earth to have temperatures in excess of 2.9×104K. In Fig. 2, the orange region is constrained by the ˜ Ts≲3×104K bound on PSR J2144–3933’s temperature [58]. As discussed in Supplemental Material [86], NS accretion of DM may be enhanced beyond the estimate given here via disk formation if these self-interacting DM subhalos are tidally disrupted, although disruption itself will depend on the binding potential of the subhalo. In addition, subhalo DM with strong self-interactions could be constrained using, e.g., bullet cluster observations—this would require a detailed simulation of subhalo collisions in the bullet cluster, and the specific self-interacting DM model [114]. The regions just discussed may be probed so long as the dynamics of the DM NS scattering process lead to capture. This condition is depicted in black in the bottom panel of Fig. 2, for velocity-independent elastic scattering, modeling the NS as a free fermion gas; in the region above the black curve all the DM flux incident on the NS is captured. As explained in Refs. [33,102] and elsewhere, the mass capture rate is independent of mχ, hence the flat capture cross section for mχ≳0.5GeV. Below this mass, Pauli blocking of final state neutrons by the Fermi-degenerate nucleon medium reduces the capture cross section in inverse proportion to mχ. Below mχ≃35 MeV the sensitivities require further study as there is a new effect: in a FIG. 1. Regions showing “frequent”vs “infrequent”NSsubhalo encounters, requiring different observational strategies. Green: neutron star cooling time as a function of surface temperature ˜ Ts, assuming superfluid nucleons. Black: The time between encounters at the solar position for various subhalo radii, as a function of the DM kinetic energy deposited in NSs per encounter, obtained from Eqs. (2) and (4). The top xaxis shows energies required to maintain corresponding steady-state ˜ Ts values in the bottom xaxis. To the left of black-green intersections, observations of a single NS with indicated temperatures would probe corresponding subhalo masses and radii. To the right, astronomical surveys are required to find NS ensembles exceeding indicated temperatures. See text for further details. PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-3 superfluid NS, DM cannot scatter elastically on nucleons, since the DM kinetic energy is insufficient to excite nucleons above the OðMeVÞsuperfluid energy gap. In principle low-mass DM can still capture in NSs via excitation of collective modes such as phonons (see, e.g., Ref. [36]) and scattering with nonsuperfluid portions of the NS crust; we leave this to future work. Infrequent NS encounters.—For the scenario of tmeet >t cool we cannot expect a steady-state temperature. Thus the strategy must be to look for NSs not long after encounter—corresponding to the right side of Fig. 1— before they cool to temperatures below telescope sensitivities. This is done by surveying the sky for anomalously hot NSs. In Fig. 2top panel (again corresponding to the black capture region in the bottom panel) we show regions where at least 100 NSs with ˜ Ts>104K are expected in the “solar vicinity,”defined as a kiloparsec-radius circle around the Sun. Note that for the range of Rsh displayed, the subhalo crossing time is less than tmeet, where our treatment is valid. FIG. 2. Subhalo masses and radii probed by our detection strategies in the top panel, corresponding to limits on DM mass and nucleon scattering cross sections shown in the bottom panel. For collisionless DM, we show regions that can be probed by the observation of a single NS with 1000 K surface temperature (black dashed lines), as well as by future telescope surveys of NSs within a kiloparsec of Earth at the temperatures indicated (blue dot dashed lines). For dissipative-fluid DM, we exclude the orange region using observations of the coldest NS, for cross sections above the “100% capture”line in the bottom panel. Also shown are regions corresponding to existing constraints from DM-cosmic ray scattering (pink) and track searches in ancient mica (brown). Above the gray dot-dashed line, the galactic DM density exceeds the density in subhalos so DM does not significantly cluster. Limits from microlensing surveys are also shown. See text for more details. PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-4 Currently, the Hubble space telescope is able to achieve point-source sensitivity as cold as 3×104K[58]. Similar directed observations of nearby pulsars with telescopes like LUVOIR [115] may provide an alternative strategy to find anomalously hot NSs. In addition, the observation of isolated >103K NSs may be undertaken (see Supplemental Material [86]) by the Roman/WFIRST survey [116] and in conjunction with deep-field observations by JWST, ELT, and TMT [33]. We note that future optical survey sensitivities may detect 104K NSs within ∼20 parsecs, as discussed in Supplemental Material [86]. To obtain these regions we first estimate the probability of ˜ TNS being at least some ˜ Ts, when heated by subhalo passage to a surface temperature ˜ Thot s,as pð˜ Thot s> ˜ TNS > ˜ TsÞ¼1−e−ts=tmeet ;ð10Þ where ts¼tcoolð˜ TsÞ−tcoolð˜ Thot sÞis the time taken for the NS to cool from ˜ Thot sto ˜ Ts; see Supplemental Material [86] for the derivation. We then multiply this probability with the number of NSs N⊙ NS in the solar vicinity obtained from the NS surface density of Model 1A* in Ref. [117]. This model assumes a Maxwell distribution of NS velocities consistent with our encounter rate, and gives N⊙ NS ≃2×105. As in Ref. [117] we assume that NSs are formed at a constant rate, so that the NS age distribution is peaked at the age of the Galaxy ≃5×109yr (We have verified that including a broader distribution of NS ages (as given in Ref. [116]) does not noticeably change our results.). This implies that the vast majority of the NSs we consider would have ˜ Ts≪1000 K in the absence of subhalo heating (cf. the NS cooling curve in Fig. 1). Other probes.—(i) Cosmic ray (CR) scattering.DM particles with stronger-than-electroweak cross sections would undergo appreciable scatters with Galactic CRs, which has been used to bound DM’s scattering cross section by studying the subsequent softening of the CR spectrum [118], and by constraining the flux of γrays [59] and up-scattered DM [60–63] produced in DM CR collisions. These limits also apply to subhalo DM, as the DM optical depth of CRs τCR over long distances is unaffected by DM clustering. In the top panel of Fig. 2, we show the M-Rsh range constrained by CR scattering, for the σNχ-mχrange shown in the bottom panel. The top panel region is bounded from below by requiring that τCR >1(setting the effective diffusion length ¼8kpc [60,119]). The upper bound is simply the requirement that DM is appreciably clustered in subhalos. (ii) Paleodetectors. Since they contain gigayear timescale archeological records of the particle interactions, ancient mineral slabs are potentially sensitive to the historic passage of DM subhalos through the Earth. Decades ago, Refs. [64,65] derived bounds on (unclumped) DM and monopole interactions by analyzing the microscopic structure of ancient mica. Reference [120] pointed out that ancient mica could be used to search for subhalo DM. Here, we recast bounds from Refs. [64,65], using the fact that for some subhalo model space the total number of DM-mica scatters is unchanged by DM clumping, discussed more in Supplemental Material [86]. We thus show existing constraints in both panels of Fig. 2similar to CR-scattering limits, with the Rsh-Mregion now bounded on the right by the requirement that the timescale of Earth encounters with DM subhalos is smaller than the age of the mica samples [determined from fission track dating [65] to be ð51Þ×108yr]. Recently, a number of proposals have suggested further mineralogic DM searches [79,80,121–124], including for subhalo DM [125]. For mχ∼0.1–1GeV (mχ∼5–105GeV), if ≳10−3 (≳10−10) of DM is diffuse, stronger DD bounds [126–128] than for CRs (mica) apply. In Fig. 2we have also shown constraints from gravitational microlensing surveys [66,67], with the reach in radii limited from below by each survey’s maximum Einstein radius, and reach in mass limited by observational cadences and finite source effects. While these constraints are agnostic to the scattering properties of DM, our probes are seen to highly complement them, reaching much larger and lighter subhalos. We have also indicated the Schwarzschild radii corresponding to large M; the region below this line is unphysical. If a small fraction of DM remains unclumped, a possibility we have not explored in this study, direct detection searches may set bounds on the σnχ-mχplane; on the other hand, our bounds on the plane apply for any distribution of DM into smooth and clumped components. Discussion.—This Letter charts the way forward for uncovering scattering interactions of subhalo DM. We have shown that cosmic ray scattering and mineral detection already provide some sensitivity, while novel methods using NS temperatures already constrain subhalo DM with strong self-interactions. In the future, subhalos over a vast range of masses and sizes may be uncovered by observations of either individual NSs or NSs in survey ensembles. Interestingly, imminent limits on the DM scattering cross section for subhalo-bound DM could become stronger than those from direct detection searches on unclumped DM (see Fig. 2bottom panel and, e.g., Ref. [128]). Several new avenues await exploration. While DM could heat old NSs to detectable brightness by imparting kinetic energy, greater NS temperatures and other interesting signals could arise if DM annihilations in the NS core are present. For instance, if DM thermalizes slowly with the NS, annihilations may turn on late, so that there are two reheating periods in our “infrequent encounter”regime. DM annihilations following capture of subhalo DM could also produce detectable signals in a variety of celestial bodies [84,129–131]. While we have focused on PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-5 DM-nucleon scattering, some of our approaches also apply to DM that scatters on leptons. The lepton content of NSs, though relatively small, provides substantial sensitivity to DM-lepton scattering [50–54]. Upscattering of DM with CR electrons should produce a detectable high-speed flux of DM, as with CR protons [61]. Exploring DM-electron scattering would therefore extend our results to even lighter DM masses than depicted in Fig. 2. While for simplicity we had assumed uniform subhalo masses, more realistic mass distributions could be investigated. Finally, our NS heating mechanism via DM subhalo encounters provides, in addition to probing the neutron portal to beyond-the-SM scenarios [132,133], further motivation to imminent astronomical missions to make exciting discoveries in fundamental physics. We thank Raghuveer Garani for discussions on neutron superfluidity and Kevin Zhou for correspondence about NS-subhalo collisions. B. J. K. thanks the Spanish Agencia Estatal de Investigación (AEI, MICIU) for the support to the Unidad de Excelencia María de Maeztu Instituto de Física de Cantabria, Ref. MDM-2017-0765. The work of J. B. and N. R. is supported by the Natural Sciences and Engineering Research Council of Canada. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities. TRIUMF receives federal funding via a contribution agreement with the National Research Council Canada. Appendix: Neutron star passive cooling.—The time for the NS to cool to ˜ Tcool, in a process that is insensitive to the starting temperature, is given by [108] tcoolð˜ T9Þ=yr ¼tenv ¼s−k 1q−γ½ð1þðs1=qÞk˜ T2−n 9Þ−γ=k −1; ˜ Tcool > ˜ Tenv; tenv þð3s2Þ−1ð˜ T−2 9− ˜ T−2 envÞ; ˜ Tcool ≤ ˜ Tenv;ðA1Þ where ˜ T9¼˜ Tcool=ð109KÞ,q¼2.75 ×10−2,s1¼8.88× 10−6,s2¼8.35 ×104,k¼ðn−2Þ=ðn−αÞ, and γ¼ð2− αÞ=ðn−αÞwith α¼2.2and n¼8. The temperature ˜ Tenv ≃4000 K corresponds to the time after which the surface and internal NS temperatures equalize due to the thinning of the crustal outer envelope. We now provide some physics background to Eq. (A1). First recall that temperatures denoted with a tilde are in the frame of a distant observer. During quiescent periods the NS temperature evolution is given by cvð˜ TÞd ˜ T dt ¼−L∞ νð˜ TÞ−L∞ γð˜ TÞ;ðA2Þ where the neutrino luminosity of the NS as measured by a distant observer of our benchmark NS is given by [134] L∞ νð˜ TÞ¼1.33 ×1039 J=yr ˜ T 109K8 ;ðA3Þ applicable for slow and modified Urca processes, which we take to be the only neutrino cooling mechanism as prescribed by the “minimal cooling”paradigm [135]. These processes dominate the NS cooling down to ˜ T¼108K before photon cooling takes over. The luminosity of photon blackbody emission from the NS surface is L∞ γð˜ TsÞ¼4πð1þzÞ2R2 NS ˜ T4 s:ðA4Þ To solve Eq. (A2) we need a relation between Tsand T, which depends on the composition of the outermost envelope, which acts as an insulating layer for temperatures above Oð103ÞK. For smaller temperatures the envelope becomes too thin for insulation, and the surface temperature reflects the internal temperature [136,137].We assume a standard iron envelope for which we have the following relation at high temperatures [138,139]: Ts¼106K MNS 1.5M⊙10 km RNS 21=4T 9.43 ×107K0.55 : ðA5Þ We identify the thin-envelope regime by solving for Ts¼T in the above equation, giving Tenv ¼3908 K. Below this temperature we model the Ts−Trelation by simply setting Ts¼T. We now have all the ingredients to solve Eq. (A2).Wedo this by following the prescription in Ref. [108], suitably improving it to account for the thin-envelope Ts−T relation, and obtain the analytic solution in Eq. (A1). The heat capacity of an NS depends sensitively on the superfluid and superconducting properties of NS nucleons. It has long been recognized that nucleons could indeed form Cooper pairs in NS cores with an MeV energy gap, corresponding to a critical temperature Tc≃1010 K [140–142]. Observational fits to cooling curves support this hypothesis [137]. The heat capacity of degenerate matter is set by particle-hole excitations close to the Fermi surface, thus the energy gap exponentially suppresses the nucleon contribution to the heat capacity for NS temperatures ≪Tc, the temperature zone we are interested in. Hence, our NS heat capacity is set by the contribution from degenerate PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-6 electrons (whose pairing-up critical temperature is <100 K), and is given by Eq. (8), obtained from analytic fits for the electronic cVprovided in Ref. [107]. In general we expect Oð10%Þvariations to the above arising from the choice of the equation of state for NS core matter and of the NS mass-radius configuration; in this work we have simply used Eq. (8). Finally, we note that while in standard NS cooling scenarios, NSs are expected to reach ∼100–103K over 1×109yr timescales [107,136], in some NS thermal evolution scenarios a temperature as high as 104K can persistent, if there is some additional source of NS heating like an anomalously large initial magnetic field strength [143] or an exceptionally fast initial rate of spin paired with rotochemical heating for certain NS equation of state models [144]. However, in both these cases there is still the expectation that many late-stage NS will have temperatures below ∼103K. *[email protected] †kav[email protected].es ‡[email protected] [1] G. Bertone, D. Hooper, and J. Silk, Phys. Rep. 405, 279 (2005). [2] L. Baudis, Phys. Dark Universe 4, 50 (2014). [3] J. M. Gaskins, Contemp. Phys. 57, 496 (2016). [4] M. Schumann, J. Phys. G 46, 103003 (2019). [5] K. P. Zybin, M. I. Vysotsky, and A. V. Gurevich, Phys. Lett. A 260, 262 (1999). [6] S. Hofmann, D. J. Schwarz, and H. Stöcker, Phys. Rev. D 64, 083507 (2001). [7] V. Berezinsky, V. Dokuchaev, and Y. Eroshenko, Phys. Rev. D 77, 083519 (2008). [8] L. Bergstrom, J. Edsjo, P. Gondolo, and P. Ullio, Phys. Rev. D 59, 043506 (1999). [9] F. C. van den Bosch, G. Ogiya, O. Hahn, and A. Burkert, Mon. Not. R. Astron. Soc. 474, 3043 (2018). [10] F. C. van den Bosch and G. Ogiya, Mon. Not. R. Astron. Soc. 475, 4066 (2018). [11] A. L. Erickcek and K. Sigurdson, Phys. Rev. D 84, 083503 (2011). [12] G. Barenboim and J. Rasero, J. High Energy Phys. 04 (2014) 138. [13] J. J. Fan, O. Özsoy, and S. Watson, Phys. Rev. D 90, 043536 (2014). [14] J. A. Dror, E. Kuflik, B. Melcher, and S. Watson, Phys. Rev. D 97, 063524 (2018). [15] P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016). [16] M. R. Buckley and A. DiFranzo, Phys. Rev. Lett. 120, 051102 (2018). [17] S. Nussinov and Y. Zhang, J. High Energy Phys. 03 (2020) 133. [18] G. Barenboim, N. Blinov, and A. Stebbins, J. Cosmol. Astropart. Phys. 12 (2021) 026. [19] A. L. Erickcek, P. Ralegankar, and J. Shelton, Phys. Rev. D 103, 103508 (2021). [20] A. L. Erickcek, P. Ralegankar, and J. Shelton, J. Cosmol. Astropart. Phys. 01 (2022) 017. [21] P. Gondolo and G. Gelmini, Nucl. Phys. B360, 145 (1991). [22] N. Hiroshima, S. Ando, and T. Ishiyama, Phys. Rev. D 97, 123002 (2018). [23] C. Blanco, M. S. Delos, A. L. Erickcek, and D. Hooper, Phys. Rev. D 100, 103010 (2019). [24] M. S. Delos, T. Linden, and A. L. Erickcek, Phys. Rev. D 100, 123546 (2019). [25] M. S. Delos, Phys. Rev. D 100, 063505 (2019). [26] K. Petraki and R. R. Volkas, Int. J. Mod. Phys. A 28, 1330028 (2013). [27] Y. Zhao, X.-J. Bi, H.-Y. Jia, P.-F. Yin, and F.-R. Zhu, Phys. Rev. D 93, 083513 (2016). [28] Y. Zhao, X.-J. Bi, P.-F. Yin, and X. Zhang, Phys. Rev. D 97, 063013 (2018). [29] K. K. Boddy, J. Kumar, A. B. Pace, J. Runburg, and L. E. Strigari, Phys. Rev. D 102, 023029 (2020). [30] A. Drukier and L. Stodolsky, Phys. Rev. D 30, 2295 (1984). [31] M. W. Goodman and E. Witten, Phys. Rev. D 31, 3059 (1985). [32] A. K. Drukier, K. Freese, and D. N. Spergel, Phys. Rev. D 33, 3495 (1986). [33] M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, Phys. Rev. Lett. 119, 131801 (2017). [34] E. Berti et al., in 2022 Snowmass Summer Study, arXiv:2203.07984. [35] N. Raj, P. Tanedo, and H.-B. Yu, Phys. Rev. D 97, 043006 (2018). [36] J. F. Acevedo, J. Bramante, R. K. Leane, and N. Raj, J. Cosmol. Astropart. Phys. 03 (2020) 038. [37] C.-S. Chen and Y.-H. Lin, J. High Energy Phys. 08 (2018) 069. [38] N. F. Bell, G. Busoni, and S. Robles, J. Cosmol. Astropart. Phys. 09 (2018) 018. [39] R. Garani, Y. Genolini, and T. Hambye, J. Cosmol. Astropart. Phys. 05 (2019) 035. [40] D. A. Camargo, F. S. Queiroz, and R. Sturani, J. Cosmol. Astropart. Phys. 09 (2019) 051. [41] K. Hamaguchi, N. Nagata, and K. Yanagi, Phys. Lett. B 795, 484 (2019). [42] W.-Y. Keung, D. Marfatia, and P.-Y. Tseng, J. High Energy Phys. 07 (2020) 181. [43] N. F. Bell, G. Busoni, S. Robles, and M. Virgato, J. Cosmol. Astropart. Phys. 09 (2020) 028. [44] B. Dasgupta, A. Gupta, and A. Ray, J. Cosmol. Astropart. Phys. 10 (2020) 023. [45] R. Garani, A. Gupta, and N. Raj, Phys. Rev. D 103, 043019 (2021). [46] T. N. Maity and F. S. Queiroz, Phys. Rev. D 104, 083019 (2021). [47] N. F. Bell, G. Busoni, T. F. Motta, S. Robles, A. W. Thomas, and M. Virgato, Phys. Rev. Lett. 127, 111803 (2021). [48] Y.-P. Zeng, X. Xiao, and W. Wang, Phys. Lett. B 824, 136822 (2022). [49] F. Anzuini, N. F. Bell, G. Busoni, T. F. Motta, S. Robles, A. W. Thomas, and M. Virgato, J. Cosmol. Astropart. Phys. 11 (2021) 056. PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-7 [50] N. F. Bell, G. Busoni, and S. Robles, J. Cosmol. Astropart. Phys. 06 (2019) 054. [51] R. Garani and J. Heeck, Phys. Rev. D 100, 035039 (2019). [52] A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Phys. Lett. B 809, 135767 (2020). [53] A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Phys. Rev. D 102, 123002 (2020). [54] N. F. Bell, G. Busoni, S. Robles, and M. Virgato, J. Cosmol. Astropart. Phys. 03 (2021) 086. [55] J. Bramante, B. J. Kavanagh, and N. Raj, arXiv: 2109.04582. [56] M. Fujiwara, K. Hamaguchi, N. Nagata, and J. Zheng, arXiv:2204.02238. [57] K. Hamaguchi, N. Nagata, and M. E. Ramirez-Quezada, arXiv:2204.02413. [58] S. Guillot, G. G. Pavlov, C. Reyes, A. Reisenegger, L. Rodriguez, B. Rangelov, and O. Kargaltsev, Astrophys. J. 874, 175 (2019). [59] R. H. Cyburt, B. D. Fields, V. Pavlidou, and B. D. Wandelt, Phys. Rev. D 65, 123503 (2002). [60] T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, 171801 (2019). [61] Y. Ema, F. Sala, and R. Sato, Phys. Rev. Lett. 122, 181802 (2019). [62] C. V. Cappiello and J. F. Beacom, Phys. Rev. D 100, 103011 (2019). [63] M. Andriamirado et al. (PROSPECT Collaboration), Phys. Rev. D 104, 012009 (2021). [64] P. B. Price and M. H. Salamon, Phys. Rev. Lett. 56, 1226 (1986). [65] D. P. Snowden-Ifft, E. S. Freeman, and P. B. Price, Phys. Rev. Lett. 74, 4133 (1995). [66] D. Croon, D. McKeen, and N. Raj, Phys. Rev. D 101, 083013 (2020). [67] D. Croon, D. McKeen, N. Raj, and Z. Wang, Phys. Rev. D 102, 083021 (2020). [68] J. A. Dror, H. Ramani, T. Trickle, and K. M. Zurek, Phys. Rev. D 100, 023003 (2019). [69] Y. Bai, A. J. Long, and S. Lu, J. Cosmol. Astropart. Phys. 09 (2020) 044. [70] J. Bramante, B. Broerman, R. F. Lang, and N. Raj, Phys. Rev. D 98, 083516 (2018). [71] J. Bramante, B. Broerman, J. Kumar, R. F. Lang, M. Pospelov, and N. Raj, Phys. Rev. D 99, 083010 (2019). [72] J. Bramante, J. Kumar, and N. Raj, Phys. Rev. D 100, 123016 (2019). [73] J. F. Acevedo, J. Bramante, A. Goodman, J. Kopp, and T. Opferkuch, J. Cosmol. Astropart. Phys. 04 (2021) 026. [74] J. F. Acevedo, J. Bramante, and A. Goodman, Phys. Rev. D 103, 123022 (2021). [75] P. Adhikari et al.,Phys. Rev. Lett. 128, 011801 (2022). [76] J. F. Acevedo, J. Bramante, and A. Goodman, Phys. Rev. D 105, 023012 (2022). [77] D. Carney et al.,arXiv:2203.06508. [78] D. M. Jacobs, G. D. Starkman, and B. W. Lynn, Mon. Not. R. Astron. Soc. 450, 3418 (2015). [79] R. Ebadi, A. Mathur, E. H. Tanin, N. D. Tailby, M. C. Marshall, A. Ravi, R. Trubko, R. R. Fu, D. F. Phillips, S. Rajendran, and R. L. Walsworth, Phys. Rev. D 104, 015041 (2021). [80] J. F. Acevedo, J. Bramante, and A. Goodman, arXiv: 2105.06473. [81] A. Das, S. A. R. Ellis, P. C. Schuster, and K. Zhou, Phys. Rev. Lett. 128, 021101 (2022). [82] D. Stiff, L. M. Widrow, and J. Frieman, Phys. Rev. D 64, 083516 (2001). [83] M. Kamionkowski and S. M. Koushiappas, Phys. Rev. D 77, 103509 (2008). [84] A. Ibarra, B. J. Kavanagh, and A. Rappelt, J. Cosmol. Astropart. Phys. 12 (2019) 013. [85] J. I. Collar, Phys. Lett. B 368, 266 (1996). [86] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.128.231801 for technicalities on neutron star observations and their encounters with subhalos, and cosmic ray and paleolithic limits, which includes Refs. [87–99]. [87] B. J. Kavanagh, T. D. P. Edwards, L. Visinelli, and C. Weniger, Phys. Rev. D 104, 063038 (2021). [88] T. D. P. Edwards, B. J. Kavanagh, L. Visinelli, and C. Weniger, Phys. Rev. Lett. 127, 131103 (2021). [89] M. J. Rees, Nature (London) 333, 523 (1988). [90] J. L. Dai, G. Lodato, and R. M. Cheng, Space Sci. Rev. 217, 12 (2021). [91] T. Piran, G. Svirski, J. Krolik, R. M. Cheng, and H. Shiokawa, Astrophys. J. 806, 164 (2015). [92] A. Ibarra, B. J. Kavanagh, and A. Rappelt, J. Cosmol. Astropart. Phys. 12 (2019) 013. [93] H. T. Diehl et al.,10.2172/1596042 (2019). [94] J. R. Peterson, 10.2172/1272167 (2016). [95] v. Ivezićet al. (LSST Collaboration), Astrophys. J. 873, 111 (2019). [96] R. E. Rutledge, D. B. Fox, M. Bogosavljevic, and A. Mahabal, Astrophys. J. 598, 458 (2003). [97] M. L. Turner, R. E. Rutledge, R. Letcavage, A. S. H. Shevchuk, and D. B. Fox, Astrophys. J. 714, 1424 (2010). [98] W. Voges et al.,Astron. Astrophys. 349, 389 (1999). [99] P. Predehl et al.,Astron. Astrophys. 647, A1 (2021). [100] J. Binney and S. Tremaine, Galactic Dynamics, 2nd ed. (Princeton University Press, Princeton, NJ, 2008), ISBN 978-0-691-13026-2. [101] A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness, Astrophys. J. 871, 120 (2019). [102] I. Goldman and S. Nussinov, Phys. Rev. D 40, 3221 (1989). [103] H. Bondi and F. Hoyle, Mon. Not. R. Astron. Soc. 104, 273 (1944). [104] H. Bondi, Mon. Not. R. Astron. Soc. 112, 195 (1952). [105] M. C. Begelman, Mon. Not. R. Astron. Soc. 181, 347 (1977). [106] E. Shima, T. Matsuda, H. Takeda, and K. Sawada, Mon. Not. R. Astron. Soc. 217, 367 (1985). [107] D. D. Ofengeim, M. Fortin, P. Haensel, D. G. Yakovlev, and J. L. Zdunik, Phys. Rev. D 96, 043002 (2017). [108] D. D. Ofengeim and D. G. Yakovlev, Mon. Not. R. Astron. Soc. 467, 3598 (2017). [109] J. H. Chang, D. Egana-Ugrinovic, R. Essig, and C. Kouvaris, J. Cosmol. Astropart. Phys. 03 (2019) 036. PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-8 [110] R. Huo, H.-B. Yu, and Y.-M. Zhong, J. Cosmol. Astropart. Phys. 06 (2020) 051. [111] J. Liu, X.-P. Wang, and W. Xue, Phys. Rev. D 100, 123012 (2019). [112] X. Shen, P. F. Hopkins, L. Necib, F. Jiang, M. BoylanKolchin, and A. Wetzel, Mon. Not. R. Astron. Soc. 506, 4421 (2021). [113] M. Ryan, J. Gurian, S. Shandera, and D. Jeong, arXiv: 2106.13245. [114] S. Tulin and H.-B. Yu, Phys. Rep. 730, 1 (2018). [115] The LUVOIR Team, arXiv:1912.06219. [116] WFIRST, arXiv:1208.4012. [117] N. Sartore, E. Ripamonti, A. Treves, and R. Turolla, Astron. Astrophys. 510, A23 (2010). [118] C. V. Cappiello, K. C. Y. Ng, and J. F. Beacom, Phys. Rev. D99, 063004 (2019). [119] T. Bringmann, F. Donato, and R. A. Lineros, J. Cosmol. Astropart. Phys. 01 (2012) 049. [120] E. A. Baltz, A. J. Westphal, and D. P. Snowden-Ifft, Phys. Rev. D 59, 023510 (1998). [121] A. K. Drukier, S. Baum, K. Freese, M. Górski, and P. Stengel, Phys. Rev. D 99, 043014 (2019). [122] T. D. P. Edwards, B. J. Kavanagh, C. Weniger, S. Baum, A. K. Drukier, K. Freese, M. Górski, and P. Stengel, Phys. Rev. D 99, 043541 (2019). [123] S. Baum, A. K. Drukier, K. Freese, M. Górski, and P. Stengel, Phys. Lett. B 803, 135325 (2020). [124] S. Baum, T. D. P. Edwards, B. J. Kavanagh, P. Stengel, A. K. Drukier, K. Freese, M. Górski, and C. Weniger, Phys. Rev. D 101, 103017 (2020). [125] S. Baum, W. DeRocco, T. D. P. Edwards, and S. Kalia, Phys. Rev. D 104, 123015 (2021). [126] A. H. Abdelhameed et al. (CRESST Collaboration), Phys. Rev. D 100, 102002 (2019). [127] E. Aprile et al. (XENON Collaboration), Phys. Rev. Lett. 123, 241803 (2019). [128] Y. Meng et al. (PandaX-4T Collaboration), Phys. Rev. Lett. 127, 261802 (2021). [129] R. K. Leane and J. Smirnov, Phys. Rev. Lett. 126, 161101 (2021). [130] R. K. Leane, T. Linden, P. Mukhopadhyay, and N. Toro, Phys. Rev. D 103, 075030 (2021). [131] R. K. Leane and T. Linden, arXiv:2104.02068. [132] D. McKeen, M. Pospelov, and N. Raj, Phys. Rev. D 103, 115002 (2021). [133] D. McKeen, M. Pospelov, and N. Raj, Phys. Rev. Lett. 127, 061805 (2021). [134] A. Cumming, E. F. Brown, F. J. Fattoyev, C. J. Horowitz, D. Page, and S. Reddy, Phys. Rev. C 95, 025806 (2017). [135] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Astrophys. J. Suppl. Ser. 155, 623 (2004). [136] D. G. Yakovlev and C. J. Pethick, Annu. Rev. Astron. Astrophys. 42, 169 (2004). [137] A. Y. Potekhin, D. A. Zyuzin, D. G. Yakovlev, M. V. Beznogov, and Y. A. Shibanov, Mon. Not. R. Astron. Soc. 496, 5052 (2020). [138] E. H. Gudmundsson, C. J. Pethick, and R. I. Epstein, Astrophys. J. 272, 286 (1983). [139] M. V. Beznogov, A. Y. Potekhin, and D. G. Yakovlev, Phys. Rep. 919, 1 (2021). [140] D. G. Yakovlev, K. P. Levenfish, and Y. A. Shibanov, Phys. Usp. 42, 737 (1999). [141] D. Page and S. Reddy, Annu. Rev. Nucl. Part. Sci. 56, 327 (2006). [142] B. Haskell and A. Sedrakian, Astrophysics and Space Science Library 457, 401 (2018). [143] D. Page, U. Geppert, and T. Zannias, Astron. Astrophys. 360, 1052 (2000). [144] K. Hamaguchi, N. Nagata, and K. Yanagi, Phys. Lett. B 795, 484 (2019). PHYSICAL REVIEW LETTERS 128, 231801 (2022) 231801-9