Halo-independent tests of dark matter annual modulation signals

I derive new halo-independent lower bounds on the product of the dark matter-nucleon scattering cross section and the local dark matter density that are valid for annual modulations of dark matter direct detection signals. They are obtained by making use of halo-independent bounds based on an expansion of the rate on the Earth's velocity that were derived in previous works. In combination with astrophysical measurements of the local energy density, an observed annual modulation implies a lower bound on the cross section that is independent of the velocity distribution and that must be fulfilled by any particle physics model. In order to illustrate the power of the bounds we apply them to DAMA/LIBRA data and obtain quite strong results when compared to the standard halo model predictions. We also extend the bounds to the case of multi-target detectors.

Recently, in ref. [40] a new framework has been established that allows to compare the constant rate of a DD signal with local energy density measurements, LHC limits, the relic abundance and indirect detection constraints (see also ref. [39] for another method to derive a lower bound on the cross section). In this work we show how the observation of an annual modulation signal in a DD detection experiment can also be used to place a lower bound on the product of the local DM density and the scattering cross section that is independent of the DM velocity distribution, in a similar fashion as was done in ref. [40] for constant rates. For this purpose, we use halo-independent bounds on the annual modulation based on an expansion of the rate on the Earth's velocity [24,25]. In order to illustrate their use we perform a numerical analysis to DAMA/LIBRA data, which we will denote by DAMA in the following. As we will see, in combination with measurements of the local energy density we obtain quite strong lower bounds on the cross section, that must be fulfilled by any particle physics model in order for the signal to be consistent with DM.
The paper is structured as follows. In section 2 we set the notation for DD and the annual modulation. In section 3 we review the case of DAMA modulation in the context of the SHM and we perform a fit to the signal. This section can be skipped by the readers familiar with the DAMA signal. Section 4 encodes the main results of this paper, where we derive lower bounds on the product of the local DM density and the scattering cross section that are independent of the DM velocity distribution and are valid for annually modulated signals. We do this in two ways: in section 4.1 by making use of the fact that constant rates need to be larger than modulations, and in section 4.2 by using an expansion of the rate on the Earth's velocity [24,25]. As an illustration we apply our halo-independent bounds to the DAMA signal in section 4.3. Finally, we give our conclusions in section 5. 1

