Improved Treatment of Dark Matter Capture in Neutron Stars II: Leptonic Targets

Neutron stars harbour matter under extreme conditions, providing a unique testing ground for fundamental interactions. We recently developed an improved treatment of dark matter (DM) capture in neutron stars that properly incorporates many of the important physical effects, and outlined useful analytic approximations that are valid when the scattering amplitude is independent of the centre of mass energy. We now extend that analysis to all interaction types. We also discuss the effect of going beyond the zero-temperature approximation, which provides a boost to the capture rate of low mass dark matter, and give approximations for the dark matter up-scattering rate and evaporation mass. We apply these results to scattering of dark matter from leptonic targets, for which a correct relativistic description is essential. We find that the potential neutron star sensitivity to DM-lepton scattering cross sections greatly exceeds electron-recoil experiments, particularly in the sub-GeV regime, with a sensitivity to sub-MeV DM well beyond the reach of future terrestrial experiments.


Introduction
The quest to identify the cosmological dark matter (DM) is one of the forefront goals of modern science. In recent years, terrestrial dark matter direct detection experiments, which search for nuclear or electron recoil signals, have provided increasingly sensitive limits on the strength of dark matter interactions with regular matter. These experiments are limited, however, by the size of the detector target mass that can be practically realized. For this reason, it makes sense to consider alternative targets with which dark matter can interact, such as stars and planets. If the interaction of dark matter with these objects could be detected, they would offer a highly sensitive probe of the interaction strength, because the drawback of having to deal with uncertain astrophysical inputs is more than compensated for by the enormous target mass.

Internal structure
NSs are primarily composed of degenerate nuclear matter. Several layers and phase transitions can be found in their interior from a thin atmosphere up to the innermost core layer. The locally homogeneous core accounts for ∼ 99% of the NS mass [50,51]. The outermost layer, the crust, even though it is just ∼ 1 km thick, is the place where thermal conduction occurs and is hence responsible for the temperature drop between the core and the surface.
The outer crust is made of ionised heavy elements in a Coulomb lattice and strongly degenerate electrons (similar to a white dwarf). At densities ρ 10 6 g cm 3 , the electron gas is already ultrarelativistic. Moving from the surface towards the star interior, the increasing density induces electron capture and the nuclei become more and more neutron rich until the neutron drip density ρ ND ∼ 4.3 × 10 11 g cm −3 is reached [52][53][54][55][56]. This density marks the onset of the inner crust, where free neutrons, dripped off neutron rich nuclei, coexist with nucleon-proton clusters and electrons. At about half the nuclear saturation density (ρ 0 = 2.8 × 10 14 g cm −3 ), the nucleon clusters dissolve into their constituents as we cross the crust-core interface.
Matter in the outer core is mainly composed of neutrons in a superfluid liquid state and an admixture of protons and electrons in beta equilibrium. Muons appear when the electron chemical potential reaches the muon mass, replacing a fraction of the electrons in order to minimise the energy of the system. This system in equilibrium is called npeµ matter and is the minimal model for the NS core. The outer core extends to densities of ∼ 2ρ 0 . The composition of NSs at higher densities is less understood; the inner core may contain meson condensates, hyperons or deconfined quark matter [50,57,58]. The appearance of these exotic species depends on the mass of the star. In this paper, we will focus only on neutron stars made of only npeµ matter and we will consider DM scattering off leptonic targets.

Equation of state
With the sole exception of the outermost crust layers (which are only a few meters thick), NS matter is mainly in a strongly degenerate state. A consequence of this is that the pressure is independent of temperature. As a result, the equation of state (EoS) of dense matter depends only on one parameter, frequently taken to be the baryon number density, n b . The EoS is the key ingredient needed to solve the NS structure equations. Nevertheless, its precise determination is still one of the key open problems in nuclear astrophysics. The EoS governing the NS core is particularly challenging, even if we assume that only npeµ matter is present, since it requires knowledge of the behaviour of strong interactions in superdense matter.
Of the several EoSs found in the literature, see e.g. refs. [58][59][60][61][62][63], we consider the unified equations of state for cold non-accreting matter with Brussels-Montreal functionals BSk19, BSk20, BSk21, BSk22, BSk24, BSk25 and BSk26 [54,61,64,65]. A unified equation of state provides a thermodynamically consistent description of a NS from the surface to the core centre. These unified EoSs assume that a NS is made of neutrons, protons, electrons and muons, neglecting the presence of exotic matter. Analytic fits for these EoSs are given in ref. [66,67] 1 . These fits not only provide us with an excellent tool for evaluating NS microscopic properties without directly performing the nuclear physics calculations, but also are easily coupled to the Tolman-Oppenheimer-Volkoff (TOV) equations [68,69] to obtain the stellar structure.
These EoS families were obtained under beta equilibrium. Inverse beta decay equilibrium and charge neutrality dictates the exact abundances Y i and chemical potentials µ F,i of the NS constituents throughout the stellar interior, . Analytic fits for these quantities in the core and the crust as a function of the baryon number density are provided in ref. [67].

