Terrestrial and Martian Heat Flow Limits on Dark Matter

If dark matter is efficiently captured by a planet, energy released in its annihilation can exceed that planet's total heat output. Building on prior work, we treat Earth's composition and dark matter capture in detail and present improved limits on dark matter-nucleon scattering cross sections for dark matter masses ranging from 0.1 to $10^{10}$ GeV. We also extend Earth limits by applying the same treatment to Mars. The scope of dark matter models considered is expanded to include spin-dependent nuclear interactions including isospin-independent, proton only, and neutron only interactions. We find that Earth and Mars heating bounds are alleviated for dark matter s-wave self-annihilation cross sections $\lesssim 10^{-44}~{\rm cm^2}$.


Introduction
Despite a preponderance of evidence for the existence of dark matter (DM), its mass and the nature of its non-gravitational interactions remain unknown. It may be the case that DM couples to Standard Model fields through a non-gravitational force. Here we will study DM's interactions with itself and with nucleons in the Earth and Mars, and by precisely modeling each, derive extended planetary heating bounds, first obtained in Reference [1].
Many underground experiments have searched for DM that scatters with nucleons through either predominantly spin-dependent or spin-independent interactions, including DEAP [2], PICO [3,4], LUX [5], PandaX [6], and XENON1T [7]. These experiments place detectors deep underground to reduce backgrounds and gain in sensitivity to weakly interacting DM. However, such experiments situated kilometers underground are less sensitive to more strongy-interacting DM, which has reduced kinetic energy after repeated scattering against the Earth's crust during its voyage underground. On the other hand, near-surface direct detection experiments [8] are more sensitive to strongly-interacting DM, along with repurposed high-altitude detectors like the XQC rocket [9]. Constraints on strongly-interacting DM from these experiments are supplemented by numerous astrophysical bounds, including analyses of the cosmic microwave background [10,11], interstellar gas cooling [12,13], neutrinos from DM annihilation [14][15][16][17][18], and a conspicuous lack of tracks in ancient mica [19][20][21]. This article will reinforce and extend existing constraints on DM's spin-independent interactions. Additionally, there are a number of simple models for which DM's non-relativistic interactions with nuclei depend sensitively on the spin of the nucleus [22,23]. For spin-dependent models, this work places bounds on some unexcluded regions of DM parameter space; to achieve this we incorporate the terrestrial distribution of nuclei with a substantial nuclear spin.
Strongly-interacting DM can appreciably raise the temperature of the Earth and Mars through capture and subsequent annihilation. Simulations indicate that galaxies like the Milky Way exist within a virialized and spherical DM halo [24,25]. From this it follows that so long as DM is a light particle (in this case lighter than a small asteroid), there will be a constant flux of DM incident upon the Earth as it orbits the Galaxy. If DM interacts sufficiently strongly with nucleons, it will scatter off elements within these planets, and if it scatters enough times, its velocity will decrease below the given planet's escape velocity, which is approximately 11 km/s for Earth, and 5 km/s for Mars. If DM is sufficiently slowed, it will stay bound to the planet and can annihilate with other similarly captured DM. Such annihilations may result in a sizable heat output. Because we know from direct measurements that no more than ∼44 TW of energy is emitted by the Earth [26][27][28][29][30][31][32][33], we are able to restrict DM candidates by requiring that they not yield more heat flow from Earth than is observed. While the origin of heat emitted by Earth is still under investigation, it is expected that around half of the observed heat emission is radiogenic [33]. Because in-situ heat measurements are not available for Mars and other planets, it is not yet possible to set a constraint using a planet other than the Earth. However, the InSight mission to Mars may soon report the first direct Martian heat flow measurements [34]. Therefore, we have also studied what constraint could be placed by observing no more than the expected ∼3.5 TW of radiogenic energy generated within Mars [35].
The most stringent limit on DM's nucleon interactions is obtained from planetary heating when the DM annihilation rate equals the DM capture rate. This occurs for DM with a sufficiently large self-annihilation cross section. This scenario is referred to as "total annihilation" [1]. Of course, a smaller DM annihilation cross section will result in less annihilation. This "partial annihilation" scenario is addressed in this article and we find that the heating bound effectively vanishes for σχ χ 10 −44 cm 2 , 10 −30 cm 2 , and 10 −38 cm 2 , for s-wave, p-wave, and impeded dark matter annihilation respectively. Intriguingly, these cross sections are around the canonical "weak" scale annihilation cross section σ weak χχ ≈ 10 −36 cm 2 [36]. Therefore, for DM with an approximately weak scale annihilation cross section, planetary heating bounds can have a non-trivial dependence on DM's self-annihilation cross section.
This work primarily considers three DM parameters: (1) the DM mass m χ , (2) the DM pernucleon scattering cross section σ χN , which determines the frequency with which the DM particle will interact with terrestrial and Martian elements, and (3) DM's self-annihilation cross section σχ χ , which determines the extent to which captured DM will result in increased heat flow out of the Earth and Mars. In Section 2 we treat the planetary capture of strongly-interacting DM, incorporating a three zone nuclear abundance model of the Earth's core, mantle, and crust, and a two zone nuclear abundance model of Mars' core and mantle. Section 4 details limits on DM from planetary heating, including both total and partial DM self-annihilation in the Earth. In Section 5 we conclude. Throughout this work, we have used natural units with = k B = c = 1.

