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

TL;DR

This study explores whether inelastic WIMP dark matter models can explain DAMA/LIBRA's annual modulation signal while remaining consistent with other experimental constraints, finding limited regions of parameter space that partially reconcile these results.

Contribution

It provides a comprehensive analysis of inelastic WIMP models with general Galilean-invariant interactions, considering all operator interferences and astrophysical uncertainties.

Findings

Inelastic scattering slightly reduces tension with DAMA and other experiments.

A small parameter region around 10 GeV WIMP mass and >20 keV mass splitting partially relieves tension.

Residual tension persists, especially with astrophysical parameter variations.

Abstract

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.

Tables

Table 1. Table 1: Non–interfering subspaces and dimensionality of the coupling constants vector space of the NR WIMP effective theory of a WIMP with j χ ≤ subscript 𝑗 𝜒 absent j_{\chi}\leq 1.

spin	couplings	dimensionality
0	$(c_1^{\tau },c_3^{\tau })$,$c_7^{\tau }$,$c_{10}^{\tau }$	2$\times \left(2+1+1\right)$=8
$\frac{1}{2}$	$(c_1^{\tau },c_3^{\tau })$,$(c_4^{\tau },c_5^{\tau },c_6^{\tau })$,	2$\times (2+3+1+2+$
	$c_7^{\tau }$, $(c_8^{\tau },c_9^{\tau })$, $c_{10}^{\tau }$,	$+1+3+1+1)$=28
	$(c_{11}^{\tau },c_{12}^{\tau },c_{15}^{\tau })$,$c_{13}^{\tau }$,$c_{14}^{\tau }$
1	$c_1^{\tau }$,$(c_4^{\tau },c_5^{\tau })$, $(c_8^{\tau },c_9^{\tau })$,	2$\times (1+2+2+1+$
	$c_{10}^{\tau }$, $c_{11}^{\tau }$, $c_{14}^{\tau }$, $c_{17}^{\tau }$, $c_{18}^{\tau }$	$1+1+1+1)$=20

Table 2. Table 2: Properties of the extreme configurations in the NR effective theory parameter space that minimize the tension N σ subscript 𝑁 𝜎 N_{\sigma} for different values of the WIMP spin j χ subscript 𝑗 𝜒 j_{\chi} and fixing the astrophysocal parameters to v 0 subscript 𝑣 0 v_{0} =220 km/s and u e ​ s ​ c subscript 𝑢 𝑒 𝑠 𝑐 u_{esc} =550 km/s.

spin	$m_{\chi ,0}$ (GeV)	$\delta \left(keV\right)$	${\sigma }_{0,max}$ (cm2)	${\widehat{𝐜}}_0$	$N_{\sigma }$
0	11.08	22.83	3.93$\times {10}^{-27}$	(${\hat{c}}_7^0$=0.68, ${\hat{c}}_7^1$=0.73)	4.0
1/2	11.64	23.74	4.68$\times {10}^{-28}$	(${\hat{c}}_4^0$=-0.0014, ${\hat{c}}_4^1$=-0.0015,	2.9
				${\hat{c}}_5^0$=-0.032, ${\hat{c}}_5^1$=-0.0166,
				${\hat{c}}_6^0$=0.692, ${\hat{c}}_6^1$=0.7217)
1	11.36	23.43	5.71$\times {10}^{-32}$	(${\hat{c}}_4^0$=0.0717, ${\hat{c}}_4^1$=0.0753,	3.2
				${\hat{c}}_5^0$=0.1892, ${\hat{c}}_5^1$=0.9764)

Table 3. Table 3: Visible energy intervals and 90% C.L. upper or lower bounds used to calculate the matrices \varmathbb ​ A j \varmathbb subscript 𝐴 𝑗 \varmathbb{A}_{j} and \varmathbb ​ B j \varmathbb subscript 𝐵 𝑗 \varmathbb{B}_{j} in Eq. ( 25 ).

Experiment	visible energy range	90% C.L. upper bound
DAMA	1 keVee1.5 keVee	0.0315 kg-1day-1keVee-1
	1.5 keVee2 keVee	0.0268 kg-1day-1keVee-1
	2 keVee2.5 keVee	0.0210 kg-1day-1keVee-1
	2.5 keVee3 keVee	0.0236 kg-1day-1keVee-1
	3 keVee3.5 keVee	0.0222 kg-1day-1keVee-1
	3.5 keVee4 keVee	0.0144 kg-1day-1keVee-1
	4 keVee4.5 keVee	0.0137 kg-1day-1keVee-1
	4.5 keVee5 keVee	0.00569 kg-1day-1keVee-1
XENON1T	3PE 70PE	11.77 events
PICO–60	$E_R>$2.45 keVnr	6.42 events
	$E_R>$3.3 keVnr	2.3 events
COSINE–100	2 keVee2.5 keVee	0.13 kg-1day-1keVee-1
	4.5 keVee5 keVee	0.13 kg-1day-1keVee-1
	7.5 keVee8 keVee	0.13 kg-1day-1keVee-1
COUPP	$E_R>$7.8 keVnr	6.68  events
	$E_R>$11 keVnr	5.32  events
	$E_R>$15.5 keVnr	11.6  events
SuperCDMS	4 keVnr100 keVnr	3.89 events
PICASSO	$E_R>$1.0 keVnr	3.45 kg-1day-1keVee-1
Experiment	visible energy range	90% C.L. lower bound
DAMA	1 keVee1.5 keVee	0.0171 kg-1day-1keVee-1
	1.5 keVee2 keVee	0.0155 kg-1day-1keVee-1
	2 keVee2.5 keVee	0.0150 kg-1day-1keVee-1
	2.5 keVee3 keVee	0.0159 kg-1day-1keVee-1
	3 keVee3.5 keVee	0.0151 kg-1day-1keVee-1
	3.5 keVee4 keVee	0.00773 kg-1day-1keVee-1
	4 keVee4.5 keVee	0.00812 kg-1day-1keVee-1
	4.5 keVee5 keVee	0.000770 kg-1day-1keVee-1