Benchmark models
In addition to QCD at high density, the NS internal structure is determined by general relativity (GR). Therefore, to obtain radial profiles of the quantities needed in our analysis, we assume a nonrotating, non-magnetized, spherically symmetric NS, and couple the EoS, P = P (n b ), ρ = ρ(n b ), to the TOV equations [68,69] 4) and the mass equation where M (r) is the mass contained within a sphere of radius r, Φ(r) is the gravitational potential, and the Schwarzschild metric is Note that the value of B(r) at the NS surface is Given an EoS, the differential equation system in Eqs. 2.3, 2.4 and 2.5 can be solved from the NS centre where ρ(0) = ρ c , with ρ c a free parameter, to the outermost layer of the crust where ρ = 10 6 g cm −3 . At that density the NS radius, R , and the gravitational mass of the star M = M (r = R ) are determined.
In Fig. 1, we show the mass radius relation for the above mentioned BSk functionals. We do not consider BSk20 and BSk21 since they yield very similar results to BSk26 and BSk24, respectively [70]. Shaded regions denote constraints on the NS maximum mass (red) [73][74][75][76][77] and radius (grey) [78] from the NS binary merger GW170817. The shaded yellow region represents the lower bound on the NS radius derived in ref. [71] from the analysis of the GW170817 event. The 2σ confidence level constraint on the radius of a 1.4 M NS is shown in black [79]. Dashed lines denote EoS families excluded by observations.
In addition, the functionals BSk19-21 were fitted to older atomic mass data than the new series of functionals BSk22, BSk24, BSk25 and BSk26. BSk19 fails to accommodate massive NSs and part of its parameter space is ruled out by the lower bound on the NS radius inferred from observations of the electromagnetic counterpart of the NS binary merger event GW170817 [71]. BSk22 is ruled out by constraints on the tidal deformability parameter also from GW170817 [70], and because it requires direct Urca processes operating in most NSs [67]. On the other hand, direct Urca processes are not allowed in stable NSs described by BSk26, which contradicts observations [67]. BSk24 and BSk25 agree with current neutron star cooling observations [72], with BSk24 giving slightly better NS mass fits to observational data [67]. Therefore, as in refs. [38,49], we select the BSk24 functional.
Coupling the BSk24 functional precision fits to the TOV equations 2.3-2.5, we solve the differential equation system from the core centre out to the outer crust, at every radial step we determine chemical potentials and particle number fractions for the different species using the appropriate functions for the core and the crust, available in ref. [67]. We also calculate GR corrections encoded in the B(r) profile. We have thus obtained radial profiles for particle number densities n i and chemical potentials µ F,i for every NS constituent. These profiles depend on the EoS choice, i.e. on the initial parameter ρ c . We have chosen the same four configurations of the functional BSk24 given in ref. [49], see Table 1, where the NS mass range is motivated by observations [80,81] and the maximum NS mass considered is limited by the GW170817 event to M 2.16M [73][74][75][76][77] 2 . In Fig. 2 we plot the corresponding lepton profiles. As mentioned above, electrons are present in the core and the crust while muons appear at baryon number densities n b 0.12 fm −3 . The kink observed in the electron chemical potential marks out the transition from the core to the inner crust. The aforementioned    Table 1. radial profiles will be used in the following section to calculate the capture rate.

Capture and Interaction Rates
In ref. [49], we derived general expressions for the capture and interaction rates of DM in NSs, valid for a broad range of DM masses, for arbitrary NS targets and DM-target cross sections. These expressions properly take into account relativistic kinematics, including the motion of the target particles in the NS; gravitational focusing; Pauli blocking (relevant at low DM masses); the effect of the star opacity and multiple scattering (important in the capture of heavy DM); and the NS internal structure. In that which follows, we summarise and extend those results.