Dark Matter Capture
As DM passes through a planet, it may scatter off its constituent particles, and will slow slightly with each scatter. For certain DM masses and cross sections, this scattering results in DM slowing below the planet's escape velocity, at which point it is gravitationally bound to the planet.
To simplify our DM capture computations, we ignore Earth and Mars' gravitational effects on the DM's trajectory and velocity. This is a reasonable approximation, given that the speed of a typical DM particle is ∼300 km/s, compared to the ∼10 km/s escape speeds in question. We also ignore the directional changes in the DM's trajectory induced by scattering. These are both conservative approximations, as both a random walk from scattering and a gravitationally curved trajectory lead the DM through more material than a straight trajectory.
For convenience, we define where m χ is the mass of the DM particle and m j is the mass of the terrestrial constituent j off of which it is scattering. The DM particle's kinetic energy after a single scatter, E f , can be defined as a function of its initial kinetic energy, E i , as where z ⊂ [0, 1] is a kinematic factor parameterizing the scattering angle of the DM-nuclear interaction [37]. On average z ≈ 1 2 , and we set it to this value in our computations. Making the substitution that kinetic energy is proportional to velocity (v) squared, and now considering τ j scatters off of each element j, the expression for DM's final velocity is When v f < 11.2 km s −1 or 5.0 km s −1 for Earth or Mars respectively, the DM is gravitationally captured, as it is no longer traveling fast enough to escape the planet's gravity. However, see later sections for discussion of DM evaporation.
To find an expression for τ j , we must first find the average number of scatters off of each element. When travelling a distance L through a medium of constant density, we define this as τ j = n j σ χj L [37], where σ χj is the nuclear cross section of the DM with element j. In a medium with non-constant elemental density like the Earth and Mars, we instead must define τ j as a function of θ, the angle between the DM particle's trajectory and the vector normal to the planet's surface. To find the mean number of scatters, we integrate along the DM's path (l) from its entry point to a distance of 2R e cos(θ), at which point it will exit the planet. A schematic diagram of the trajectory through the Earth is given in Figure 1. Explicitly, the expectation value for the number of scatters is given by the integral Because number density, n j , is a function of radius from the planet's center, the substitution r = l 2 + R 2 e − 2lR e cos θ must be made. Scattering events are independent and discrete, so τ j (θ) actually follows a Poisson distribution with a mean of τ j . In computations we required the number of scatters off of element j to follow the probability distribution It is also important to distinguish between DM-nucleus cross section (σ χj ) and the DM-nucleon cross section (σ χN ). We follow the conventions of [1,37,38]. For a spin-independent cross section we used the conversion σ (SI) where m N is the mass of a nucleon (∼1 GeV), µ(m j or m N ) is the reduced mass of the DM particle and atom j or a single nucleon respectively, and A j is the number of nucleons in atom j. In the spin-dependent case, Here, we define J j as nuclear spin of atom j, S p j and S n j as its average proton and neutron spins, and a p and a n as proton and neutron coupling constants. In this work, we consider three cases: (1) isospin-independent scattering (a p = a n = 1), (2) proton-only scattering (a p = 1, a n = 0), and (3) neutron only scattering (a p = 0, a n = 1).
Finally, we note that in a number of recent publications [20,37,39], it has been pointed out that the spin-independent DM-nucleon scattering cross section as presented in Eq. (6), reaches a theoretical transition point at σ χN ≈ 10 −26 cm 2 . At larger cross sections, the implied DMnucleus cross section is larger than the physical area of nuclei, implying either long-range forces (e.g. mediation by light dark photons [13,39]) or composite dark matter [40].