Table 4. Table 4: The operating thresholds with corresponding exposures and measured events for COUPP  Behnke et al. ( 2012 ) .

$E_{th}$ (keV)	exposure (kg day)	measured events
7.8	55.8	2
11	70	3
15.5	311.7	8

Table 5. Table 5: Observed number of events and 1–sigma statistical fluctuations (extracted from Fig. 4 of Ref. Behnke et al. ( 2017 ) ) for each operating threshold used in PICASSO.

$E_{th}$ (keV)	Event rate (events/kg/day)	Fluctuation
1.0	-1.5	3.8
1.5	-0.2	1.0
2.7	0.3	0.8
6.6	-0.8	1.8
15.7	-1.4	2.3
36.8	0.3	1.0

Equations

$$\mathcal{H}(\mathbf{r})$$

$$\text{event rate}\propto {\mathbf{c}}^T\text{varmathbb}X\mathbf{c}.$$

$$v_{min}^{\ast }=\sqrt{\frac{2\delta }{{\mu }_{\chi N}}},$$

$$v_{min}=\frac{1}{\sqrt{2m_N{E}_R}}\left|\frac{m_N{E}_R}{{\mu }_{\chi N}}+\delta \right|,$$

$$v_{min}^{\ast Na}<v_{esc}^{lab}<v_{min}^{\ast F},$$

$${\mathcal{O}}_1$$

$${\mathcal{O}}_3$$

$${\mathcal{O}}_5$$

$${\mathcal{O}}_7$$

$${\mathcal{O}}_9$$

$${\mathcal{O}}_{11}$$

$${\mathcal{O}}_{13}$$

$${\mathcal{O}}_{15}$$

$${\mathcal{O}}_{16}$$

$${\mathcal{O}}_{17}$$

$$(v_T^{\perp }{)}^2=v_T^2-v_{min}^2.$$

$$R_{[E_1^{\prime },E_2^{\prime }]}$$

$$\frac{dR}{d{E}^{\prime }}$$

$$E_{ee}$$

$$\frac{d{R}_{\chi T}}{d{E}_R}(t)=\sum _T{N}_T\frac{{\rho }_{\text{mbox}WIMP}}{m_{\chi }}{\int }_{v_{min}}{d}^3{v}_{T}f({\vec{v}}_T,t)v_T\frac{d{\sigma }_T}{d{E}_R},$$

$$\frac{d{\sigma }_T}{d{E}_R}=\frac{2m_T}{4\pi v_T^2}\left[\frac{1}{2j_{\chi }+1}\frac{1}{2j_T+1}|{\mathcal{M}}_T{|}^2\right],$$

$$\frac{1}{2j_{\chi }+1}\frac{1}{2j_T+1}|{\mathcal{M}}_T{|}^2=$$

$$\frac{4\pi }{2j_T+1}\sum _{\tau =0,1}\sum _{{\tau }^{\prime }=0,1}\sum _k{R}_k^{\tau {\tau }^{\prime }}\left[c_j^{\tau },(v_T^{\perp }{)}^2,\frac{q^2}{m_N^2}\right]W_{Tk}^{\tau {\tau }^{\prime }}(y).$$

$$f({\vec{v}}_T,t)$$

$$N$$

$$|{\vec{v}}_E(t)|=v_{\odot }+v_{orb}\cos \gamma \cos \left[\frac{2\pi }{T_0}(t-t_0)\right],$$

$$S_{m,[E_1^{\prime },E_2^{\prime }]}\equiv \frac{2}{T_0}{\int }_0^{T_0}\cos \left[\frac{2\pi }{T_0}(t-t_0)\right]R_{[E_1^{\prime },E_2^{\prime }]}(t)dt,$$

$$S_{0,[E_1^{\prime },E_2^{\prime }]}\equiv \frac{1}{T_0}{\int }_0^{T_0}{R}_{[E_1^{\prime },E_2^{\prime }]}(t)dt.$$

$$S_{m,n}^{DAMA}(m_{\chi },\delta )$$

$$S_{0,n}^{exp}(m_{\chi },\delta )$$

$${\mathbf{c}}^T\text{varmathbb}{A}_n^{exp}(m_{\chi },\delta )\mathbf{c}\equiv {\mathbf{c}}^T\frac{\text{varmathbb}{S}_{0,n}^{exp}(m_{\chi },\delta )}{N_{0,n}^{exp}}\mathbf{c}<1.$$

$${\mathbf{c}}^T\text{varmathbb}{A}_n^{DAMA}(m_{\chi },\delta )\mathbf{c}\equiv {\mathbf{c}}^T\frac{\text{varmathbb}{S}_{m,n}^{DAMA}(m_{\chi },\delta )}{S_{m,n}^{DAMA,max}}\mathbf{c}<1,$$

$${\mathbf{c}}^T\text{varmathbb}{A}_j(m_{\chi },\delta )\mathbf{c}<1,j\in \text{varmathbb}E.$$

$${\mathbf{c}}^T\text{varmathbb}{B}_n(m_{\chi },\delta )\mathbf{c}\equiv {\mathbf{c}}^T\frac{\text{varmathbb}{S}_{m,n}^{DAMA}(m_{\chi },\delta )}{S_{m,n}^{DAMA,min}}\mathbf{c}>1.$$

$$\sum _{i\in \text{varmathbb}E}{\xi }_i<1,$$

$$\sum _{i\in \text{varmathbb}E}{\xi }_i\text{varmathbb}{A}_k-\text{varmathbb}{B}_n\text{mbox}isapositivematrix.$$

$$\left[\sum _{i\in \text{varmathbb}E}{\xi }_i\text{varmathbb}{A}_k-\frac{\text{varmathbb}{S}_{m,n}^{DAMA}}{{\hat{S}}_{m,n}^{DAMA,min}}\right]{\hat{\mathbf{c}}}_0\equiv \text{varmathbb}{M}_{max}{\hat{\mathbf{c}}}_0=0,$$

