Thermalization and annihilation of dark matter in neutron stars

The capture of dark matter, and its subsequent annihilation, can heat old, isolated neutron stars. In order for kinetic heating to be achieved, the captured dark matter must undergo sufficient scattering to deposit its kinetic energy in the star. We find that this energy deposit typically occurs quickly, for most of the relevant parameter space. In order for appreciable annihilation heating to also be achieved, the dark matter must reach a state of capture-annihilation equilibrium in the star. We show that this can be fulfilled for all types of dark matter-baryon interactions. This includes cases where the scattering or annihilation cross sections are momentum or velocity suppressed in the non-relativistic limit. Importantly, we find that capture-annihilation equilibrium, and hence maximal annihilation heating, can be achieved without complete thermalization of the captured dark matter. For scattering cross sections that saturate the capture rate, we find that capture-annihilation equilibrium is typically reached on a timescale of less than 1 year for vector interactions and 104 years for scalar interactions.

F Capture-annihilation equilibrium for the EFT operators 30 1 Introduction There has been much recent work on the capture of dark matter (DM) in neutron stars (NSs) , as a sensitive probe of dark matter interactions with ordinary matter.This can potentially be used to test dark matter interactions in a way that is highly complementary to experiments on Earth, especially since DM is accelerated to relativistic speeds during infall to a NS.In some cases, neutron star techniques could potentially probe interactions that would be difficult or impossible to ever observe in dark matter direct detection experiments.This includes dark matter that is too light to leave a detectable signal in nuclear-recoil experiments, or interactions for which the non-relativistic scattering cross section is momentum suppressed.It was recently pointed out that old, isolated, NSs in the Solar neighborhood could be heated by DM capture [13], leading to a temperature increase of ∼ 2000 K.At ages greater than ∼ 10 Myr, isolated NSs are expected to cool to temperatures below this, provided they are not reheated by accretion of standard matter or by internal heating mechanisms [40].As a result, the observation of a local NS with a temperature O(1000 K) could provide stringent constraints on DM interactions.Importantly, NS temperatures in this range would result in near-infrared emission, potentially detectable by future telescopes.
There are two distinct contributions to this heating: 1. Kinetic heating, where the DM kinetic energy is deposited in the NS medium.
2. Annihilation heating, where the DM rest mass energy is deposited through the annihilation of DM to particles that are trapped in the star.
The kinetic heating occurs as follows: The initial scattering interaction, which leads to gravitational capture of the DM particle, transfers only a small portion of the DM kinetic energy to the star.The rest of the kinetic energy is transferred through subsequent scattering interactions of the gravitationally bound DM until, eventually, the DM thermalizes with constituents in the centre of the star.In general, a large number of collisions are typically required for full thermalization.Moreover, if the scattering cross section is momentum suppressed in the non-relativistic limit, the time between collisions, and hence the time required to achieve full thermalization, can become very large, larger than the age of the Universe.Importantly, we shall see that even in cases where the full thermalization process is slow, the majority of the kinetic energy is deposited very quickly.
DM annihilation occurs in a region very close to the centre of the star, where the thermalized DM accumulates.The annihilation rate will increase as the DM abundance in the star grows until, eventually, a state of equilibrium between capture and annihilation is reached.When this occurs, the annihilation heating is maximized.The complete thermalization of the DM will result in a smaller, denser sphere of thermalized DM in the centre of the star, assisting in annihilation.However, we shall find that capture-annihilation equilibrium, and hence maximal annihilation heating, can be achieved without full thermalization.
The thermalization process was previously examined in refs.[41,42] for a subset of interaction types.These previous studies did not consider the importance of annihilation heating from partially thermalized DM.In this paper, we present a detailed calculation of the timescales required for thermalization and for capture-annihilation equilibrium, utilising our improved treatment of DM capture in NSs [22,24,25,27].We perform these calculations for a full set of DM-nucleon interaction types, parameterized by a set of Effective Field Theory (EFT) operators for fermionic DM.This includes operators for which either the scattering cross section or the annihilation cross section is suppressed in the non-relativistic limit.By properly accounting for the annihilation of partially thermalized DM, we show that full kinetic plus annihilation heating can be achieved for most of the interesting parameter space, on a short timescale.
This paper is organized as follows: We briefly review NS structure and composition in Section 2, and outline the calculation of the NS temperature due to DM energy deposition in Section 3.
Section 4 reviews the DM capture and interaction rates, while Section 5 examines the thermalization process.The annihilation of captured DM is discussed in Section 6, where we explain how to modify the usual capture-annihilation criterion to account for the annihilation of partially thermalized DM.The timescales for kinetic and annihilation heating are discussed in Section 7 and concluding remarks given in Section 8.

Neutron Star Composition
Massive stars, with masses above ∼ 8 − 10M ⊙ and below ∼ 20M ⊙ , end their life in core-collapse supernova explosions, leaving behind the densest known compact stellar object, a neutron star.The strongly degenerate matter found in a NS provides a unique means to test fundamental interactions and to look for new physics.Despite recent theoretical and observational breakthroughs, many uncertainties remain regarding NS composition and their equation of state (EoS).
Several phase transitions are expected to be found towards the interior of a NS.Below a thin atmosphere, the stellar interior can be broadly divided into a crust and a core.The outer layers of the former insulate the surface from the hot core.The core accounts for ∼ 99% of the stellar mass and is composed mostly of strongly degenerate neutrons.However, as the density increases towards the centre of these stellar remnants, beta equilibrium allows for the existence of protons, electrons and muons; the so-called npeµ matter.In addition, the extremely dense inner core of massive NSs may harbor exotic species, such as hyperons, which also appear under beta equilibrium.
For a consistent estimation of the DM capture rate, thermalization times and annihilation rate, an equation of state (EoS) that relates pressure and density is needed in order to solve the NS structure equations.Of the several EoSs of dense matter available in the literature [43][44][45][46][47][48][49], we have chosen the quark-meson coupling (QMC) model [50,51] as in refs.[25,27].This relativistic EoS is consistent with current constraints and permits the existence of ∼ 2M ⊙ NSs that contain hyperonic matter [44].We also assume a non-rotating, non-magnetized, spherically symmetric NS.We then solve the coupled system of Tolman-Oppenheimer-Volkoff (TOV) equations [52,53] and QMC EoS for given central baryon number densities n c B (see Table 1), from the NS centre to the surface.In this way, the radius, R ⋆ , and the gravitational mass, M ⋆ , of the star are determined, as well as radial profiles for the number density and effective mass m eff i of each particle species, and general relativity corrections embedded in B(r) (the coefficient that appears in the time part of the Schwarzschild metric).For an example of these radial profiles see Fig. 1 of refs.[22,27] and Fig. 2 of ref. [24].
In Table 1, we show the relevant NS quantities needed for our calculations, for three of the four NS configurations obtained in ref. [27].These three benchmark stars will be used to illustrate the dependence of our results on the compactness of the NS.