Elemental Makeup
The Earth's material composition is divided into three parts: the crust, the mantle, and the core. Taking the Earth's center as r = 0, the crust begins at the Earth's surface at r = R e = 6371 km and ends at r = 6346 km [41]. From that radius until r = 3480 km is the mantle, [42] and the remainder of the Earth is the core in our model [43]. The relative sizes of these regions are shown in Figure 1, and the material composition of these regions is given below in Table 1.
Also shown in Table 1 are the compositions of Mars' mantle and core. We assume Mars' crust to have a thickness of 50 km [44]. Its composition is taken as the same as that of the mantle for calculating DM drift times in section 4.4, and is conservatively omitted for capture and annihilation computations in sections 4.1 and 4.3. We give the core and mantle thicknesses of 2000 and 1340 km respectively, for a total Martian radius of R m = 3390 km. All radial thickness values given here have been chosen among values presented in the above references, to minimize DM capture on Earth and Mars.
The Earth's density also varies as a function of distance from the center. The preliminary reference Earth model [47] is a reasonable approximation of Earth's mass density. To be conservative, we used the minimum possible Martian density at all radii from the models given in [44,46]. We henceforth reference these densities as ρ(r), which is plotted in Figure 2. To convert ρ(r) to n j (r), the number density of element j, we divide by the mass of that element m j and multiply by the mass fraction of that element for either the core, mantle, or crust, as given in Table 1.
It is also necessary to determine the planets' temperature profiles. There have been multiple Earth temperature models proposed; there is not strong consensus on the matter [28,48,49]. To be conservative, we take the highest reasonable proposed temperature at all radii to construct a "maximum temperature profile" of the Earth. This is shown in Figure 2. The atmospheric temperatures used for Earth can be found in [50]. We use the highest possible temperature profile in [46] as our model for Mars, and used the Martian atmospheric temperatures given in [51]. Many geological studies and models have been performed that estimate the total heat flux from within the Earth, all of which find the total heat flux to be around 44 TW. Some of this flux has been attributed to known processes, such as emission from radiogenic sources like uranium and thorium present in the Earth [28,52]. However, to be conservative we will attribute all of the observed 44 TW to DM annihilation, when setting bounds on DM parameters. We similarly take a total heat flux of 3.5 TW for Mars, which is the maximum value of the projected range given in [35].
To calculate spin-dependent DM cross sections, we must determine the spin properties of the various elements in the Earth. The natural abundances of several elements with non-zero nuclear spins (J) have been tabulated in [53]. The values used in this work can be found in Table 2, which gives the fraction of the corresponding element in Table 1 that exists as the non-zero nuclear spin isotope listed. Later references to number density, n j (r), will carry different meanings in the spin-independent and spin-dependent cases. For the former, it will refer to the number density of a given element. For the latter, it will refer to the number density of a given elemental isotope with non-zero nuclear spin.
In general, a nucleus is able to have spin if one or more of its nucleons is unpaired. Paired nucleons' spins will cancel, leading to a net-zero nuclear spin. However, even if all nucleons of a given kind are paired, the expectation value of the spin of paired protons or neutrons ( S p and S n respectively) may be non-zero, which can result in spin-spin interactions. These terms are required to calculate spin-dependent per-nucleon cross sections (see Equation 7), but their exact values remain unknown, and are slightly model dependent. The values used in this work are tabulated for the elements of interest in Table 3.  For the purpose of this work, we have taken both the Earth and Mars to be modeled as perfect spheres with isotropic densities, temperatures, and compositions. We also ignore atmospheric scattering of DM, as the bulk of atoms within these planets is far larger and therefore dominates all scattering.