Dark matter direct detection and the annual modulation signal
In this section we review the relevant expressions for DD. We focus on elastic scattering of DM particles χ off a nucleus with mass number A depositing a nuclear recoil energy E R . The differential rate for a detector with n different target nuclei is given by: where ρ χ is the local DM mass density and v A m is the minimal velocity of the DM particle required for a recoil energy E R , with µ χA being the reduced mass of the DM-nucleus system. f A is the mass fraction of the detector nucleus labeled by A. For singletarget detectors, there is just one contribution and thus the sum over A is absent. f det ( v, t) describes the DM velocity distribution in the detector rest frame, with the normalization d 3 v f det ( v, t) = 1. The velocity distributions in the rest frames of the detector, the Sun and the galaxy are related by where v e (t) is the velocity vector of the Earth relative to the Sun and v s is the velocity of the Sun relative to the galactic halo.
We will concentrate on elastic spin-independent (SI) and spin-dependent (SD) contact interactions, where the differential scattering cross section dσ A (v)/dE R scales as 1/v 2 , and we will take equal couplings of the DM to neutrons and protons. Then the SI cross section becomes 2 where σ SI is the total DM-proton scattering cross section at zero momentum transfer, µ χp is the DM-proton reduced mass, and F A (E R ) is a nuclear form factor. For SD interactions there is no A 2 enhancement, the form factor is different, and σ SD will denote the zero-momentum DM-proton scattering cross section. Using eqs. 2.1 and 2.2 the event rate can be written in a very compact way as where we defined the detector-independent quantity Due to the time-dependence of the Earth's velocity and its orientation with respect to the Sun's velocity DD signals exhibit annual modulations [2,3]. The time dependence in the event rate is introduced through where we have defined the constant contribution as 6) and the modulated one comes from an expansion to first order in v e = 29.8 km/s v sun 230 km/s: , where the phase t * (E R ) is the time of year at which the event rate is maximum, and thus A η (v A m ) > 0. Notice that observing a sign flip in the modulation amplitude would be a strong evidence that the signal is due to DM. However, the phase and the energy at which the flip may occur depends on the velocity distribution. Moreover, there could also be several sign flips at different energies. For the SHM, t SHM * = 0.42 y, corresponding to June 2nd, and there is a sign flip at v A m ≈ 210 km/s, see for instance ref. [12]. Similarly, the amplitude of the modulated rate is: which is always positive. For a specific detector the number of DM induced events in an energy range [E 1 , E 2 ] is given by where M and T are the detector mass and exposure time, respectively, and we define which is an integration over energy of the weighted sum of the target nuclei for the quantity is the detector response function describing the probability that a DM event with true recoil energy E R is reconstructed in the energy interval [E 1 , E 2 ], including efficiencies, energy resolution, and quenching factors.
3 The DAMA annual modulation signal within the SHM In section 4 we will derive halo-independent bounds valid for annual modulation signals, and in order to illustrate their use, we will apply them to the modulation observed by DAMA [41]. In the following we review and analyze the interpretation of this modulation within the SHM, and therefore this section can be skipped by those readers familiar with this signal.
The DAMA experiment reports an annual modulation of the rate in their NaI scintillator detector, with a period of one year and a maximum around June 2nd. With an exposure of 1.33 ton y the modulation is claimed to be compatible with DM at a very high statistical significance (9.3σ) [41]. The amplitude of the annual modulation observed in the [2,6] keVee range is M [2,6] = (0.0112 ± 0.0012) counts/keVee/kg/day [41], while above this range the data is consistent with no annual modulation. As it is well known, the signal can be explained under standard assumptions regarding the DM halo (SHM) and WIMP interactions. However, the DM interpretation of the modulation is strongly disfavored by other experiments independently of the halo, both for elastic SI and SD interactions [25], and for inelastic interactions [28]. Despite this, we will use the DAMA signal in the following to illustrate the use of the lower bounds derived in section 4. Notice that if DAMA were not due to DM, current limits (for instance by LUX [42]) imply that observing annual modulations will be very hard in the future (see for instance ref. [43]).
First let us discuss which is the dark mass range consistent with the data for the SHM. There is no observation of a sign flip and all the observed modulations in the [2,6] keVee range are positive, which means that for the SHM they are above the sign flip, implying that v m 210 km/s. Therefore, in order for the modulation of the first bin to be positive, if the scattering is on sodium, from eq. (2.1) the DM mass needs to be smaller than ∼ 30 GeV, while if the scattering is on iodine, it needs to be smaller than ∼ 90 GeV. In addition, there are also lower bounds on the DM masses in order for the DM particles to be energetic enough to deposit recoil energies above threshold. This means that if the scattering is on iodine, DM masses below ∼ 20 GeV cannot produce observable recoils. On the contrary, for m χ 20 GeV and SI interactions, iodine typically dominates the rates due to the A 2 enhancement. 3 Furthermore, by demanding that there are observable recoils in the 6 keVee range, we obtain that m χ 8 (30) GeV for scattering on Na (I).
This is illustrated in figure 1, where we plot the minimum velocity v m (E R ), eq. (2.1), versus the DM mass for sodium (red) and iodine (blue) for recoil energies relevant for the DAMA experiment: E R = 2, 4 , 6 keVee (from bottom to top as solid, dotted and dashed curves). We also show as dotted black the minimum velocity below which the modulation flips sign for the SHM, and as a dashed black line the typical escape velocity in the detector rest frame. The shaded area encoded between the dotted vertical red (blue) lines shows the allowed mass range 8 m χ (GeV) 30 (30 m χ (GeV) 90) for scattering on sodium (iodine) within the SHM. For more complicated haloes with streams or debris flows the phase of the modulation will be different (see for instance ref. [12]), and so will the allowed mass range.
In order to analyze more quantitatively DAMA's modulation within the SHM, we perform a χ 2 fit to the data in the [2,6] keVee energy range, both for SD and SI elastic interactions. For our numerical analysis we take equal couplings to protons and neutrons, both for SI and SD, 4 and quenching factors q Na = 0.3 for sodium and 3 For this reason lowering the threshold of the DAMA experiment would in principle allow to distinguish between the two SI solutions, the ∼ 10 GeV DM mass, for which a lower threshold would imply an enhanced rate due to scatterings on iodine, and the ∼ 80 GeV one, which would show a phase flip in the lowest energy bin, see ref. [44]. 4 Regarding SD, both 23 Na and 127 I are dominated by the spin of the protons. Therefore, suppressed couplings to neutrons help to reduce the inconsistency with other null-result experiments that have  Table 1. Results of the DAMA fit to the SHM, done for SI and SD interactions, assuming equal couplings to protons and neutrons. We use v esc = 550 km/s, ρ χ = 0.4 GeV cm −3 and the quenching factors q Na = 0.3 for sodium and q I = 0.09 for iodine.
q I = 0.09 for iodine. For the SI form factor we use the Helm parameterization, For SD we take the structure functions from ref. [46]. For the SHM we use a Maxwellian velocity distribution with mean velocityv = 220 km/s truncated at the galactic escape velocity v esc = 550 km/s, and a local energy density ρ χ = 0.4 GeV cm −3 , see ref. [47] for a recent review.
We use the latest results from figure 8 of ref. [41]. As expected from the literature, there are two minima: for the large DM mass solution the scattering is predominantly on iodine, while for the low DM mass one it is on sodium. The best-fit values are shown in table 1. In all cases we obtain a good fit, with χ 2 min /dof < 1.5, where dof stands for the number of degrees of freedom, which is 6 (8 data points minus 2 parameters). Notice that for SI interactions, for the large DM mass solution the scattering is predominantly on iodine and therefore the cross section is much lower than in the case of low DM masses (scattering on sodium), due to the larger A 2 enhancement of the former. For SD both solutions imply roughly the same cross sections. 5 Our results agree with those present in the literature, see for instance refs. [45,48,49].
In figure 2, we show the allowed 90% CL parameter space for the SHM as solid blue lines, together with the ∆χ 2 = 2.3, 5.99, 9.21 contours (in blue, green, light blue respectively) corresponding in two dimensions to a CL of 68.27, 95, 99% respectively. This is done for SI (top) and SD (bottom), and we show the low (high) mass solution in the left (right), which correspond to scattering on Na (predominantly on I). We also plot the best-fit points with blue marks.