$${\mathbf{c}}_{0,max}=\sqrt{\frac{{\hat{S}}_{m,n}^{DAMA,max}}{{\hat{\mathbf{c}}}_0^T\text{varmathbb}{S}_{m,n}^{DAMA}{\hat{\mathbf{c}}}_0}}{\hat{\mathbf{c}}}_0$$

Full text

I Introduction

For more than 15 years the DAMA collaboration Bernabei et al. (2008a, 2010, 2013, 2018) 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 Ade et al. (2014) 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 Drukier et al. (1986). In particular, with the release of the latest DAMA/LIBRA-phase2 data Bernabei et al. (2018) the statistical significance of DAMA effect has reached almost $12\sigma$. However, in the most popular WIMP scenarios used to explain the DAMA signal as due to 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 Kang et al. (2019) also after the bound from the COSINE–100 collaboration Adhikari et al. (2018a), 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 Kurylov and Kamionkowski (2004); Savage et al. (2004); Foot (2008); Feng et al. (2008); Bottino et al. (2008); Kouvaris (2008); Tucker-Smith and Weiner (2005); Savage et al. (2009a); Bai and Fox (2009); Savage et al. (2009b); Fox and Poppitz (2009); Feldstein et al. (2010); Kopp et al. (2010); Del Nobile et al. (2015); Kang et al. (2018a, 2019). 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 Fitzpatrick et al. (2013); Anand et al. (2014), of the form:

[TABLE]

where $t^{0}=\mathbbold{1}$ is the identity in isospin space, $t^{1}=\tau_{3}$ is the third Pauli matrix, and ${\mathbf{r}}$ denotes the dark matter-nucleon relative distance and the operators $\mathcal{O}_{j}$ depend on the exchanged momentum $\vec{q}$, the WIMP incoming velocity $\vec{v}$, the WIMP spin $\vec{S}_{\chi}$ and nuclear spin $\vec{S}_{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\equiv|\vec{q}|$ and $v\equiv|\vec{v}|$ the ensuing Hamiltonian contains 8 independent couplings $c_{j}^{\tau}$ for a scalar WIMP, 28 independent couplings for a spin–1/2 WIMP Fitzpatrick et al. (2013); Anand et al. (2014) and 20 couplings for spin–1 Dent et al. (2015). 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 orders of magnitude. While, due to the large dimensionality, its direct scanning appears to be challenging this has been achieved by using matricial techniques Catena et al. (2016), exploiting the fact that in terms of the couplings vector ${\bf c}=\left(c_{1}^{(0)},c_{1}^{(1)},...,c_{n}^{(0)},c_{n}^{(1)}\right)^{T}$ for all direct detection experiments the expected event rate can be written in the form:

[TABLE]

with $\varmathbb{X}$ is a real symmetric $2n\times 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 effective field theory parameter space the surfaces of constant signal in different detectors are ellipsoids, and, as discussed in Refs. Catena et al. (2016, 2018), 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 non-relativistic effective theory, without making any a priori assumptions regarding the relative size of the various Wilson coefficients $c_{k}^{\tau}$. In this way for a standard halo model the DAMA/LIBRA–phase1 result Bernabei et al. (2008a, 2010, 2013) was shown to be incompatible to the constraints from other experiments in the case of elastic interactions of a WIMP particle with spin $\leq$ 1/2 and the interaction Hamiltonian of Eq.( 1) with arbitrary couplings combinations Catena et al. (2016).

As pointed out by the authors, the analysis in Catena et al. (2016) 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) Scopel and Yoon (2016); Scopel and Yu (2017); Kang et al. (2019). In such model the WIMP particle interacts with nuclear targets through 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) Tucker-Smith and Weiner (2001); Chang et al. (2009); March-Russell et al. (2009) reconciles the above scenario to fluorine detectors. In IDM a DM particle $\chi_{1}$ of mass $m_{\chi_{1}}=m_{\chi}$ interacts with atomic nuclei exclusively by up–scattering to a second heavier state $\chi_{2}$ with mass $m_{\chi_{2}}=m_{\chi}+\delta$. 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:

[TABLE]

with $\mu_{\chi 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:

[TABLE]

with $m_{N}$ the nuclear mass.

In particular, indicating with $v_{min}^{*Na}$ and $v_{min}^{*F}$ the values of $v_{min}^{*}$ for sodium and fluorine, and with $v_{esc}$ the WIMP escape velocity, in Refs. Scopel and Yoon (2016); Scopel and Yu (2017); Kang et al. (2019) constraints from WIMP–fluorine scattering events in droplet detectors and bubble chambers were shown to be evaded when the WIMP mass $m_{\chi}$ and the mass gap $\delta$ are chosen in such a way that the hierarchy:

[TABLE]

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. Catena et al. (2016) to extend the analyses of Scopel and Yoon (2016); Scopel and Yu (2017); Kang et al. (2019) to the general interaction Hamiltonian of Eq. (1), or, conversely, we wish to extend the analysis of Catena et al. (2016) 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. Catena et al. (2016). 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_{\chi}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 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_{\chi}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 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) dam (2018); Baum et al. (2018), 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 Kang et al. (2018b), with the exception of $\mathcal{O}_{1}$ and $\mathcal{O}_{4}$, which in the notation of Fitzpatrick et al. (2013); Anand et al. (2014) correspond respectively to the standard SI or Spin–Dependent (SD) interactions, this problem is not present for all the other operators $\mathcal{O}_{j}$ of the Hamiltonian (Eq. 1). This represents an additional motivation to update the result of Ref. Catena et al. (2016) 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 Thompson (2017) 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. Catena et al. (2016) 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 $\leq$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 $\varmathbb{X}$ introduced in Eq. (2), needed to evaluate the expected rate to compare to the experimental data.

The full list of operators $\mathcal{O}_{j}$ entering the Hamiltonian of Eq. 1 for the nuclear scattering process of a WIMP particle of spin $J{\chi}\leq$ 1 is given by:

[TABLE]

In the equation above $1_{\chi N}$ is the identity operator, $\vec{q}$ is the transferred momentum, $\vec{S}_{\chi}$ and $\vec{S}_{N}$ are the WIMP and nucleon spins, respectively, while $\mathcal{S}=\frac{1}{2}(\epsilon^{\dagger}_{i}\epsilon_{j}+\epsilon^{\dagger}_{j}\epsilon_{i})$ is the symmetric combination of polarization vectors in the case of a spin–1 DM particle and $\vec{v}^{\perp}=\vec{v}+\frac{\vec{q}}{2\mu_{\chi{\cal N}}}$ (with $\mu_{\chi{\cal N}}$ the WIMP–nucleon reduced mass) is the relative transverse velocity operator satisfying $\vec{v}^{\perp}\cdot\vec{q}=0$. For a nuclear target $T$ the quantity $(v^{\perp}_{T})^{2}\equiv|\vec{v}^{\perp}_{T}|^{2}$ can also be written as Barello et al. (2014):

[TABLE]

where $v_{min}$ is given by Eq.(4).

Operator ${\cal O}_{2}$ is of higher order in $v$ compared to all the others, implying a cross section suppression of order ${\cal O}(v/c)^{4})\simeq 10^{-12}$ for the non–relativistic WIMPs in the halo of our Galaxy. Moreover it cannot be obtained from the leading-order non relativistic reduction of a manifestly relativistic operator Fitzpatrick et al. (2013). So, following Refs.Fitzpatrick et al. (2013); Anand et al. (2014), we will not include it in our analysis. Moreover, operator $\mathcal{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 $\delta>$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}^{\tau}$ as independent on the transferred momentum $q$ and neglect propagator effects. Moreover, we will consider real $c_{j}^{\tau}$’s, although in general for an inelastic process they can be complex Barello et al. (2014). Finally, we will not consider the possibility of inelastic scattering among states of different spins Barello et al. (2014).

The expected rate in a given visible energy bin $E_{1}^{\prime}\leq E^{\prime}\leq E_{2}^{\prime}$ of a direct detection experiment is given by:

[TABLE]

with $\epsilon(E^{\prime})\leq 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, ${\cal G}_{T}(E^{\prime},E_{ee}=q(E_{R})E_{R})$ is the probability that the visible energy $E^{\prime}$ 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:

[TABLE]