Neutron Star Temperature from Dark Matter Heating
We now discuss the potential NS heating that can be achieved by DM scattering and annihilation.We assume a nearby NS, located in the Solar neighborhood, and thus take ρ χ = 0.4 GeV cm −3 , v ⋆ = 230 km s −1 and v d = 270 km s −1 as the DM density, NS velocity, and DM dispersion velocity respectively.DM deposits energy into the NS via two mechanisms: (i) kinetic heating due to scattering with the constituents of the NS and (ii) annihilation of DM to SM particles that do not escape the star.We define the DM contribution to the temperature, measured by an observer far From left to right, these quantities are NS mass; radius; neutron effective mass m eff n (0) and Fermi energy ε F,n (0) at the NS centre; and B(0) and B(R ⋆ ), which are related to the escape velocity.from the star, as T ∞ χ = B(R ⋆ ) T χ , where, B(r) is the time component of the Schwarzchild metric.Assuming black body radiation, this temperature will be given by

EoS
where Ėχ = Ėχ,kin + Ėχ,ann is the rate of energy deposition and we have assumed the absence of any other source of heating.
The DM kinetic energy is deposited at the rate Ėχ,kin ≃ m χ 1 where C geom is the maximum DM capture rate (defined later in Eq. 4.10) and f quantifies how efficiently DM is captured, where we sum over the capture rates C i for scattering on all baryonic species i in the star.Note that we have used B(0) in Eq. 3.2, instead of B(R ⋆ ), which was previously used in the literature.This is because gravitational potential energy is converted to kinetic energy as the DM falls deeper into the NS.Therefore, the total energy the DM can deposit is equal to the kinetic energy it gains when moving from infinity to the centre of the star.If this were the only source of heating, the observed temperature would be T ∞ χ,kin ∼ 1870 K f 1/4 for the QMC-2 (1.5M ⊙ ) benchmark NS.For the 1M ⊙ and 1.9M ⊙ NSs, we find ∼ 1510 K f 1/4 and ∼ 2240 K f 1/4 , respectively.
Annihilation of DM in the centre of the NS causes further heating.The annihilation rate Γ ann , and hence the annihilation heating, is maximized when capture-annihilation equilibrium has been achieved.In this limit, the DM annihilation rate is given by Γ ann = C/2.Then, the rate at which DM deposits all of its energy, both kinetic and rest-mass, can be expressed as This rate implies a temperature of T ∞ χ,kin+ann ∼ 2410 K f 1/4 for the 1.5M ⊙ NS.For the lightest NS considered (1M ⊙ ) this value decreases to ∼ 2160 K, while for the heaviest NS (1.9M ⊙ ), this temperature reaches ∼ 2640 K f 1/4 .Therefore, annihilation heating contributes an additional ∼ 400 − 650 K to the NS temperature compared to kinetic heating alone, depending on the NS configuration.

Scattering and Capture Rates
In this paper, we restrict our attention to DM scattering with baryons, the most abundant particle species in NSs. 1 To compute the DM scattering and capture rates in neutron stars, we must account for a number of effects including relativistic kinematics, general relativistic corrections, the motion of the target particles in the NS, and Pauli blocking (relevant to the scattering of light DM) [22,24].In addition, the extreme conditions in the NS interior give rise to two additional important effects: (i) baryons experience strong interactions, which result in an effective mass different from their in vacuum value, and (ii) the momentum transfer in each collision is large enough to resolve the internal structure of the nucleon [25,27]. 2 A detailed account of the scattering and capture rate calculations can be found in refs.[22,24,25,27].Here, we briefly outline the aspects of the scattering calculation relevant to this paper.

Scattering Rate
The DM-baryon scattering rate is a key ingredient in both the capture and thermalization processes.It is used in constructing the DM energy loss probability distribution function, which, in turn, determines the probability that the DM is captured [22], and the average energy lost in each collision during the thermalization process.
The most general expression for the DM down-scattering rate, expressed in terms of the DMtarget response function S(q 0 , q) [22,41], is where S(q 0 , q) = 2 ) are the DM initial and final momenta, p µ = (E i , ⃗ p) and p ′µ = (E ′ i , ⃗ p ′ ) are the target particle initial and final momenta, q 0 = E ′ i − E i is the DM energy loss, ⃗ q = ⃗ p − ⃗ p ′ is the three-momentum exchanged in the scattering, m i is the mass of the target, |M (s, t, m i )| 2 is the spin-averaged squared matrix element, and f FD is the Fermi Dirac distribution.Note that the purpose of Θ(q 0 ) is to select only the down-scattering interactions.
When the energy transfer is small compared to the Fermi energy, Pauli blocking strongly suppresses the scattering rate.Since the targets are in a highly degenerate Fermi plasma, we can take the target energy levels to be either full or empty in the zero temperature limit.Therefore, the interaction rate depends on the number of targets with energy E i in the initial state, and on the number of empty final states with energy E i + q 0 , where q 0 is the dark matter energy loss.As a result, the interaction rate necessarily vanishes in the limit that q 0 → 0.More generally, Pauli blocking suppresses the differential interaction rate when E i + q 0 ≤ ε F,i .
To calculate the squared matrix elements, |M | 2 , we take the DM-quark couplings to be described by the four-fermion effective field theory (EFT) operators listed in Table 2, where the strength of the coupling is parameterized by a cutoff scale, Λ q .For each of these operators, the corresponding |M | 2 is expressed in terms of the Mandelstam variables s and t, and the DM-target mass ratio In the following sections, we shall derive approximations for interaction rates and timescales that depend on the form of the matrix element.We therefore define the parameter n to denote the t-dependence of the squared matrix element as with n = 0, 1, 2.
For DM interactions with baryons, the squared matrix elements contain hadronic coefficients, c I i (t), which depend on the transferred momentum t.In the case of scalar and pseudoscalar interactions, these coefficients also depend on the baryon mass m i .In general, these coefficients can be defined as a function of their values at zero momentum transfer as c I i (0) [54] where S, P, V, A and T denote scalar, pseudoscalar, vector, axial and tensor interactions, respectively.The coefficients c I i (0) are given in appendix A of ref. [27], while F (t) is the square of the dipole formfactor, Note that the energy scale Q 0 depends on the specific hadronic form factor.As in refs.[25,27], we conservatively assume Q 0 = 1 GeV for all operators and baryonic target species.The effect of the strong interactions between the baryons is incorporated in the proper calculation of the Fermi energies of the baryonic species, ε F,i , and in the use of baryon effective masses for the target masses.These depends on the microphysics embedded in the EoS [25,27].
For the collisions that result in DM capture, analytic scattering rate expressions can be obtained [22,24] because the DM kinetic energy is always significantly higher than the NS temperature and hence the zero temperature T ⋆ → 0 approximation holds.During the thermalization process, however, the DM kinetic energy and NS temperature become comparable, and so finite temperature effects become important.It is numerically intensive to calculate the interaction rate directly from Eq. 4.1 for these low temperatures, hence we keep only the lowest order thermal corrections.This amounts to a modification of the differential interaction rate, i.e. the q 0 integrand of Eq. 4.1, such that As a result, we can no longer obtain fully analytic expressions for the interaction rate.Therefore, all the results we present below have been calculated numerically.