Interaction Rate
The DM scattering rate as a function of arbitrary DM energy, and the corresponding differential interaction rate, are the key elements of the capture calculations. In addition, they are a necessary input in constructing the probability density function of the DM energy loss. This is required to define the capture probability after N scatterings [49], which is relevant to the capture of heavy DM via multiple scattering.
Following refs. [29] and [49], we define the DM scattering rate as where k µ = (E χ , k), k µ = (E χ , k ) are the DM initial and final momenta, p µ = (E , p) and p µ = (E , p ) are the target particle initial and final momenta, m is the lepton target mass, q 0 = E −E is the DM energy loss and f FD is the Fermi Dirac distribution. We define the quantity S(q 0 , q) to be the response function, which contains the dependence on the squared matrix element, |M | 2 . To calculate the interaction rate, we consider the interaction of Dirac DM with SM leptons, described by the dimension 6 effective operators listed in Table 2, where the strength of the coupling is parametrised by the cutoff scale Λ and µ = m χ /m . In ref. [49], we showed that Eq. 3.1 can be solved analytically for DM-nucleon differential cross sections that depend on powers of the Mandelstam variable t (but not on s) as a function of the DM energy E χ . In that case, assuming that T → 0, Eq. 3.1 reduces to where t E = −t = q 2 − q 2 0 and µ F, is the target chemical potential after subtracting the target rest mass energy as in Fig. 2. From here on and as in ref. [49], we denote by µ F, the Fermi energy (without the rest mass) of the leptonic species . The quantity Table 2: EFT operators [82] and squared matrix elements for the scattering of Dirac DM from leptons. The coefficient of each operator is given as a function of the lepton Yukawa coupling, y , and the cutoff scale, Λ. The fourth column shows the squared matrix elements at high energy as a function of the Mandelstam variables s and t.
is the minimum energy of the target before the collision, which is obtained from kinematics, and g 0 (x) is a step function with a smooth transition, (3.5) The explicit integrals over t E are given in appendix B of ref. [49]. The differential interaction rate dΓ dq 0 (E χ , q 0 ) is the integrand of Eq. 3.3. These expressions can be applied to squared matrix elements that depend on linear combinations of t n , with n = 0, 1, 2. As seen in Table 2, this is applicable to the D1-D4 operators. For the remaining operators D5-D10, we require either a numerical computation or an analytical approach that generalises that of ref. [49] to now handle s-dependent interaction rates. We derive such analytical expressions for s-dependent interaction rates for the first time, with our results presented in appendix A. It is worth noting that these expressions are valid only in the zero temperature approximation. With these results, the interaction rates for operators D5-D10 can be obtained as linear combinations of those for simple power laws |M | 2 ∝ t n s m . There are 6 possible power laws in total, namely 1, t, t 2 , s, st, s 2 . The methodology to calculate the full expressions for |M | 2 ∝ t n s m is similar to that adopted for s-independent matrix elements in ref. [49], with a few additions that are outlined in appendix A. We do not report the full expressions for |M | 2 ∝ t n s m due to their length. Table 3: Typical values of m * and σ th χ for lepton targets. The exact value of σ th χ depends on the DM mass, and the operator. We show here the simplest case of constant matrix element; other operators give similar results. The threshold cross section is approximately constant in the range 1 GeV m χ m * , and takes larger values outside that range with a 1/m χ or m χ scaling for small and large masses, respectively.

Capture Rate
Below we provide a summary of the various expressions for the capture rate and the regimes for which they apply, the details of which can be found in ref. [49].
1. Optically thin, single scatter: σ σ th χ and m χ m * We begin by defining c 1 to be the probability that a single scattering interaction will result in capture of the DM particle. The simplest regime occurs when the cross section is much smaller than the threshold cross section, σ σ th χ , and the DM mass is much smaller than m * , where m * is the DM mass for which multiple scattering becomes relevant. Both quantities σ th χ and m * depend on the specific target i. To calculate m * for operators D1-D10, we use the differential interaction rates dΓ dq 0 computed in section 3.1 and follow the approach outlined in ref. [49]. In Table 3, we show typical values of m * and σ th χ for electron and muon targets. Note that the exact value of m * depends on B(r), µ F, (r) and the type of interaction.
When the above mentioned conditions are met, the capture probability is of order one, c 1 ∼ 1, and the neutron star can be treated as optically thin. In this limit, the capture rate is given by where ρ χ is the local DM density, v is the NS velocity, v d is the DM velocity dispersion, ζ(r) = n (r) n f ree (r) is a correction factor due to the use of a realistic target number density n (r), E and E are the initial and final energy respectively of the target , and The integration intervals for s, t and E i are given in ref. [49]. Note that Eq. 3.7 correctly accounts for Pauli blocking, given by the 1 − f FD term, which, for muons, is relevant for m χ m µ . Electrons, on the other hand, are ultra-relativistic throughout the inner crust and the core, with Pauli suppression effective for m χ µ F,e .