Planetary Heating Limits on Annihilating Dark Matter
With our planetary compositions and DM interaction models established, we will now bound DM's couplings using anomalous heating of the Earth and Mars. Results are given in terms of m χ , σ χN , and the DM-DM annihilation cross section σ χχ . Full exclusion limits obtained by requiring s-wave DM annihilations to not exceed the energy emitted from Earth's surface are shown in Figures 3, 4, 5, and 6. A full set of plots for s-wave, p-wave, and impeded [55] dark matter annihilation for both Earth and Mars can be found in the appendix.

Total Annihilation
To set a lower bound on the DM-nucleon cross section, we first consider the case that all captured DM annihilates, so that the rate at which DM is captured by the Earth or Mars is also the rate at which it annihilates. As we will see, this capture-annihilation equilibrium is reached in all parameter space of interest for a DM s-wave annihilation cross section σ χχ 10 −34 cm 2 .
To compute a limit for the case of DM capture-annihilation equilibrium, we began by running one thousand Monte Carlo simulations for each point on a grid of m χ -σ χN values. We chose to run one thousand simulations, for the following reasons. The total potential dark matter energy flux through the Earth, as described by equation 10, is around 3300 TW. Because our exclusion condition requires 44 TW of heating, only 44/3300 ≈ 0.01 of the velocity normalized dark matter flux must be captured to set a bound. By running a thousand simulations, we expect to have sampled ∼ 10 times more of the DM distribution than is strictly necessary to set a bound. For each of the thousand simulations, the DM particle in question was given a random initial velocity, distributed by a three dimensional Maxwell-Boltzmann distribution where this expression is the rate-normalized Maxwell-Boltzmann distribution [56]. We took standard values of v 0 = 220 km s −1 , with σ v = v 0 3/2 as the velocity dispersion, and with the velocity of DM in the planet's rest frame defined as v d = v i + v e where | v e | ≈ 230 km s −1 is the Earth's velocity in the galactic rest frame [56,57]. We normalize N e in this Maxwell Boltzmann distribution to match a conservatively low background DM density of ρ χ = 0.3 GeV/cm 3 [58][59][60] and truncate the distribution at a conservatively low galactic escape velocity of v esc = 528 km s −1 [61]. For each simulated DM particle, an entry angle θ into the planet in question was randomly chosen, distributed according to the probability density function which yields the average chord length traveled through the Earth, as described by Dirac's formula [62]. Using the generated θ and v i values in Equations 3, 4, and 6 or 7 for the spin-independent or spin-dependent cases respectively, and using Poisson-distributed τ j (see Equation 5), we found v f for all angles generated.
In each simulation, we found P cap , the probability of capturing a single DM particle entering the Earth, for a given DM mass and cross section. Defining v i as the average initial velocity of a captured DM particle, the mass capture rate (Γ m ) was calculated using the following equation for the total mass capture rate where we note that this is half the total flux of DM expected through the Earth's surface, since we are only interested in ingoing (and not outgoing) DM particles. By simulating one thousand DM particles per parameter space point tested, we identified the minimum nucleon scattering cross section for which more than ∼44 TW of DM would be captured by the Earth. We continued simulation iterations, until the range of cross section values converged upon varied by less than one percent, and among these values selected a cross section that implied slightly more than 44 TW of heating by DM. Assuming that all captured DM annihilates, this amount of darkogenic Earth heating can be safely excluded, since the total heat output of the Earth is 44 TW. This method was then repeated for Mars, this time using a maximum heat flux of 3.5 TW. The masses and cross sections that resulted inṀ χ ≈ 44 TW in the Earth are shown as the lower solid blue limits in Figures 3, 4, 5, and 6. The values generated for Mars are shown as solid red limits. The limit has been truncated at m χ = 10 10 GeV, a mass cutoff advocated in Ref. [1], chosen by requiring that the non-relativistic DM s-wave self-annihilation cross section required for DM to reach capture-annihilation equilibrium not exceed a unitarity limit [64]. The planetary heating bounds in this paper can be trivially extended to higher masses, though for such large masses and self-annihilation cross sections, it would be preferable to consider an explicit model.  Figure 3: Spin-independent DM-nucleon Earth heating limit assuming all DM annihilates (blue) and for DM self-annihilation cross sections σ χχ given in cm 2 as labelled, and the prospective Mars heating exclusion limit for σ χχ = 10 −36 cm 2 (red). The edges of these exclusion limits are determined as follows: the bottom edge is set by requiring a large enough cross section that enough DM is captured and anomalously heats the planet, the top edge is set by requiring that, despite having a large scattering cross section, DM can drift to the center of the planet within a billion years, and the left edge requires that DM is not so light that it evaporates out of the planet's interior. While the right edge is fixed at 10 10 GeV by s-wave annihilation unitarity considerations, the bounds shown could be consistently extended to higher masses. Juxtaposed are limits from [1, 7-13, 19-21, 63]. Underlaid, the dark grey line shows the Earth heating bound set by Mack et al. [1]. A number of recent references [20,39] have pointed out that only a few dark matter models will consistently provide an effective DM-nucleon cross section in excess of σ χN ≈ 10 −26 cm 2 . For a brief discussion of dark matter models which validly imply a DM-nucleon cross section greater than σ χN ≈ 10 −26 cm 2 , see the last paragraph of Section 2.