Name
Operator [55] for the coupling of Dirac DM to quarks (column 2), together with the squared matrix elements for DM-hadron scattering (column 5), where s and t are Mandelstam variables, µ = m χ /m i , and m i is the hadron mass.The quark-level coefficients g q are expressed in terms of the Yukawa couplings, y q , and the cutoff scale, Λ q .The squared hadron-level coefficients, g 2 i , depend on the momentum-dependent couplings c I i (t), as in Eq. 4.5.The final column indicates the dominant term in the squared matrix element, in the low momentum-transfer limit relevant for the late stage of the thermalization process.

Capture Rate
For DM-baryon cross sections much smaller than a threshold value, which we denote as σ th iχ , a neutron star can be considered as an optically thin medium.In this regime, DM passes through the star and capture can occur anywhere in the NS interior.The capture rate can then be calculated using the following expression where ρ χ is the local DM density, u χ is the DM velocity far away from the NS, and f MB is the distribution of relative velocities between NS targets and DM particles far away from the NS, which we assume to be of Maxwell-Boltzmann form.We define Ω − (r) to be the interaction rate of Eq. 4.1 in the special case where the initial DM kinetic energy is fixed to , as relevant for a DM particle infalling to the NS.This can be expressed as where the integration intervals for are given in refs.[22,27], together with further details about the capture rate calculations.
If the DM-target cross section reaches or surpasses the threshold value σ th iχ , the NS optical depth τ χ can no longer be taken as close to zero.We must then introduce an optical factor, η(r) = exp[−τ χ (r)], in Eq. 4.8 to suppress the capture rate; see ref. [22] for further details.In the optically thick limit, the NS is seen as a rigid sphere by the DM particles, and all of the capture occurs on the surface of the star.This is the so-called geometric limit.In this case, the capture rate is given by [15]

Thermalization
After becoming gravitationally bound to the NS, the DM particles continue to scatter with NS targets, losing energy in each collision until reaching thermal equilibrium at the centre of the star.
We outline the calculation of the thermalization time in terms of the average DM energy lost in a single collision, and use first-order approximations to derive scaling relations that allow us to understand the qualitative features of our numerical results.

Average DM energy loss
The average energy a DM particle loses per collision can be calculated by weighting the DM energy loss, q 0 , with differential interaction rate.We thus obtain where q MAX 0 is the maximum energy lost in a single scatter.Figure 1 shows q MAX 0 as a function of the DM kinetic energy, K χ = E χ − m χ , for DM-neutron collisions.We see that heavier DM particles lose a smaller fraction of their kinetic energy per collision than lighter DM3 .Nevertheless, as K χ approaches the Pauli blocked region, K χ ≪ m i ε F,n /m χ (dashed blue line), the maximum energy loss per collision becomes independent of the DM mass.
For the initial collision that results in capture, Pauli blocking represents, at most, a sub-leading correction to the capture rate for DM masses above the Fermi energy of the targets.Following capture, however, the DM energy will continue to decrease as a result of subsequent scattering, eventually reaching kinetic energies where Pauli blocking is an important effect.Consequently, Pauli blocking will strongly impact the rate at which dark matter is thermalized, for a wide DM mass range that extends well above ε F,i .
It is useful to define a critical DM mass, above which Pauli blocking is never in effect throughout the entire thermalization process.We do this by analysing the regions of the interaction rate phase space that are suppressed by Pauli blocking, arriving at 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 K χ (GeV) q MAX 0 (GeV) For neutron targets with ε F,n = 200 MeV, and assuming an equilibrium temperature of 10 3 K, i.e.K χ ≳ 10 3 K, we find (neglecting the nucleon form factors) m crit χ ∼ 9.65 × 10 9 GeV. 4 Pauli blocking will then suppress at least some part of the thermalization process for all DM masses below this value.
In either regime, we can obtain first-order approximations for the average fraction of energy that a DM particle loses in a single collision by making use of the zero temperature approximation, a constant target mass, and nucleonic form factors at zero momentum transfer, i.e., F (t) ∼ 1.First, we consider the regime in which Pauli blocking is negligible, m χ ≳ m crit χ .For a constant cross section (i.e.|M | 2 ∝ t 0 ) we find at first order in m crit χ /m χ .(See appendix A for the corresponding approximation of the interaction rate, Eq.A.16.)For cross sections proportional to t 1 and t 2 , the average energy losses can be obtained in the same way, starting from the relevant expressions for dΓ dq 0 .Figure 2 shows the average energy loss fraction per collision.We see the ⟨∆K χ phase of the evolution, where the kinetic energy is driven down toward values where Pauli blocking eventually becomes active.The latter Pauli blocked phase is indicated by the horizontal arrows in Fig. 2.
Moving to the case where Pauli blocking suppresses the scattering rate, m χ ≲ m crit χ , the average energy loss per collision for the case of a constant cross section is (5.4) The average energy loss now scales linearly with K χ (the flat regions in Fig. 2) in contrast to the K 1/2 χ dependence of Eq. 5.3.As the DM kinetic energy decreases, the average fraction of energy transferred to the targets progressively increases until K χ no longer satisfies Eq. 5.2 and consequently 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 K χ (GeV)  the interaction rate becomes Pauli blocked.As expected from Eq. 5.2, the Pauli-suppressed region starts at higher kinetic energy for lower DM masses.
For interactions with t-dependent matrix elements, the average energy loss per collision also scales linearly with K χ in the Pauli blocked regime.For |M | 2 ∝ t n , with n = 1, 2, we find respectively, where we have used Eqs.A.11 and A.12 for the interaction rates.Note that the average energy loss fraction per collision exhibits similar behavior for all the interaction types considered, as seen by comparing the left and right panels of Fig. 2.This is true in both the Pauli-blocked and non-blocked regimes.

