Probing DAMA/LIBRA in the full parameter space of WIMP effective models of inelastic scattering

We discuss the compatibility of the combined annual modulation effect measured by DAMA/LIBRA-phase1 and DAMA/LIBRA-phase2 with an explanation in terms of inelastic scattering events induced by the most general Galilean-invariant effective contact interaction of a Weakly Interacting Massive Particle (WIMP) dark matter particle of spin 0, 1/2 or 1. We take into account all the possible interferences among operators by studying the intersections among the ellipsoidal surfaces of constant signal of DAMA and other experiments in the space of the coupling constants of the effective theory. In our analysis we assume a standard Maxwellian velocity distribution in the Galaxy. We find that, compared to the elastic case, inelastic scattering partially relieves but does not eliminate the existing tension between the DAMA effect and the constraints from the null results of other experiments. Such tension is very large in all the parameter space with the exception of a small region for WIMP mass $m_{\chi}\simeq$ 10 GeV and mass splitting $\delta\gtrsim$ 20 keV, where it is partially, but not completely relieved. In such region the bounds from fluorine targets are evaded in a kinematic way because the minimal WIMP incoming speed required to trigger upscatters off fluorine exceeds the maximal WIMP velocity in the Galaxy, or is very close to it. As a consequence, we also find that the residual tension between DAMA and other results is more sensitive on the astrophysical parameters compared to the elastic case. We find that the configurations with the smallest tension can produce enough yearly modulation in some of the DAMA bins in compliance with the constraints from other experiments, but the ensuing shape of the modulation spectrum is too steep compared to the measured one. For such configurations the recent COSINE-100 bound is evaded in a natural way due to their large expected modulation fractions.