Thermal Evaporation
A gravitationally bound DM particle will not necessarily stay captured. Thermally vibrating nuclei in the Earth or Mars may scatter with captured DM and increase the DM's kinetic energy. For a sufficiently low DM mass, this kinetic energy increase can result in DM escaping its planet's gravity. This process is known as evaporation (see e.g. [66]). Evaporation implies a minimum mass for which Earth heating exclusions are valid, since evaporated DM will not annihilate within the planet.
First we consider the thermal radius, which is the radius of containment for DM that has reached thermal equilibrium with the planet -where DM's average thermal energy will be equal to its average potential energy in the Earth's gravitational well. The thermal radius is found using the virial theorem DM masses less than ∼100 MeV yield thermal radii greater than the radius of the Earth, meaning that the sphere of containment for DM at this mass extends outside the Earth. For Mars, this mass : Spin-dependent DM-nucleon Earth heating limit (blue) for DM self-annihilation cross sections σ χχ as labelled, and the Mars heating exclusion limit for σ χχ = 10 −36 cm 2 (red). Juxtaposed with limits given by [3,4,7,11,65]. This exclusion assumes equal spin-dependent coupling to neutrons and protons, aka isospin-independent scattering with a n = a p = 1 in Eq. (7).  : Spin-dependent DM-proton Earth heating limit (blue) for DM self-annihilation cross sections σ χχ as labelled, and the Mars heating exclusion limit for σ χχ = 10 −36 cm 2 (red). Juxtaposed with limits given by [3,4,11]. This exclusion assumes only spin-dependent coupling to protons, with a p = 1 and a n = 0 in Eq. (7).  Figure 6: Spin-dependent DM-neutron Earth heating limit (blue) for DM self-annihilation cross sections σ χχ as labelled, and the Mars heating exclusion limit for σ χχ = 10 −36 cm 2 (red). Juxtaposed with limits given by [7,11]. This exclusion assumes only spin-dependent coupling to neutrons, with a p = 0 and a n = 1 in Eq. (7).
is ∼420 MeV. To be conservative, we end our limit at this mass, as the upper exclusion limit loses its meaning when there is no defined thermal radius to which a DM particle would drift.
A caveat to the evaporation process as detailed above is that, for moderately large cross sections, low mass particles can be re-captured as they scatter against terrestrial or Martian constituent particles along their exit trajectory. Using the root mean square thermal velocity of a particle in equilibrium with the Earth or Mars, v th (r) = 3T (r) m , one could determine for what cross section light DM would be re-captured, and the new effective volume within which DM annihilates. We discuss what future work might be done along these lines in Section 5.