Thermalization timescale
Once a DM particle is captured, it becomes gravitationally bound to the NS and follows an orbit that may or may not lie completely within the NS.If the orbit lies partly outside the NS, subsequent scatterings will be required in order for the DM particle to lose enough energy so that the complete orbit lies within the NS.This is the first stage in the thermalization process.When estimating the amount of time needed for the DM orbit to lie completely within the star, we find that this time is always much shorter than the full time required for DM to reach thermal equilibrium with the neutron targets.Consequently, this first step in the thermalization process can be safely neglected.This finding is in agreement with ref. [16].
We shall also assume that up-scattering of the DM to larger kinetic energy does not play an important role 5 .These effects will become relevant as the DM approaches thermal equilibrium, increasing the thermalization time.We estimate that up-scattering will, at most, increase the thermal equilibrium time by O(10%), and thus we neglect this correction.
For DM of mass much larger than the target mass, m χ ≫ m eff i , there is an additional stage in the thermalization process where either Pauli blocking plays no role, or the interaction rate has a different power law relationship with the temperature than those identified in Section 5.1.These initial scatterings make a negligible contribution to the thermalization time, as Γ − is a sharply decreasing function of the DM kinetic energy K χ .
Let us denote the number of initial collisions prior to reaching the Pauli blocked regime as N 1 , and the number of additional collisions required for complete thermalization as N 2 .For light DM, m χ ≲ m eff i , Pauli blocking affects the entire thermalization process, i.e.N 1 = 0. Let K N be the kinetic energy after N scatterings.After N 1 + N 2 collisions, the DM will reach the equilibrium temperature T eq , which can be written as where we have used the fact that the average fractional energy loss is the same in each collision.The thermalization time can then be defined as the sum of the average time between collisions, up until the final energy transfer is equal to T eq [41] t For m χ ≲ m crit χ , the fraction of energy lost in the last few scatters is still a considerable fraction of the DM kinetic energy before the collision.Furthermore, these scatterings may take a considerably long time to occur, indicating that the process is discrete.As an example, consider thermalization to a temperature of 10 3 K, for a DM particle of mass m χ = 1 TeV and constant cross section σ nχ ∼ 10 −45 cm 2 . 6Hundreds of collisions are required to fully thermalize; the last 10 or so are spaced longer than a second apart; and the last couple are longer than 10 kyr apart.
To compute the thermalization time, we numerically integrate Eq. 4.7 to obtain the interaction rate and the average energy lost in each collision for each of the EFT operators in Table 2.We find that the thermalization time for a particular interaction type scales according to the dominant power of the Mandelstam variable t in the corresponding matrix element; see the last column of Table 2. Thus, to understand how the thermalization times scale, it is enough to consider differential cross sections that are proportional to a given power of t, i.e. dσ ∝ t n , with n = 0, 1, 2. Below, we present results for operators that depend only on t (D1-4 in Table 2), and not on the centre of mass energy, s.In Appendix C, we outline the procedure used to obtain analytic expressions for those operators with an explicit dependence on s (operators D5-10).
In Figure 3, we show the full numerical results for the thermalization time as a function of the DM mass for different equilibrium temperatures.It is clear that the power law scaling of the thermalization time with DM mass depends on whether m χ is larger or smaller than the nucleon mass.In order to understand these results, we make use of the analytic approximations for the average DM energy loss derived in Section 5.1 and valid in the zero temperature limit.
We begin by studying the Pauli blocked regime, m χ ≲ m crit χ .In this case, the majority of the thermalization time is dictated by the final few scatters, for which the form factors are close to their value at zero momentum transfer.These last collisions occur close to the NS centre, so we can take the target mass as constant and equal to the value at the centre of the star, m eff i (0).For the case of a constant DM-neutron cross section, dσ ∝ t 0 , the thermalization time can thus be obtained by using Eqs.A.6, 5.6 and 5.4 in Eq. 5.7.This leads to where σ n=0 iχ is the DM-baryon cross section.The scaling of this expression with T eq , DM mass, and DM-target cross sections agrees with ref. [41].Numerically we obtain where we have set m eff i (0) = 0.5 m n .The m χ dependence of these expressions explains the features observed in the top left panel of Fig. 3.For m χ ≪ m eff n (0), the thermalization time scales with the DM mass as the number of scatterings needed for thermalization increases.Conversely, for therm is inversely proportional to m χ due to the reduced Pauli blocking, explaining the change of slope around the value of the neutron effective mass in the NS centre.
We repeat the same analysis for cross section proportional to higher powers of t, i.e., dσ ∝ t n with n = 1, 2. Using Eqs.A.11, A.12 and 5.5, we find therm ∼ 17 where the factors involving B(R ⋆ ) arise from fixing the cross section to its value at the surface; see Appendix A for details.
To gain insight into the order of magnitude of these thermalization times, we set B(R ⋆ ) = 0.5 and m eff i (0) = 0.5 m n , yielding (5.12) , m χ ≫ m eff i (0). (5.13) As anticipated, we see that the momentum-suppressed cross sections translate into significantly longer thermalization times than for the case of a constant (unsuppressed) cross section.These expressions also allow us to understand the dependence of t therm on the DM mass.For dσ ∝ t n , the thermalization time scales as m n+1 χ for m χ ≪ m eff n (0), and as the inverse of this quantity for m χ ≫ m eff n (0).The choice of EoS has a small but non-negligible impact on the thermalization time, as indicated by the widths of the shaded regions in Fig. 3.For a constant cross section, we observe almost no variation in t therm with the NS configuration, except for the m χ ≲ m n region.This is due to the dependence of m eff n (0) on the NS model; see Table 1 and Eq.5.8.For cross sections dσ ∝ t n , with n = 1, 2, the dependence on B(R ⋆ ) in Eqs.5.10 and 5.11 adds an extra dependence on the choice of NS model.For these momentum-suppressed interactions, DM requires more time to reach an equilibrium temperature in heavier NSs.This is due to the combination of two effects: the effective mass of the targets in the centre of the NS is smaller in more massive NS configurations, while B(R ⋆ ) increases.Nonetheless, the dependence on NS configuration remains relatively mild.
We now turn to the m χ ≳ m crit χ regime, which is observed only for temperatures above 10 4 K in Fig. 3.This regime change is indicated by the change of slope that occurs at large DM masses, clearly evident for T eq = 10 5 K (orange) at a DM of mass m χ ≳ 5 × 107 GeV.In this regime, the energy lost in each collision is a tiny fraction of the initial DM kinetic energy, and the time between scatterings is of order a fraction of a second.This indicates that a continuous approximation becomes more appropriate than a discrete sum to estimate t therm .In this case, the momentum-dependent part of the form factor, F (t), will be relevant only at the beginning of the thermalization process and become less and less relevant as the average momentum transfer decreases in each subsequent scatter. 7It is these low momentum-transfer collisions that dominate the thermalization time.For a constant cross section (n = 0), in the zero temperature approximation, we obtain (see Appendix B for details) In this super heavy DM mass regime, we see that the thermalization time is an increasing function of m χ .
It is worth remarking that for a constant DM-neutron cross section (top left panel) thermalization will always occur within the age of the Universe.However, this is not true for momentumsuppressed cross sections, for a range of DM masses.Specifically, for T eq = 10 3 K and the assumed reference cross section, DM of mass m χ ≲ 10 TeV (m χ ≲ 1 PeV) will not have enough time to thermalize for dσ ∝ t (dσ ∝ t 2 ).Importantly, however, we shall see below that even when full thermalization takes longer than the age of the Universe, the majority of the kinetic energy is deposited on a much shorter timescale.
Finally, we must incorporate the fact that DM will scatter with various baryonic species in the NS rather than just the neutrons.To do this, we combine the thermalization times for scattering from each individual species in an appropriate way.Specifically, we sum the inverse single-species thermalization times, weighted by their relative abundance at the centre of the NS, such that where Y i (0) is the abundance of the species in the centre of the NS, and the sum runs over all possible baryons.For the case of the heaviest NS we consider, 1.9 M ⊙ , this includes the Λ 0 , Ξ 0 and Ξ − hyperons.The resulting thermalization time then lies between the fastest and slowest single-species times, as is expected.
6 Capture-Annihilation Equilibrium The captured DM will accumulate in the centre of the NS, where it will begin to annihilate.The annihilation rate will grow until sufficient time has elapsed for the capture and annihilation processes to reach equilibrium.In this limit, the total amount of DM in the NS is maximized, and will remain constant.Once this occurs, annihilation efficiently deposits the DM mass-energy into the star.Let us begin by assuming that the dark matter has fully thermalized.After reviewing the standard capture-annihilation equilibrium calculation, we will relax this assumption to consider the more general case of partially thermalized dark matter and derive new expressions that hold in that scenario.Importantly, we shall see that capture-annihilation equilibrium, and hence efficient annihilation can occur without full thermalization.