I. INTRODUCTION
For more than 15 years the DAMA collaboration [1][2][3][4] has been measuring a yearly modulation effect in a large-mass low-background sodium iodide target compatible to the signal of the Dark Matter (DM) particles that are believed to make up 27% of the total mass density of the Universe [5] and more than 90% of the halo of our Galaxy. Indeed, Weakly Interacting Massive Particles (WIMPs), which provide one of the most popular DM explanations, are expected to have feeble interactions with nuclear targets in a terrestrial detector with a scattering rate that presents a modulation with a period of one year due to the Earth revolution around the Sun [6]. In particular, with the release of the latest DAMA/LIBRA-phase2 data [4] the statistical significance of DAMA effect has reached almost 12σ. However, in the most popular WIMP scenarios used to explain the DAMA signal as due to DM WIMPs, the DAMA modulation appears incompatible with the results from many other DM experiments that have failed to observe any signal so far. Nevertheless, until recently none of the experiments ruling out the DAMA effect used the same target nuclei as DAMA/LIBRA, so that such incompatibility relied on both Particle-Physics and Astrophysics assumptions. Such model dependence has been shown to persist [7] * Electronic address: francis735@naver.com † Electronic address: scopel@sogang.ac.kr ‡ Electronic address: tomar@sogang.ac.kr also after the bound from the COSINE-100 collaboration [8], that has recently published an exclusion plot for a standard elastic, spin-independent isoscalar WIMP nucleus interaction and a WIMP Maxwellian velocity distribution that for the first time rules out the DAMA effect at low WIMP masses using 106 kg of NaI, the same target of DAMA. Given the strong statistical significance of the DAMA/LIBRA signal, and the scientific implications, this prompted the need to extend the class of WIMP models. Indeed, several scenarios have been introduced trying to reconcile the DAMA effect with other null results [7,9]. A more systematic approach is to compare DAMA and other null results exploiting the non-relativistic (NR) nature of the WIMP-scattering process, that allows to express the interaction in terms of the more general effective Hamiltonian allowed by Galilean invariance [10,11], of the form: where t 0 = 1 is the identity in isospin space, t 1 = τ 3 is the third Pauli matrix, and r denotes the dark matter-nucleon relative distance and the operators O j depend on the exchanged momentum q, the WIMP incoming velocity v, the WIMP spin J χ and nuclear spin J N . In Eq. (1) the isoscalar and isovector (dimension -2) coupling constants c 0 j and c 1 j , are related to those to protons and neutrons c p j and c n j by c p j = (c 0 j + c 1 j ) and c n j = (c 0 j − c 1 j ). Truncating the effective theory expansion to operators at most quadratic in q ≡ | q| and v ≡ | v| the ensuing Hamiltonian contains 8 independent couplings c τ j for a scalar WIMP, 28 independent couplings for a spin-1/2 WIMP [10,11] and 20 couplings for spin-1 [12]. The ultimate assessment of the compatibility of DAMA with other constraints requires a full exploration of such large parameter space, where the relative sensitivities of different nuclear targets to DM scattering events may vary by order of magnitudes. While, due to the large dimensionality, its direct scanning appears to be challenging this has been achieved by using matricial techniques [13], exploiting the fact that in terms of the couplings vector c = c (0) 1 , c (1) 1 , ..., c (0) n , c (1) n T for all direct detection experiments the expected event rate can be written in the form: with is a real symmetric 2n × 2n matrix, which encodes all the information about nuclear responses, the dark matter velocity distribution, experimental efficiencies, etc., but which is independent of the underlying particle physics model (for a given dark matter mass). Due to this factorization in the EFT parameter space the surfaces of constant signal in different detectors are ellipsoids, and, as discussed in Refs. [13,14], the determination of their geometrical intersections allows to efficiently compare the results of various direct detection experiments in the high-dimensional parameter space of the nonrelativistic effective theory, without making any a priori assumptions regarding the relative size of the various Wilson coefficients c τ k . In this way for a standard halo model the DAMA/LIBRA-phase1 result [1][2][3] was shown to be incompatible to the constraints from other experiments in the case of elastic interactions of a WIMP particle with spin ≤ 1/2 and the interaction Hamiltonian of Eq.( 1) with arbitrary couplings combinations [13].
As pointed out by the authors, the analysis in [13] did not cover several alternative scenarios. One of them is inelastic scattering. Indeed, one of the very few scenarios that reconcile DAMA with the constraints of other experiments is proton-philic Spin-dependent Inelastic Dark Matter (pSIDM) [7,15,16]. In such model the WIMP particle interacts with nuclear targets trough a spin-dependent coupling that is suppressed on neutrons, in order to comply to constraints using neutron-odd scattering targets (germanium and xenon). Moreover inelastic scattering (IDM) [17] reconciles the above scenario to fluorine detectors. In IDM a DM particle χ 1 of mass m χ 1 = m χ interacts with atomic nuclei exclusively by up-scattering to a second heavier state χ 2 with mass m χ 2 = m χ + δ. A peculiar feature of IDM is that there is a minimal WIMP incoming speed in the lab frame matching the kinematic threshold for inelastic upscatters and given by: with µ χN the WIMP-nucleus reduced mass. This quantity corresponds to the lower bound of the minimal velocity v min (also defined in the lab frame) required to deposit a given recoil energy E R in the detector: with m N the nuclear mass.
In particular, indicating with v * Na min and v * F min the values of v * min for sodium and fluorine, and with v esc the WIMP escape velocity, in Ref. [7] constraints from WIMP-fluorine scattering events in droplet detectors and bubble chambers were shown to be evaded when the WIMP mass m χ and the mass gap δ are chosen in such a way that the hierarchy: is achieved. In fact, in such case WIMP scatterings off fluorine turn kinematically forbidden while those off sodium can still serve as an explanation to the DAMA effect. So the pSIDM mechanism rests on the trivial observation that the velocity v * min for fluorine is larger than that for sodium.
In the present paper we wish to apply the technique introduced in Ref. [13] to extend the analyses of [7,15,16] to the general interaction Hamiltonian of Eq. (1), or, conversely, we wish to extend the analysis of [13] to the case of inelastic scattering, updating it to the present experimental situation and including the DAMA/LIBRA-phase2 data release and the state of the art of all the constraints from other experiments.
The two most important improvements of DAMA/LIBRA-phase2 compared to the previous phases are that now the exposure has almost doubled and the energy threshold has been lowered from 2 keV electron-equivalent (keVee) to 1 keVee. In particular this latter feature introduces an important difference between the present analysis and that of Ref. [13]. In the latter, for kinematic reasons and irrespective of the effective interaction the DAMA/LIBRA-phase1 data were only sensitive to scattering events off Sodium for a WIMP mass m χ < ∼ 20 GeV, implying that in such range of mass the scaling law induced by the effective Hamiltonian (1) only entered in the comparison of the scattering rate off Sodium and that off the targets of other experiments. However, due to the lower threshold, now for m χ < ∼ 20 GeV DAMA/LIBRA-phase2 is sensitive to both the target nuclei, with WIMP-Iodine scattering events contributing to the expected rate in the new low-energy range below 2 keVee and Sodium at higher energy. This implies that the scaling law among different targets is now also relevant in explaining the energetic spectrum of the modulation amplitudes measured by DAMA alone. Indeed, due to this reason for a standard Maxellian velocity distribution the the goodness-of-fit of a WIMP explanation of the DAMA/LIBRA-phase2 data has already been showed to worsen compared to DAMA/LIBRA-phase1 for a standard Spin-Independent interaction (SI) [18,19], requiring to tune the ratio between the WIMP-proton and the WIMP-neutron couplings in order to suppress WIMP-Iodine scattering events below 2 keVee. As shown in [20], with the exception of O 1 and O 4 , which in the notation of [10,11] correspond respectively to the standard SI or Spin-Dependent (SD) interactions, this problem is not present for all the other operators O j of the Hamiltonian (Eq. 1). This represents an additional motivation to update the result of Ref. [13] to the DAMA/LIBRA-phase2 data. Finally, the recent COSINE-100 bound is potentially relevant to our analysis because any probe of the DAMA effect using NaI is expected to lead to conclusions independent on the WIMP-nucleus cross section scaling law and so on the particular choice of the couplings of the Hamiltonian (Eq. 1). However, as we will show, an important dependence on the Hamiltonian (Eq. 1) is still present when comparing COSINE-100 and DAMA. This is due to the fact that, although an initial modulation analysis of COSINE-100 with two-year data is forthcoming and an additional low-threshold analysis is also actively under development, COSINE-100 needs to collect several years of data [21] in order to reach the sensitivity required to probe the DAMA signal, and until then COSINE-100 will only exploit the average count rate. So the results of the two experiments are presently based on two different observables, the yearly modulation expected from the rotation of the Earth around the Sun and the time-averaged rate, and their relative size does depend on the specific model of WIMP-nucleus interaction besides a standard SI or SD interaction with nuclei.
The paper is organized as follows: in Section II we summarize how we calculate WIMP direct detection rates in NR effective theory; in Section III we outline the geometrical method of Ref. [13] that we use to study the intersections among the ellipsoidal surfaces of constant signal of DAMA and other experiments in the space of the coupling constants of the effective theory; our quantitative analysis is contained in Section IV. We devote Section V to our conclusions. For completeness we summarize the response functions of a WIMP of spin ≤1 in Appendix A and we provide the details of our treatment of experimental constraints in Appendix B.