Using that the constant rate is larger than the annual modulation
An upper bound on the halo integral η(v m ) 6 can be easily obtained using the normalization condition of the velocity distribution as well as the definition of η(v m ), eq. (2.6), as was done in refs. [32,36,40]: where from the first to the second line we used that below the threshold v 1 the minimum is obtained for constant η(v) = η(v 1 ), as it is a monotonously decreasing function, and in the last line we dropped the second term. From these expressions an upper bound on the constant rate can be derived [40]. In order to have a bound valid for annual modulation signals, one can simply use that the constant rate is larger than the amplitude of the modulation: If the spectrum is not measured, one can derive a bound in terms of the measured number of modulated events in an energy interval [E 1 , E 2 ] by combining eq. (2.9) with eqs. (4.3) and (4.4). Using the definition of C, eq. (2.4), one gets: Notice that this bound can be used for multitarget experiments.
If the spectrum is measured, one can obtain a stronger bound using eq. (4.2) instead of eq. (4.3). In this case, for a single-target experiment with perfect energy resolution, for fixed DM mass, a measurement of the spectrum M(E R ) allows a determination of the halo integral via eqs. (2.3) and (2.8): .
Then, in combination with eq. (4.2), the following bound can be obtained: If a DD experiment observes an annual modulation spectrum, eq. (4.7) provides a lower bound on the product ρ χ σ SI which is independent of the local DM velocity distribution. In ref. [40] bounds valid for constant rates analogous to those of eqs. (4.5) and (4.7) were derived. Notice that if a sign flip is observed, a better approach is to apply the bound for the region below and above the sign flip separately, in order to obtain stronger bounds.
If the detector has different elements, in order to use eq. (4.7) one needs to assume that one of them gives the dominant contribution. The bound will be stronger in general for the lowest mass nuclei in the case of SI, due to the A 2 suppression, but so will the SHM prediction, and therefore we expect the ratio of the bound and the SHM to remain approximately the same for different nuclei. One also needs to take into account that, although for SI interactions the heaviest element typically dominates the rate due to the A 2 enhancement, for low DM masses the scattering on the heavy nuclei may not be possible, as there are no DM particles with speeds larger than the escape velocity in the detector rest frame, ∼ 750 km s −1 . We will illustrate these features for the DAMA data in section 4.3. In appendix C we generalize eq. (4.7) to the case of multi-target detectors.