where $\rho_{\mbox{\tiny 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

[TABLE]

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:

[TABLE]

In the above expression $j_{\chi}$ and $j_{T}$ are the WIMP and the target nucleus spins, respectively, $q=|\vec{q}|$ while the $R_{k}^{\tau\tau^{\prime}}$’s are WIMP response functions (that we report for completeness in Eq.(29)) which depend on the couplings $c^{\tau}_{j}$ as well as the transferred momentum $\vec{q}$ and $(v^{\perp}_{T})^{2}$. In equation (13) the $W^{\tau\tau^{\prime}}_{Tk}(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 Fitzpatrick et al. (2013); Anand et al. (2014), $k$=$M$, $\Phi^{\prime\prime}$, $\Phi^{\prime\prime}M$, $\tilde{\Phi}^{\prime}$, $\Sigma^{\prime\prime}$, $\Sigma^{\prime}$, $\Delta$, $\Delta\Sigma^{\prime}$. The $W^{\tau\tau^{\prime}}_{Tk}(y)$’s are function of $y\equiv(qb/2)^{2}$, where $b$ is the size of the nucleus. For the target nuclei $T$ used in most direct detection experiments the functions $W^{\tau\tau^{\prime}}_{Tk}(y)$, calculated using nuclear shell models, have been provided in Refs. Anand et al. (2014); Catena and Schwabe (2015) 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 two (2013, 2015) which are only available for a few isotopes and can be important when the one–body contribution is suppressed. Finally, $f(\vec{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:

[TABLE]

with $z=3u_{esc}^{2}/(2v_{rms}^{2})$. In the isothermal sphere model hydrothermal equilibrium between the WIMP gas pressure and gravity is assumed, leading to $v_{rms}$=$\sqrt{3/2}v_{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_{\odot}$ and $v_{\odot}$=$v_{0}$+12 km/sec (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:

[TABLE]

where $\cos\gamma\simeq$0.49 accounts for the inclination of the ecliptic plane with respect to the Galactic plane, $T_{0}$=1 year and $v_{orb}$=2$\pi r_{\oplus}/(T_{0})\simeq$ 29 km/sec ($r_{\oplus}$=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 Koposov et al. (2010) and $u_{esc}$=550 km/sec Piffl et al. (2014) as reference values, although in Section IV we will also discuss the dependence of the results when the same parameters are varied in the ranges $v_{0}$=(220$\pm$ 20) km/s Koposov et al. (2010) and $u_{esc}$=(550$\pm$ 30) km/s Piffl et al. (2014). Our reference choice of parameters corresponds to the WIMP escape velocity in the lab rest frame $v_{esc}^{lab}\simeq$ 782 km/s. To make contact with other analyses, for the dark matter density in the neighborhood of the Sun we use $\rho_{\mbox{\tiny WIMP}}$=0.3 GeV/cm3, which is a standard value commonly adopted by experimental collaborations, although observations point to the slightly higher value $\rho_{\mbox{\tiny WIMP}}$=0.43 GeV/cm3 Salucci et al. (2010); Nesti and Salucci (2013). Notice that direct detection experiments are only sensitive to the product $\rho_{\mbox{\tiny WIMP}}\sigma_{p}$, so the results of the next Section can be easily rescaled with $\rho_{\mbox{\tiny WIMP}}$.

In particular, in each visible energy bin DAMA is sensitive to the yearly modulation amplitude $S_{m}$, defined as the cosine transform of $R_{[E_{1}^{\prime},E_{2}^{\prime}]}(t)$:

[TABLE]

with $T_{0}$=1 year and $t_{0}$=2nd June, while other experiments put upper bounds on the time average $S_{0}$:

[TABLE]

Using the ingredients listed above, for a given value of the two parameters $m_{\chi}$ and $\delta$ 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$:

[TABLE]

III Maximal DAMA signals compatible to null results

Following the analysis in Catena et al. (2016), in this section we will use the property that, for a fixed value of the WIMP mass $m_{\chi}$ and of the mass splitting $\delta$, 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:

[TABLE]

As far as the DAMA modulation amplitudes are concerned, the experimentally observed intervals $[S^{DAMA,min}_{m,k},S^{DAMA,max}_{m,k}]$ in energy bins $k=1...N$ imply the additional upper bounds:

[TABLE]

which add to the previous constraints. From now on we will indicate all upper bound matrices as $\varmathbb{A}_{j}$, for $j\in\varmathbb{E}$, with $\varmathbb{E}$ the full set of experimental upper constraints including the upper bounds on the DAMA modulation amplitudes, so that the following conditions must be verified:

[TABLE]

In our analysis we will include the 8 DAMA modulation amplitudes for 1 keVee $\leq E^{\prime}\leq$ 5 keVee Bernabei et al. (2018), and selected energy bins from XENON1T Aprile et al. (2018a), PICO–60 ($C_{3}F_{8}$ target) Amole et al. (2017, 2019), COSINE–100 Adhikari et al. (2018a), COUPP Behnke et al. (2012), SuperCDMS Agnese et al. (2018a) and PICASSO Behnke et al. (2017). The details of how we implemented the DAMA effect and the bounds are provided in Appendix B and in Table 3. Given the DAMA modulation amplitudes, an explanation of the effect in terms of WIMPs implies also the lower bounds:

[TABLE]

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 ${\bf c}^{T}\varmathbb{B}_{n}(m_{\chi},\delta){\bf c}$=1, $n=1,...N$ and the volume inside the ellipsoids ${\bf c}^{T}\varmathbb{A}_{k}(m_{\chi},\delta){\bf c}$=1, $k\in\varmathbb{E}$ is non–vanishing. To prove this it is sufficient to find a set of real parameters $\xi_{k}\leq 0$ that, for each DAMA energy bin $n$ satisfy Pólik and Terlaky (2007):

[TABLE]

In particular, if the binary test above is verified no set of couplings ${\bf 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 $\varmathbb{A}_{k}$ is fully contained in the DAMA ellipsoid of $\varmathbb{B}_{n}$. On the other hand, when the matrix of Eq. (25) is not positive–defined such set of couplings exists. Notice that, in such case, since the $\varmathbb{A}_{k}$ matrices include the upper bounds on the DAMA modulation amplitudes, for that choice of $(m_{\chi},\delta)$ the condition $S^{DAMA,min}_{m,n}\in[S^{DAMA,min}_{m,n},S^{DAMA,max}_{m,n}]$ is automatically satisfied in the energy bin $n$ in compliance to all other existing constraints, although this is not guaranteed for the modulation 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 $\varmathbb{B}_{n}$ as a free parameter and find the minimal value $\hat{S}^{DAMA,max}_{m,n}$ for which the condition (25) is verified Catena et al. (2016, 2018). A schematic view of the intersection between the ${\bf c}^{T}\varmathbb{S}_{m,n}^{DAMA}(m_{\chi},\delta)/\hat{S}^{DAMA,max}_{m,n}\;{\bf 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_{\sigma}$ away from the measurement (in absolute value). The tension between DAMA and the other experiments can then be quantified as the maximum of $N_{\sigma}\equiv\max(n_{\sigma})$ among the DAMA energy bins calculated in the following way: i) we fix one target energy bin $n$; ii) we maximize the modulation signal in $n$; iii) for the corresponding set of couplings we calculate the modulation signal also in the other bins; iv) we take $n_{\sigma}$ as the maximum tension among all the bins; v) we loop over the target bin $n$ and take the minimum (since each target bin yields a different model). Notice that the procedure described above may not yield the model which better reproduces the data in the bins where the maximal allowed modulation exceeds the central value of the measurement.

In total, we have solved Eq. (25) using 27 matrices (8+8 DAMA matrices plus 11 matrices for null results).

IV Analysis

In this section we discuss the $N_{\sigma}$ parameter solving Eq. (25) using PICOS pic , a Python interface to conic optimization, together with the CVXOPT solver cvx .

The main results of our analysis is shown in Fig. 2. In such figure, for different value of the WIMP spin $j_{\chi}$, the upper band shows $N_{\sigma}$ as a function of $m_{\chi}$ in the elastic case ($\delta$=0), while the lower bands represent $N_{\sigma}$ minimized in terms of $\delta$ at fixed $m_{\chi}$. The bands indicate the variation of $N_{\sigma}$ when the parameters $v_{0}$ and $u_{esc}$ of the Maxwellian distribution of Eq. (15) are varied in the ranges $v_{0}$=(220$\pm$ 20) km/s Koposov et al. (2010) and $u_{esc}$=(550$\pm$ 30) km/s Piffl et al. (2014), 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_{\chi}$=0, the red bands $j_{\chi}$=1/2 and the purple ones represent $j_{\chi}$=1. As far as the $\delta$=0 case is concerned, a DAMA explanation is excluded at more than $\simeq$ 7 sigmas for all WIMP masses below 200 GeV. Compared to the elastic case, inelastic scattering partially relieves this tension with values as low as $\simeq$ 4.0$\sigma$ for $j_{\chi}$=0, $\simeq$ 2.9$\sigma$ for $j_{\chi}$=1/2 and $\simeq$ 3.2$\sigma$ for $j_{\chi}$=1. However $N_{\sigma}$ 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 effective field theory nor the modified kinematics due to inelasticity can eliminate completely the tension between DAMA and experimental constraints.

When the quantity $N_{\sigma}$ is plotted in the $m_{\chi}$–$\delta$ plane a general result common to all $j_{\chi}$ values is that the region of parameter space where the tension is relieved is localized in a narrow region with \delta\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 20 keV. In Figs. 3,4 and 5 we provide contour plots of $N_{\sigma}$ centered on such localized regions of the $m_{\chi}$–$\delta$ parameter space for $j_{\chi}$=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 $\delta$ beyond which the minimal speed $v^{*Na}_{min}$ introduced in Eq. (3) exceeds the escape velocity in the Lab frame, $v_{esc}^{lab}$, while the (blue) short dashes show the minimal value of $\delta$ beyond which $v^{*F}_{min}<v_{esc}^{lab}$. This implies that between the two lines the condition of Eq. (5) is verified. The closed contour where the tension $N_{\sigma}$ 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 Scopel and Yoon (2016); Scopel and Yu (2017); Kang et al. (2019) summarized in the Introduction. This also explains why, as observed in Fig. 2, $N_{\sigma}$ 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 (Scopel and Yoon (2016); Scopel and Yu (2017); Kang et al. (2019)) 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_{\sigma}$. 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 ${\bf\hat{c}}_{0}$ ($|{\bf\hat{c}}_{0}|$=1):

[TABLE]

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 ${\bf c}_{0,max}$, joining the origin to the intersection points and given by:

[TABLE]

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 ${\bf 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. (29,30). 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 1 for different values of the WIMP spin $j_{\chi}$. This implies that also the matrix $\varmathbb{M}_{max}$ is block–diagonal, so that ${\bf c}_{0,max}$ must belong to one of the subspaces of Table 1.

The properties of the extreme configurations found in this way are given in Table 2. Interestingly, the configuration with the smallest tensions corresponds to a $c_{7}$ coupling for $j_{\chi}$=0 (a spin–dependent interaction with explicit velocity dependence and momentum suppression $q^{2}$)) and to approximately a $c_{6}$ coupling for $j_{\chi}$=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 Scopel and Yu (2017). On the other hand, for $j_{\chi}$=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 Scopel et al. (2015).

Our procedure does not correspond to the minimization of a $\chi$–square, since we minimize the tension in one target bin at a time. However, once a minimal tension configuration is obtained, the quantity $\chi^{2}$=$\sum[S^{DAMA}_{m,n}-S^{DAMA}_{m,n,exp}]^{2}/\sigma_{exp}^{2}$ (with $S^{DAMA}_{m,n,exp}$ and $\sigma_{exp}$ the measured modulation amplitudes and standard deviations) can be calculated. In this way, for the configurations of Table 2 we get $\chi^{2}$=60.6, 25.7 and 33.5 for $j_{\chi}$=0,1/2 and 1, respectively.

The corresponding predictions for the DAMA modulation amplitudes are shown in Figs. 6,7 and 8 for the different values of $j_{\chi}$ and compared to the measured ones Bernabei et al. (2018). 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_{\sigma}$. 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 ${\bf c}_{0}\equiv c_{0}{\bf\hat{c}}_{0}$ with $c_{0}$ a free normalization, $c_{0}$ is common to all signals inside the eigenspace of ${\bf\hat{c}}_{0}$. A convenient parameterization is through the introduction of the reference cross section:

[TABLE]

with $\mu_{\chi{\cal N}}$ the WIMP–nucleon reduced mass. The direction in coupling space singled out by the unit vector ${\bf\hat{c}}_{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 Feng et al. (2011). Once ${\bf\hat{c}}_{0}$ and $\delta$ are fixed the DAMA signal and the bounds can be discussed in a familiar $m_{\chi}$–$\sigma_{0}$ plane. This is done in Figs. 9, 10 and 11 for the extreme configurations summarized in Table 2 and for $j_{\chi}$=1, 1/2 and 1, respectively. In the same figure the (red) circle represents the point $(m_{\chi,0},\sigma_{0,max})$ with $m_{\chi,0}$ the value of the WIMP mass for the extreme configuration and $\sigma_{0,max}=c^{2}_{0,max}\frac{\mu_{\chi{\cal N}}^{2}}{\pi}$, $c_{0,max}\equiv|{\bf c}_{0,max}|$. The point $(m_{\chi,0},\sigma_{0,max})$ intersects one or more of the most constraining boundaries on $\sigma_{0}$ at $m_{\chi}$=$m_{\chi,0}$, providing a nice confirmation of the numerical solution of Eq. (25). Actually, this can be directly observed in Figs. 9, 10 and 11 for for $j_{\chi}$=0, 1/2 and 1, respectively, where the point $(m_{\chi,0},\sigma_{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. 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_{\sigma}$ of the tension between DAMA and the constraints from null results. Indeed, as already observed in Kang et al. (2019), the condition (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}^{*}\equiv E_{R}(v^{*Na}_{min})$=$|\delta|\mu_{\chi 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 shape111The DAMA/LIBRA-phase1 data showed instead a maximum, and for this reason inelastic scattering could provide a good fit to the data Scopel and Yoon (2016). closer to elastic scattering. As a consequence, the DAMA data pull to low values of $\delta$. However, the solutions of Eq. (25) with smallest tension with the constraints require sizeable values of $\delta$ (\delta\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 20keV) in order to verify Eq. (5). As a consequence, when $\delta$ is fixed to such values the DAMA data pull to higher values of the WIMP mass $m_{\chi}$ to dilute the effect of $\delta$. This explains why, systematically, the DAMA regions in Figs. 9,  10 and  11 are at higher WIMP masses compared to the values of $m_{\chi,0}$ in Table 2. 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 2 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. Kang et al. (2019), obtained 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 3, such as background subtraction or the optimal-interval method Yellin (2002). In Figs. 9, 10 and 11 all this standard machinery Kang et al. (2018c) can instead be applied to the full set of existing experiments, providing an a posteriori confirmation that the set of bounds $\varmathbb{E}$ 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 Agnese et al. (2018b), CRESST-II Angloher et al. (2016); Angloher et al. (2017a), the upper bound from the average count rate of DAMA Bernabei et al. (2008b)), DarkSide–50 Agnes et al. (2018) and the $CF_{3}I$ target run of PICO–60 Amole et al. (2015) (the details of such bounds implementation are provided in Appendix B). None of these additional null results further constrains the extreme configurations ($m_{\chi,0},\sigma_{0,max}$).

We conclude by noting that the recent bound from COSINE–100 Adhikari et al. (2018a), 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 is that the bound of Ref. Adhikari et al. (2018a) is on the time–averaged signal $S_{0}^{COSINE}$, while DAMA measures the yearly modulation amplitudes expected from the rotation of the Earth around the Sun $S_{m}^{DAMA}$. In appendix B we estimate in COSINE–100 a residual count rate $b\simeq$0.13 events/kg/day/keV after background subtraction while modulation fractions in DAMA are of the order of 0.02 events/kg/day/keV. The bound $S_{0}^{COSINE}<b$ implies $S_{m}^{DAMA}/S_{0}^{DAMA}$ = $S_{m}^{DAMA}/S_{0}^{COSINE}\times S_{0}^{COSINE}/S_{0}^{DAMA}$ \mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 0.12, including a factor $S_{0}^{COSINE}/S_{0}^{DAMA}\simeq$ 0.8 due to a difference between the energy resolutions and efficiencies in the two experiments. For a standard Maxwellian WIMP velocity distribution in the SI elastic case the predicted modulation fractions $S_{m}^{DAMA}/S_{0}^{DAMA}$ are below such bound (for instance, for $m_{\chi}$=10 GeV $S_{m}^{DAMA}/S_{0}^{DAMA}$ is between $\simeq$0.05 and $\simeq$0.12 for $E_{ee}<$ 3 keVee) explaining why in Ref. Adhikari et al. (2018a) 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 Tucker-Smith and Weiner (2001); Scopel and Yoon (2016), and this is particularly true when the condition (5) is verified (for instance, in the extreme configurations of Table 2 we find S_{m}^{DAMA}/S_{0}^{DAMA}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}0.8 for $E^{\prime}\leq$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 Bernabei et al. (2008a, 2010, 2013, 2018) with an explanation in terms of inelastic scattering events induced 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. Catena et al. (2016). 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 Aprile et al. (2018a), PICO–60 ($C_{3}F_{8}$ target) Amole et al. (2017, 2019), COSINE–100Adhikari et al. (2018a), COUPP Behnke et al. (2012), SuperCDMS Agnese et al. (2018a) and PICASSO Behnke et al. (2017). The tension, quantified as the maximum of $N_{\sigma}\equiv\max(n_{\sigma})$ among the DAMA energy bins below 5 keVee, exceeds 7$\sigma$ in all the parameter space $m_{\chi}<$ 200 GeV with the exception of a small region of parameter space for $m_{\chi}\simeq$ 10 GeV and \delta\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}20, where it drops to values as low as $\simeq$ 4$\sigma$ for $j_{\chi}$=0, $\simeq$ 2.9$\sigma$ for $j_{\chi}$=1/2 and $\simeq$ 3.2$\sigma$ for $j_{\chi}$=1, and that overlaps to the proton–philic spin–dependent inelastic Dark Matter (pSIDM) scenario Scopel and Yoon (2016); Scopel and Yu (2017); Kang et al. (2019) 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_{\sigma}$ is partially relieved can easily produce enough yearly modulation in the lowest–energy 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.

Acknowledgements.

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, grant number 2016R1D1A1A09917964.

Note added

After the submission of the present paper COSINE–100 has released its first annual modulation analysis, consistent at 68.3% C.L. with both a null hypothesis and DAMA/LIBRA’s 2–6 keVee best–fit value Adhikari et al. (2019).

Appendix A WIMP response functions

We collect here the WIMP particle–physics response functions introduced in Eq.(13) and for the general case of complex couplings Barello et al. (2014) (although in the present analysis real couplings have been assumed). For a WIMP particle of spin $J_{\chi}\leq\frac{1}{2}$ they are given by Fitzpatrick et al. (2013); Anand et al. (2014):

[TABLE]

On the other hand, for a WIMP particle with spin $J_{\chi}$=1 Dent et al. (2015):

[TABLE]

Appendix B Experiments

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. In the solution of Eq. (24) we have used the 8 DAMA modulation amplitudes for 1 keVee $\leq E^{\prime}\leq$ 5 keVee Bernabei et al. (2018), and selected energy bins from XENON1T Aprile et al. (2018a), PICO–60 ($C_{3}F_{8}$ target) Amole et al. (2017, 2019), COSINE–100Adhikari et al. (2018a), COUPP Behnke et al. (2012), SuperCDMS Agnese et al. (2018a) and PICASSO Behnke et al. (2017), as shown in Table 3.

On the other hand, in the exclusion plots in Figs. 9, 10 and 11 we have included an extensive set of constraints that are representative of the different techniques used to search for DM: XENON1T Aprile et al. (2018a), CDMSlite Agnese et al. (2018b), SuperCDMS Agnese et al. (2018a), PICASSO Behnke et al. (2017), PICO–60 (using a $CF_{3}I$ target  Amole et al. (2015) and a $C_{3}F_{8}$ one Amole et al. (2017, 2019)), CRESST-II Angloher et al. (2016); Angloher et al. (2017a), DAMA (modulation data Bernabei et al. (1998, 2008a, 2010, 2018) and average count rate Bernabei et al. (2008b)), DarkSide–50 Agnes et al. (2018). providing an a posteriori confirmation that the limited set of bounds in Table 3 used to solve Eq. (25) did not miss any relevant constraint.

In the following, if not specified otherwise we adopt for the energy resolution a Gaussian form, ${\cal G}(E^{\prime},E_{ee})=Gauss(E^{\prime}|E_{ee},\sigma)=1/(\sqrt{2\pi}\sigma)exp(-(E^{\prime}-E_{ee})/2\sigma^{2})$. The quenching factor of bolometers (SuperCDMS, CRESST-II) is assumed to be equal to 1.

B.1 Xenon: XENON1T

For XENON1T we have assumed 7 WIMP candidate events in the range of 3PE $\leq S_{1}\leq$ 70PE, as shown in Fig. 3 of Ref. Aprile et al. (2018a) 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 Aprile et al. (2018a) 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 $\langle n_{ph}\rangle/E$ from Fig. 7 in Aprile et al. (2018b) with an electric field of $90~{}{\rm 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 $\sigma_{PMT}=0.5$, so that:

[TABLE]

with $Poiss(n,\lambda)=\lambda^{n}/n!exp(-\lambda)$.

B.2 Argon: DarkSide-50

The analysis of DarkSide-50 Agnes et al. (2018) 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 Agnes et al. (2018), 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.Agnes et al. (2018) we have subtracted the background minimizing the likelihood function:

[TABLE]

where $i$ represents the energy bin, $x_{i}$ the measured spectrum with error $\sigma_{i}$, while $\sigma S_{i}$ and $\rho b_{i}$ are the DM signal and the background, respectively, with $\sigma$ and $\rho$ arbitrary normalization factors ($\sigma$ is identified with the effective WIMP-proton cross section $\sigma_{p}$). In particular we obtain the 90% C.L. upper bound on $\sigma_{p}$ by taking its profile likelihood with $-2{\cal L}-[-2{\cal L}]_{min}=n^{2}$ and $n$=1.28. We take $x_{i}$, $\sigma_{i}$ and $b_{i}$ from Fig.7 of Agnes et al. (2018). The ionization yield of argon has been measured only down to \mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 10 keVnr, while DS50 uses a model fit to calibration data. We use the latter as taken from Fig. 6 of Agnes et al. (2018) with a hard cut at 0.15 keVnr, the lowest energy for which it is provided. We take the efficiency from Fig. 1 of Agnes et al. (2018).

B.3 Germanium: SuperCDMS and CDMSlite

The latest SuperCDMS analysis Agnese et al. (2018a) observed 1 event between 4 and 100 keVnr with an exposure of 1690 kg days. We have taken the efficiency from Fig.1 of Agnese et al. (2018a) and the energy resolution $\sigma=\sqrt{0.293^{2}+0.056^{2}E_{ee}}$ from Ahmed et al. (2010). To analyze the observed spectrum we apply the optimal interval method Yellin (2002).

For CDMSlite we considered the energy bin of 0.056 keV$<E^{\prime}<$ 1.1 keV with a measured count rate of 1.1$\pm$0.2 [keV kg day]-1 (Full Run 2 rate, Table II of Ref. Agnese et al. (2018b)). We have taken the efficiency from Fig.4 of Agnese et al. (2018b) and the energy resolution $\sigma=\sqrt{\sigma_{E}^{2}+BE_{R}+(AE_{R})^{2}}$, with $\sigma_{E}$=9.26 eV, $A$=5.68$\times 10^{-3}$ and $B$=0.64 eV from Section IV.A of Agnese et al. (2018b).

B.4 Fluorine: COUPP, PICASSO and PICO–60

Bubble chambers are threshold experiments for which we employ the nucleation probability:

[TABLE]

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 4. For fluorine and carbon we use $\alpha$=0.15 in Eq.(33). For iodine we adopt instead a step function with nucleation probability equal to 1 above the energy threshold.

The PICASSO experiment Behnke et al. (2017) 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 5 (extracted from Fig. 4 of Ref. Behnke et al. (2017)). For the nucleation probability we used Eq.(33) with $\alpha_{C}$=$\alpha_{F}$=5.

One of the target materials used by PICO–60 is $C_{3}F_{8}$, for which we used the complete exposure Amole et al. (2019) consisting in 1404 kg day at threshold $E_{th}$=2.45 (with 3 observed candidate events and 1 event from the expected background, implying an upper bound of 6.42 events at 90%C.L. Feldman and Cousins (1998)) and 1167 kg day keV at threshold $E_{th}$=3.3 keV (with zero observed candidate events and negligible expected background, implying a 90% C.L. upper bound of 2.3 events). For the two runs we have assumed the nucleation probabilities in Fig. 3 of Amole et al. (2019).

B.5 Fluorine+Iodine: PICO–60

PICO–60 can also employ a $CF_{3}I$ target. For the analysis of Ref.Amole et al. (2015) 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 Amole et al. (2015).

B.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 Bernabei et al. (2008b) (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*-1day-1keVee-1* in the energy range 1 keVee $<E^{\prime}<$ 8 keVee from Ref.Bernabei et al. (2018). In both cases we assume a constant quenching factors $q$=0.3 for sodium and $q$=0.09 for iodine, and the energy resolution $\sigma$ = 0.0091 (Eee/keVee) + 0.448 $\sqrt{E_{ee}/{\rm keVee}}$ in keV.

The exclusion plot for COSINE–100 Adhikari et al. (2018a) relies on a Montecarlo Adhikari et al. (2018b) to subtract the different backgrounds of each of the eight crystals used in the analysis. In Ref. Adhikari et al. (2018a) 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 reproduce the exclusion plot in Fig.4 of Ref. Adhikari et al. (2018a) for the isoscalar spin-independent elastic case. The result of our procedure yields $b\simeq$0.13 events/kg/day/keVee, which implies a subtraction of about 95% of the background. We take the energy resolution $\sigma/\mbox{keV}=0.3171\sqrt{E_{ee}/\mbox{keVee}}+0.008189E_{ee}/\mbox{keVee}$ averaged over the COSINE–100 crystals cos and the efficiency for nuclear recoils from Fig.1 of Ref. Adhikari et al. (2018a). Quenching factors for sodium and iodine are assumed to be equal to 0.3 and 0.09 respectively, the same values used by DAMA.

B.7 $CaWO_{4}$: CRESST-II

CRESST-II measures heat and scintillation using $CaWO_{4}$ crystals. We considered the Lise module analysis from Angloher et al. (2016) with energy resolution $\sigma$=0.062 keV and detector efficiency from Fig. 4 of Angloher et al. (2017b). 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.