II. WIMP INELASTIC SCATTERING IN NON-RELATIVISTIC EFFECTIVE MODELS
In the present Section we briefly summarize the ingredients that we use to calculate for each experiment and for each energy bin used in our analysis the matrix introduced in Eq. (2), needed to evaluate the expected rate to compare to the experimental data.
The full list of operators O j entering the Hamiltonian of Eq. 1 for the nuclear scattering process of a WIMP particle of spin Jχ ≤ 1 is given by: In the equation above 1 χN is the identity operator, q is the transferred momentum, S χ and S N are the WIMP and nucleon spins, respectively, while S = 1 2 (ǫ † i ǫ j + ǫ † j ǫ i ) is the symmetric combination of polarization vectors in the case of a spin-1 DM particle and v ⊥ = v + q 2µ χN (with µ χN the WIMPnucleon reduced mass) is the relative transverse velocity operator satisfying v ⊥ · q = 0. For a nuclear target T the quantity (v ⊥ T ) 2 ≡ | v ⊥ T | 2 can also be written as [22]: where v min is given by Eq.(4).
Operator O 2 is of higher order in v compared to all the others, implying a cross section suppression of order O(v/c) 4 ) ≃ 10 −12 for the non-relativistic WIMPs in the halo of our Galaxy. Moreover it cannot be obtained from the leadingorder non relativistic reduction of a manifestly relativistic operator [10]. So, following Ref. [10,11], we will not include it in our analysis. Moreover, operator O 16 is a linear combination of other operators, so can be omitted. This implies a maximal number of 16 operators in the effective Hamiltonian in Eq. (1), namely 4 operators for a spin-0 DM particle, 14 operators for spin 1/2 and 10 operators for spin 1.
To reduce the parameter space of the effective interaction of Eq. (1) in the following we will make a few simplifying assumptions. First, we will only consider the case δ >0, i.e. upscatters of a light state to a heavier one; then we will assume a contact effective interaction between the WIMP and the nucleus, i.e., we will assume the coefficients c τ j as independent on the transferred momentum q and neglect propagator effects. Moreover, we will consider real c τ j 's, although in general for an inelastic process they can be complex [22]. Finally, we will not consider the possibility of inelastic scattering among states of different spins [22].
The expected rate in a given visible energy bin E ′ 1 ≤ E ′ ≤ E ′ 2 of a direct detection experiment is given by: with ǫ(E ′ ) ≤ 1 the experimental efficiency/acceptance. In the equations above E R is the recoil energy deposited in the scattering process (indicated in keVnr), while E ee (indicated in keVee) is the fraction of E R that goes into the experimentally detected process (ionization, scintillation, heat) and q(E R ) is the quenching factor, G T (E ′ , E ee = q(E R )E R ) is the probability that the visible energy E ′ is detected when a WIMP has scattered off an isotope T in the detector target with recoil energy E R , M is the fiducial mass of the detector and T the live-time of the data taking. For a given recoil energy imparted to the target the differential rate for the WIMP-nucleus scattering process is given by: where ρ WIMP is the local WIMP mass density in the neighborhood of the Sun, N T the number of the nuclear targets of species T in the detector (the sum over T applies in the case of more than one nuclear isotope), while with m T the nuclear target mass and, assuming that the nuclear interaction is the sum of the interactions of the WIMPs with the individual nucleons in the nucleus: In the above expression j χ and j T are the WIMP and the target nucleus spins, respectively, q = | q| while the R ττ ′ k 's are WIMP response functions (that we report for completeness in Eq.(A1)) which depend on the couplings c τ j as well as the transferred momentum q and (v ⊥ T ) 2 . In equation (13) the W ττ ′ T k (y)'s are nuclear response functions and the index k represents different effective nuclear operators, which, crucially, under the assumption that the nuclear ground state is an approximate eigenstate of P and CP, can be at most eight: following the notation in [10,11] where b is the size of the nucleus. For the target nuclei T used in most direct detection experiments the functions W ττ ′ T k (y), calculated using nuclear shell models, have been provided in Refs. [11,23] under the assumption that the dark matter particle couples to the nucleus through local one-body interactions with the nucleons. In our analysis we do not include two-body effects [24,25] which are only available for a few isotopes and can be important when the one-body contribution is suppressed. Finally, f ( v T ) is the WIMP velocity distribution, for which we assume a standard isotropic Maxwellian at rest in the Galactic rest frame truncated at the escape velocity u esc , and boosted to the Lab frame by the velocity of the Earth. So for the former we assume: with z = 3u 2 esc /(2v 2 rms ). In the isothermal sphere model hydrothermal equilibrium between the WIMP gas pressure and gravity is assumed, leading to v rms = √ 3/2v 0 with v 0 the galactic rotational velocity.
With the exception of DAMA, all the experiments included in our analysis are sensitive to the time average of the expected rate for which < v E >=v ⊙ and v ⊙ =v 0 +12 (accounting for a peculiar component of the solar system with respect to the galactic rotation). In the case of DAMA, the yearly modulation effect is due to the time dependence of the Earth's speed with respect to the Galactic frame, given by: where cos γ ≃0.49 accounts for the inclination of the ecliptic plane with respect to the Galactic plane, T 0 =1 year and v orb =2πr ⊕ /(T 0 ) ≃ 29 km/sec (r ⊕ =1 AU neglecting the small eccentricity of the Earth's orbit around the Sun). In our analysis for the two parameters v 0 and u esc we take v 0 =220 km/sec [26] and u esc =550 km/sec [27] as reference values, although in Section IV we will also discuss the de-pendence of the results when the same parameters are varied in the ranges v 0 =(220± 20) km/s [26] and u esc =(550± 30) km/s [27]. Our reference choice of parameters corresponds to the WIMP escape velocity in the lab rest frame v lab esc ≃ 782 km/s. To make contact with other analyses, for the dark matter density in the neighborhood of the Sun we use ρ WIMP =0.3 GeV/cm 3 , which is a standard value commonly adopted by experimental collaborations, although observations point to the slightly higher value ρ WIMP =0.43 GeV/cm 3 [28,29]. Notice that direct detection experiments are only sensitive to the product ρ WIMP σ p , so the results of the next Section can be easily rescaled with ρ WIMP .
In particular, in each visible energy bin DAMA is sensitive to the yearly modulation amplitude S m , defined as the cosine with T 0 =1 year and t 0 =2 nd June, while other experiments put upper bounds on the time average S 0 : Using the ingredients listed above, for a given value of the two parameters m χ and δ both S 0 and S m can be expressed as quadratic forms like Eq. (2), i.e. for each of the experimental observable considered in our analysis a real symmetric matrix can be obtained. Schematically, for each energy bin n and both for DAMA and for each of the other experiments exp:

III. MAXIMAL DAMA SIGNALS COMPATIBLE TO NULL RESULTS
Following the analysis in [13], in this section we will use the property that, for a fixed value of the WIMP mass m χ and of the mass splitting δ, constant-rate surfaces in the couplings vector space are ellipsoids. In particular, given the experimental upper bound N exp 0,n for experiment exp and energy bin n, the condition S exp 0,n < N exp 0,n implies that allowed configurations must lie inside the ellipsoid: which add to the previous constraints. From now on we will indicate all upper bound matrices as º j , for j ∈ ¾, with ¾ the full set of experimental upper constraints including the upper bounds on the DAMA modulation amplitudes, so that the following conditions must be verified: In our analysis we will include the 8 DAMA modulation amplitudes for 1 keVee ≤ E ′ ≤ 5 keVee [4], and selected energy bins from XENON1T [30], PICO-60 (C 3 F 8 target) [31], COSINE-100 [8], COUPP [32], SuperCDMS [33] and PI-CASSO [34]. The details of how we implemented the DAMA effect and the bounds are provided in Appendix B and in Table III. Given the DAMA modulation amplitudes, an explanation of the effect in terms of WIMPs implies also the lower bounds: In total, we have solved Eq. (25) using 26 matrices (8+8 DAMA matrices plus 10 matrices for null results). Compatibility between DAMA and the other experiments is achieved only if in the coupling constants parameter space the intersection between the volumes outside the ellipsoids c T » n (m χ , δ)c=1, n = 1, ...N and the volume inside the ellipsoids c T º k (m χ , δ)c=1, k ∈ ¾ is non-vanishing. To prove this it is sufficient to find a set of real parameters ξ k ≤ 0 that, for each DAMA energy bin n satisfy [35]: In particular, if the binary test above is verified no set of couplings c exists for which the two conditions (23) and (24) are satisfied at the same time. Geometrically, this implies that the volume of intersection among the ellipsoids º k is fully contained in the DAMA ellipsoid of » n . On the other hand, when the matrix of Eq. amplitudes in the other energy bins. An alternative way to show the result of the test above is to calculate the maximal value of the modulation amplitude allowed by present constraints, i.e. to take S DAMA,max m,n in » n as a free parameter and find the minimal valueŜ DAMA,max m,n for which the condition (25) is verified [13,14]. A schematic view of the intersection between the c T DAMA m,n (m χ , δ)/Ŝ DAMA,max m,n c = 1 ellipse and the edge of the experimentally allowed volume in the coupling constants space is provided in Fig. 1 for the case of a single coupling and two upper bounds. Such value can then be converted in a number of standard deviations n σ away from the measured value. The tension between DAMA and the other experiments can then be quantified as the maximum of N σ ≡ max(n σ ) among the DAMA energy bins.

