Direct Detection Anomalies in light of $Gaia$ Data

Measurements from the Gaia satellite have greatly increased our knowledge of the dark matter velocity distributions in the Solar neighborhood. There is evidence for multiple cold structures nearby, including a high-velocity stream counterrotating relative to the Sun. This stream could significantly alter the spectrum of recoil energies and increase the annual modulation of dark matter in direct detection experiments such as DAMA/Libra. We reanalyze the experimental limits from Xenon1T, CDMSlite, PICO-60 and COSINE-100, and compare them to the results of the DAMA/Libra experiment. While we find that this new component of the dark matter velocity distribution can greatly improve the fit to the DAMA/Libra data, both spin-independent and spin-dependent interpretations of the DAMA/Libra signal with elastic and inelastic scattering continue to be ruled out by the null results of other experiments, in particular Xenon1T.


I. INTRODUCTION
There is strong gravitational evidence that 26.5% of the Universe's energy budget and the majority of the mass within galaxy clusters is composed of dark matter; an unknown, invisible, non-relativistic material. No particle in the Standard Model (SM) has the appropriate properties to be dark matter, thus evidence of dark matter is evidence of new physics. If dark matter consists of a particle with a mass GeV and has small but non-zero couplings to the SM baryons (as suggested by the thermal relic freeze-out scenario), then direct detection experiments can be a powerful probe of the dark sector.
Direct detection searches for a dark matter particle passing through a low-background detector undergoing a nuclear recoil with one of the detector's atoms. This recoil can then be seen in the experiment, assuming that the recoil momentum is above the detector threshold. The number of events in a given direct detection experiment is then set by a confluence of factors: the parameters in the dark matter theory space (e.g., mass, scattering mechanism, and cross section), the detector energy threshold and efficiency function, and the astrophysical distribution of dark matter (the local density and local velocity distribution). From a particle theorist's point of view, these latter parameters are confounding variables which prevent a straightforward mapping from the experimental results to exclusion or discovery within the theory space.
In particular, if dark matter near the Sun is moving with velocities significantly different from the baseline expectations, then the relative sensitivity between different experiments can also be drastically altered. Indeed, as has been well-studied in N -body simulations [1][2][3][4][5][6][7][8][9], the dark matter distribution in the Milky Way is not expected to trace the Maxwell-Boltzmann Standard Halo Model (SHM). Since the Milky Way halo was constructed through the hierarchical merger of smaller subhalos [10], we should expect streams and tidal debris [4] of dark matter, giving additional velocity structure on top of the smooth halo. These results have significant implications for dark matter direct detection, even if the stream contributes only a small component of the dark matter density [11,12].
Given this background, the results of the Gaia mission [13] are of timely interest. Gaia measures the position and motion of the nearest ∼ 1.4 billion stars, mostly within ∼ 2 kpc of the Sun's location. Combined with metallicity information, this enormous sample allows for kinematic studies of the stellar halo with a resolution and coverage not possible before. Given that the stars in the halo can be used to trace the dark matter [14], Gaia can infer a great deal about the dark matter velocity distribution in the Solar neighborhood. As expected from N -body simulations, the smooth halo component of dark matter is suppressed at high-velocities [15][16][17], leading to weaker direct detection sensitivity for low mass dark matter as compared to the predictions of the SHM. Perhaps more surprisingly, there is evidence that the Sun is in the path of a number of high velocity dark matter streams [18]. One such stream (the S1 stream) is counter-rotating relative to the Sun's motion through the Galaxy, increasing the relative motion of the dark matter within the stream and the Earth. This stream is consistent with being part of a remnant of a ∼ 2 × 10 10 M dwarf galaxy, which was tidally stripped by the Milky Way several billion years ago. While the dark matter density of S1 has not yet been measured, it could be a significant component locally, perhaps O(10%) [19].
The impact of both the Gaia-derived halo velocity distribution and that of the S1 stream on the future direct detection limits have been studied in prior literature [15,19]. However, until this work, their effect on existing limits and potential signals have not been fully calculated. Of particular interest is the effect of S1 on direct detection experiments that search for a yearly modulation of scattering events, rather than the overall rate averaged over a year.
As the Earth moves around the Sun, the fraction of dark matter particles with a relative velocity capable of surpassing the detector threshold increases when the Earth's motion is into the local dark matter "wind," and decreases as the Earth moves with the flow of dark matter. As a result, the direct detection rate should modulate over a year.
As the S1 stream is itself fast-moving and counter-rotating to the Sun's motion, it could induce a very large yearly modulation -in some cases nearly an order of magnitude more than one would expect from the dark matter halo itself. As the stream is kinematically cold, the velocity profile is narrow, leading to a sharp peak in the nuclear recoil spectrum, as would be seen by a direct detection experiment (at a recoil energy set by a combination of the dark matter mass and the target nuclei). Additionally, the S1 stream happens to be oriented in such a way that the modulation peak occurs on a day in early June, near the day one would predict from the non-rotating halo.
As is well-known, the DAMA/Libra experiment [20] observes a yearly modulation in scattering events with a peak date in June [21]. This signal is also compatible with the initial modulation analysis of the COSINE-100 experiment [22], though the results have nearly equal preference for the null (i.e., no dark matter signal) hypothesis. Assuming the SHM or a smooth halo profile derived from simulation, it is difficult to interpret the DAMA/Libra signal as dark matter, given the negative results from other experiments, for both spin-independent (SI) and spin-dependent (SD) interactions [23][24][25][26][27][28]. However, in light of the Gaia data and the S1 stream we should revisit this conclusion, as these could greatly increase the modulation signal without as significant a change in the limits set by other experiments, as well as alter the observed recoil spectrum in a way that might better fit the data.
Given this motivation, using the velocity distributions extracted from Gaia data, we reanalyze the direct detection limits from the DAMA/Libra experiment along with the experiments that set the strongest current limits for spin-dependent and spin-independent scattering: Xenon1T [29], CDMSlite [30], PICO-60 [31], and the COSINE-100 rate measurement [32], considering a range of possibilities for the couplings of dark matter to protons and neutrons. We show that Gaia-derived halo models mildly weaken the bounds on dark matter interactions with baryons for lower mass dark matter when compared with the predictions from the SHM. Addition of the S1 stream can improve the statistical fit of the DAMA/Libra signal to a dark matter interpretation (compared to the fit to a dark matter velocity distribution without a stream), and shifts the best-fit region to lower masses and cross sections. Despite these significant changes, we find that for all tested scenarios of elastic and inelastic scattering (assuming either SD or SI coupling) the best-fit regions of the DAMA/Libra signal continue to be excluded by null results of other experiments (most importantly, Xenon1T), unless the S1 stream is the dominant contributor to the local dark matter density ( 80%).
After a review of dark matter scattering physics in Section II, in Section III we examine the local velocity distribution as determined by the Gaia satellite in the context of dark matter direct detection experiments, including both the background halo distribution and the S1 stream. We then examine the impact of these velocity distributions on the existing experimental limits of Xenon1T, CDMSlite, COSINE-100, and PICO-60, and the best-fit regions of DAMA/Libra, for both SD and SI couplings, varying degrees of isospin violation, and varying S1 dark matter densities. Details of our recasting of experimental results are described in Appendix A.