Capture-annihilation equilibrium of thermalized dark matter
The thermalized DM will collect within an isothermal sphere at the centre of the NS where it will begin to annihilate.The efficiency of the annihilation will depend on the volume of this sphere, which is expected to be very small for the DM masses we consider.Very close to the centre of the NS, the density does not vary significantly and can be taken to be constant.Then, working in the weak field approximation such that B(r) ∼ 1 + 2Φ(r), we can obtain the gravitational potential inside the NS, where ρ c is the central density of the NS.The number density of DM particles that have thermalized to a temperature T eq as a function of radius will then be given by a Maxwell-Boltzmann distribution where N χ is the total number of DM particles within the isothermal sphere, and r iso is the radius of the DM isothermal sphere.Applying the viral theorem leads to the following expression for r iso , The total number of DM particles enclosed in this sphere is then and the velocity distribution of the thermalized DM is given by In the absence of evaporation, which can safely be neglected for m χ ≳ 1.17 × 10 −8 GeV for a Gyr old NS with core temperature ∼ 10 3 K [24], the time evolution of the total number of DM particles present inside the NS is governed by Table 3: Thermally averaged annihilation cross sections ⟨σ ann v χ ⟩ for the dimension 6 EFT operators, expanded to second order in v χ .The y q factors are the fermion Yukawa couplings [56].
Here C is the capture rate and A is related to the DM annihilation rate, Γ ann , through where and ⟨σ ann v χ ⟩ is the thermally averaged DM annihilation cross section.These are given in Table 3 for the EFT interactions we consider, where ⟨v 2 χ ⟩ = v 2 th = 6T eq /m χ .We note that the cross sections shown in Table 3 are quark-level expressions.For most of the mass range of interest, these provide excellent approximations to the hadron-level annihilation cross sections, provided we impose a lower bound on the DM mass for which an annihilation channel is open, taken to be the pion mass.See Appendix E for details.
The solution to Eq. 6.6 in terms of the capture and annihilation rates is Ultimately, we are interested in the behavior of Eq. 6.9 at late stages in the NS evolution, i.e., for t → t ⋆ , where t ⋆ is the age of the NS, which we take to be ∼ 1 Gyr.In this limit, the hydrostatic NS structure, and hence C, are not expected to change with time.Of particular interest is whether or not an equilibrium is reached between the capture and annihilation rates.Such a state is reached for timescales greater than For earlier times, t < t eq , one can neglect the loss of DM particles due to annihilation, leaving N χ ∼ Ct.

Capture-annihilation equilibrium of partially-thermalized dark matter
The standard calculation of the annihilation rate, using Eq.6.8, assumes that the DM has thermalized, i.e., t > t therm .If thermalization has not been achieved by a time t ∼ t ⋆ < t therm , the DM kinetic energy distribution will peak around the lowest temperature that DM has had enough time to reach.This is given by , (6.11) where n is the exponent of the dominant t n term in the differential cross section, dσ ∝ t n , as given in the last column of Table 2. See Appendix D for details.We can then find the radius of the DM distribution (which is no longer isothermal) and the ⟨σ ann v χ ⟩ corresponding to the peak of the energy distribution K χ .We obtain A via the replacement where α = 3/2 for s-wave annihilation, and α = 1/2 for p-wave.Making this replacement in Eq. 6.10 leads to a capture-annihilation equilibrium time of No previous estimate of the capture-annihilation equilibrium time has considered the case of partially thermalized DM.If thermalization has not been achieved, the additional factor in Eq. 6.13, compared to Eq. 6.10, increases the equilibrium time.(Or, equivalently, increases the cross sections required to reach equilibrium within a specified time.)However, it is critical to realize that t eq can be shorter than t therm .In fact, it is possible for annihilation to occur efficiently, even if complete thermalization never occurs.In this scenario, we must use Eq.6.13.
Assuming the DM is captured at the geometric limit, we arrive at the following result for our benchmark NS QMC-2 T eq 10 3 K

.14)
Comparing this expression with the thermalization times in the previous section, we anticipate that t eq will typically be shorter than t therm , often by many orders of magnitude.

