Reevaluation of spin-dependent WIMP-proton interactions as an explanation of the DAMA data

We reexamine the interpretation of the annual modulation signal observed by the DAMA experiment as due to WIMPs with a spin-dependent coupling mostly to protons. We consider both axial-vector and pseudo-scalar couplings, and elastic as well as endothermic and exothermic inelastic scattering. We conclude that the DAMA signal is in strong tension with null results of other direct detection experiments, particularly PICASSO and KIMS.


Introduction
Direct Dark Matter (DM) searches aim at detecting Weakly Interacting Massive Particles (WIMPs) scattering with nuclei in a target material. The most stringent limits on the DM mass and cross section parameter space are currently set by LUX [1] and SuperCDMS [2] for WIMPs with spin-independent interactions and spin-dependent interactions with neutrons, and by PICASSO [3], SIMPLE [4], COUPP [5], and KIMS [6] for those with protons. While none of these experiments have detected a DM signal so far, three other experiments have signals that can be interpreted as due to WIMP scattering: DAMA and DAMA/LIBRA (called here DAMA from now on) [7], CoGeNT [8,9], and CDMS-II with silicon detectors [10]. CRESST-II [11] has not confirmed a previous DM hint found by the same collaboration [12]. Among the potential DM signals, DAMA's observation of an annually modulated rate has the highest statistical significance. However, these potential signals are challenged by the null results of other direct detection experiments which exclude the possibility of WIMP scattering in a large number of particle models. In particular, the scattering cross section fitting the DAMA data for WIMPs with spin-independent interactions is several orders of magnitude above the 90% CL LUX limit [13].
Over the years, several particle candidates have been proposed with properties which enhance a potential signal in DAMA and weaken the main limits imposed by direct DM searches with negative results, among them WIMPs with spin-dependent coupling mostly to protons, which we reevaluate in this paper.
Inelastic DM [14][15][16][17][18][19][20] scatters to another particle state, either heavier (endothermic scattering) or lighter (exothermic scattering, see e.g. [21,22]), when colliding with a nucleus. Endothermic scattering favors heavier targets, thus enhancing scattering off I in DAMA while reducing scattering off lighter targets such as Ge. Moreover, this type of interaction enhances the annual modulation amplitude, thus pushing the cross section needed to fit the DAMA data to lower values. However, experiments employing Xe as target material, which is heavier than I, rule out endothermic scattering of DM as an explanation to the DAMA data unless there is an additional feature of the interaction that favors a signal in DAMA. Two types of WIMP couplings favor Na and I (DAMA) over Xe (LUX, XENON10) and Ge (CDMS, SuperCDMS): a spin-dependent coupling mostly to protons and a magnetic dipole moment coupling [19,23]. The reason for the first is that the spin of a nucleus is mostly due to an unpaired nucleon and Na and I have an unpaired proton, while Xe and Ge have an unpaired neutron. The reason for the second is the large magnetic moment of both Na and I. We do not consider here inelastic magnetic DM, which Ref. [20] recently found still marginally compatible with all negative results of direct searches, for a small value of the I quenching factor (without the aid of inelasticity, instead, this candidate has been shown to be ruled out [24]).
The possibility that a DM candidate with spin-dependent interactions mostly with protons would explain the DAMA signal was, to the best of our knowledge, first studied in Ref. [25]. Compared to spin-independent interactions, spin-dependent couplings reduce the bounds from experiments with heavy targets, most notably LUX, due to the lack of the usual A 2 T enhancement factor proper to the spin-independent interaction (A T being the mass number of the target nucleus). This interaction might explain why the DAMA signal is not seen by LUX and SuperCDMS. In this case, bounds from PICASSO, SIMPLE, COUPP, and KIMS become relevant, since they contain F, I and Cs, all nuclei with an unpaired proton. This candidate was further studied in Ref. [26], in the context of both elastic and inelastic endothermic scatterings. The inelastic endothermic kinematics reduces the expected rate in experiments employing F (PICASSO, SIMPLE) because it is light, thus making the COUPP (CF 3 I) and KIMS (CsI) bounds the most relevant constraints on WIMP scatterings off I in DAMA. Ref. [20] found that a small portion of the parameter space favored by DAMA for inelastic spin-dependent couplings with protons can still escape all bounds from null experiments.
Inelastic exothermic scattering [21,27] favors lighter targets, so it favors Na in DAMA over heavier nuclei (Ge and Xe). In this case the most important limits are set by experiments containing F (PICASSO and SIMPLE).
[28] studied a Dirac WIMP candidate coupled to standard model (SM) fermions through a light pseudo-scalar mediator, and claimed that with a contact interaction and elastic scattering it reconciles the DAMA data with the null results of other experiments at the 99% credible level. The model produces a non-standard spin-dependent interaction, with the noteworthy feature that, for universal flavor-diagonal quark couplings to the pseudoscalar mediator, the WIMP couples mainly to protons. The couplings of pseudo-scalar light bosons (m φ < 7 GeV) with quarks are strongly constrained by rare meson decays [29][30][31], and unless the pseudo-scalar coupling to the DM (called g DM below) is very large, g DM 10 3 , the scattering cross section required in Ref. [28] is rejected [31]. The flavor physics bounds on pseudo-scalar couplings to quarks proportional to the quark mass are less stringent [31], but in this case Ref. [28] founds that the resulting proton to neutron coupling ratio is not large enough to reconcile a DM signal in DAMA with the results of other direct detection experiments. Leaving aside the limits from other types of experiments we concentrate here on direct detection.
Ref. [28] employed a Bayesian analysis, where a number of uncertain parameters like quenching factors and background levels, as well as astrophysical quantities, are marginalized over. While the process of marginalization (i.e. integrating over nuisance parameters with assumed prior probability distributions) is the proper treatment of uncertain and uninteresting parameters in the context of Bayesian statistics, it makes it unclear whether there exists at least one set of values of the uncertain parameters, in particular one halo model, that produces the same result of the analysis.
In this paper we reconsider the viability of a signal due to WIMPs with spin-dependent coupling mostly to protons as an explanation of the DAMA data. We study both axial-vector and pseudo-scalar couplings, which lead respectively to s χ · s p and ( s χ · q)( s p · q) couplings in the non-relativistic limit ( s χ and s p are the spins of the WIMP χ and the proton respectively, and q is the momentum transfer). We assume the mediator to be either heavy enough for the contact interaction limit to be valid, or otherwise to be massless. The possibilities of elastic and inelastic scattering, both endothermic and exothermic, are considered.
In Section 2, we present the differential cross sections for axial-vector and pseudo-scalar couplings, which can be used in the direct detection rate formula in Section 3. The analysis methods we adopt for experimental data are described in Section 4. Our results are presented in Section 5, followed by conclusions in Section 6. couplings to nucleons, mediated by a vector boson of mass m φ , is Here q µ is the momentum transfer four-vector, and N is a nucleon, p or n. g DM and a N are the mediator coupling constants to χ and N , respectively, and they are real. The scattering amplitude is We now follow Ref. [32] because we will largely use the nuclear form factors given in this reference. We first take the non-relativistic limit of the Dirac spinors, in the chiral representation, for both χ and N : u s ( p) 1/4m (2m − p · σ) ξ s , (2m + p · σ) ξ s T , where σ are the Pauli matrices. This limit is justified by the fact that the DM initial speed and the exchanged momentum are small. The matrix element for scattering off a single nucleon then reads Eqs. (44), (47d), and (49) of Ref. [33]). Notice that this matrix element assumes the usual form for the Dirac spinors with normalizationū s ( p) u s ( p) = 2mδ ss ; Quantum Mechanical amplitudes usually assume a different state normalization, which differs by a factor of 2 m 2 + p 2 . With this normalization Eq.
For a model of inelastic DM, one could introduce two Dirac fields, χ 1 and χ 2 , with slightly different masses and the DM-nucleon effective Lagrangian The g DM coupling can now be complex. We use the same symbol, g DM , for the couplings in Eqs. (2.1) and (2.4), because then the expression of σ AV p in Eq. (2.6) is valid both for elastic and inelastic scattering. χ 1 is the DM particle entering the scattering process, with mass m, while χ 2 is the DM particle in the final state, with mass m = m + δ. The sign of the mass splitting δ determines the different kinematic regimes: δ > 0 implies endothermic scattering, δ < 0 implies exothermic scattering, while δ = 0 implies elastic scattering.
One may attempt to build an inelastic DM model without introducing additional degrees of freedom by assuming the interaction in Eq. (2.1) and adding a small Majorana mass term which produces two almost degenerate Majorana fermions, χ 1 and χ 2 , in which case χ = χ 1 + iχ 2 becomes a quasi-Dirac fermion. However, as noted in Ref. [26], this interaction produces diagonal terms which result in elastic scattering rather than inelastic. In this case, one can instead write an effective tensor interactionχσ µν χN σ µν N , which produces inelastic scattering since the diagonal interaction terms vanish identically. The non-relativistic limit of this operator is also s χ · s N .
The differential cross section for DM-nucleus scattering, for both the elastic and inelastic interactions introduced above in Eqs. (2.1) and (2.4), is where E R = q 2 /2m T is the nuclear recoil energy, v is the incoming WIMP speed, and µ p is the DM-proton reduced mass. F 2 AV (q 2 ) is a nuclear form factor including spin dependence, and will be defined in Section 2.3. σ AV p is the total DM-proton cross section in the limit of The term in parenthesis in Eq. (2.5) accounts for long-range interactions, when m φ is smaller than or comparable to q = √ 2m T E R . For typical target masses of a few tens of GeV and recoil energies around few to tens of keV, the interaction becomes effectively long-range if the mediator mass is smaller than several MeV: m φ q 20 MeV (E R /10 keV)(m T /20 GeV). In order to plot our results for long-range interactions, we express the differential cross section in terms of a reference total cross section σ AV, ref corresponding to a reference mediator mass m ref φ , which we set equal to 1 GeV: The massless mediator limit thus corresponds to setting m φ = 0 in the equation above. In the following we will refer to any scenario with m 2 φ q 2 as massless mediator limit or long-range limit.
Given that a large value of a p /a n is needed to suppress the strong LUX and SuperCDMS constraints, we will assume the maximally isospin-violating coupling a n = 0 for the AV interaction.