2.
Optically thin, large mass -multiple scattering: σ σ th χ and m χ m * , For m χ m * and σ σ th χ , the assumption that the capture probability is c 1 ∼ 1 no longer holds. In fact, it is significantly smaller than 1. We calculate the capture rate using the following approximation: , (3.10) where the capture probability is given by where n * represents the average number of interactions with target species required to remove a DM particle from the incoming flux. In this way, a suitable approximation that accounts for multiple scattering is obtained. Typical values are reported in Table 3 for lepton targets.
3. Optical depth: σ ∼ σ th χ . If σ ∼ σ th χ , the optically thin limit is not valid and hence we must modify the capture rate expressions above (Eqs. 3.6 and 3.10) to include an optical factor η(r) [49]. This is an extinction factor that accounts for the star opacity. We then have 4. Geometric limit: σ σ th χ . In this case we can safely estimate the capture rate using the geometric limit calculated in ref. [36],

Capture Rate
In this section, we present results for the capture rate CΛ 4 for each of the EFT operators in Table 2, calculated in the optically thin limit using Eq. 3.6 for m χ m * and Eq. 3.10 for m χ m * 3 . Figs. 3 and 4 show the results for electron (light blue) and muon (magenta) targets, considering three NS benchmark models: BSk24-1 (dashed, 1M ), BSk24-2 (solid, 1.5M ) and BSk24-4 (dot-dashed, 2.16M ). In addition, we assume a nearby NS, located in the Solar neighbourhood, and thus take In these figures, we observe that the capture rate is suppressed due to Pauli blocking when m χ m µ . The change of slope at m χ ∼ m * ∼ 10 5 GeV, observed for both targets, is due to multiple scattering. Note that the slope of the capture rate for the three distinctive regions, low mass, intermediate (m µ m χ m * ) and large mass (multiple scattering) is very similar for operators m χ (GeV) We require Λ to be sufficiently large that the capture rates are smaller than the geometric limit, C geom .
D5-D10 (Fig. 4), while for D1-D4 (Fig. 3), the shape of C is controlled by the power of t that dominates the interaction, which in general is the lowest power [49]. The sole exception to this are the capture rates for operators D1 and D2 with electron targets, which show a distinctive feature in the region m e m χ 100 MeV that does not occur for the other operators. The C rate for D1 and D2 is more suppressed in that particular region, similarly to D3 and D4, respectively. This is due to the form of the corresponding matrix elements together with the smallness of the electron mass. Namely, D1 and D2 are the only two operators that contain a factor (t − 4m 2 ) in their scattering amplitudes, for electrons this means that the lowest power of t in |M | 2 is multiplied by m 2 e , i.e. these terms are suppressed in the m e m χ 100 MeV interval. Consequently, the capture rate in that DM mass region is dominated by the unsuppressed t-terms in |M | 2 , t for D1 (as for D3) and t 2 for D2 (see Table 2), while below m e this additional suppression disappears and the capture rate follows the lowest power of t as for muon targets.
From Fig. 3, we note that for the same cutoff scale Λ, the muon contribution to the total capture rate for operators D1-D4 surpasses that of the electron by approximately 4 orders of magnitude for most of the DM mass range, and by about 8 orders of magnitude at very low mass for operators D1-D2 (because of the additional suppression described above). This is due to the large hierarchy between DM couplings to electrons and muons, which is of order ( mµ me ) 2 . Conversely, for operators D5-D10, electrons and muons have the same couplings (see Table 2). However, despite similar couplings and a lower abundance, muons are able to capture DM at a rate comparable to electrons (see light blue regions in Fig. 4), thanks to their larger mass and lower chemical potential (see Fig. 2, right panels), i.e., their interactions with DM are less Pauli suppressed. The small difference between the rates at which electron and muon are able to capture DM particles reduces for heavier 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 m χ (GeV) BSk26-1 Geometric Opt. thin Figure 5: Capture rate in the optically thin limit for muon targets (magenta) and geometric (orange) limit as a function of the DM mass for constant cross section σ µχ = 10 −45 cm 2 , ρ χ = 1 GeV cm −3 and BSk26 functional for M 1.52M and R 11.6 km denoted as BSk26-1. Capture rate calculations from ref. [30] for a NS configuration with EoS BSk20-1 [66] equivalent to BSk26-1, are shown for comparison. NS configurations, e.g. from a factor ∼ 5 (BSk24-1) to ∼ 1.5 (BSk24-4) for D6 and D10; see the light blue shaded regions in Fig. 4. Recall that muons are expected to be found in larger fractions in massive NSs (see Fig. 2, left panels).
It is also worth noting that different EoS assumptions can lead to variations in the capture rate for electron targets of at least two orders of magnitude in the Pauli suppressed region and ∼ 2. The DM capture rate for muon targets was calculated in ref. [30], for constant cross section and light DM, m χ ≤ 10 GeV. That calculation accounts for the NS internal structure and Pauli blocking, but neglects general relativity (GR) corrections and assumes that muons are non-relativistic. In order to compare our capture rate calculation with that of ref. [30], as in ref. [49], we have selected a NS model that matches that of Fig. 12 of ref. [30], namely Model A (BSk20-1): M 1.52M , R 11.6 km. This new benchmark model is denoted as BSk26-1. Note that there are no public fits for chemical potentials and particle abundances for BSk20; however, as mentioned in section 2.1, BSk26 yields NS configurations that are almost indistinguishable from those obtained with BSk20 [70].
In Fig. 5, we compare both capture rate calculations for σ µχ = 10 −45 cm 2 and the same assumptions about ρ χ , v and v d as in ref. [30]. Comparing the geometric limit, Eq. 3.13 (solid orange), which properly accounts for gravitational focusing in NSs, with the non-relativistic computation in ref. [30] (dot-dashed brown), we observe a ∼ 67% enhancement, due to the 1/B(R ) factor that encodes GR corrections [22,86]. In the region not affected by Pauli blocking, m χ m µ , our calcu- lation in the optical thin limit (solid magenta) exceeds that of ref. [30] (dot-dashed blue) by a factor of ∼ 4, which increases as we move to the Pauli suppressed region where our computation is more than one order of magnitude higher. Unlike ref. [30], our formalism incorporates GR corrections and made use of relativistic kinematics (recall that muons in NSs are mildly relativistic). We also show in dashed light blue, an estimation of the capture rate using the approximation δp/p F,µ ∼ m χ v esc /p F,µ for m χ < m µ [25], where p F,µ is the muon Fermi momentum and v esc is the escape velocity. This approximation overestimates the capture rate by a factor of approximately 2 in the Pauli blocked region below 10 MeV and underestimates it in the region of larger DM masses.