Using upper bounds on the annual modulation
One a can derive stronger tests for the case that the spectrum is measured by using the results of refs. [24,25], where upper bounds on the annual modulation signal in terms of the constant rate were obtained by doing a first order expansion of the rate on the Earth's velocity.
Lets consider a spectral measurement of both the modulation M(E R ) and the constant rate R(E R ) in the energy range [E 1 , E 2 ], which for a given DM mass can be related to a velocity interval [v 1 , v 2 ] via eq. (2.1). By doing an expansion on the Earth's velocity and some mild assumptions about the halo (see below) one can derive the following upper bounds on the modulation amplitude [24,25]: For the "General halo" it is assumed that the velocity distribution is constant in time on the scale of 1 year (i.e., that all the time dependence of the rate comes from the velocity of the Earth) and that it is constant in space on the scale of the size of the Sun-Earth distance. Notice that halo substructures like streams with velocities larger than v e are covered by the bounds. The "Symmetric halo" further assumes that there is a preferred direction for the DM flow, specified by α, the angle between the DM flow and the orthogonal to the ecliptic. The most conservative bound is obtained for sin α = 1, which corresponds to a DM stream parallel to the ecliptic. However, it may happen that the DM velocity is aligned with the motion of the Sun, and therefore sin α 0.5. This is the case of isotropic velocity distributions and, up to a small correction due to the peculiar velocity of the Sun, it also holds for triaxial halos or a dark-disc. For the "General halo" the bound is independent of the phase, which is free. However, for the "Symmetric halo" the phase is independent of v m (and therefore independent of E R ), and fixed. The SHM is included in this last case, with sin α 0.5 and the phase being equal to June 2nd. As for the case of the "Spectrum bound", eq. (4.7), if a sign flip in the modulation is observed in the velocity range considered, one may apply the bound for the region below and above the sign flip separately, in order to get stronger bounds. More details regarding the different assumptions of the bounds can be found in the original ref. [24] (and also in refs. [25,28]).
From the above discussion it is clear that in order to decide whether to apply the "General halo" or the "Symmetric halo" one can use the information encoded in the phase. First of all, one needs to do a fit keeping as free parameters the modulation amplitude, the period and the phase. If data disfavours a constant phase across the different energy bins, this invalidates the assumption of a "Symmetric halo" and one is forced to use the "General halo". Regarding the polar angle α, unfortunately, the phase of the modulation does not carry information about it, and different values of α will give the same phase. It only depends on the azimuthal angle within the ecliptic, φ, on which however the bounds do not depend. This is easy to visualize if one imagines a stream parallel to the ecliptic, with α = π/2. In this case it is clear that the phase is completely determined by the angle φ at which the Earth's velocity and the stream direction are aligned. Notice that the amplitude of the modulation is proportional to sin α, but, unfortunately, its value cannot be disentangled without specifying the halo.
In order to check the consistency of a DD signal one can apply these bounds if both the modulation and the constant rate have been measured [24]. In some cases one may not be able to separate the background from the DM contribution in the measured constant rate, as is the case of DAMA, but it is clear that if the inequalities are violated for an unknown constant background plus DM rate, they will also be for just the DM component, and the conclusion (the exclusion of the DM as an interpretation of the modulation) will still hold. However, if the bounds are fulfilled for some unknown constant background plus DM rate, they may well be violated for just the DM rate. In these cases, in order to check if the modulation observed in one experiment is consistent, one can compare it with upper limits on the rate from a different experiment [25,28]. Interestingly, as we will show now, one can obtain upper bounds that do not depend in any way on the constant rate nor on its upper limits. In combination with other measurements, these set halo-independent constraints on the DM models [40].
By applying eq. (4.2) in the "General halo", eq. (4.8), we can obtain an upper bound on the modulation (see appendix A for the derivation): Similarly, for the "Symmetric halo", eq. (4.9), we obtain: Both bounds are written completely in terms of velocities, with no factors containing the constant rate. Therefore, even if just the annual modulations are measured, one can get an upper bound. This is crucial, since, as we argued, some experiments are not able to disentangle constant DM signals from constant backgrounds.
If we now consider a spectral measurement of the modulation such that we can extract A η via eq. (4.6), then, using eq. (4.10), we can get an upper bound on the integrated modulation over the velocity range: from which we can derive a lower bound in the product of the energy density and the cross-section which does not depend on f (v), and is valid for very generic haloes: "General bound". (4.13) An analogous bound can be obtained for symmetric haloes, eq. (4.11): 7 "Symmetric bound". (4.14) If a DD experiment observes an annual modulation, eq. (4.5) ("Events bound"), eq. (4.7) ("Spectrum bound"), eq. (4.13) ("General bound") and eq. (4.14) ("Symmetric bound") provide lower bounds on the product ρ χ σ SI which are independent of the local DM velocity distribution. The first two just use that modulations are smaller than constant rates, while the last two are based on an expansion of the rate on the Earth's velocity, and are clearly a factor ∼ v/v e (with v > v m ) stronger than the former. These are the main results of this paper. Notice if the detector has different nuclei, in order to use eqs. (4.7), (4.13) and (4.14), one needs to assume that one of them gives the dominant contribution. In appendix C we generalize these bounds to the case of multi-target detectors.