IV. ANALYSIS
In this section we discuss the N σ parameter solving Eq. (25) using PICOS [36], a Python interface to conic optimization, together with the CVXOPT solver [37].
The main results of our analysis is shown in Fig. 2. In such figure, for different value of the WIMP spin j χ , the upper band shows N σ as a function of m χ in the elastic case (δ=0), while the lower bands represent N σ minimized in terms of δ at fixed  [26] and u esc =(550± 30) km/s [27]. The gray shaded regions represent j χ =0, the red bands j χ =1/2 and the purple ones represents j χ =1. For each value of j χ the upper band represents the elastic case (δ=0), while the lower one the inelastic case, when N σ is minimized in terms of δ. at fixed m χ . varied in the ranges v 0 =(220± 20) km/s [26] and u esc =(550± 30) km/s [27], while the solid line indicates the result for the reference values v 0 =(220 km/s and u esc =(550 km/s. The gray shaded regions represent j χ =0, the red bands j χ =1/2 and the purple ones represent j χ =1. As far as the δ=0 case is concerned, a DAMA explanation is excluded at more than ≃ 7 sigmas, with the exception of the case j χ =1/2 at m χ ≃ 10 GeV, were the tension drops just below that value. Compared to the elastic case, inelastic scattering partially relieves this tension with values as low as ≃ 3.6σ for j χ =0, ≃ 3.3σ for j χ =1/2 and ≃ 3.5σ for j χ =1/2. However N σ is considerably more sensitive on the astrophysical parameters v 0 and u esc compared to the elastic case. From this figure one can conclude that neither the large range of different interactions provided by the EFT nor the modified kinematics due to inelasticity can eliminate completely the tension between DAMA and experimental constraints. When the quantity N σ is plotted in the m χ -δ plane a general result common to all j χ values is that the region of parameter space where the tension is relieved is localized in a narrow region with δ > ∼ 20 keV. In Figs. 3,4 and 5 we provide contour plots of N σ centered on such localized regions of the m χ -δ parameter space for j χ =0, 1/2, 1, respectively, and for the three combinations (v 0 ,u esc )=(220,550),(200,520),(240,580) km/s. In all the plots the dotted (red) line represents the maximal value of δ beyond which the minimal speed v * Na min introduced in Eq. (3) exceeds the escape velocity in the Lab frame, v lab esc , while the (blue) short dashes show the minimal value of δ beyond which v * F min < v lab esc . This implies that between the two lines the condition (5) is verified. The closed contour where the tension N σ drops tracks for different astrophysical parameters the region between the two lines, an unequivocal indication that the same kinematic mechanism is at work as in the pSIDM scenario [7,15,16] summarized in the Introduction. This also explains why, as observed in Fig. 2, N σ is considerably more sensitive on the astrophysical parameters v 0 and u esc compared to the elastic case. So we conclude that the pSIDM scenario described in ([7, 15, 16]) emerges as the unique mechanism to ease the tension between DAMA and other constraints from a general scan of the inelastic DM parameter space.
In spite of the residual tension between DAMA and other constraints it is interesting to discuss in detail the solutions corresponding to the minimum values of N σ . In order to do so one needs to go back to Eq. (25). At the boundary of its positivity the smallest eigenvalue of the matrix in that equation is vanishing and the corresponding eigenvectorĉ 0 (|ĉ 0 |=1): individuates the line joining the origin to the points of intersections between the extreme DAMA ellipsoid and those of the constraints, as shown schematically by the solid line in Fig. 1. On the other hand, the vector c 0,max , joining the origin to the intersection points and given by: , c τ 10 , +1 + 3 + 1 + 1)=28 (c τ 11 , c τ 12 , c τ 15 ),c τ 13 ,c τ 14 1 c τ 1 ,(c τ 4 , c τ 5 ), (c τ 9 , c τ 9 ), 2×(1 + 2 + 2 + 1+ c τ 10 , c τ 11 , c τ 14 , c τ 17 , c τ 18 1 + 1 + 1 + 1)=20 contains the set of couplings for the configuration of maximal modulation amplitude in the bin n, and in Fig. 1 is represented by the black arrow. The properties of the specific point c 0,max in the space of couplings can be further elucidated if one takes a closer look at the matrices in Eqs. (19,20). Depending on the spin of the WIMP particle they can have a different dimensionality, as can be simply read-off from the WIMP response functions in Eqs. (A1,A2). Moreover, not all couplings interfere, so that the matrices can be decomposed into block-diagonal form. The dimensionalities and non-interfering subspaces are indicated in Table I for different values of the WIMP spin j χ . This implies that also the matrix max is block-diagonal, so that c 0,max must belong to one of the subspaces of Table I.
The properties of the extreme configurations found in this way are given in table II. Interestingly, the configuration with the smallest tensions corresponds to a c 7 coupling for j χ =0 (a spin-dependent interaction with explicit velocity dependence and momentum suppression q 2 )) and to approximately a c 6 coupling for j χ =1/2 (a spin-dependent interaction with momentum suppression q 4 ). These two couplings combinations correspond to two of the possible generalizations of the pSIDM scenario already discussed in [16]. On the other hand,  for j χ =1 the extreme configuration corresponds to a dominant c 5 coupling (associated to the WIMP coupling to the orbital angular momentum operator) also with momentum suppression q 4 , and a non-negligible c 4 contribution. The role of momentum suppression in relieving the tension between the DAMA result and other constraints has already been pointed out in [38].
The corresponding predictions for the DAMA modulation amplitudes are shown in Figs. 6,7 and 8 for the different values of j χ and compared to the measured ones [4]. In such figures the experimental points marked with a (red) circle correspond to the energy bins included in the solution of Eq. (25), while the point marked with an additional (blue) inner circle corresponds to the DAMA energy bin where the maximal tension with the bounds arises and that drives the determination of N σ . As one can see in the lower energy bins the allowed modulation signal is large enough to explain the DAMA signal, although the amplitudes spectrum decays faster with energy.
Introducing the vector c 0 ≡ c 0ĉ0 with c 0 a free normalization, c 0 is common to all signals inside the eigenspace ofĉ 0 . A convenient parameterization is through the introduction of the reference cross section: with µ χN the WIMP-nucleon reduced mass. The direction in coupling space singled out by the unit vectorĉ 0 , and that individuates a specific set of coupling ratios that eases the tension between DAMA and the constraints, can be seen as the generalization in an arbitrary number of dimensions of the concept of isospin-violating DM [39]. Onceĉ 0 and δ are fixed the DAMA signal and the bounds can be discussed in a familiar m χ -σ 0 plane. This is done in Figs. 9, 10 and 11 for the extreme configurations summarized in Table II and Table II, while the effective cross section σ 0 is defined in Eq. (28).The (red) circle represents the point in parameter space with minimal N σ . more of the most constraining boundaries on σ 0 at m χ =m χ,0 , providing a nice confirmation of the numerical solution of Eq. (25). Actually, this can be directly observed in Fig. 10 for j χ =1/2 and Fig. 11 for j χ =1 where the point (m χ,0 , σ 0,max ) lies on the intersection between the bounds from XENON1T and PICO-60 (in a realization of the mechanism shown schematically in Fig. 1) and in Fig. 9 where it is constrained by PICO-60 only. In all three cases, the extreme configuration hits also the 90% C.L. upper bound on the modulation fraction in the first bin, as shown in Figs. 6, 7 and 8. In Figs. 9, 10 and 11 the distance from the DAMA region and the extreme point provides an additional visual indicator besides N σ of the tension between DAMA and the constraints from null results. Indeed, as already observed in [7], the condition (Eq. 5) implies that inside the energy range of the DAMA effect the spectrum of the predicted modulation amplitudes has a maximum corresponding to the recoil energy E * R ≡ E R (v * Na min )=|δ|µ χN /m N for scattering events off sodium. On the other hand, the data from DAMA/LIBRA-phase2 are more compatible to a monotonically decreasing shape 1 closer to elastic scattering. As a consequence, the DAMA data pull to low values of δ. However, the solutions of Eq. (25) with smallest tension with the constraints require sizeable values of δ (δ > ∼ 20keV) in order to verify Eq. (5). As a consequence, when δ is fixed to such values the DAMA data pull to higher values of the WIMP mass m χ to dilute the effect of δ. This explains why, systematically, the DAMA regions in Figs. 9, 10 and 11 are at higher WIMP masses compared to the values of m χ,0 in Table II. Moreover, Figs. 6, 7 and 8 show that configurations allowed by constraints from null results can produce enough yearly modulation in some of the DAMA bins, but the ensuing shape of the modulation spectrum is too steep, so that the maximal modulation at high energy is constrained by the bins at low energy. In light of this observation, we interpret the fact that all the smallest tension configurations of Table II have an interaction with explicit momentum suppression as a way to alleviate this problem by suppressing the DAMA modulation amplitudes in the lowest-energy bins.
These findings are in agreements to those of Ref. [7], ob-tained for the specific case of a standard spin-dependent interaction. Equation (25) can only be solved for a limited selection of experimental bounds both because of computing time limits, and because some of the constraints require more refined treatments beside a simple comparison between theoretical predictions and upper bounds as in Eq. (23) and Table III, such as background subtraction or the optimal-interval method [40]. In Figs. 9, 10 and 11 all this standard machinery [41] can instead be applied to the full set of existing experiments, providing an a posteriori confirmation that the set of bounds ¾ used to solve Eq. (25) did not miss any relevant constraint. In particular, besides the experiments included in the solution of Eq. (25), in such figures we have added CDMSlite [42], CRESST-II [43,44], the upper bound from the average count rate of DAMA [45]), DarkSide-50 [46] and the CF 3 I target run of PICO-60 [47] (the details of such bounds implementation are provided in Appendix B). None of these additional null results further constrains the extreme configurations (m χ,0 , σ 0,max ).
We conclude by noting that the recent bound from COSINE-100 [8], obtained with the same NaI target material as DAMA, is not particularly binding in our analysis, as can be seen again in Figs. 9, 10 and 11. The reason of this its that the bound of Ref. is between ≃0.05 and ≃0.12 for E ee < 3 keVee) explaining why in Ref. [8] the DAMA effect is ruled out. However expected rates for inelastic scattering are sensitive to the high-speed tail of the WIMP velocity distribution for which the modulation fractions are sizeably higher [15], and this is particularly true when the condition (Eq. 5) is verified (for instance, in the extreme configurations of Table II we find S DAMA m /S DAMA 0 > ∼ 0.8 for E ′ ≤3.5 keVee).