Pseudo-scalar (PS) interaction
Here the DM particle is a Dirac fermion χ, coupled to a real PS boson φ with mass m φ , (as in Refs. [28,34]), with a real coupling constant g DM . The PS field couples also to the SM quarks with real coupling g q , While PS couplings to quarks are usually taken to be proportional to the fermion mass (see e.g. [34]), we will assume instead a flavor-universal coupling g q = g, which introduces a larger |a p /a n | ratio, a p /a n −16.4 [28] (see below).
To model inelastic scattering we assume a non-diagonal coupling of two Dirac DM fields χ 1 and χ 2 with φ, Again, the g DM coupling can now be complex. With this definition of g DM , Eqs. (2.8) and (2.10) yield the same expression for σ PS p in Eq. (2.16). Diagonal interaction terms as well as non-diagonal mass terms can be forbidden by assuming a Z 2 symmetry under which both φ, χ 1 (or χ 2 ), and all the SM electroweak doublets have charge −1.
The DM-nucleon effective Lagrangian, for elastic scattering (for inelastic scattering we should haveχ 2 γ 5 χ 1 instead ofχγ 5 χ), is a Nχ γ 5 χN γ 5 N , (2.11) and yields in the non-relativistic limit where s χ and s N were defined after Eq. (2.3). This is a different type of spin-dependent interaction than in Eq. (2.3). Due to the extra factors of q, the PS cross section receives a large q 4 /m 2 N m 2 suppression with respect to the AV cross section. Therefore, the normalization of the signal and its spectrum, and also the nuclear form factors are different in the two cases [32] (see Section 2.3). Given the large momentum suppression, one needs to check the existence of unsuppressed radiative corrections to this tree-level cross section, that would spoil the setup. The PS interaction in Eq. (2.11) with a Dirac fermion χ has been proven not to produce such corrections [35], while this would not be the case if χ were a Majorana fermion.
The proton and neutron couplings appearing in Eq. (2.11) are given by factors parametrize the quark spin content of the nucleon, and are usually determined experimentally or computed with lattice calculations. As in Ref. [28], we adopt the following values from [36]: with which a p −0.4g. As a natural feature of this model, the proton coupling a p is larger (in modulus) than the neutron coupling a n , by an amount that depends on the choice of the ∆ (N ) f 's. As noted in [28], the values in Eq. (2.14) are conservative in the sense that they minimize the ratio a p /a n with respect to other values encountered in the literature (see e.g. Table 4 in [33]). In this case a p /a n = −16.4.
The differential cross section for the PS interaction, for both elastic and inelastic scattering, is , to be defined in Section 2.3, is the nuclear form factor including spin dependence. σ PS p is the total DM-proton cross section in the limit of contact interaction, (2. 16) In this case, the total DM-proton cross section has a v 4 dependence, and, for the purpose of plotting our results in terms of the reference cross section, in Eq. (2.16) we evaluate σ PS p at a reference speed v ref . We set v ref equal to the rotational speed of our Local Standard of Rest, 220 km/s, which is representative of the WIMP speeds with respect to Earth.
For long-range PS interactions we proceed in the same manner as for the long-range AV interactions, by writing the differential cross section in terms of a reference total cross section σ PS, ref