Kinetic heating timescale
It has commonly been assumed that the DM kinetic energy deposition occurs instantaneously.However, it is not immediately obvious that this is true.In particular, for scattering interactions suppressed by powers of the momentum transfer, t, full thermalization can take longer than the age of the Universe.We must therefore determine whether a significant fraction of the initial kinetic energy can be deposited on a shorter timescale.To quantify the timescale on which kinetic heating takes place, we define t kin to be the time required for a DM particle to lose 99% of its maximum kinetic energy, K χ = m χ (1/ B(0)−1).This time is calculated in the same way as the thermalization time, while also keeping track of the time the DM spends outside the star along its orbit, i.e., the initial stage of the thermalization process that was neglected in Section 5.2.For simplicity, the DM particle orbits are taken to be linear, passing through the centre of the star, with the radial extent calculated using the geodesic equations.We expect an O(1) correction to our results when considering circular orbits.Additionally, we randomize the radial position in the NS where the DM interacts with a target.
Figure 4 shows the time required for kinetic heating to be achieved, assuming DM-neutron interactions of the form dσ ∝ t n , with cross sections normalized to σ nχ = 10 −45 cm 2 at the surface of the star.As the location of each interaction is randomized, these results are obtained by averaging over several simulations for each DM mass.For light DM, t kin decreases with increasing DM mass, due to the decreased effects of Pauli blocking, with the change of slope at m χ ∼ 0.1 GeV indicating the point where Pauli blocking affects only a fraction of the total process.For masses ≳ 10 GeV, Pauli blocking is not relevant this part of the thermalization process, and hence the t kin monotonically increases with the DM mass, as was seen in the thermalization of super-heavy DM.
Fig. 4 illustrates two key facts.First, t kin differs by orders of magnitude for the different cross section types, dσ ∝ t n , with larger values of t kin for the most highly momentum-suppressed interactions, as expected.Second, and importantly, the kinetic heating occurs relatively quickly for 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 m χ (GeV) p−wave p−wave t kin t therm t eq 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 m χ (GeV)  The operator D1 has an unsuppressed scattering cross section and a p-wave suppressed annihilation cross section, while D4 has a q 4 tr suppressed scattering cross section and an unsuppressed (s-wave) annihilation cross section.The interaction strength has been chosen to give maximal capture.(Specifically, we used Λ q values corresponding to a capture cross section at the geometric limit, assuming scattering with the neutron targets in the QMC-2 NS.) all interaction types, on timescales much shorter than a typical NS age.Indeed, for the case of a constant cross section, t kin is much shorter than a year.

Annihilation heating timescale
Figure 5 shows all relevant timescales for DM-induced heating of old neutron stars.These timescales have been calculated considering DM scattering off the neutron targets in the benchmark NS QMC-2, for the case of maximal capture, f = 1 (i.e., we have set the EFT parameter Λ q as required to achieve capture at the geometric limit).We show these results for two indicative operators: The scalar-scalar interaction D1 (left), which has a p-wave suppressed annihilation cross section, and the pseudoscalar-pseudoscalar operator D4 (right), which has an s-wave annihilation cross section.
As anticipated, capture-annihilation equilibrium takes longer to achieve for the velocity-suppressed p-wave annihilation cross section than for the s-wave.Nonetheless, equilibrium (and hence maximal annihilation heating) is reached relatively quickly in both cases, on timescales of 10 4 years for the scalar interaction, and even quicker for the pseudoscalar, well within the age of a typical NS.
For both interaction types, the kinetic and annihilation heating contributions are both realized on timescales much shorter than that required for full thermalization.If the scattering cross section is momentum suppressed (as with the dσ ∝ t 2 = q 4 tr dependence for D4), the thermalization time is increased; if the annihilation cross section is velocity suppressed (as with the p-wave annihilation cross section for D1) the capture-annihilation equilibrium time is increased.
Finally, note that there are parameters for which the annihilation timescale t eq is shorter than the kinetic heating timescale t kin .In this case, the annihilation process deposits the DM mass energy and any remaining kinetic energy, which is carried by the annihilation products.Therefore, the minimum time required for DM to deposit all of its energy, both kinetic and rest-mass, is t eq .