V. CONCLUSIONS
In the present paper we have discussed the compatibility of the combined annual modulation effect measured by DAMA/LIBRA-phase1 and DAMA/LIBRA-phase2 [1][2][3][4] with an explanation in terms of inelastic scattering events in-duced by the most general Galilean-invariant effective contact interaction of a spin 0, 1/2 or 1 WIMP dark matter particle taking into account all the possible interferences among operators by studying the intersections among the ellipsoidal surfaces of constant signal of DAMA and other experiments in the space of the coupling constants of the effective theory, following the approach introduced in Ref. [13]. In our analysis we have assumed a standard Maxwellian velocity distribution in the Galaxy. Compared to the elastic case, inelastic scattering partially relieves but does not eliminate the existing tension between the DAMA effect and the constraints from the null results of other experiments. We have determined the ellipsoids using 90% C.L. upper bounds from selected energy bins from XENON1T [30], PICO-60 (C 3 F 8 target) [31], COSINE-100 [8], COUPP [32], SuperCDMS [33] and PICASSO [34]. The tension, quantified as the maximum of N σ ≡ max(n σ ) among the DAMA energy bins below 5 keVee, exceeds 7σ in all the parameter space m χ < 200 GeV with the exception of a small region of parameter space for m χ ≃ 10 GeV and δ > ∼ 20, where it drops to values as low as ≃ 4σ for j χ =0, ≃ 2.9σ for j χ =1/2 and ≃ 3.2σ for j χ =1/2, and that overlaps to the proton-philic spin-dependent inelastic Dark Matter (pSIDM) scenario [7,15,16] already discussed in the literature for the specific case of a standard spin-dependent interaction, where the bounds from fluorine targets are evaded in a kinematic way because the minimal WIMP incoming speed required to trigger upscatters off fluorine exceeds the maximal WIMP velocity in the Galaxy, or is very close to it. In particular, from a general scan of the inelastic DM parameter space such kinematic feature, together with momentum suppression in the effective operator, emerge as instrumental in easing the tension between DAMA and other constraints. As a consequence, the latter is considerably more sensitive on the astrophysical parameters compared to the elastic case. The configurations for which the tension N σ is partially relieved can easily produce enough yearly modulation in the lowestenergy bins of the modulation spectrum measured by DAMA in compliance with the constraints from other experiments. However, the ensuing shape of the modulation spectrum is too steep, so that, when not excluded by other constraints, the maximal allowed modulation at higher energies is constrained by the modulation measured in the lowest energy bins.
The present analysis extends the scope of previous ones in the task to explore a DAMA explanation in the full WIMP direct detection parameter space, but is still not the most general one. Possible extensions include: i) long-range interactions; ii) allowing for complex couplings; iii) assuming a WIMP velocity distribution that departs from a standard Maxwellian. In particular, given the large dependence on the astrophysical parameters that we observed in our results we expect the latter generalization as very promising in order to find effective models that reconcile the DAMA result with the null observations of other experiments.
We collect here the WIMP particle-physics response functions introduced in Eq. (13) and for the general case of complex couplings (although in the present analysis real couplings have been assumed). For a WIMP particle of spin J χ ≤ 1 2 they are given by [10,11]: On the other hand, for a WIMP particle with spin J χ =1 [12]: (A2)