Partial Annihilation
For a small enough self-annihilation cross section, DM that is efficiently captured by the Earth or Mars may not lead to anomalous heating. Including the effect of DM annihilation, the differential equation describing the number of DM particles accumulated is where C χ =Ṁ χ /m χ is the number of DM particles captured per second and V th = 4πr 3 th /3 is the thermalized volume. Note that V th will decrease with increasing m χ . As is customary in the literature, we define the s-wave, p-wave, and impeded [55] dark matter annihilation rates such that the thermal velocity dependence of the DM annihilation rate in Eq. (12) scales as σ χχ v (χ) th = σ χχ , σ χχ v 2 th , σ χχ v th , respectively. Solving Equation 12 gives an expression for the number of DM particles captured over a period of time. For convenience, we define c ann = σ χχ v (χ) th /V th , which gives a compact expression for the number of DM particles in the planet, accounting for self-annihilation N χ , DM will have reached capture-annihilation equilibrium once t ∼ 1/ C χ c ann .
The equation for the annihilation rate of DM is then C ann = cannN 2 , and each DM-DM annihilation will have an energy of approximately 2m χ . Thus, we define a new equation for the heating rate induced by DMQ χ that accounts for the DM-DM annihilation cross section [67]: Again simulating many DM capture events as described in Section 4.1, but now also using the self-annihilation cross section to determine the total heating rate as described above, we determined values of m χ , σ χχ , and σ χN corresponding to a C χ that yielded a heat flux of ∼44 TW. By running simulations with multiple combinations of these three parameters, we found exclusion limits for values of σ χχ , shown as dashed lines in Figures 3, 4, 5, 6. Note that for a sufficiently large DM-DM s-wave annihilation cross section, the white partial annihilation lines converge on the total annihilation limit, at σ χχ 10 −44 cm 2 for both Earth Mars.

Drift Time
DM with a large enough nuclear cross section will be captured efficiently by the Earth and Mars. However, too large a nuclear cross section will result in DM stopping near the surface of the planet. As a result, it would be contained in a larger volume than DM that is free to descend to its thermalization radius. With its annihilation rate volumetrically suppressed, the DM will not necessarily annihilate quickly enough to heat the planet in which it is captured. Here we set an upper cross section limit by ensuring that the captured DM is able to drift through the planet in question on a timescale less than the planet's age.
To compute this upper limit on our cross section exclusion, we used a treatment similar to [1,68,69]. Assuming DM-nuclear interactions are frequent enough, the planet's gravitational force will balance against viscous drag. We require DM to drift from the surface of the planet to its thermal radius (described in section 4.2) within 4.5 Gyr, the approximate age of both the Earth and Mars. Balancing these two forces yields where G is the gravitational constant, v j (r) = 3T (r)/m j is the thermal velocity of molecule j, M (r) is the mass enclosed in radius r (equal to the volume integral over ρ(r)), and v drif t = ∂r/∂t is the drift velocity of the DM [69]. Making these substitutions, Equation 15 becomes where we conservatively set t = 1 Gyr in calculations. Our upper exclusion limits, shown as the upper solid lines in Figures 3, 4, 5, 6, were found by solving Equation 16 for σ χN at a given m χ . In the spin-independent case, we used Equation 6 to convert from per-nucleon to nuclear cross section.
In the spin-dependent case, we used Equation 7, and values given in Table 3.