Finite Temperature Effects and Evaporation
In section 4.1, we have restricted our computation of the capture rates to the DM mass range m χ ∈ [1 keV, 10 8 GeV]. It is worth noting that this calculation can also be performed for smaller or larger DM masses. Note, however, that the analytic expressions for the DM interaction rate in ref. [49] and appendix A were derived in the zero temperature approximation. Therefore, they can be used safely only for m χ T , where T is the NS temperature. For m χ O(10)T , thermal effects play an important role and increase the capture rate of very light DM [30]. Consequently, the complete Fermi Dirac distribution should be used in Eqs. 3.2 and 3.7. To illustrate the effect of the NS temperature, we show in Fig. 6 the ratio of the capture rate in a NS with T = 10 5 K 8.6 eV to the corresponding C rate in the T → 0 limit, assuming scattering on electrons, the targets for which this effect is most relevant. From this figure, we immediately notice that the ratio starts to depart from 1 at m χ ∼ 100 eV ∼ 10T for all operators. Operators whose matrix element depends on higher powers of the exchanged momentum t feature a larger increment in the capture rate due to finite temperature. In fact, the operator D4 (|M | 2 ∝ t 2 ) receives the largest correction, followed by D2-D3 (whose |M | 2 is a linear combination of t 1 , t 2 ), D1 (|M | 2 is a linear combination of t 0 , t 1 , t 2 ) and finally by D5-D10 (whose |M | 2 include all powers of the kind t n s m ).
In the very light DM regime, there is another process we should be aware of: evaporation. This occurs when the dark matter up-scatters to a state where the final DM kinetic energy is larger than the energy required to escape the star, and hence DM particles are expelled. Thus, opposite to capture, evaporation drains energy from the star. To estimate the evaporation rate, we convolve the DM distribution within the star, with the interaction rate for up-scattering, Γ − + , retaining the temperature dependence. Assuming the DM distribution to be isothermal with temperature T χ = T , we have where n c is the DM number density at the centre of the star, while the interaction rate for up- where dΓ − dq 0 is the differential interaction rate in the T → 0 approximation derived in ref. [49] and appendix A (for details of the derivation of Eq. 4.2, see appendix B). The evaporation rate then reads When the DM distribution is concentrated very close to the centre of the star, this expression can be approximated by The rate at which DM particles accumulate in NSs is then given by assuming that DM annihilation is negligible. The solution of this equation is where t is the age of the NS. The term in brackets quantifies the negative contribution of the evaporation process to the total number of accumulated DM particles. Note that this factor is of order 1 except when E(m χ ) t O(1). Therefore, we define the evaporation mass as the DM mass for which the previous relation holds, i.e. E(m evap )t ∼ 1. For DM masses below this threshold, m χ m evap , the capture and evaporation processes are in equilibrium with each other. In that limit, the net energy exchange in the star due to combined effects of DM capture and evaporation would be negligible, and hence we would be unable to constrain DM interactions using the NS temperature as a probe.
Using Eq. 4.3, we find the evaporation mass to be of order m evap ∼ O(100T ) for all scattering targets in old NSs with t ∼ O(10 Gyr). For instance, for T = 10 5 K and electron targets, we obtain m χ (MeV) The solid blue (electron) and magenta (muon) lines represent σ th , computed assuming the NS model BSk24-2, while the shaded bands represent the expected range due to variation of the EoS. For comparison we show leading electron recoil bounds for heavy mediators from SENSEI [87], DAMIC [88], Xenon10 [89], Xenon1T [90], projected sensitivities from DAMIC-M [91] as well as the neutrino floor for silicon detectors [92].
m evap 300 eV. From Fig. 6, we note that the evaporation mass (dashed brown line) is larger than the mass at which temperature effects on the capture rate are important. In other words, evaporation comes into play before finite temperature effects become relevant so that these effects can be safely neglected when calculating the capture rate with the aim to constrain DM interactions.