Xenon: XENON1T
For XENON1T we have assumed 7 WIMP candidate events in the range of 3PE ≤ S 1 ≤ 70PE, as shown in Fig. 3 of Ref. [30] for the primary scintillation signal S1 (directly in Photo Electrons, PE), with an exposure of 278.8 days and a fiducial volume of 1.3 ton of xenon. We have used the efficiency taken from Fig. 1 of [30] and employed a light collection efficiency g 1 =0.055; for the light yield L y we have extracted the best estimation curve for photon yields n ph /E from Fig. 7 in [49] with an electric field of 90 V/cm.
For XENON1T experiment we have modeled the energy resolution combining a Poisson fluctuation of the observed primary signal S 1 compared to < S 1 > and a Gaussian response of the photomultiplier with σ PMT = 0.5, so that: with Poiss(n, λ) = λ n /n!exp(−λ).

Argon: DarkSide-50
The analysis of DarkSide-50 [46] is based on the ionization signal extracted from liquid argon with an exposure of 6786.0 kg days. The measured spectrum for N e − < 50 (with N e − the number of extracted electrons) is shown in Fig. 7 of [46], and shows an excess for 4 < N e − <7 N e − compared to a simulation of the background components from known radioactive contaminants. Following Ref. [46] we have subtracted the background minimizing the likelihood function: where i represents the energy bin, x i the measured spectrum with error σ i , while σS i and ρb i are the DM signal and the background, respectively, with σ and ρ arbitrary normalization factors (σ is identified with the effective WIMPproton cross section σ p ). In particular we obtain the 90% C.L. upper bound on σ p by taking its profile likelihood with −2L − [−2L] min = n 2 and n=1.28. We take x i , σ i and b i from Fig.7 of [46]. The ionization yield of argon has been measured only down to < ∼ 10 keVnr, while DS50 uses a model fit to calibration data. We use the latter as taken from Fig. 6 of [46] with a hard cut at 0.15 keVnr, the lowest energy for which it is provided. We take the efficiency from Fig. 1 of [46].