Applying the halo-independent bounds: the case of DAMA
Now we will apply the halo-independent bounds derived in the previous section to the annual modulation observed by DAMA [41] (see section 3 for a review and analysis of DAMA in the context of the SHM). As discussed in section 4.2, in order to know if one can apply the "Symmetric bound" or only the "General bound" one should perform a fit to the modulation data leaving the amplitude, the period and the phase free in order to determine whether the phase is constant or not across the energy range of the signal. For DAMA these fits have been performed [41], obtaining that the phase is compatible with being constant across different energy bins within 3σ, and equal to June 2nd (see table 4 of ref. [41]). It should be emphasized that the DAMA modulation is compatible with the SHM (and thus with sin α 0.5), see section 3, but it can also be caused for instance, as discussed in section 4.2, by a stream with a different polar angle, sin α = 0.5, but with the same azimuthal angle φ as in the SHM (where it is given by the Sun's velocity) such that the phase is June 2nd, and with a stream speed such that the final modulation amplitude has the same value as the one observed by DAMA. Therefore, for DAMA one can use the "Symmetric halo" for different values of α. In this case, the most conservative option is sin α = 1, and if one is interested in analysing a DM flow aligned with the velocity of the Sun one can use sin α 0.5. In order to illustrate their different strengths we will also show the "General halo", which does not require a constant phase. 8 Notice that the differences in the modulation amplitudes and their errors in the case of fixed (free) phase, c.f. table 3 and 4 of ref. [41], are negligible, and therefore we can use the same data when applying both bounds. However, in general this may not be the case, see for instance the CoGENT analysis of ref. [24], which uses the extracted modulation amplitudes for the different assumptions of ref. [50].
In order to apply the bounds to real data we need binned expressions, which are given explicitly in appendix B. To extract A η (v m ) from the observed modulation amplitude M(E R ) via eq. (4.6), and therefore to apply our bounds, one needs to assume that the scattering occurs on a particular nucleus. We will assume in the following that the scattering is dominated by either sodium or iodine, which, as we explained before, is a very good approximation, given the large mass hierarchy between them. In appendix C we generalize the bounds to multi-target experiments and apply them to both nuclei at the same time.
In figure 3 we show the 90% CL lower bounds on the scattering cross section, together with the SHM preferred regions, already shown in figure 2. For definiteness, we use the mass range allowed for the SHM (plotted in blue), see section 3 for details. From bottom to top, we plot the "Spectrum bound", eq. (4.7), which just uses that A η ≤ η, as solid red; we also show the bounds that use the expansion on the Earth's velocity: the "General bound", eq. (4.13), as dashed red, and the "Symmetric bound", eq. (4.14), as dotted (dotted-dashed) red for sin α = 1 (sin α = 0.5, i.e., for the DM flow aligned with the velocity of the Sun). As expected, in all cases the latter bounds are much stronger than the "Spectrum bound" (solid red), due the v/v e enhancement (c.f. eqs. (4.13) and (4.14) to eq. (4.7)).
One can use these lower bounds on the cross section to compare with LHC limits, with relic abundance constraints or with indirect detection limits [40]. This has to be done within the context of a particular particle physics model and can be used to check its validity. Another useful and general comparison can be done by combining the lower bound on ρ χ for a given DM mass with local density measurements [40]. This combination provides a lower bound on the cross section of any given particle physics model. In the analysis, we use values for ρ χ in the range between 0.2 and 0.6 GeV/cm 3 (see ref. [47] for a recent review), and we fix the DM mass to be equal to the best-fit values obtained in the χ 2 fit, see table 1. In this case one would need to obtain the DM mass by other means, for instance from LHC measurements, from an indirect detection signal (i.e., a gamma ray line) or from more than one DD signal, see for instance refs. [19,37]. Notice that results for different DM masses can be obtained by a simple rescaling of figure 3. This is shown in figure 4, where we plot the local energy density versus the scattering cross section. We show in blue the 90% CL SHM predictions and in red the 90% CL lower bounds, for SI (top) and SD (bottom), and for scattering on Na (left) and I (right). We plot from left to right the "Spectrum bound" of eq. (4.7) (solid red), the "General bound" of eq. (4.13) (dashed red) and the "Symmetric bound" of eq. (4.14) (dotted red for sin α = 1, dotted-dashed red for sin α = 0.5). As expected, the last three give much stronger constraints than the former.
For the very conservative "General bound" (dashed red) we obtain that SI cross sections σ SI 8 · 10 −42 (3 · 10 −43 ) cm 2 are disfavoured for m χ = 12 GeV, scattering on Na (m χ = 79 GeV, scattering on I). For SD σ SD 3·10 −38 (2×10 −38 ) cm 2 , for m χ = 12 GeV, scattering on Na (m χ = 63 GeV, scattering on I). A bit stronger constraints are obtained for the "Symmetric bound", eq. (4.14) (dotted red), for sin α = 1 (if the DM flow is aligned with the velocity of the Sun, i.e., sin α = 0.5, the limits, shown with a Figure 3. Results for the DAMA data for SI (top) and SD (bottom), and scattering on Na (left) and I (right). The SHM allowed parameter space at the 90% CL is encoded between the solid blue lines. The 90% CL lower bounds are shown in red from bottom to top: the "Spectrum bound" of eq. (4.7) (solid red), the "General bound" of eq. (4.13) (dashed red) and the "Symmetric bound" of eq. (4.14) (dotted red for sin α = 1, dotted-dashed red for sin α = 0.5. Notice that the latter is a factor of 2 stronger than the former.). The SHM ∆χ 2 = 2.3, 5.99, 9.21 contours corresponding in two dimensions to a CL of 68.27, 95, 99% are shown in blue, green and light blue, respectively, together with the best-fit points, depicted with blue marks.
dotted-dashed red curve, are a factor of 2 stronger than in this last case). Compared to the SHM, these are roughly between ∼ 0.5 and ∼ 1.5 orders of magnitude weaker.
Notice that these cross sections are excluded independently of the unknown velocity distribution. In addition, the limits are conservative in the sense that they still hold if the particle that gives rise to the DD signal constitutes only part of the total DM, since then the lower bound on the total local DM density can only be larger.

Conclusions
In this paper we have derived new halo-independent lower bounds on the product of the DM-nucleon scattering cross section and the local DM density that are valid for annually-modulated signals, extending the work done in ref. [40] for constant rates. First, we have used that the amplitude of the annual modulation should be smaller than the rate. This allows to obtain the "Events bound", eq. (4.5), if only the number Figure 4. Lower bounds on the local energy density applied to DAMA modulation SI (top) and SD (bottom), and for scattering on Na (left) and I (right). The DM masses in each case are taken to be close to the best-fit points. We show the SHM prediction (solid blue), together with the bounds in red, from left to right: the "Spectrum bound" of eq. (4.7) (solid red), the "General bound" of eq. (4.13) (dashed red) and the "Symmetric bound" of eq. (4.14) (dotted red for sin α = 1, dotted-dashed red for sin α = 0.5. Notice that the latter is a factor of 2 stronger than the former.). All are shown at 90% CL.
of events is measured, and the "Spectrum bound", eq. (4.7), if the spectrum is observed. Then we have derived new tests by making use of halo-independent bounds obtained in refs. [24,25] by doing an expansion of the rate on the Earth's velocity. These are the main results of the paper: the "General bound", eq. (4.13), valid for general haloes, and the "Symmetric bound", eq. (4.14), which assumes some preferred direction for the DM velocity. These bounds can also be extended to multi-target detectors, see appendix C.
In order to illustrate their use, we have applied them to DAMA data. In combination with local energy measurements, we are able to exclude cross sections that are roughly between one and one and a half orders of magnitude smaller than the SHM ones, see figure 4, but with the reward of being independent of the unknown velocity distribution. It is important to stress that these limits must be fulfilled by any particle physics model that gives rise to elastic SI/SD interactions in order for the DAMA signal to be consistent with DM.
In the spirit of being as conservative as possible, we would like to encourage the community to show, for positive DM direct detection signals, in addition to the typical SHM best-fit regions, the halo-independent lower bounds derived in this work, for instance the "Symmetric bound" of eq. (4.14) (dotted red) for the cases where the phase is constant, or the "General bound" of eq. (4.13) (dashed red) when the phase varies, in a similar fashion as was done in figures 3 and 5 for the DAMA signal.
Acknowledgements: I am very grateful to Thomas Schwetz for fruitful discussions regarding not only this work but also our previous related projects. I would also like to thank Mattias Blennow and Stefan Vogl for useful discussions on our recent related project.

A Derivation of the upper bounds on the modulation
We derive here eq. (4.10) valid for the "General halo". By applying eq. (4.2) in the "General halo", eq. (4.8), we get: Integrating the third term, and by parts the last term, we obtain for eq. (A.1): Now we can use that the last term of eq. (A.2) obeys and inserting eq. (A.3) back into eq. (A.2), we get: The last term of eq. (A.4) depends on η(v) and we can conservatively drop it. We finally obtain: which is an upper bound on the modulation amplitude expressed solely in terms of velocities. In the end the result is the same as if we had used in eq. (4.8) the weaker bound of eq. (4.3) instead of the stronger one of eq. (4.2). An analogous derivation can be performed for eq. (4.11).
Notice that strictly speaking all these inequalities should read "less than" (<), instead of "less or equal than" (≤), because: 1) eq. (4.2) (and similarly eq. (4.3)) is very conservative, assuming just that η is a non-increasing function of velocity and, 2) eq. (4.8), and therefore eq. (4.10) (and similarly eq. (4.9) and (4.11) for symmetric haloes), are based on a first order expansion on the Earth's velocity. Moreover, in the last step we dropped a positive term. 9 However, we decide to keep "less or equal than" in all inequalities to make explicit that for very special (and probably unrealistic) haloes the inequalities could be (at least approximately) saturated, see ref. [24] for examples. In practice this means that any modulation signal that even saturates the bounds is very unlikely have a DM origin.