Threshold Cross Section
In ref. [49], we defined the threshold cross section, σ th χ , as the cross section for which the capture rate C(σ(Λ), m χ ), calculated in the optically thin regime (i.e. without the optical factor η), is ∼ C geom . This definition is general and is applicable to both relativistic and non-relativistic targets. The threshold cross section restricts the NS sensitivity to DM-target interactions, since for σ ≥ σ th χ the capture rate saturates to the geometric limit C geom .
In Fig. 7, we show the threshold cross sections for lepton targets, electrons and muons, and compare them with existing limits and expected sensitivities of future experiments. The neutrino floor for electron recoil experiments for silicon targets [92] is shown as a shaded yellow region. The solid light blue and magenta lines correspond to the value of σ th for electrons and muons respectively, calculated using the NS model BSk24-2 (1.5M ), while the shaded bands in light blue and magenta denote the expected range for σ th for the two different targets, obtained by varying the NS configuration along the BSk24 family. BSk24-1 (1M ) gives the upper bound on σ th and BSk24-4 (2.16M ) the lower bound. Note that the variation in σ th due to the NS EoS increases with the DM mass and for muons goes from about one order of magnitude in the low mas range to two orders of magnitude in the multiple scattering region. For electrons, this effect is slightly less pronounced. All the limits for existing experiments are orders of magnitude weaker than the expected NS reach, with only the future DAMIC-M [91] (dashed brown line) expected to overcome NS electron scattering sensitivity and approach that of muons, in the DM mass range 3 MeV m χ 30 MeV. Moreover, NS sensitivity to DM interactions with lepton targets is expected to be well below the neutrino m χ (GeV)  Figure 8: Comparison of the reach in Λ for D1 with the approach of ref. [42]. In these references, a NS with constant chemical potentials and particles abundances averaged over the core volume was assumed, these quantities were taken from ref. [38] and are consistent with a NS with EoS BSk24-2. The shaded regions denote the difference in Λ (top) and CΛ 4 (middle) between the two approaches and their ratio is shown in the bottom panels.

Ratio
floor for m χ 100 MeV and, in the case of muons, even for m χ 1 MeV. Note that NSs have a better sensitivity to vector-vector interactions (operator D5, see right panel) than scalar-scalar interactions (operator D1, see left panel) in the low DM regime for both leptonic targets, especially for electrons, since as mentioned in section 4.1 there is an additional suppression in the capture rate of scalar operators that stems from a m 2 e t term in their scattering amplitudes. Similar threshold cross sections can be estimated for the remaining operators. Operators with s-dependent matrix elements (D6-D10) have σ th that behaves like that of D5 for both electrons and muons. D2 presents the same features as D1 in the sub-GeV regime for electrons, due to the similar shape of their capture rates (see Fig. 3) and D3-D4 show a steeper slope in the m χ m e region with respect to D1-D2, due to the capture rate dependence on higher powers of t (see Table 2 and Fig. 3).  Figure 9: Comparison of the reach in Λ for D5 with the approach of refs. [41,42]. The shaded regions denote the difference in Λ (top) and CΛ 4 (middle) between the two approaches and their ratio is shown in the bottom panels.