Neutron star heating sensitivity for various interaction types
We now examine the regions of parameter space where maximal heating can be achieved, for DMhadron interactions described by the four-fermion operators of Table 2.As the extent of DM-induced heating depends on how efficiently the DM is captured, it is clear that maximal heating corresponds to the case of f = 1 in Eqs.3.2 and 3.4, i.e., when the dark matter scattering cross section is at or above the geometric limit.
In Figure 6, we show the parameters for which the DM deposits its entire kinetic and rest mass energy (yellow region) or only its full kinetic energy (light blue region).Above these shaded regions, the value of Λ q is too large (and hence the scattering cross section is too small) for maximal capture.The overall shape of the shaded regions is dictated by the behaviour of the capture rate.For sub-GeV DM, Pauli blocking suppresses the capture rate, and so smaller Λ q values are required to reach the geometric limit.At large DM mass, ≳ 10 5 GeV, the capture rate is suppressed because a single collision does not transfer enough energy to result in capture.
We see from the t eq contours (magenta) that capture-annihilation equilibrium, and hence full annihilation heating, is achieved on a timescale much shorter that the NS lifetime, which we take to be t ⋆ ∼ 1 Gyr.We do not show the contours for the kinetic heating timescale, t kin , as this is significantly shorter than the age of the star.Moreover, for masses where DM annihilation to hadrons is possible, m χ > m π , it is not necessary for kinetic heating to occur before capture-annihilation equilibrium is established, as the total kinetic plus mass energy will be deposited when the DM annihilates.For completeness, we include the contours of the thermalization time (grey).We stress again that the DM does not need to fully thermalize to achieve maximal heating.
To highlight how the sensitivity varies for different interaction types, we show results for operators D1, D4 and D5 in Figure 6, with the remaining operators presented in Figure 8.These operators were chosen because they allow us to compare interactions with and without momentum or velocity suppressed scattering or annihilation cross sections.Specifically: • D1 (scalar) -unsuppressed scattering cross section; p-wave suppressed annihilation.
Comparing the projected NS heating sensitivity with limits from terrestrial direct detection experiments (shown as green and orange curves in Fig. 6) we find similar behaviour for D1 and D5.This is expected, as both of these operators give rise to unsuppressed spin-independent DM-nucleon scattering cross sections.The p-wave suppression of the D1 annihilation cross section increases t eq with respect to that for D5; nonetheless, equilibrium is reached relatively quickly compared to t ⋆ .The D4 (pseudoscalar) interaction has dismal prospects of being observed in direct detection experiments, due to the severe q 4  tr suppression for the scattering of non-relativistic DM.In contrast, NS heating has much greater sensitivity.
Note that the time required for complete thermalisation (grey contours) is much longer for D4 (momentum suppressed scattering) than for D1 and D5.In fact, for operator D4, full thermalisation is not achieved for most of the interesting parameter space.This illustrates the importance of correctly identifying t eq as the timescale on which full heating is achieved, rather than the much longer t therm .10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 m χ (GeV)  We have used the QMC-2 (1.5M ⊙ ) benchmark NS configuration.We show the regions where DM kinetic and annihilation heating both contribute to the NS luminosity (yellow) and where kinetic heating alone contributes (light blue).Contour lines for capture-annihilation equilibrium (magenta) and thermalization (grey) are shown for indicative timescales.Lower limits on Λ q from leading direct detection experiments [57][58][59][60] are also shown.
A Interaction rate in the zero temperature approximation In this section, we calculate the interaction rate in the zero temperature approximation for |M | 2 = α t n , where n = 0, 1, 2 and α is a constant, in the low energy, Pauli suppressed regime where K χ = E χ − m χ < ε F,i .We assume the simplest scenario of constant target mass and point-like targets, as justified in Section 5.2.
In this approximation, the interaction rate in the energy regime relevant for thermalization is given by [22] Γ where t E = −t, and The D n functions can be found in Appendix B of ref. [22].As we shall see, the integration intervals in Eq.A.1 depend on whether or not Pauli blocking suppresses any part of the thermalization process.
In both cases, we can find simple analytic approximations to these integrals.The minimal DM mass for which Pauli blocking is never in effect is denoted by m crit χ .We first consider the case where m χ ≲ m crit χ .For the cases of µ ≪ ε F,i /K χ or µ ≫ K χ /ε F,i , Γ − at first order in K χ is given by where t ± E are defined in ref. [22].For matrix elements independent of t (n = 0), we have D 0 (q 2 0 , t ± E ) = 1 and this result simplifies to We can rewrite the previous expression in terms of the DM-baryon scattering cross-section using the following expression giving the interaction rate at first order in K χ This result has the same K χ and µ scaling as that of ref. [41].
Performing a similar analysis for |M | 2 = α(−t) n , n = 1, 2, we find The expressions for the cross sections for n = 1, 2 are These cross sections must be normalized to sensible momentum transfer.We take this reference point to be the surface of the star, such that The interaction rates for n = 1, 2 can then be written as (A.12) We now look at the interaction rate in the super-heavy DM mass regime, m χ ≳ m crit χ .The exact value of m crit χ will depend on the NS configuration.However, we can take some typical values relevant to thermalization to give an estimate of its value.Taking K χ = 10 3 K, ε F,i = 200 MeV, we see that The maximum energy transfer in this regime will always be q MAX 0 < K χ , with Performing a similar analysis as the m χ ≲ m crit χ regime leads to the following expression for Γ − , where t + µ − is defined in ref. [22].For the simplest case of constant |M | 2 this results in B Thermalization of super-heavy DM For DM that is heavier than the critical mass m χ ≳ m crit χ , the energy lost in each scatter is a tiny fraction of the total DM kinetic energy.Moreover, the average time between collisions is typically on the order of fractions of a second.This warrants the use of a continuous approximation in this regime rather than performing the discrete summation.The thermalization time is then found by integrating the rate at which the DM kinetic energy changes, from the initial kinetic energy, − 1 , to the final value T eq ≪ m χ .For a constant cross-section (n = 0), we substitute Eqs.A.16 and 5.3 into the expression above leading to Taking the final temperature to be T eq = 10 3 K and B(R ⋆ ) = 0.5, this yields Repeating for dσ ∝ t n (n = 1, 2), we calculate the thermalization time for n = 1 to be and that for n = 2 to be C Thermalization time for s-and t-dependent interactions In Section 5.2, we assumed |M | 2 ∝ t n when deriving analytical approximations for the thermalization timescale.To understand the behavior of the thermalization time for the operators in Table .2, we can make use of the results for t n dependent interactions.For cross sections that are linear combinations of different powers of t, we can approximate the thermalization time using the previous results in the following way Hence, the inverse of the thermalization time will be given by a weighted linear combination of the inverse times for each contribution.As higher powers of t require significantly longer thermalization times, for coefficients of similar size, the resulting sum will be dominated by the lowest power of t appearing in |M | 2 .We can thus identify the dominant terms for operators D1-D4 based on power counting, which we have listed in Table 2.
For s-dependent amplitudes, we can in principle use the interaction rates calculated in Appendix A of ref. [24], perform a series expansion in K χ and repeat the same procedure outlined in Section 5.2 for s-independent matrix elements.Interestingly, we find that for the purpose of calculating the thermalization time, there is an easier way to obtain the correct result.One can indeed check that, at zero order in ε F,i /m eff i , the resulting time for s 1 , s 2 is equivalent to the constant case, with the matrix element calculated by setting while the st case has a result equivalent to the t case, with the matrix element calculated using the same substitution.This is, in practice, equivalent to setting both the DM and neutron targets at rest.There is, however, an important exception, when it comes to calculating the thermalization time of a linear combination of these terms.In particular, when the amplitudes at order O(t 0 ), are proportional to combinations of 1, s, s 2 such as All these combinations give a null result after applying substitution C.4.In such a case, one may think that the dominant term is given by some remaining t n term.It is worth noting that the expressions in Eq.C.5 appear in operators that, at low energy, are known as velocity-dependent, because their matrix elements are proportional to positive even powers of the DM-target relative speed.Consequently, it is important not to neglect the motion of the targets in the neutron star, moving at relativistic speeds that are of the order of the Fermi velocity v 2 F = 2ε F,i /m eff i .In those cases, one should instead set s to In summary, operators D5, D8 and D9 can be safely expanded using C.4, while operators D6, D7 and D10 have velocity-dependent amplitudes and require Eq.C.6.The dominant terms for each operator can be found in Table 2.For equal values of the leading term in |M | 2 , the thermalization time for each operator will be the same as the relevant t n power law.

D Temperature distribution of captured dark matter
As seen in Fig. 3, interactions that depend on the momentum transfer, namely dσ ∝ t n with n = 1, 2, there are regions of the DM mass parameter space where thermalization does not occur within the age of the star.For the DM masses and NS temperatures of interest, this region of non-thermalization always occurs in the m χ ≪ m crit χ regime.From Eqs. 5.8, 5.10 and 5.11, we can estimate the time required for the DM to reach a kinetic energy K If the DM does not thermalize within the age of the star, it will instead reach a minimum temperature, K min χ .Comparing the time required to achieve this temperature to the thermalization time, t therm i.e. to have reached the equilibrium temperature T eq , we find Accounting for the case where the DM reaches thermalization, we can write The population of captured DM will have a distribution of energies at any given time, with this distribution being peaked at this minimum energy.As the orbital periods of the DM will be much shorter than the average time between interactions, the DM will be able to virialize between each interaction.Therefore, we can treat the DM as being contained within an isothermal sphere with temperature K min χ > T eq .Finally, it is worth noting that at times t > t therm , even though the thermalization condition has been reached, the captured DM would consist of two components: a fraction of it (whose amount depends on time) would be in thermal equilibrium with the NS at temperature T eq ; and another component still in the cooling down process.Assuming a capture rate constant over time, the fraction of thermalized DM is f therm (t) = t − t therm t . (D.5)