II. REVIEW OF DARK MATTER DIRECT DETECTION
The differential rate dR/dE R for dark matter scattering off a target nucleus of atomic mass M in a particular direct detection experiment is given by expressed in counts/day/kg/keV with E R the nuclear recoil energy. The quantity ρ 0 is the local dark matter density (to which we assign a canonical value 1 of 0.3 GeV/cm 3 ), µ = M m χ /(m χ + M ) is the nucleus-dark matter reduced mass, and F (q) is a form-factor which depends on the transferred momentum q = √ 2M E R . For spin-independent scattering we adopt the Helm form factor [40] F (q) = 3(sin qr − qr cos qr) where (following Ref. [41]), r = c 2 + 7/3 π 2 a 2 − 5s 2 is the nuclear radius, c = 1.23A 1/3 − 0.6, a = 0.52, s = 0.9, with A the atomic mass of a specific nuclear target. The spin-dependent form-factor goes as where S(q) = a 2 0 S 00 (q) + a 0 a 1 S 01 (q) + a 2 1 S 11 (q) [41][42][43][44] and S(0) is the normalization of the structure functions in the zero momentum transfer limit. The S ij contain information on the protons, neutrons, and their interference. For the different target nuclei we use the nuclear shell model fits derived in the appendix of Ref. [43] and for the nucleon coefficients a i , we use the values in Refs. [45,46].
The quantityf (v) in Eq. (1) is the lab frame dark matter velocity distribution and is integrated from the detector dependent velocity v min (E R ) to the maximum velocity where the density distribution has support. This is set by the escape velocity of the Galaxy, v esc . The velocity v min is the minimum dark matter velocity capable of inducing a nuclear scattering event with recoil energy E R . For elastic collisions, it is For inelastic collisions, a dark matter ground state χ may up-scatter off nuclei in the detector to an excited dark sector state χ * , with mass difference given by δ = m χ * − m χ . In this case the minimum velocity required to induce a nuclear recoil of energy E R is [47] v idm For either elastic or inelastic scattering, the observed number of events in a detector is given by the differential rate integrated over some range of recoil energies times a detector-dependent efficiency factor (E R ) and the detector exposure. As a result, the rate depends on experimental effects (exposure, thresholds, target mass, etc), particle physics parameters (dark matter mass and scattering cross section), and astrophysical parameters (local dark matter density and velocities). The particle physics information is encapsulated in the dark matter-nucleus scattering crosssection, σ N , in Eq. (1). To make contact with theoretical models, and to be able to compare experimental results between target nuclei, it is more convenient to write this in terms of the scattering cross section between dark matter and individual nucleons.
In this paper, we consider two options for the dark matter-nucleon scattering: spin-independent and spin-dependent interactions. Though momentum-dependence can be introduced into the cross section by an appropriate choice of the particle physics interaction [48][49][50][51][52], in this paper we assume the cross-section is independent of q 2 , with the momentum dependence only in the form factor.
For spin-independent couplings, σ N can be expressed in terms of the dark matter-proton and -neutron couplings f p and f n as where Z is the atomic number, A the atomic mass of the nucleus, µ p the reduced proton-dark matter mass, and σ SI p the dark matter-proton scattering cross section. Dark matter direct detection limits are canonically presented in terms of σ p assuming isospin-conservation. We will also consider isospin-violating interactions, as such interactions can change the relative signal rate between experiments with different target nuclei (with different ratios of Z and A). For isospin-violating couplings, we parametrize the violation with the ratio and report the overall cross section in terms of the equivalent cross section σ SI p . For dark matter scattering off a target that is composed of multiple isotopes (of one or more elements), the scattering rate measured would be the sum over the rates from Eq. (1) for each isotope, weighted by the isotope abundance. In principle, this means that certain direct detection target materials could have greatly suppressed rates, if the isospin-violating coupling f is tuned such that  (9), for spin-independent couplings which would apply to experiments using sodium, iodine, xenon, or germanium targets, assuming natural abundances of isotopes. Vertical grey lines at f ∼ −0.700 and −0.785 denote the minimum sensitivity of xenon and germanium detectors, respectively. Right: Average suppression factor S for spin-dependent couplings which apply to sodium, iodine, xenon, germanium, or fluorine targets, and with the constraint |ap| + |an| = 1. Vertical grey lines as ap ∼ −0.098 and −0.88 correspond to minimum sensitivity for fluorine and xenon, respectively, see Eqs. (11) and (12).
for the target isotope [53]. Even ignoring the coincidence this would imply, for many of the experiments setting the strongest current constraints the scattering cross section cannot be tuned to zero, because the natural abundance of the target elements consists of multiple isotopes. For such targets, which include xenon-and germanium-based detectors, a single value of f cannot cancel away interactions with all of these nuclei. In Figure 1, we show the abundanceaveraged suppression factor -defined as the reduction in scattering cross section relative to the isospin-conserving case -as a function of the isospin violating parameter f : where the sum runs over isotopes i with abundance P i , atomic numbers Z i , and mass numbers A i . This factor S(f ) gives an estimate in the reduction in sensitivity for a given experiment. For spin-dependent interactions, the equivalent of Eq. (6) can be written as where a p and a n are the proton-and neutron-couplings (the spin-dependent equivalents of f p and f n ). The total nuclear spin is J, and the spin-expectation values for the proton and neutrons are S p and S n . The spin parameters for nuclei relevant to direct detection are listed in Table I. As nuclear spins are typically small, and Eq. (10) lacks the equivalent of the A 2 enhancement of Eq. (6), the experimental reach on spin-dependent scattering (expressed in terms of σ SD p ) are much weaker than the spin-independent ones. Somewhat confusingly, in the direct detection literature, the "proton" and "neutron" couplings a p and a n are defined in terms of the tree-level Lagrangian couplings. Pion exchange will mix these couplings, introducing a coupling to neutrons even in the limit a n → 0, or a coupling to protons when a p is zero. We follow the calculations of Ref. [43] for the size of this mixing effect: a p → a p + 1 2 (a p − a n )δa, a n → a n − 1 2 (a p − a n )δa, where the momentum-dependent mixing parameter δa is generically ∼ −20%. As in the spin-independent case, isospinviolation for spin-dependent interactions corresponds to a ratio a n /a p = 1. In the right-hand panel of Figure 1, we  We use the proton and neutron expectation values compiled from nuclear shell models in Ref. [43] show the equivalent of Eq. (9) for spin-dependent scattering, defined as the suppression (or enhancement) of the cross section relative to the isospin-conserving assumption of a p = a n : where we have imposed the constraint |a p | + |a n | = 1 (this normalization is chosen as several nuclei relevant for spindependent scattering have either S p S n or S p S n , making a n /a p 1 and a n /a p 1 both potentially interesting). Following Ref. [43], at zero momentum exchange, the mixing induced by pion exchange in Eq. (11) is taken to be δa ≈ −0.2.
Importantly, in Eq. (1) the dependence of the rate on the dark matter velocity distribution can be factorized from the particle physics and experimental effects up to the dependence on v min (E R ). Different experiments can then be compared in an "astrophysics-independent" manner [54,55] by comparing the rates observed as functions of v min . A particular dark matter direct detection experiment is sensitive to those dark matter particles moving fast enough (in the lab-frame) to scatter with a nucleus in the target, imparting a recoil energy above the detector threshold. The rate can then be expressed as a convolution between the experimental response and a piece depending on the dark matter velocity distribution [15,27], in particular the average of the inverse velocity, integrating over all velocities above an experiment-dependent v min : The velocity distribution in the lab-framef must be related to the distribution in the Galactic rest frame (to which the lab is moving with a relative velocity v lab ) viaf ( v, t) = f ( v + v lab (t)). However, this "astrophysics independent" comparison is made more difficult for the DAMA/Libra results, as this experiment does not attempt to measure the overall rate of dark matter scattering, but rather the modulation in the event rate due to the motion of the Earth around the Sun, which causesf (v, t) to vary over the year. In order to compare this experiment with those that only measure the time-integrated rate, some model of the dark matter velocity distribution must be adopted. In the next section, we review the standard choice for this distribution, as well as the discoveries about the local velocity distribution made possible by the Gaia mission.