B Binned versions of the bounds
In order to apply the bounds to real data, we need binned expressions, which necessarily involve some approximations. We will discuss the validity of the binning prescriptions here used at the end of the section. Let us define the bin average of a quantity X(E) as where E i1 and E i2 are the boundaries of bin i and ∆E i = E i2 − E i1 . We denote by M i ≡ M i the observed modulation in bin i, where i = 1, ..., N . For a fixed DM mass, A η can be obtained from the amplitude of the modulation in bin i, i.e., M i , by: where q is a possible quenching factor which we assume to be energy-independent (as is the case of DAMA), v i ≡ v m (E i ) is the minimum velocity centered in the corresponding energy bin using eq. (2.1), F 2 A i is the nucleus form factor averaged over the bin width and f A is the nucleus mass fraction if different elements or isotopes contribute to the rate. Similarly, η i ≡ η(v i ) can be expressed in terms of the constant rate R i . Now we have to perform a bin average of the bounds. This involves quantities like g(v) M(E R ) i , which we replace by g(v) i M i , which is a good approximation whenever g(v) does not vary drastically within a bin. The integral is approximated by a sum over bins. Using eq. (B.2), the "Spectrum bound", eq. (4.7), becomes: where we denote by ∆v i the width in velocity space of bin i. Similarly, for the "General bound", eq. (4.13), we get: And finally the "Symmetric bound", eq. (4.14), reads in its binned version: The inevitable inaccuracies involved in the binning procedure are small whenever the quantities within a bin do not show drastic changes. By using different prescriptions for the bounds, for instance 1/v 1 , vs 1/v 1 , vs 1/ v 1 , we have checked that inaccuracies due to the binning procedure are negligible. Notice also that the averaging of the velocity v can be performed analytically, as the relation with recoil energy E R is known, c.f., eq. (2.1). In addition we need to extract information from the measured modulation amplitude in a particular bin (of size 0.5 keVee for DAMA), which necessarily implies unfolding form factors and the detector response, which involves energy resolution, efficiencies and quenching factors. We have also checked that by including energy resolution in a similar fashion as in ref. [44] (see also refs. [51,52]) the quantities involved in our bounds do not vary much, and therefore energy resolution is not expected to change the results significantly.