Discussion
We have derived Earth and Mars heating bounds on both spin-dependent and spin-independent DM scattering. Of course, these bounds will apply to DM which annihilates to Standard Model particles. DM that does not annihilate, or which annihilates to particles which freely stream out of the Earth or Mars are not excluded by this result. We have also determined how Earth and Mars heating bounds change as DM's self-annihilation cross section is varied. The limits we have found on dark matter annihilation in the Earth are especially interesting, considering that the canonical WIMP self-annihilation cross section is σ χχ ∼ 10 −36 cm 2 [70]. While the Earth heating limit appears to exclude the canonical WIMP self-annihilation cross section for s-wave dark matter annihilation, p-wave and impeded annihilating dark matter are not entirely excluded. As such, a major result presented in this work for the first time, is that DM planetary heating bounds do not necessarily constrain p-wave and impeded DM models.
In future work, there are additional improvements that could be made to the treatment of DM heating planets. For lighter DM, heating is limited by dark matter evaporation. The treatment of DM evaporation presented here can be improved using Monte Carlo simulations that account for incoming light DM scattering repeatedly not only during capture, but also during evaporation. Captured DM that interacts strongly enough with nuclei, may avoid evaporation by being trapped via back-scattering. Such a future analysis might result in better limits on light DM.
The potential effectiveness of this analysis for other planets in our solar system should also be considered. While the Earth and Martian interior compositions are known with most certainty, based on seismic data and surface mineral sampling, relatively thermally inactive bodies such as Earth's moon, which has a maximum DM capture rate of ∼250 TW, and a projected internal heat flow of ∼0.75 TW [71], could be studied to set additional limits on DM. Analysis of lunar chemical composition indicates a significant presence of heavier elements with non-zero nuclear spin [72], meaning that improved spin-dependent DM scattering limits may be possible. On the other hand, temperatures very near the Moon's center are still largely unknown, and would require further study.
Finally, for the case of DM with predominantly spin-dependent nucleon interactions, our results have shown that substantial parameter space exists for low dark matter masses, where experiments like PICO might be able to set new limits on or discover strongly-interacting DM by placing detectors above-ground. We leave a detailed study of the prospects for above-ground spin-dependent dark matter searches to future work.   Figure 9: Spin-dependent DM-nucleon Earth (left) and prospective Mars (right) heating exclusion limits for DM s-wave self-annihilation cross sections σ χχ as labelled, assuming spin-dependent coupling to protons, with a p = 1 and a n = 0 in Eq. (7).  Figure 12: Spin-dependent DM-nucleon Earth (left) and prospective Mars (right) heating exclusion limits for DM p-wave self-annihilation cross sections σ χχ as labelled, assuming equal spin-dependent coupling to neutrons and protons, aka isospin-independent scattering with a n = a p = 1 in Eq. (7).  Figure 13: Spin-dependent DM-nucleon Earth (left) and prospective Mars (right) heating exclusion limits for DM p-wave self-annihilation cross sections σ χχ as labelled, assuming spin-dependent coupling to protons, with a p = 1 and a n = 0 in Eq. (7).  Figure 14: Spin-dependent DM-nucleon Earth (left) and prospective Mars (right) heating exclusion limits for DM p-wave self-annihilation cross sections σ χχ as labelled, assuming spin-dependent coupling to neutrons, with a p = 0 and a n = 1 in Eq. (7).   Figure 16: Spin-dependent DM-nucleon Earth (left) and prospective Mars (right) heating exclusion limits for DM impeded self-annihilation cross sections σ χχ as labelled, assuming equal spindependent coupling to neutrons and protons, aka isospin-independent scattering with a n = a p = 1 in Eq. (7).  Figure 18: Spin-dependent DM-nucleon Earth (left) and prospective Mars (right) heating exclusion limits for DM impeded self-annihilation cross sections σ χχ as labelled, assuming spin-dependent coupling to neutrons, with a p = 0 and a n = 1 in Eq. (7).