Ratio
In Figs. 8 and 9, we compare our results for D1 and D5 with those of refs. [41,42] 4 . The formalism in refs. [41,42] is valid for relativistic and non-relativistic targets in a broad mass range, but neglects the DM velocity distribution and the NS internal structure. Instead, constant chemical potentials and particle abundances, averaged over the core volume, are assumed. These quantities correspond to the NS model BSk24-2 and were calculated in ref. [38]. In the top panels, we compare the reach in Λ for DM-lepton scattering cross sections in refs. [41,42] with the cutoff scale we obtain for the maximum capture rate C(Λ, m χ ) = C geom . Our results differ the most for electron targets in the Pauli suppressed region by a factor of ∼ 2.5 and we find Pauli blocking is active at a slightly lighter DM mass. Recall that we have obtained the radial profiles for chemical potentials and number densities with a unified EoS, i.e. for the core and the crust, and note that most light DM particles whose interactions are subject to Pauli blocking are captured close to the surface [49]. The difference between our approach and that of refs. [41,42] is reduced to a factor of ∼ 1.25 in the intermediate mass region and there is almost no difference in the large mass regime, except for the DM mass at which multiple scattering becomes relevant, which in our case is once again slightly lighter. For muons, we find a Λ that is, on average, a factor ∼ 1.33 (D5) greater than that of refs. [41,42] along the whole DM mass range and is in almost perfect agreement in the m χ m µ region for D1.
In the middle panels of Figs. 8 and 9, for operators D1 and D5 respectively, we show how these apparently small differences in the two approaches translate to differences in the capture rate. To that end, we compare CΛ 4 = Λ 4 C geom obtained with the two formalisms. Since the geometric limit of the capture rate is not defined in refs. [41,42], we use a definition similar to Eq. 3.13 and compliant with assumptions made by these authors. For electron scattering, we see that the formalism that does not account for the NS internal structure underestimates the capture rate in the region affected by Pauli blocking by a factor ∼ 40 (bottom LH panels of Figs. 8 and 9). This difference is slightly larger in the region where the Pauli suppression is stronger, becoming almost a factor of ∼ 100 for D1 in the range m e m χ 100 MeV. For muons, the difference between the two approaches is less pronounced, with a maximum ratio of ∼ 3.5 for both operators.

Conclusions
Neutron stars (NSs) are potential cosmic laboratories to study dark matter (DM) interactions with ordinary matter under extreme conditions. Gravitational focusing enhances the rate at which DM particles can accumulate in these stars. Thus, NSs emerge as potential DM probes, complementary to direct detection experiments which are restricted by recoil thresholds and small momentum transfers. DM scattering off NS targets is, however, not free of limitations; in the sub-GeV regime DM scattering off strongly degenerate targets is suppressed by Pauli blocking and in the large mass region multiple collisions are required to capture heavy DM. Moreover, there is a natural threshold for the maximum cross section that can be probed in NSs, above which the capture rate saturates to its geometric limit.
NSs are systems in beta equilibrium such that, even though they are primarily composed of degenerate neutrons, protons and electrons are present throughout the star with abundances of order a few percent. Muons are also found in NS cores at higher densities. Unlike nucleons, the leptonic species in the NS are relativistic. In this paper, we have examined the reach of NSs to probe the interactions of fermionic DM with the leptonic NS constituents, in the context of an effective field theory (EFT). To that end, we generalised our formalism to calculate the DM interaction rate, presented in ref. [49], to enable it to handle any differential cross section parametrized in terms of the Mandelstam variables s and t. With this extended formalism, we calculated the capture rate for the full list of dimension 6 effective operators for a broad DM mass range, properly including Pauli blocking, multiple scattering, NS internal structure and general relativity (GR) corrections.
To be consistent, the aforementioned calculation requires knowledge of the microscopic properties of the target species, such as chemical potential, number density and abundance, as well as GR corrections. These quantities have a radial dependence and hence require the assumption of a NS equation of state (EoS). This is particularly relevant for leptons, as their particle fractions are heavily dependent on the NS mass. To account for that uncertainty, we have assumed the unified EoS with Brussels-Montreal functional BSk24, which is well motivated by observations. We find that scattering off muons dominates the leptonic contribution to the capture rate for scalar and pseudoscalar DM-lepton interactions, despite the muon abundance being lower than that of electrons. This is due to the fact that the couplings for these interactions scale with the lepton mass. For other interaction types, electrons and muons have the same coupling scale. In spite of that, and the fact that electrons are ultra-relativistic while muons are only mildly relativistic, the capture rates for electron and muon targets are comparable, a consequence of the large muon mass and lower muon chemical potential. This effect is enhanced in heavy NSs where muon and electron number densities are very similar.
The NS sensitivity to DM interactions with leptons greatly surpasses that of any current direct detection (DD) experiment, particularly in the sub-GeV DM regime. Only future DD experiments such as DAMIC-M could be competitive and, even then, only in a narrow mass range 3 MeV m χ 30 MeV. Finally, note that the evaporation mass, the DM mass below which capture and evaporation processes are expected to be in equilibrium, is much smaller for NSs than for other stars or planets. This provides sensitivity to sub-MeV DM with scattering cross sections even below the DD neutrino floor. These findings are particularly relevant for leptophilic DM for which couplings to nucleons arise at loop level.

Acknowledgements
NFB and SR were supported by the Australian Research Council and MV by the Commonwealth of Australia. We thank Filippo Anzuini and Tony Thomas for helpful discussions.