C Bounds for multi-target detectors
With the exception of the "Events bound", eq. (4.5), in order to apply the rest of the bounds here derived to multi-target experiments one needs to assume that scattering is dominated by a particular nucleus. Let us now generalize the bounds to experiments with different n elements, labeled by A. For a similar procedure for multi-target bounds see ref. [24]. In this case the modulation amplitude in an energy bin i receives contributions from each element: M i = n A=1 M A i . We use the notation v A i ≡ v A m (E i ), see eq. (2.1). Summing eq. (4.2) in its binned version for all the n elements we get: i is the total rate in bin i, and min A (X i ) means that we need to take the minimum X between all the nuclei present in the target, for each bin i. From the first to the second line we used the equivalent of eq. (B.2) for η i (in terms of constant rates), and from the second to the third line we changed the order of the sums and used the fact that C is detector independent, see eq. (2.4). Using that constant rates are larger than modulations, R i ≥ M i , and the definition of C, eq. (2.4), we get: which is the "Spectrum bound" valid for multi-target experiments. Now we can do a similar thing for the bounds based on an expansion on the Earth's velocity, the "General bound", eq. (4.13), and the "Symmetric bound", eq. (4.14). Summing eq. (4.10) in its binned version for all the n elements we get: (C.5) By changing the order of the sums, we finally obtain: (C.6) which is the "General bound", eq. (4.13), for multi-target experiments. Doing the same procedure in eq. (4.11), we get for the "Symmetric bound", eq. (4.14): which is now valid for multi-target experiments. It is easy to check how these multi-target bounds reduce to the single-target ones for n = 1, as it has to be. Notice however that the bounds here derived are weak due to the crude approximations needed in going from eq. (C.2) to eq. (C.3). For instance, for two isotopes of the same element, the multi-target bounds should be as strong as the mono-target ones, but due to the simplifications made they are clearly weaker by more than a factor of 2.
In figure 5 we show in black from bottom to top the "Spectrum bound" of eq. (C.4) (solid), the "General bound" of eq. (C.6) (dashed) and the "Symmetric bound" of eq. (C.7) (dotted for sin α = 1, dotted-dashed for sin α = 0.5) for DAMA data for SI (top) and SD (bottom), using both sodium and iodine. The "Events bound", eq. (4.5), which is valid for multi-target, is plotted as solid purple. We also show the SHM prediction in blue (dotted for scattering on Na, dashed on I, and solid on both). The SHM ∆χ 2 = 2.3, 5.99, 9.21 contours (already shown in figure 2) corresponding in two dimensions to a CL of 68.27, 95, 99% are shown in blue, green and light blue, respectively, together with the best-fit values, depicted with blue marks. As expected, in all cases the multi-target bounds are weaker than the ones obtained for the singletarget bounds, cf. figure 3. In particular, for SI interactions in the low DM mass region (m χ 20 GeV) the difference between the multi-target bounds and the SHM is greater than two orders of magnitude, as the former ones are suppressed by the iodine A 2 factor with respect to the latter, which are dominated by scatterings on sodium. Notice that the "Events bound", eq. (4.5) (solid purple), is stronger than the "Spectrum bound" for multi-target detectors (solid black) due to the 1/n suppression of the latter (n = 2 in this case), see eq. (C.4). shown at 90% CL. The 90% CL SHM range is plotted in blue (dotted for scattering on Na, dashed on I, and solid on both) together with the "Spectrum bound" of eq. (C.4) (solid), the "General bound" of eq. (C.6) (dashed) and the "Symmetric bound" of eq. (C.7) (dotted for sin α = 1, dotted-dashed for sin α = 0.5), using both Na and I. The "Events bound" of eq. (4.5), which is valid for multi-target, is plotted as solid purple. The SHM ∆χ 2 = 2.3, 5.99, 9.21 contours corresponding in two dimensions to a CL of 68.27, 95, 99% are shown in blue, green and light blue, respectively, together with the best-fit points, depicted with blue marks.