Nuclear form factors
We adopt the form factors computed in Ref. [32] using standard shell model techniques, for the nuclides for which they are available, namely the main stable isotopes of Ge, Xe, Na, I, and F. In these cases, we define for the AV interaction, and for the PS interaction. The (squared) nuclear form factors F Σ and F Σ are tabulated in [32] for the nuclides mentioned above. These form factors can be employed unmodified also for inelastic scattering [20]. We include a factor of 1/3 in the definition of F 2 AV in order to normalize F 2 AV and F 2 PS to be 1 in the limit of zero momentum transfer when the target is an isolated proton. This factor traces back to F Σ being twice as large as F Σ at q = 0, which is consistent with the fact that F Σ corresponds to the component of the nucleon spin along the direction of the momentum transfer, while F Σ corresponds to the transverse component. F 2 AV (q 2 ) can be expressed (see Eqs. (59), (60) and (77c) in Ref. [32]) in terms of the usual nuclear spin structure function S(q 2 ) = a 2 0 S 00 (q 2 ) + a 0 a 1 S 01 (q 2 ) + a 2 1 S 11 (q 2 ) [37] (with a 0 = a p + a n and a 1 = a p − a n the isoscalar and isovector parameters): with J T the spin of the target nucleus. At zero momentum transfer and where S z p is the component of S p ≡ protons s p along the z-axis [38] ( S n is defined analogously). Notice that S p and S n are often denoted with boldface style in the literature, although they are not vector quantities. F 2 AV can then be expressed in terms of the usually called spin-dependent form factor F 2 SD (q 2 ) = S(q 2 )/S(0) as For the nuclides for which no form factors have been computed in Ref. [32] (Cl, C and Cs), we define F 2 AV (q 2 ) by means of Eq. (2.22), with the spin-dependent form factor in Gaussian form with A T the mass number of the target nucleus [39]. In this case we also assume F 2 PS = F 2 AV , which we expect to be approximately valid at low q 2 . For Cs, a component of KIMS's target material, we take S p = −0.370, S n = 0.003 [40,41]. For SIMPLE, we use S p = −0.051, S n = −0.0088 for both 35 Cl and 37 Cl [38], and S p = −0.026, S n = −0.155 for 13 C [42]. Notice that there are large uncertainties in the hadronic matrix elements S p and S n and the nuclear form factors, which differ in different nuclear models (see e.g. Fig. 1 of Ref. [43] for the Xe nuclear structure functions). Factors of 2 difference in different calculations are not uncommon [44].