A Interaction rate for s-dependent amplitudes
In ref. [49], we obtained analytic expressions for the DM interaction rate, for squared matrix elements that depend on t, but not on the centre of mass energy s. In the following, we generalise our previous result to |M | 2 that can be written as a polynomial function of the variables s and t, i.e. |M | 2 = αs m t n , where n are m are integers and α is a constant.
Following refs. [29] and [49], we define the DM scattering rate as where k µ = (E χ , k), k µ = (E χ , k ) are the DM initial and final momenta, p µ = (E i , p) and p µ = (E i , p ) are the target particle initial and final momenta, and q 0 = E i − E i is the DM energy loss.
Note that now |M | 2 depends on p, so we leave it inside the response function S(q 0 , q). Integrating the response function over d 3 p using the delta function leaves After that the final target energy is fixed to where θ is the angle between p and q. To perform the integral over d 3 p we change it to d 3 p = pE i dE i d cos θdφ and use the delta function to integrate over θ [29,93]. Note that this gives rise to Θ(1 − cos 2 θ) [93]. Then, we obtain we can determine the integration interval for E i . For the q 2 > q 2 0 case, q µ is expected to be space-like, t = q µ q µ < 0, and the response function becomes where E t − i is the minimum energy of the neutron before the collision, which is obtained from kinematics and given by To perform the integration over the azimuth angle φ we rewrite s in terms of the other kinematic variables, q, q 0 , where the value of the scalar product of the two momenta is We are mostly interested in values of m = 1, 2. For instance, for m = 1 From now on, we do not give explicit expressions for the integrals, due to their length, but just sketch the procedure to easily obtain the solutions using any symbolic calculation language. The previous expression A. 11 can be written as where U m (q 2 , q 0 , E χ , E i ) is a polynomial of degree m in E i that can be rewritten as where V m,i (q 2 , q 0 , E χ ) are also polynomials. We therefore need to calculate integrals of the form Using we solve the integrals, e.g. for m = 1 where F 1 (x, z) is obtained by integrating by parts and P L is the PolyLog function. Similarly, we can solve the m = 2 integral over E i for For F 0 , three distinct regimes can be identified as noted in ref. [49] (s-independent case) Next, we use the following results for each of the above E i intervals, which are valid in the T → 0 limit The expressions above, which resemble a step function with a smooth transition, can be recast in terms of the function Then, the contribution of F 0 to the response function in the T → 0 limit is Note that Eq. A.29 is the result we found in ref. [49]. We proceed in a similar way for F 1 , F 2 and obtain lim T →0 (A.33) Then, the general expression for the response function is (A. 35) Comparing to the case where the amplitude does not depend on s, there are two additional transition functions.
We now return to the scattering rate Γ − = α d cos θk 2 dk 64π 2 E χ E χ m 2 i t n Θ(E χ − q 0 − m χ )Θ(q 0 )S − (q 0 , q), (A. 36) change variables from k , cos θ to q 0 , q, q 0 = E χ − k 2 + m 2 χ , (A.37) q 2 = k 2 + k 2 − 2kk cos θ, (A. 38) and substitute in the result for S − , Eq. A.34, to obtain, (A.39) To simplify the integration over t, we define t E = −t = q 2 − q 2 0 , which leads to (A.40) As V m,i (t E + q 2 0 , q 0 , E χ ) are polynomials, we need to calculate integrals of the form which can be solved by decomposing the integration interval using the primitives of Using the operator defined in ref. [49] that encodes the t E integral over the correct intervals, we obtain Γ − as a linear combination of terms of the kind for further details on the calculation of t ± E and t ± µ ± see appendix B of ref. [49].

B Interaction rate for up-scattering
We perform a similar calculation to that in the previous section, assuming that q 0 is now negative. The derivation of the up-scattering interaction rate is essentially the same as for down scattering until we arrive at the point of the identification of the three different regimes in the response function. For q 0 < 0 these read We consider finite values of T , and take the leading contribution, i.e. the terms of order e −|q 0 |/T . For matrix elements independent of s we have and the response function is now given by which is equal to Eq. B.23 of ref. [49] times a factor that depends on the temperature and the momentum transfer q 0 (and implicitly on µ F,i and B), S − (q 0 , q, T ) = − e q 0 /T 1 − e q 0 /T S − (q 0 , q).
(B.9) Therefore, the differential interaction rate for up-scattering can be estimated with the following expression dΓ − + dq 0 (q 0 , T ) = − e q 0 /T 1 − e q 0 /T dΓ − dq 0 (q 0 ) , q 0 < 0. (B.10) One can then integrate over the variable t E , as in the previous section, to obtain the interaction rate.