E Quark-level vs hadron-level annihilation cross sections
The annihilation cross sections shown in Table 3 are for DM annihilation to quark final states.More properly, we should consider the hadron-level annihilation cross section.However, we are primarily concerned with the capture-annihilation equilibrium timescale, and not the details of the annihilation process.Therefore, if the annihilation rate to hadrons is not significantly different from the quark level result, this subtlety can be avoided.
For DM masses in the range m π < m χ ≲ m charm = 1.27GeV, we find that the cross section for annihilation to hadrons differs by less than an order of magnitude than that for annihilation to quarks.For larger DM mass, the difference is negligible.Therefore, to simplify the discussion, we consider DM annihilation to quark final states for DM masses above the pion mass.Below the pion mass, the only kinematically allowed DM annihilation channels would be to leptons or photons.The size of the DM couplings to these states would, in general, be unrelated to the DM-quark couplings we have assumed.(They are expected to be non-zero, because they would be induced at loop level [17], even if absent at tree level.)However, due to the considerable Fermi energies of the electrons and muons near the centre of the NS, these channels will be Pauli blocked for the whole DM mass range below m π , forbidding these annihilations from occurring.To remain as model-independent as possible, we will not consider lepton and photon annihilation channels.
Figure 7 shows t eq contour lines in the Λ q − m χ plane for the NS benchmark model QMC-2, T eq = 1000 K and two representative operators D7 (left) and D8 (right).Operators whose thermally averaged annihilation cross section ⟨σ ann v χ ⟩ has a m q /m χ leading order term, namely D1-D4 and D8 (see Table 3), exhibit a sudden change in the slope wherever a new annihilation channel opens (see dotted lines on the right panel of Fig. 7).Note that the higher the cutoff scale Λ q , the lower the scattering and annihilation cross sections, resulting in a larger t eq timescale.For lower T eq temperatures, DM requires more time to reach both equilibrium conditions, thermalization and capture-annihilation.The variation of these results with respect to the NS configuration amounts at most to a factor of ∼ 2 in the t eq contours (see shaded regions in the left panels) from the lightest configuration (QMC-1, 1M ⊙ ) to the heaviest (QMC-3, 1.9M ⊙ ) for most operators, with the sole exception of D4 for which this factor rises up to ∼ 2.4.

F Capture-annihilation equilibrium for the EFT operators
In Fig. 8, we show isocontours of maximal capture (yellow and light blue lines) and captureannihilation equilibrium (magenta lines) timescale in the Λ q − m χ plane for all EFT operators.Values of Λ q below the t eq lines result in smaller capture-annihilation equilibrium timescales.We can see in Fig. 6 that for all operators the t eq timescale is always smaller than the time required for captured DM to thermalize.Captured DM achieves the steady state condition in a timescale as short as ∼ 1 yr (D1, D6-D10) or as long as 10 5 yr (D3).Note that the displayed values of t eq are the values for which the entire parameter space relevant for capture reach equilibrium with annihilation.
For D1 and D6-D10, we observe that captured DM achieves thermal equilibrium in less than ∼ 1 Myr (grey lines).On the other hand, as expected from Fig. 3 captured DM whose interactions are momentum suppressed, namely operators D2-D4, would never thermalize to the temperature expected from DM induced heating in 1 Gyr, or even in less than the age of the Universe, with the sole exception of the corner region of very light DM m χ ≲ 2 MeV and an even narrower corner of the parameter space for D4.
Therefore, for all EFT operators, the energy released in the annihilation process adds up to the energy deposited via capture increasing the DM induced heating for m χ ≳ m π (yellow shaded area).Recall that we have made no assumptions about the energy scale that controls DM interactions with leptons.For DM of mass below m π (light blue lines), at least kinetic heating is expected to contribute to the star luminosity.

Figure 1 :
Figure 1: Maximum energy loss per collision with neutron targets, as a function of the DM kinetic energy.We have assumed ε F,n = 200 MeV.

Figure 2 :
Figure 2: Average fraction of energy loss per DM-neutron collision for constant cross section (left) and dσ ∝ t (right) as a function of the DM kinetic energy.Horizontal arrows indicate the Pauli blocked (PB) regime, m χ ≲ m crit χ .

Figure 4 :
Figure4: Timescale on which the DM deposits 99% of its initial kinetic energy in the NS.We have assumed a NS with configuration QMC-2, and a DM-neutron scattering cross section of σ nχ = 10 −45 cm 2 at the surface of the star.

Figure 5 :
Figure5: Timescales for kinetic heating (blue), thermalization (orange) and capture-annihilation equilibrium (magenta), for operators D1 (left) and D4 (right).The operator D1 has an unsuppressed scattering cross section and a p-wave suppressed annihilation cross section, while D4 has a q 4 . + A n n .H e a t i n g t eq = 100 yr t eq = 100 yr tthe rm = 1 G y r tthe rm = 1 G y r

Figure 6 :
Figure6: Projected NS heating sensitivity for maximal capture efficiency, for DM-baryon interactions described by operators D1, D4 and D5.We have used the QMC-2 (1.5M ⊙ ) benchmark NS configuration.We show the regions where DM kinetic and annihilation heating both contribute to the NS luminosity (yellow) and where kinetic heating alone contributes (light blue).Contour lines for capture-annihilation equilibrium (magenta) and thermalization (grey) are shown for indicative timescales.Lower limits on Λ q from leading direct detection experiments[57][58][59][60] are also shown.

Figure 7 :
Figure 7: Contours of the capture annihilation timescale, t eq , in the Λ q −m χ plane for operators D7 (left) and D8 (right) and T eq = 1000 K. Solid lines represent the calculation for the NS benchmark model QMC-2, and shaded regions denote the variation with the NS choice for the QMC EoS family.Dotted lines in the right panel indicate the mass thresholds for various annihilation channels.

Table 1 :
[25,27]ies of the benchmark NS configurations as determined in refs.[25,27], using the QMC equation of state and varying the central baryon number density n c B .