Direct Detection Rate
The DM-nucleus scattering rate for a target nuclide T with mass m T is where m is the WIMP mass, ρ is the local DM density, f ( v, t) is the DM velocity distribution in Earth's reference frame, and dσ T /dE R is the differential cross section for a WIMP scattering off the target T . v min (E R ) is the minimum WIMP speed needed to impart to the target nucleus a recoil energy E R . The detectors do not measure directly the recoil energy, but they measure a related energy E (sometimes expressed in 'keV electron equivalent' or keVee, or else in number of photoelectrons). The rate in Eq. (3.1) is related to the event rate measured by an experiment within a detected energy interval [E 1 , E 2 ] as where C T is the mass fraction of the nuclide T in the detector, G T (E R , E ) is the targetdependent resolution function of the detector and (E , E R ) is the experimental acceptance.
The resolution function is defined as the probability distribution for a nuclear recoil energy E R to be measured as E , and incorporates the mean value E = Q T (E R )E R , with Q T the target's quenching factor, and the detector energy resolution. The detector acceptance is a function of both E R and E in general, but it is often given to be E R independent, or the dependencies on E R and E are factorizable.
In this paper we assume the standard halo model (SHM) for the dark halo of our galaxy, where the DM local density is ρ = 0.3 GeV/cm 3 and the velocity distribution of WIMPs in the galactic frame is a truncated Maxwell-Boltzmann distribution Here v 0 is the velocity dispersion and v esc is the escape speed from our galaxy. The normalization factor We consider v 0 to be the same as the velocity of the Local Standard of Rest v 0 = 220 km/s, and we take v esc = 533 km/s, according to recent Radial Velocity Experiment (RAVE) 2013 results [45]. The velocity distribution in Earth's frame, f ( v, t) in Eq. (3.1), can be obtained with the Galilean transformation where v and v ⊕ (t) are the velocity of the Sun with respect to the galaxy and the time dependent velocity of Earth with respect to the Sun, respectively. We take v = 232 km/s, and v ⊕ = 30 km/s in an orbit inclined at 60 • with respect to the galactic plane [46]. Earth's revolution around the Sun causes the velocity distribution given in Eq. (3.5), and therefore causes the scattering rate in Eq. (3.1) to modulate in time. Both the timeaverage rate and the modulation amplitude of the rate can be measured. In the SHM, the velocity of Earth with respect to the Galaxy is maximum at the end of May or the beginning of June. The rate has a maximum at this moment if v min > 200 km/s and it has instead a minimum at this moment if v min < 200 km/s. Thus, choosing the modulation phase so that the modulation amplitude is positive for v min > 200 km/s, the amplitude is negative for v min < 200 km/s. Only v min > 200 km/s are compatible with the phase of the DAMA modulation data such that the rate has the maximum on the 2nd of June [7]. The value of v min = 200 km/s is shown as a vertical line in Fig. 1.
The minimum speed the DM particle must have in the rest frame of the target nuclide in order to impart a nuclear recoil energy where µ T is the DM-nucleus reduced mass. The mass splitting δ = m − m can be either positive for endothermic scattering [14], negative for exothermic scattering [21,27,47], or zero for elastic scattering. By inverting this equation one obtains the minimum and maximum recoil energies E ± R (v) that are kinematically allowed for a fixed DM speed v, (3.7) The minimum possible value of v for the interaction to be kinematically allowed is v δ = 2δ/µ T for endothermic scattering, and v δ = 0 for exothermic. This speed value corresponds to the point of intersection of the two E ± R branches, which occurs at The maximum value of the DM speed allowed for a given halo model is the sum of the escape speed v esc and the modulus of Earth's velocity in the galactic rest frame, which is equal on average to the Sun's velocity v ; therefore, for our choice of parameter values, v max = 765 km/s. This is indicated with a vertical line in Fig. 1. As a result, endothermic scattering with a target nucleus will be kinematically forbidden for δ larger than 3.3 keV(µ T /GeV). The left panel in Fig. 1 shows the two branches for the Na component of DAMA, for δ = −30 and −50 keV, for WIMP masses that correspond to the best fit regions we will present later. In these ranges for m and δ, endothermic scattering with Na in DAMA is kinematically forbidden. On the other hand, scattering off I in DAMA is allowed for elastic or endothermic scattering, as can be seen in the right panel of Fig. 1 for δ = 50 and 100 keV.