A. Dark Matter in the Halo
The simplest assumption of the local motion of dark matter is that the local dark matter velocity in the Galactic rest frame is an isotropic Maxwell-Boltzmann distribution -the Standard Halo Model (SHM): where the normalization factor is  [16] from Gaia data for the halo (red dotted), substructure (red dashed) and sum (red solid lines, with shaded region indicating uncertainties). The blue line is the normalized velocity distribution for the S1 stream, with shading indicating uncertainties. Recall that, though we normalize all distributions to unity in these plots, the stream will not contain 100% of the local dark matter density.
We adopt the SHM parameters v 0 = 220 km/s for the dark matter mean speed, and a Galactic escape velocity of v esc = 544 km/s [56]. The SHM distribution f (v) × 4πv 2 in the Galactic rest-frame is shown in the left panel of Figure 2.
To make contact with experiments, which are sensitive to the lab-framef rather than the rest-frame f , we must also specify the relative velocity of the Earth and Sun with respect to the Galactic standard of rest. In Galactic coordinates (x radial from the Galactic Center, z out of the disk), the Sun's velocity is v = (10, 233, 7) km/s. (16) To this, we must add the time-dependent Earth velocity [57] v where ω = 2π/365.25 days, with the phase shift t 1 = 79.62 (relative to January 1 st ). The yearly modulation of the total velocity of the Earth+Sun system will modify the integral overf (v) for a fixed v min , thus causing a yearly modulation in a direct detection experiment's sensitivity and resulting number of events, as reported by DAMA/Libra. For a dark matter distribution that can be modeled as an isotropic Maxwell-Boltzmann in the Galactic rest frame, such as the SHM, the η function has an analytic form (see e.g. Ref. [57]) The η function of the SHM with the benchmark parameters is shown in the right panel of Figure 2, while the modulation of η over the year is shown in Figure 3. Notice (in the left panel of Figure 3) that the modulation ∆η(t) ≡ η(t) −η switches sign as v min is increased, from a minimum around June 2 nd to a maximum on that same day. However, it has long been known [1-3, 5-8, 58, 59] that the true distribution of dark matter in the Milky Way-like halos should deviate significantly from this idealized distribution. Dark matter direct detection experiments often take this uncertainty into account when reporting limits [30,[60][61][62].
With the Gaia DR2, our expectation that dark matter should deviate from the SHM can be directly tested, and a halo distribution inferred using low-metallicity stars in the Galactic halo [15,16]. Thus, the astrophysical uncertainty in direct detection experiments can be significantly reduced. We adopt the distributions of Ref. [16] for the model of the dark matter halo. We also note the recent work of Ref. [62], which also derived modifications to the dark matter velocity distribution and the resulting implications on direct detection from Gaia data; in that case the primary effect coming from the changes in the escape velocity of the halo.
The Gaia-derived distributions used in this work have two components: dark matter from a smooth halo which has been integrated into the Galaxy, and a substructure of dark matter from tidal debris [4,16,63], which is likely the remnant in velocity-space of long-ago major merger events. The total dark matter distribution in the local volume of the Milky Way is then c halo f halo + c sub f sub with [16] In Figure 2 we show the SHM velocity distribution and η functions with these data-driven distributions. The datadriven model of the velocity distribution has fewer high-velocity particles, which lead to weaker bounds for light dark matter (since as the dark matter mass is lowered, v min correspondingly increases). Similarly, as seen in Figure 2, the amplitude of ∆η decreases for v min 350 km/s for this Gaia-derived halo model (though the peak day is nearly unchanged). Given the four to six orders of magnitude between the DAMA/Libra best fit region and the Xenon1T limits (assuming spin-independent isospin-conserving couplings), the DAMA/Libra region is safely ruled out for modifications of the smooth dark matter velocity distribution (assuming velocity-independent cross sections), even with these slightly weakened limits (see, e.g., Ref. [8]).

B. The S1 Stream
The possibility that the Earth lies in the path of multiple kinematically cold streams of dark matter has long been known. Such streams will modify the dark matter velocity distribution in ways not captured by simulations of the overall halo and can have large effects on the modulation measured in direct detection [11]. The stream S1, identified in Ref. [18], provides a concrete example and requires us to revisit the conclusion that the null results of other experiments definitively exclude the DAMA/Libra best-fit regions. S1 is a very high velocity stream, with a velocity of ∼ 300 km/s relative to the Galactic rest frame. These stars are coherent in velocity space, and though fairly spread out in position space, do overlap with the position of the Sun within the Milky Way. Such "tidal debris" is expected from simulations [4,16,63], and the stream S1 can provisionally be identified as the stellar remnants of a ∼ 2 × 10 10 M dwarf galaxy which was absorbed by the Milky Way some 10 billion years ago. The contribution of the S1 stream to the local dark matter density is not known at this point, and in this paper we will treat it as a free parameter. Note that the S1 stream is within the Galactic disk, and so was not part of the Gaia stellar sample that Ref. [16] used to construct their model of the overall halo distribution.
Relevant for the DAMA/Libra experiment, the stream is essentially anti-parallel to the Sun's motion through the Milky Way. This is important for two reasons. First, as a result of this antiparallel motion, the local dark matter wind from the stream has a very large relative velocity, which can increase the magnitude of a modulation signal. Second, as the Earth moves around the Sun, the velocity distribution of the dark matter stream peaks very close to the date one would predict from the SHM alone.
The DAMA/Libra search (and claimed positive signal) is reliant on the yearly modulation of the dark matter η function, due to the Earth's motion around the Sun. In order to match the observed signal, the modulation signal must peak around June 2 nd (t = 152.5). This is close to the expected signal peak for a SHM dark matter distribution, as it is close to the day that Earth's velocity relative to Galactic rest is maximized (as noted previously, the improved dark matter distributions derived from Gaia data have the same peak day).
Though the density is not known, the velocities of the stars within S1 are known, and from this, a velocity distribution can be modeled. We assume in this paper that the stars are good tracers of the dark matter; this conclusion has been supported by N -body simulation [16]. Certainly one should expect streams of dark matter to exist without accompanying stars, but not necessarily the reverse: stars from tidally disrupted dwarf galaxies should be accompanied by dark matter.
The stars in S1 are counter-rotating relative to the Sun, with mean velocity v S1 = (8.6, −286.7, −67.9) km/s in Galactic cylindrical coordinates. Assuming that the velocity distribution of the stars in the Galactic rest frame is drawn from a three-dimensional Maxwell-Boltzmann distribution with diagonal velocity dispersion matrix v 0 : The f (v) and η distribution for the S1 stream are shown in Figure 2, and the yearly modulation of η is shown in Figure 3. Notice that for large v min 450 km/s, the magnitude of ∆η can be nearly an order of magnitude higher than in either the SHM or the data-driven model. While η itself also increases, for very high v min ∼ 500 − 600 km/s, the relative increase in ∆η as compared to η means that it is possible for experiments which measure modulation to get very large boosts in sensitivity while those that measure only average rate would be comparatively unchanged. Furthermore, even if the stream contributes only O(10%) to the local dark matter density, the sharp feature in the recoil spectrum it could induce could be visible in a direct detection experiment, and change the particle physics parameters of a best-fit point. It is these relative changes in signal sensitivity due to the S1 stream that we are interested in here.

IV. EXPERIMENTAL LIMITS
We are now in the position to reassess the current experimental limits in light of our new understanding of the local dark matter velocity distribution.
We are interested in determining the exclusion reach and signal regions of the existing direct detection experiments in light of the Gaia-derived halo mode and in the presence of the high-velocity S1 stream. In particular, does dark matter moving in this stream allow the DAMA/Libra best-fit region to evade the other experiment's null constraints?
The presence of the S1 stream induces features in the recoil spectrum that would not be expected in a smooth halo model, and so may be of experimental interest. As astrometric surveys of the local neighborhood are still in their early days, it may be expected that more streams and tidal debris will be identified in the future, in which case this study of the S1 stream can provide an example of what effects on the results of direct detection experiments can be expected. As we do not know the fraction of the local dark matter density which results from the S1 stream we will treat this as a free parameter; the density of the S1 stream may be better constrained with future astrometric measurements and comparison with simulations.
We consider the current strongest constraints on dark matter spin-independent and spin-dependent scattering, for dark matter masses between 1 and 10 4 GeV. Under the assumption of the SHM, for spin-independent searches, the strongest constraints are set by Xenon1T [29] for most of the dark matter mass range, surpassed by the germanium-based CDMSlite detector [30] in the low mass region. The COSINE-100 [32] experiment does not set the world-leading limits for any mass range, however, we include it in this analysis as it is composed of the same target material as DAMA/Libra: sodium-iodide crystals. This is important as we consider isospin-violating couplings, which can change the relative strength of the scattering of dark matter against the different target materials. The strongest limits on spin-dependent scattering are set by the PICO-60 experiment [31] for scattering against the proton, while the Xenon1T limits can be reinterpreted in terms of spin-dependent scattering to give the strongest limits for spin-dependent scattering against the neutron [64].
In order to extract limits (or for DAMA/Libra, best-fit signal regions) for velocity profiles other than the SHM, we must recast each of these experimental results. Our procedure for each experiment is explained in Appendix A. Our methods do not recover the exact experimental results: for the low mass region our exclusion regions are somewhat weaker than the official limits. As this is the mass range that will be of interest in comparing with DAMA/Libra, our results are in that sense conservative. We show the 90% CL upper limits from each experiment under various assumptions of scattering interactions and velocity distribution.
For the DAMA/Libra fits, we are interested in both the spectrum of the recoil energies and the yearly modulation. For the former, data in narrow energy bins have been made available by DAMA/Libra, but these assume a yearly modulation peaking on the day predicted by the SHM. Data binned in time-series are only available for a few overlapping energy bins. We therefore provide two sets of fits to DAMA/Libra: "Amplitude" fits to the binned recoil spectrum (where the binning in time assumes a yearly modulation), and "modulation" fits to the time series data, binning recoil events between 1 and 6 keV ee . As will be seen, the modulation fits demonstrate that the S1 stream is consistent with the same peak date as the smooth halo, and so the amplitude fits can be used for more fine-grained distinction. For both modulation and amplitude fits, we fit the signal hypothesis by minimizing a χ 2 fit to the available data. Unless otherwise noted, we will show the region of parameter space that is within 2σ of the minimum χ 2 , taking into account the appropriate number of degrees of freedom.

A. Spin-Independent Elastic Scattering
Combining the limits and best-fit regions for all the experiments, we show in the left panel of Figure 4 the isospinconserving limits assuming the data-driven dark matter velocity distribution derived from Gaia data in Ref. [16]. For each limit curve, the substructure ratio c sub /c halo is varied over the 1σ range of Eq. (21), while for the fits to DAMA/Libra, this is treated as a free parameter and is allowed to float with the interaction cross-section and dark matter mass. The limits from the SHM are also displayed. As can be seen, Xenon1T sets the strongest limits over the majority of the mass range, with CDMSlite taking over below ∼ 5 GeV. DAMA/Libra is decisively excluded by Xenon1T, with COSINE-100 excluding the claimed signal for the entire 2σ region as well. For the fit to the binned amplitude data, the best-fit DAMA/Libra point with f = +1 has χ 2 /d.o.f. = 3.5 (treating m χ , σ SI p , and c sub /c halo as free parameters), and the 2σ regions are relative to this minimum.
In the right panel of Figure 4, we show the limits assuming the isospin violating parameter f = −0.7 (the value that minimizes the sensitivity of xenon-based detectors, illustrated in Figure 1). As can be seen, changing this parameter also has the effect of suppressing scattering on iodine and enhancing sodium scattering, moving the DAMA/Libra best fit region of the amplitude data from m χ ∼ 60 GeV to m χ ∼ 15 GeV. The good-fit regions to the daily modulation data are consistent with those of the amplitude data, but much less constraining. This is not surprising, given that the modulation data is binned much less finely in recoil energy, which is where most of the discriminating power in mass comes from. Given this disparity, we will mostly refer to the good-fit regions for the amplitude data, using the modulation fits to demonstrate the important point that the non-SHM velocity distributions continue to have a peak day which can generally fit the DAMA/Libra observations.
Plotting the DAMA/Libra 2σ regions relative to the minimum χ 2 assuming a given value of f obscures how good the overall fit is for the assumed parameter of each plot. To investigate this, we show in Figure 5 both the best overall fit with f = +1 as well as the fit with m χ = 15 GeV and f = −0.7. This point is a good fit to the daily modulation data, but a significantly worse fit to the recoil energy spectrum, with χ 2 /d.o.f. = 7.4. The poor fit of the low mass DAMA/Libra region to the Phase-2 spectrum assuming the SHM has been noted previously [27]. The resulting predictions for the daily modulation of the event rate, compared to the published Phase-2 data, is shown in Figure 6, which demonstrates that the Gaia halo models have the appropriate phase shift to match the observations. We next consider the effect of the S1 stream, treating its contribution to the local dark matter density as a free parameter. To demonstrate the effect of the new velocity profile on the direct detection experiments, on the right panel of Figure 6 we plot the daily modulation event rate on top of the DAMA/Libra Phase-2 data and show that the S1 stream peaks on nearly the same day as the data for our best fit parameters. We also show in the left panel of Figure 7 the extrapolated limits if the S1 stream was 100% of the local density, assuming isospin conservation (f = +1). We vary stream parameters from the 1σ limits of Eq. (23).
As expected from the behavior of the η function shown in Figure 2 the exclusion limits for all experiments strengthen somewhat at low dark matter mass and weaken at high mass. Further, the DAMA/Libra regions move to lower mass and lower cross section; the best-fit parameter point occurs at either ∼ 26 GeV or ∼ 7 GeV, depending on the stream velocity parameters. However, the χ 2 /d.o.f. for these parameter points is between 8 and 15 when fit to the binned recoil spectrum (counting the stream parameters as additional degrees of freedom). This indicates that the stream, by itself, is not a particularly good fit to the measured DAMA/Libra recoil spectrum (though it is still statistically preferred over the no-signal hypothesis). An example of one of these best-fit parameters for 100% stream density and f = +1 is shown in the right panel of Figure 5, while the modulation fit is shown in Figure 6. Critically, in the latter fit, it is apparent that the phase shift of the stream is largely consistent with the observations, with the peak modulation occurring at t = 151 days.
Of course, it is not reasonable to expect that the S1 stream will constitute 100% of the local dark matter density. To display limits in the multi-parameter space of dark matter mass m χ , cross section σ SI p , isospin-violating parameter f , and stream fraction, we take two approaches. First, as an example, we show in the center and right panels of Figure 7 the extrapolated limits and best-fit DAMA/Libra regions assuming the stream is 25% of the local density -this assumption is perhaps on the high end of the O(10%) estimate, but not implausible. Again, we show limits for two choices of the isospin parameter: f = +1 and f = −0.7. The former has a minimum χ 2 /d.o.f. = 3.9, while the latter has a minimum goodness-of-fit of 6.7. Again, we see that the stream can shift the DAMA/Libra best-fit to lower masses and therefore bring it closer to the edge of the various exclusion limits, but in the representative parameter choices, the DAMA/Libra region continues to be ruled out.
However, this set of plots does not completely guarantee that there is not some combination of dark matter mass, cross section, isospin-violation, and stream fraction that would not allow the DAMA/Libra region to evade the  existing limits. We therefore show in Figure 8 two scans: the first over dark matter mass m χ and f (left panel) and the second over m χ and the fraction of the local density in the S1 stream (right panel). As before, we select stream parameterizations within the 1σ errors of Eq. (23). In these plots we display the regions of parameter space that are within 1σ, 2σ, and 3σ of the best-fit. The best fit point corresponds to m χ = 38 GeV, f = −0.60, σ SI p = 2.5×10 −39 cm 2 , and a stream fraction of 52%, with a χ 2 /d.o.f. = 0.79. A slightly worse fit (χ 2 /d.o.f = 1.2) exists at m χ = 27 GeV, f = −0.67, σ SI p = 8.4 × 10 −39 cm 2 , and a stream fraction of 34%. We show the recoil spectrum and daily modulation for this latter point in Figures 5 and 6.
From this analysis, we can see that the addition of the stream allows for a markedly improved fit to the observed spectrum, for reasonable values of the stream density. We can therefore conclude that the addition of streams can FIG. 7: 90%CL upper limits from Xenon1T (black), CDMSlite (blue) and COSINE-100 (red) on the spin-independent dark matter-proton scattering cross section σ SI p assuming the S1 stream is 100% of the local dark matter density and an isospin parameter of f = +1 (left), 25% of the local density and f = +1 (center), and 25% of the density and f = −0.7 (right). The 2σ regions around the best-fit to DAMA/Libra data are shown in orange for fits to the reported modulation amplitudes (with the lighter shaded region indicating the 1σ variation of the S1 stream velocity distribution), and in yellow for fits to the reported daily modulation rates. Corresponding regions assuming the SHM with f = +1 are shown in grey shading.
noticeably alter the expected spectrum of dark matter recoils in direct detection, when compared to the SHM or smooth halo models predicted from simulation. This is perhaps interesting in light of the DAMA/Libra Phase-2 results, which have been noted [27,65] to be poor fits to low-mass dark matter scattering, due to the observed recoil spectrum at low energies. As seen, even relatively small admixtures of dark matter streams such as S1 can improve the quality of fit.
However, in Figure 8, along with the fits to the DAMA/Libra results, we also show shaded regions which indicate those points where the scattering cross section that corresponds to the minimum χ 2 fit to the DAMA/Libra data is itself ruled out by one of the other direct detection experiments, typically Xenon1T, hence for brevity and simplification of our results we only show the exclusion by Xenon1T. As can also be seen, for all choices of spinindependent elastic scattering, the DAMA/Libra modulation is ruled out for all possible S1 stream parameters (for cross sections up to 3σ from the best-fit point), with the exception of a very small region around m χ ∼ 30 GeV and stream fractions of more than 80%, which is marginally within the 3σ contour around the global best-fit point. Note that no equivalent 3σ region appears in the scan over mass and f ; this indicates that a non-excluded region requires some of the model parameters to be perturbed from the values that would minimize the χ 2 fit at fixed m χ and f .

B. Spin-Dependent Elastic Scattering
We now turn to spin-dependent elastic scattering. In Figure 9 we show the extrapolated limits on the proton scattering cross section σ SD p , assuming the isospin-conserving couplings a p = a n = 1/2 and the Gaia halo model without the S1 stream. We also display two isospin-violating scenarios, one where a p = −0.098 and a n = +0.902 (chosen to minimize the scattering off fluorine and thus relax the PICO-60 bound,) and the second assuming a p = −0.88 and a n = +0.12 (chosen to minimize the Xenon1T limit). 3 As with the spin-independent scattering, the DAMA/Libra regions continue to be excluded assuming a smooth halo model as extracted from the Gaia data. Assuming a p = a n = 1/2, we find a best fit point of m χ = 48 GeV, σ SD p = 3.6 × 10 −37 cm 2 and χ 2 /do.f. = 1.9. Profiling over a p , we find little improvement in the halo-only model, with the best fit occurring for a p = −0.67, but at essentially the same mass, cross-section, and χ 2 value. We can now introduce the S1 stream. As in the spin-independent case, we demonstrate the effect of the S1 stream in Figure 10, showing the extrapolated results assuming 100% of the local density comes from S1, and isospin-conserving couplings. We also show the more reasonable 25% stream density, with a p = a n = 1/2 and a p = −0.88 (chosen to minimize the xenon couplings). The minimum χ 2 /d.o.f. under these assumptions is 6.9 for the 100% stream density, and 1.4 (1.45) for the 25% density with a p = 1/2 (a p = −0.88). Thus, as with elastic spin-independent scattering, the stream by itself does not improve the fit to the observed recoil spectrum, but a reasonable admixture of the smooth halo and the stream results in a much improved goodness-of-fit parameter. As the resulting recoil spectrum and daily modulation strongly resembles that of the spin-independent scattering, we omit the spin-dependent equivalents of Figures 5 and 6.
As with spin-independent scattering, given the breadth of possible parameter space, we cannot ensure that all regions of DAMA/Libra parameter space are excluded by other null results by taking slices through a p space. We , CDMSlite (blue) and COSINE-100 (red) on the spin-dependent dark matter-proton scattering cross section σ SI p assuming the S1 stream is 100% of the local dark matter density and an isospin parameter of ap = an = +1/2 (left), 25% of the local density and ap = an = +1/2 (center), and 25% of the density and ap = −0.88 (right). The 2σ regions around the best-fit to DAMA/Libra data are shown in orange for fits to the reported modulation amplitudes (with the lighter shaded region indicating the 1σ variation of the S1 stream velocity distribution), and in yellow for fits to the reported daily modulation rates. Corresponding regions assuming the Gaia halo model with ap = an = +1/2 are shown in grey shading. again perform two scans over mass and a p space, fitting to the DAMA/Libra amplitude data. In the first scan, we assume the best fit stream fraction and σ SD p . For the second, we scan over mass and stream fraction, for the best fit a p and cross section. For each parameter point, we check to see if the required cross section for a 3σ fit to the DAMA/Libra data is ruled out by any other experiments at 90%CL. The results are shown in Figure 11. As can be seen, no region of parameter space survives this test.
It should also be noted that the modulation data alone can have good fits in regions of m χ − σ SD p parameter space which are not excluded by any null result. This occurs without the S1 stream for a small region of low-mass parameter space with a p = −0.88 (right panel, Figure 9), and becomes more pronounced as the stream is added (Figure 10). These regions are not within in the 2σ good-fit regions to the amplitude data. So, if the DAMA/Libra signal from the modulation data alone were to be interpreted as a signal of dark matter scattering through an elastic spin-dependent interaction, then the observed energy recoil spectrum as measured by DAMA/Libra must be in significant tension with the actual recoil spectrum.
For example, consider a parameter point which is a good fit to the modulation data and is not ruled out by other experiments: a p = 1/2, a stream fraction of 25%, m χ = 4 GeV and σ SD p = 4 × 10 −36 cm 2 . This parameter point results in an improvement of ∆χ 2 ∼ 24 over the null hypothesis for the modulation data. However, this point is only a χ 2 /d.o.f. ∼ 1.5 improvement over the null hypothesis in the recoil spectrum data, and would be χ 2 /d.o.f. ∼ 30 away from the best possible fit in the m χ /σ SD p plane. Therefore, while including the stream at reasonable density fractions can fit the observed DAMA/Libra modulation data while evading all other null constraints, some significant systematic errors would have to be present to reconcile the reported DAMA/Libra recoil spectrum with that caused by dark matter.

C. Inelastic Dark Matter
Having reanalyzed the elastic scattering constraints in the presence of the S1 stream, we now turn to the constraints assuming the inelastic dark matter scenario. Here, scattering proceeds through a ground dark matter state which up-scatters into an excited state, with mass difference δ. In elastic scattering, the dark matter mass sets the recoil spectrum through both the implicit dependence of v min on m χ and the explicit mass dependence in Eq. (1). Inelastic scattering on the other hand allows an additional handle by modifying the minimum dark matter velocity required for a scattering of recoil energy E R , as described in Eq. (5). However one new degree of freedom: the mass splitting δ is added to the model parameters (cross-section, mass, isospin-violating parameter f and stream density fraction) and to study the viable parameter space we must scan over these. Since we are dealing with a multi-dimensional parameter space, similar to the elastic scattering case above, we test whether the DAMA/Libra best fit cross-sections are excluded by the other experiments, when varying the extra inelastic degree of freedom. Hence in Figure 12 we fit the DAMA/Libra amplitude data, scanning over dark matter mass m χ and mass spitting δ (for both spin-independent and spin-dependent scattering). We show the 1, 2 and 3σ regions around the best fit point, portrayed in yellow, green and orange respectively. For spin-independent scattering (left panel), the best-fit point for spin-independent scattering is m χ = 230 GeV, δ = 20 keV, f = −0.79, σ SI p = 7.7 × 10 −37 cm 2 with a stream density fraction of 47% and χ 2 /d.o.f. = 0.9. There are nearly as good fits (χ 2 /d.o.f. ∼ 1) at lower mass points as well. However, for the entire range of m χ and δ, the best fit parameter points are ruled out by other experimental results, in particular Xenon1T (as shown by the grey shaded region).
For spin-dependent scattering, the best-fit point occurs when the inelastic scattering is turned off (δ = 0), and m χ = 28.5 GeV, a p = −0.15, σ SD p = 9 × 10 −34 cm 2 , and benchmark stream density of 24%. An island of similarly good fits occurs around these parameters for δ 20 keV. As seen in Figure 12, the DAMA/Libra 3σ and Xenon1T exclusion regions nearly coincide for m χ ∼ 35 GeV. This is due to the fact that, as we move away from the 3σ region of parameter space, the worsening fit to the DAMA/Libra data drives the cross section of the best possible fit lower and hence closer to the Xenon1T limit.

V. DISCUSSION AND CONCLUSIONS
The dark matter velocity distribution in the local Galactic neighborhood has long been postulated to follow a Gaussian Standard Halo Model. Astrophysical data however, points to the presence of new dark matter substructure that contributes significantly to the velocity distribution on top of the smooth background halo model. Recently, the Gaia space telescope observed a number of dark matter streams and clumps moving at high velocities near the Solar System. One of these streams, S1, is moving on a low inclination, rapidly counter-rotating orbit through the local Solar neighborhood. Due to its trajectory, the relative velocity between dark matter within S1 and the Earth is increased, enhancing the modulation rate of scattering events in experiments on the Earth.
In this work we reanalyzed the limits from current direct detection experiments, including not only the S1 stream, but also data driven local velocity distributions derived from the Gaia low metallicity stellar population. We considered various scenarios for dark matter scattering, including elastic and inelastic scattering together with spinindependent and spin-dependent dark matter-nucleus interactions in both scattering circumstances. In each case, we also considered various isospin violating couplings, which could result in the relative suppression or enhancement of scattering rates between different direct detection experiments. We directly compared the SHM, the Gaia derived velocity distributions, and the S1 stream and illustrated the effect each of these distributions has on the number of dark matter particles with sufficient kinetic energy to initiate a scattering event above threshold for a direct detection experiment. We showed that the S1 stream has a distinctive recoil spectrum which could be observed in direct detection experiments as compared to the SHM and the Gaia distribution, while maintaining a peak modulation date close to the SHM expectation.
We then studied how the existing experimental limits placed on particle physics parameters (m χ and σ) change in the presence of S1 and the Gaia distribution. To do this we first focused on the DAMA/Libra experiment which claimed to have seen a modulation signal of dark matter recoils at 13σ CL. We fit the DAMA/Libra Phase-2 amplitude and annual modulation data for different parameters in the dark matter theory and astrophysics space (e.g. interaction type, isospin ratio, stream fraction, etc). In addition to DAMA/Libra we considered those experiments which report a null signal for dark matter and place the strongest current limits on spin-dependent or spin-independent scattering (CDMSlite, Xenon1T and PICO-60). We also took into account the recent limits from the COSINE-100 experiment which use the same target material as DAMA/Libra.
As S1 increases the modulation rate and introduces new structures in the recoil spectrum, we find that even for relatively modest contributions of the stream to the local density (∼ 20−30%), the statistical fit to the DAMA/Libra data is dramatically improved for both SI and SD interactions with different isospin ratios. We find that fits including S1 prefer generally lower dark matter masses compared to the SHM or the Gaia model for the smooth halo. However, despite these changes in the best fit regions, we find that the null experiments continue to exclude the DAMA/Libra allowed region even in the presence of S1 and Gaia distribution, for our chosen benchmark parameters.
To ensure we were not missing any crucial model points by simply taking slices in parameter space, we scanned over the stream fraction and isospin factors for the best fitting DAMA/Libra dark matter-nucleon cross-section and compared with that for Xenon1T. We found that for either elastic or inelastic scattering (both SI and SD), the DAMA/Libra (modulation amplitude) best fit point continues to be excluded by Xenon1T, with the exception of spin-independent elastic isospin-violating scattering with large contributions (> 80%) to the local density from the stream. We can therefore conclude that the anomalous DAMA/Libra results continue to be excluded by the null results from other experiments unless one (or both) of the following conditions are met: 1. The S1 stream (or a collection of other streams with similar kinematic parameters) composes the vast majority of the local density of dark matter (greater than ∼ 80%). This would not be the standard expectation, but would allow for a recoil spectrum matching that of DAMA/Libra while not being excluded at the 90%CL by any other experiment.
2. The dark matter signal is present, but the dark matter recoil spectrum from the modulating signal is in significant ( 3σ) tension with the DAMA/Libra results. As we have shown, the best-fits from the yearly modulation data alone allows for low-mass fits when including the S1 stream at reasonable density that are not excluded by any other experiment. However, while the overall modulation rate matches the data, the resulting recoil spectrum would be in significant tension with DAMA/Libra.
Measuring the contribution of the S1 stream to the local density would go far in testing the first option, in addition to being extremely interesting in its own right as a probe of dark matter substructure and galaxy evolution, as well as the impact on other direct detection experiments of a non-SHM recoil spectrum. The second possibility will be tested by increased data from NaI experiments such as COSINE-100 or SABRE [66]. We note also that it is highly likely that the effects of S1 or similar streams would also affect lower mass dark matter scattering off electrons, we leave this for future considerations.
DAMA/Libra uses a sodium-iodide (NaI) crystal as the detector target. Critically, unlike other direct detection experiments, DAMA/Libra does not aim for zero-(or at least very low)-background. Rather, the experiment looks for the annual modulation of scattering events, resulting from the yearly modulation of f (v) as the Earth orbits the Sun. For many years now, the DAMA/Libra experiment has reported a positive signal, observing a yearly modulation with 12.9σ significance in the full (Phase-2) data-set [21]. The phase of the modulation matches the date expected from the Earth's motion through the Milky Way's dark matter halo, peaking at day ∼ 150 after January 1 st , depending on the range of recoil energies considered.
The DAMA/Libra Phase-2 results (corresponding to 1.13 ton×year) [21] provides the modulation data in two ways. First, they provide the modulation amplitude S m (fitting to a sinusoidal modulation) as a function of electron recoil energy from 1 to 20 keV. The signal is limited to recoil energies from 1 to 6 keV. Secondly, they provide the residual rate of scattering events as a function of time for three recoil energy ranges: 1-3 keV, 2-6 keV, and 1-6 keV. We consider fits to both slicings of the data separately. Importantly, though the modulation of a single dark matter component (i.e., the Galactic halo or the S1 stream) is typically sinusoidal, the sum of two or more components does not need to be. While the S1 stream has modulation peaks which are nearly in phase with the halo itself (indeed, this coincidence is required for S1 to possibly fit the DAMA/Libra signal), the daily modulation data provides important complimentary information which is somewhat obscured by the more finely-binned amplitude results.
To convert between the measured electron recoil energy E e and the dark matter-induced nuclear recoil energy, we convolute the response function [20,27] with the differential nuclear recoil distribution dR dE R and a detector efficiency function . Here α = (0.448±0.035) keV 1/2 and β = (9.1 ± 5.1) × 10 −3 . We use quenching factors Q I = 0.09 and Q Na = 0.3 [67].
The electron equivalent recoil distribution is given by For the purposes of this work we assume = 1. The total recoil rate in units of [1/keV/kg/day] is obtained by integrating the electron equivalent recoil distribution per energy interval per time: where the E min e and E max e are the lower and upper bounds for the energy bin under consideration, while t min and t max are the lower and upper bounds for the time interval, with ∆E e and ∆t the energy and time bin widths. The DAMA/Libra annual modulation data is usually presented in terms of a residual recoil rate defined as R −R, wherē R is the the mean rate averaged over the total number of time bins in an energy interval. We perform a χ 2 fit to the DAMA/Libra Phase-2 residual data using where O j is the observed residual rate in bin j, S j is the predicted signal residual rate (dependent on the input dark matter mass and scattering cross-section), and σ j is the uncertainty per bin.
For the modulation amplitude we assume a period of one year for the modulation, ω = 2π/(365.25 days), and extract the amplitude of this mode from Eq. (A3) per energy bin: Using a similar form as Eq. (A5) above, we perform a chi-squared fit to the DAMA/Libra phase-2 amplitude data [21] binned in 0.5 keV bins from E e = 1−4 keV, 1 keV bins from 4−7 keV, and one bin from 6−20 keV. We fit the mass m χ and dark matter-proton scattering cross section σ p , assuming a SHM velocity distribution and isospin-conserving interaction. The resulting best-fit regions as a function of m χ and σ SI p are shown in Figure 13, (spin-independent on the left and spin-dependent on the right) and broadly match the results of Ref. [27]. Both amplitude and annual modulation confidence regions are overlaid in Fig. 13 with the 1-3 keV, 1-6 keV, and 2-6 keV energy range modulation fits represented by the green, red and yellow shaded contours respectively and the amplitude data represented by the blue shaded region.

Xenon-1T
Xenon1T is a liquid xenon detector which is sensitive to the scintillation light of dark matter-nucleon interactions. With a total mass of 1.3 tonnes and 280 days of exposure, Xenon1T sets, at the time of writing, the strongest limits on dark matter-nucleon spin-independent scattering. We adopt the limits from Ref. [29] for spin-independent scattering and from [64] for spin dependent scattering. Xenon1T, like most liquid xenon-based detectors, uses two scintillation signals (S1 and S2) to reject backgrounds. We use the published efficiency curves from Refs. [29,64] for validation of our calculations. In the Xenon1T signal region we set our 90% CL upper limit (solid black) to correspond to cross sections that give 3.7 events. Our SHM validation curves are shown in Figure 14 (on the left are spin-independent (dashed black) and on the right are the spin-dependent limits for dark matter -proton (dashed blue) and -neutron (dashed black) interactions). As can be seen, at low recoil energy (corresponding to low dark matter mass) it is difficult to fully model the detector response, which is rapidly changing in this regime, resulting in recast limits that are off by O(1) from the official collaboration results.

CDMSlite
CDMSlite is part of the SuperCDMS experiment [68,69], run in a low-threshold mode to maximize sensitivity to low mass ( 10 GeV) dark matter. It uses germanium targets that measure both ionization and phonons to search for dark matter recoils, rejecting background using the ratio of ionization to phonon energy (the ionization yield).  [21]. The blue region is the best-fit to the reported modulation amplitude as a function of electron recoil energy, the green, red, and yellow regions are the fits to the daily modulation data in the [1,3], [2,6], and [1,6] keVee recoil energy bins respectively. Left: 2σ spin-independent limits and right: 2σ spin-dependent limits. All fits assume the SHM velocity distribution and isospin-conserving couplings.
FIG. 14: Left: 90% CL upper limits (solid black) on isospin-conserving spin-independent proton cross section σ SI p set by Xenon1T [29] with our validation at 90% CL shown as the black dashed line. Right: 90% CL upper limits on the isospin conserving spin-dependent proton (blue solid) and neutron (black solid) cross-sections with our validation represented as the dashed lines. The y-axis label σ SD N is intended to represent the respective nucleon. Our validations both assume the Standard Halo Model velocity distribution.
We adopt the limits set by CDMSlite in Ref. [30], using that reference's parameters for the Lindhard model for the relation between electron equivalent and nuclear recoil energies. We recast the limits from CDMSlite's Run 2, using the Ref. [30] reported energy-dependent efficiency and an exposure of 70.1 kg×day. We set our recast limits using the published background spectrum and the maximum gap method [70] in the energy range of [0.3, 1] keV ee . The comparison between our recast limits assuming the SHM and the CDMSlite results are shown in Figure 15 as the solid blue line.

PICO-60
PICO-60 is a superheated bubble chamber filled with 52.2 kg of C 3 F 8 target. Nuclear recoils that deposit energy above the 3.3 keV threshold cause a bubble nucleation. The acoustic properties of these bubbles can be used to discriminate between alpha decays and the dark matter-signal nuclear recoils. Due to the nuclear high spin of fluorine, PICO-60 provides the strongest constraints on spin-dependent, dark matter-proton interactions. We use the results of Ref. [31], corresponding to 1167 kg×days of exposure.
We use the bubble nucleation efficiency curve for fluorine as reported by PICO-60 in Ref. [71]. With an exposure of 1167 kg×days PICO-60 observed zero single bubble events, and expected a background of 0.25 events. Assuming a Poisson distribution, the 90% C.L. upper limit on the number of signal events is 2.05. Our validation of the upper limit is shown in the blue solid line with the PICO-60 reported limit shown in black. Our validation assumes the SHM velocity distribution and we show only the spin-dependent limits since they are currently the most competitive.

Cosine-100
COSINE-100 is a direct detection experiment using NaI crystal targets, recording nuclear recoil signals using the emitted light, as measured by photomultiplier tubes. This experiment, using the same target as DAMA/Libra, is an attempt to directly measure that experiment's claimed signal in an identical material. At present, the sensitivity of COSINE-100 is not sufficient to see the claimed modulation signal from DAMA/Libra, but a limit can be set on the dark matter scattering rate averaged over the detector live-time (59.5 days between Oct. 20 and Dec. 19,2016). Total target mass is 105 kg, though only 79 kg was used in the analysis [32].
We use the efficiency curve provided in Ref. [32], along with their background model in the recoil energy range 2 − 6 keV. Using the provided event spectrum, background rate, and estimated errors, we set a 90% CL limit on the scattering cross section using a simple χ 2 measure. The comparison between our estimate for the limits on cross section assuming the SHM model and the results of Ref. [32] are shown in Figure 17. Our validation is shown in blue, while the reported COSINE-100 limit is shown in black with 1 and 2σ uncertainties.