Germanium: SuperCDMS and CDMSlite
The latest SuperCDMS analysis [33] observed 1 event between 4 and 100 keVnr with an exposure of 1690 kg days. We have taken the efficiency from Fig.1 of [33] and the energy resolution σ = 0.293 2 + 0.056 2 E ee from [50]. To analyze the observed spectrum we apply the optimal interval method [40].
For CDMSlite we considered the energy bin of 0.056 keV< E ′ < 1.1 keV with a measured count rate of 1.1±0.2 [keV kg day] −1 (Full Run 2 rate, Table II of Ref. [42]). We have taken the efficiency from Fig.4 COUPP is bubble chamber using a CF 3 I target. For each operating threshold used in COUPP the corresponding exposure and number of measured events are summarized in Table  IV. For fluorine and carbon we use α=0.15 in Eq.(B3). For iodine we adopt instead a step function with nucleation probability equal to 1 above the energy threshold.
The PICASSO experiment [34] uses C 4 F 10 as a target and operated its runs with six energy thresholds. For each threshold we provide the corresponding number of observed events and statistical fluctuations in Table V (extracted from Fig. 4 of Ref. [34]). For the nucleation probability we used Eq.(B3) with α C =α F =5.
The target material of PICO-60 is C 3 F 8 . For PICO-60 [31] only the threshold E th =3.3 keV is analyzed with a total exposure of 1167.0 kg days and no event detected. We have assumed the nucleation probability in Fig. 4 of [51].

Fluorine+Iodine: PICO-60
PICO-60 can also employ a CF 3 I target. For the analysis of Ref. [47] we adopt an energy threshold of 13.6 keV and an exposure of 1335 kg days. The nucleation probabilities for each target element are taken from Fig.4 in [47]. 6. Sodium Iodide: DAMA and COSINE-100 For DAMA we consider both the upper bound from the average count rate (DAMA0) and the latest result for the annual modulation amplitudes. For DAMA0 we have taken the average count rates from [45] (rebinned from 0.25-keVee-to 0.5-keVee-width bins) from 2 keVee to 8 keVee. We use the DAMA modulation amplitudes normalized to kg −1 day −1 keVee −1 in the energy range 1 keVee < E ′ < 8 keVee from Ref. [4]. In both cases we assume a constant quenching factors q=0.3 for sodium and q=0.09 for iodine, and the energy resolution σ = 0.0091 (E ee /keVee) + 0.448 √ E ee /keVee in keV. The exclusion plot for COSINE-100 [8] relies on a Montecarlo [52] to subtract the different backgrounds of each of the eight crystals used in the analysis. In Ref. [8] the amount of residual background after subtraction is not provided, so we have assumed a constant background b at low energy (2 keVee< E ee < 8 keVee), and estimated b by tuning it to repro-duce the exclusion plot in Fig.4 of Ref. [8] for the isoscalar spin-independent elastic case. The result of our procedure yields b ≃0.13 events/kg/day/keVee, which implies a subtraction of about 95% of the background. We take the energy resolution σ/keV = 0.3171 √ E ee /keVee+0.008189E ee /keVee averaged over the COSINE-100 crystals [53] and the efficiency for nuclear recoils from Fig.1 of Ref. [8]. Quenching factors for sodium and iodine are assumed to be equal to 0.3 and 0.09 respectively, the same values used by DAMA.

CaWO 4 : CRESST-II
CRESST-II measures heat and scintillation using CaWO 4 crystals. We considered the Lise module analysis from [43] with energy resolution σ=0.062 keV and detector efficiency from Fig. 4 of [54]. For our analysis we have selected 15 events for 0.3 keVnr< E R < 0.49 keVnr with an exposure of 52.15 kg days.