Data Analysis
In this section we examine the compatibility of the WIMP interpretation of the DAMA annual modulation signal with various null results for the AV and PS models described above. We consider elastic, endothermic, and exothermic scattering. We start by describing our data analysis, which follows the procedure already presented in [24,[48][49][50].
For the DAMA annual modulation signal, we take the data plotted in Fig. 8 of Ref. [7]. We determine the DAMA favored regions in the DM parameter space by performing a Maximum Likelihood analysis, assuming the data are Gaussian distributed. We consider scattering off Na and I separately. Due to the uncertainties residing in the quenching factors of Na and I, which play an important role in the analysis, we choose two values for each target, namely Q Na = 0.4 and 0.3 for Na, and Q I = 0.09 and 0.06 [51] for I (see e.g. Ref. [52] and references therein).
We also compute an upper limit on the WIMP cross section using the total rate measured by DAMA, employing the data points plotted in Fig. 1 of Ref. [53]. We restrict our analysis to energies above the experimental threshold of 2 keVee. Given the very large number of observed events in each bin, and the resulting small statistical fluctuations, we compute an upper bound on the cross section by requiring that the predicted rate does not exceed the observed rate in any energy bin. This limit is particularly important for exothermic scattering, which reduces the modulation amplitude with respect to the average rate.
To compute the LUX bound, following Ref. [50], we apply the Maximum Gap method [54] to the variable S 1 in the range 2-30 photoelectrons. We choose several numbers of observed events, i.e. 0, 1, 3, 5 and 24 as described in Ref. [50]. In our SHM model the maximum WIMP speed is v max = 765 km/s, thus the maximum recoil energy for a WIMP lighter than 11.5 GeV in an elastic scattering with Xe is ∼ 12 keV. Using the approximated recoil energy contours in Fig. 4 of Ref. [1], and dropping all observed events in and above the electron-recoil band (plotted at 1.28σ) in the same figure, only five observed events remain below ∼ 12 keV. This means that for elastic scattering, choosing 5 events provides a safe upper limit for WIMP masses m < 11.5 GeV. However, if m > 11.5 GeV, only using all the 24 events lying outside the electron-recoil band provides a reliable upper limit. For inelastic scattering with δ = −50, −30, 50 and 100 keV, the maximum WIMP masses for which using the 5 events bound is reliable are 8.2, 9.3, 19.2 and 56.2 GeV, respectively. Since our procedure does not depend on the WIMP distribution in the S 1 − log 10 (S 2 /S 1 ) plane [1], our Maximum Gap upper limits are conservative and safe to be applied to any WIMP-nucleus interactions.
For the SuperCDMS bounds, we use the data set collected by the seven Ge detectors between October 2012 and June 2013, corresponding to an effective exposure of 577 kg-days. We use the Maximum Gap method with the eleven observed events listed in Table 1 of Ref. [2], which passed the selection criteria introduced by the collaboration to discriminate signal from background events. We do not incorporate the uncertainty in recoil energy in our analysis. For the detector acceptance we take the red curve in Fig. 1 of Ref. [2].
To compute the SIMPLE limits, we consider only the Stage 2 [4], a C 2 ClF 5 detector with an exposure of 6.71 kg-day, with one observed event above 8 keV compatible with the expected background of 2.2 events. We use the Feldman-Cousins method [55] to place a 90% CL upper limit of 2.39 signal events for 2.2 expected background events and 1 observed event.
For PICASSO we perform a Maximum Likelihood analysis using the data in Fig. 5 of Ref. [3]. The target material in PICASSO is C 4 F 10 , but the collaboration only considers scattering off F in their analysis [3]; we do the same, noting that the contribution of C for DM spin-dependent interactions with protons is anyway negligible. We construct our Gaussian likelihood using the expected rate above each one of the eight energy thresholds adopted by the collaboration (1.7, 2.9, 4.1, 5.8, 6.9, 16.3, 38.8, and 54.8 keV), and the measured rate with its uncertainty, which are already background subtracted.
For KIMS we perform again a Maximum Likelihood analysis using the data points with their 68% CL intervals from Fig. 4 of Ref. [6], assuming Gaussian distributed data. Because Cs and I have similar atomic masses, their quenching factors are not measured separately [56]. As for DAMA, we perform our analysis of the KIMS data adopting two values for Q I = Q Cs in CsI: 0.05 and 0.1 (see Fig. 2

Results
The plots in Figs. 2-4, 6-8, and 11-13 show 90% CL upper bounds and 68% and 90% CL allowed regions in the m − σ p plane. The regions shaded in green and labeled 'DAMA 1 ' are the allowed regions compatible with the DAMA annual modulation, for quenching factors Q Na = 0.4 and Q I = 0.09 in dark green, and Q Na = 0.3 and Q I = 0.06 in light green. The lower the quenching factor, the higher is the DM mass needed to fit the data. The low and high WIMP mass regions correspond to the interpretation of the DAMA data as the WIMP scattering off Na and I in the detector, respectively. The upper limit due to DAMA total rate (black, and labeled 'DAMA 0 Na') is shown for scattering off Na assuming Q Na = 0.4. 90% CL upper limits from LUX data are shown as various magenta curves. As in Ref. [50] the different dashing styles of the lines indicate different selections of candidate events used in the Maximum Gap analysis: dotted (0 events), double-dot-dashed (1 event), dot-dashed (3 events), dashed (5 events) and solid (24 events) curves. Two purple lines show the 90% CL upper limits from KIMS data with quenching factors Q I = Q Cs = 0.1 (solid) and 0.05 (dashed). 90% CL upper limit from SIMPLE (brown), PICASSO (cyan), and SuperCDMS (dark yellow) are also drawn.  Figure 2. 90% CL bounds and 68% and 90% CL allowed regions in the WIMP-proton reference cross section σ p vs WIMP mass plane, assuming the SHM, for elastic proton-only contact AV interactions. The unmodulated DAMA rate limit (black) corresponds to a Na quenching factor of 0.4.  Fig. 3 is the same as Fig. 2, but for PS interactions with flavor-universal coupling a n /a p = −1/16.4 (left panel), and with proton-only coupling a n = 0 (right panel). As expected, the only limits that change from one case to the other are those of LUX and SuperCDMS, due to their enhanced sensitivity to DM-neutron couplings. The DAMA region for WIMP scattering off Na is entirely excluded by SIMPLE and PICASSO, and the region for scattering off I is excluded by KIMS when assuming similar values for the I quenching factor in both experiments. This result is different from what was found in Ref. [28], where some portion of the Na and I DAMA regions are compatible with all null experiments for the PS flavor-universal coupling.

Elastic contact interactions
Ref. [28] uses Bayesian statistics to infer 99% credible level exclusion limits, and 90% and 99% credible regions for DAMA, marginalizing over the SHM parameters using Gaussian priors (taking central values for the velocities v 0 = 230 km/s and v esc = 544 km/s, and for the local WIMP density ρ = 0.3 GeV/cm 3 , with standard deviations ∆v 0 = 24.4 km/s, ∆v esc = 39 km/s and ∆ρ = 0.13 GeV/cm 3 ). As a result, regions and limits at a specific point in parameter space do not necessarily correspond to a fixed set of values for the SHM parameters. In our analysis instead we assumed the same set of SHM parameter values across all experimental results. We found that the regions and limits move approximately in the same manner in the parameter space as we vary the DM velocities, and the DAMA regions fail to escape the upper bounds at the 90% CL. This can be seen in Fig. 4 (left panel), which shows the results for elastic PS interactions with flavor-universal coupling where both v esc and v 0 are taken 3σ below their central values in Ref. [28] (however, we keep v ref = 220 km/s to plot σ p ). Our choice for the SHM velocities roughly matches the low mass SIMPLE limit in Fig. 1 of Ref. [28]. The right panel of Fig. 4 shows also the 99% CL upper bounds (dotted lines) and 99% CL regions for the same set of parameters. In this case, the I component of the DAMA region corresponding to a quenching factor of 0.9 escapes the KIMS upper limit for quenching factor 0.5. The Na component of the DAMA region is still rejected by PICASSO at the 99% CL. In Fig. 5 we present our results from the left panels of Figs. 3 and 4 in the same plane as Fig. 1 of Ref. [28], namely in the m − Λ φ plane, where Λ φ ≡ m φ / √ g DM g. The allowed DAMA regions shown in Ref. [28] are much larger than the regions we found. We believe that this is due to their marginalization over the SHM parameters and experimental parameters including quenching factors.
Ref. [31] found for PS interactions that the LUX bound excludes the DAMA region (see Fig. 9 in Ref. [31], where the I region in DAMA is completely excluded by LUX in both the contact and long-range limits). This is in disagreement with our conclusions, possibly because of the different analysis of the LUX data.  Fig. 2 and the right panel of Fig. 6). This is expected from the E R dependence of the differential cross sections given in Eqs. (2.5) and (2.17): disregarding the form factors, the differential cross section for the long-range PS and contact AV interactions is proportional to E 0 R , for contact PS it is proportional to E 2 R , and for long-range AV it is proportional to E −2 R . As it can be seen in Fig. 6, considering long-range elastic interactions does not help to bring compatibility between the DAMA regions and the upper limits from the experiments with null results.

Exothermic contact interactions
Figs. 7 and 8 show the results for exothermic inelastic proton-only AV and PS interactions, respectively, with δ = −30 keV (left panels) and δ = −50 keV (right panels). As |δ| increases, the DAMA regions move to lower masses, for the following reason. The lowest reach in DM mass for a direct detection experiment is obtained when E + R (v max ) = E th (see Fig. 1), where the threshold energy E th is the lowest detectable recoil energy. The mass reach can be found by extracting m as a function of δ and v from and evaluating it at v = v max . The DAMA region will be therefore located at WIMP masses higher thanm(δ, v max ), taking Na as the target element. For DAMA, E th = 5 keV for scattering off Na with quenching factor 0.4. The lowest reach of DAMA is the lower green line plotted in Fig. 9. Also shown in the same figure are the mass reaches of PICASSO and SIMPLE, for which we used E th = 1.7 keV and 8 keV, respectively. We only considered scattering off F in both experiments.
An estimated upper limit on the DM mass for the DAMA region comes from requiring that DM particles with speeds below 200 km/s always scatter below threshold and are therefore undetectable, because otherwise DAMA should have observed a sign change in the modulation amplitude. In other words, scatterings of DM particles slower than 200 km/s would yield a different phase for the modulated signal with respect to that measured by DAMA, and therefore an acceptable fit requires these scatterings to occur below threshold. Since a fixed nuclear recoil energy can be imparted by heavier DM particles traveling at lower velocities, the condition v min (E th ) > 200 km/s implies an upper limit on the DM mass in DAMA, given bym(δ, v = 200 km/s). This upper limit is the higher green line plotted in Fig. 9. The other possible condition to avoid scatterings of DM particles slower than 200 km/s, i.e. having a large enough E − R (v = 200 km/s), would imply a very odd spectrum in DAMA, with more events at higher energy instead of the observed spectrum vanishing at high energy.
Since exothermic scattering decreases the value of v min for a given recoil energy, the modulation amplitude becomes smaller with respect to the time-average rate. For large enough |δ| the DAMA modulation signal becomes inconsistent with the DAMA time-average rate. For values of δ lower than about −30 keV (for AV interactions) and −50 keV (for PS interactions), the DAMA total rate limit rules out the modulation signal in Na, as indicated by the black curve excluding the DAMA Na region in the right panels of Figs. 7 and 8.
For the values of δ allowed by the DAMA rate, the limit by PICASSO (and also the SIMPLE limit in most instances) rejects the allowed regions. For each value of δ on the horizontal axis of Fig. 9, the DAMA region spans a mass range enclosed within the green belt, while the SIMPLE and PICASSO lines indicate the mass value where the limits in the m − σ p plane become vertical. From the plot it becomes clear that exothermic scattering brings compatibility between SIMPLE and DAMA for large enough |δ|, as suggested by Figs. 7 and 8, however the region is rejected by the DAMA average rate measurements.  While this is true for Q Na = 0.4, smaller quenching factors move the DAMA region to larger DM masses, thus potentially compromising this compatibility with SIMPLE. In any event, the DAMA region does not escape the PICASSO limit.

Endothermic contact interactions
Figs. 11 and 12 show the result for the proton-only spin dependent endothermic scattering with AV and PS interactions, respectively, in the contact limit, with δ = 50 keV (left panels) and δ = 100 keV (right panels). As δ increases, scattering off light targets becomes kinematically forbidden since v δ becomes larger than v max . For δ = 100 keV, the only remaining limits are from KIMS and LUX. We can see that the DAMA region for I scattering moves towards the left compared to the KIMS upper bound as δ increases. For PS interactions with δ = 50 keV, it is only the combination of larger quenching factor Q I = 0.09 for I in DAMA and smaller quenching factor Q I = 0.05 in KIMS that allows the DAMA signal to be compatible with all upper limits. For δ = 100 keV, the I region corresponding to Q I = 0.06 barely escapes the limit with Q I = 0.05 from KIMS, and the situation remains tense for quenching factors 0.09 (DAMA) and 0.1 (KIMS). Raising δ further makes it progressively more difficult to find a region of the DAMA signal that is kinematically allowed. For AV interactions the DAMA regions are even more severely constrained: only the larger quenching factor Q I = 0.09 DAMA region is allowed by the KIMS upper bound with smaller quenching factor Q I = 0.05 for δ = 100 keV.
These results are largely consistent with those of Ref. [20], where the framework of non-relativistic operators introduced in Ref. [32] was generalized to inelastic scattering and a model-independent analysis was performed on a series of effective operators. In Table V of Ref. [20] the AV and PS interactions correspond to fermion operators 15 and 4, respectively. For those interactions, Ref. [20] quotes best fit parameters for DAMA (corresponding to δ = 106 keV for AV and δ = 57 keV for PS) that are consistent with the KIMS data only for DAMA quenching factor 0.09 and KIMS quenching factor 0.05. In Ref. [20], however, scattering off Cs in KIMS was neglected due to the lack of the form factor in Ref. [32]. Since the contribution of Cs to the scattering rate is sizable for interaction with protons, we adopt an approximated form factor for Cs as discussed in Sec. 2.3, resulting in stronger KIMS bounds.
At this point it is important to recall the flavor physics bounds we mentioned in the introduction. Fig. 9 of Ref. [31] shows that the quark couplings needed for PS inelastic scattering on I to fit the DAMA data, where κ 0.2 and κ 0.1 for δ = 50 keV, and κ 1.6 and κ 0.7 for δ = 100 keV, is rejected [31] for any reasonable value of g DM (unless g DM > 10 5 ). Note that the quark coupling used in Ref. [31] is g/ √ 2.

Endothermic long-range interactions
For the AV interaction via a massless mediator shown in Fig. 13, only a very small portion of the I DAMA region with quenching factor 0.09 escapes the KIMS limit with quenching factor 0.05. Therefore, quite different quenching factors are needed for I in DAMA and KIMS to have the DAMA region escape the KIMS limit. Although quenching factors of a given element in different crystals can have different values in general, large differences may be questionable. Again, we do not plot the results for long-range PS interactions since these are qualitatively similar to those for contact AV interactions, as commented above.

Conclusions
We investigated the possibility of interpreting the annual modulation signal observed in the DAMA experiment as due to WIMPs with spin-dependent coupling mostly to protons. We considered both an axial-vector (AV) interaction, which is what is usually referred to as 'spin-dependent interaction', and a pseudo-scalar (PS) interaction, proposed in Ref. [28] to reconcile DAMA with the null experiments. We also extended our analysis to inelastic scattering, and considered both contact and long-range interactions. Due to the similar E R dependence of the differential cross sections, we find for the long-range PS interaction the same results as for contact AV interactions, up to a shift in σ p .
Spin-dependent WIMP couplings mostly to protons effectively weaken the bounds from experiments using Xe and Ge as target elements, whose spin is due mostly to neutrons. However, the bounds from experiments with F and I targets, such as PICASSO, SIMPLE and KIMS, remain relevant since their spin is due mostly to protons. For elastic scattering (see Figs. 2 to 6) we found that, in all the cases investigated here, the DAMA regions for Na are entirely excluded by SIMPLE and PICASSO, while the regions for I are excluded by KIMS.
For exothermic scattering (see Figs. 7 and 8), the DAMA regions move progressively to smaller WIMP masses with respect to the upper limits, because the modulation phase observed by DAMA forces v min > 200 km/s, and this is possible only for progressively lighter WIMPs as |δ| increases (see Fig. 9). Thus, exothermic scattering brings compatibility between the DAMA region for Na and the upper bounds from SIMPLE. However, it does not suppress the PICASSO limit which continues to rule out the DAMA region. Furthermore, exothermic scattering reduces the modulation amplitude with respect to the time-average rate. Thus the upper limit derived from the DAMA average rate measurement rejects the interpretation of the signal as due to scattering off Na for values of δ < −30 keV. The DAMA region for scattering off I is excluded by the SIMPLE and KIMS upper bounds.
For endothermic scattering (see Figs. 11 to 13), only KIMS provides relevant bounds. Scattering in all detectors besides KIMS and LUX becomes kinematically forbidden for large enough δ. We showed results for δ = 50 and 100 keV, because as δ increases further, scattering off I becomes kinematically forbidden as well. For δ = 50 keV, only assuming a larger quenching factor Q I = 0.09 for I in DAMA, and a smaller quenching factor Q I = 0.05 in KIMS, the DAMA region is compatible with all present limits for PS couplings, although the possibility that the same nuclide has such different quenching factors in different crystals is questionable. The same holds for contact and long-range AV and long-range PS interactions for δ = 100 keV. For contact PS interactions and δ = 100 keV, a small sleeve of the 90% CL DAMA region for scattering off I escapes the 90% CL KIMS limit with similar Q I for both experiments. These results are largely consistent with the results of Ref. [20]. However, for PS interactions the DAMA regions are rejected by flavor physics bounds on the PS coupling to quarks [31] (unless g DM can be very large g DM > 10 5 ).