From direct detection to relic abundance: the case of proton-philic spin-dependent inelastic Dark Matter

TL;DR

This paper investigates how direct detection data can inform the relic abundance of proton-philic spin-dependent inelastic dark matter, revealing that only certain interactions align with thermal relic conditions, especially under non-Maxwellian velocity distributions.

Contribution

It extends the analysis of pSIDM to general spin-dependent interactions, assesses the impact of velocity distributions, and introduces a consistency test for thermal relic scenarios based on direct detection data.

Findings

Only ${ m O}_{spin}$ interaction can produce a thermal relic WIMP.

Degeneracy between local density and scattering cross section limits direct relic abundance inference.

Non-Maxwellian velocity distributions are favored for the thermal relic interpretation.

Abstract

We discuss strategies to make inferences on the thermal relic abundance of a Weakly Interacting Massive Particle (WIMP) when the same effective dimension-six operator that explains an experimental excess in direct detection is assumed to drive decoupling at freeze-out, and apply them to the proton-philic Spin-dependent Inelastic Dark Matter (pSIDM) scenario, a phenomenological set-up containing two states $\chi_1$ and $\chi_2$ with $m_{\chi_2}>m_{\chi_1}$ that we have shown in a previous paper to explain the DAMA effect in compliance with the constraints from other detectors. We update experimental constraints on pSIDM, extend the analysis to the most general spin-dependent momentum-dependent interactions allowed by non-relativistic Effective Field Theory (EFT), and consider for the WIMP velocity distribution in our Galaxy both a halo-independent approach and a standard Maxwellian. The problem of calculating the relic abundance by using direct detection data to fix the model parameters is affected by a strong sensitivity on $f(v)$ and by the degeneracy between the WIMP local density and the WIMP-nucleon scattering cross section. As a consequence, a DM direct detection experiment is not directly sensitive to the physical cut-off scale of the EFT, but on some dimensional combination that does not depend on the actual value of the relic abundance. However, such degeneracy can be used to develop a consistency test on the possibility that the WIMP is a thermal relic in the first place. When we apply it to the pSIDM scenario we find that only a WIMP with a standard spin-dependent interaction ${\cal O}_{spin}$ with quarks can be a thermal relic, for a galactic velocity distribution that departs from a Maxwellian. However all the $\chi_2$ states must have already decayed today, and this requires some additional mechanism besides that provided by the ${\cal O}_{spin}$ operator.

Tables

Table 1. Table 1 : Relativistic Effective Field Theories for the nucleon–IDM interaction having as a low–energy limit a generalized spin–dependent elastic scattering. We do not include models leading to an explicit velocity dependence. The quantities m 𝒩 subscript 𝑚 𝒩 m_{\cal N} , S → 𝒩 subscript → 𝑆 𝒩 \vec{S}_{\cal N} and q 𝑞 q represent the nucleon mass, the nucleon spin and the transferred momentum, respectively. Some of the interaction terms in the second column contain an arbitrary scale m M subscript 𝑚 𝑀 m_{M} to ensure correct dimensionality: in our analysis we will fix it to the nucleon mass m 𝒩 subscript 𝑚 𝒩 m_{\cal N} . Adapted from Table 1 of Ref. [ 37 ] .

	Relativistic EFT	Non-relativistic limit	${\sum }_i{𝒪}_i$	cross section scaling
${𝒪}_4^{AA}$	${\overline{\chi }}_1{\gamma }^{\mu }{\gamma }^5{\chi }_2\overline{𝒩}{\gamma }_{\mu }{\gamma }^5𝒩$+h.c.	$-4{\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{S}}_{𝒩}$	-4${𝒪}_4$	$W_{{\Sigma }^{\prime \prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)+W_{{\Sigma }^{\prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$
${𝒪}_9^{VA}$	${\overline{\chi }}_1{\gamma }^{\mu }{\chi }_2\overline{𝒩}{\gamma }_{\mu }{\gamma }^5𝒩$+h.c.	$+\frac{2}{m_{\chi }}i{\overrightarrow{S}}_{\chi }\cdot \left({\overrightarrow{S}}_{𝒩}\times \overrightarrow{q}\right)$	$\simeq 2\frac{m_{𝒩}}{m_{\chi }}{𝒪}_9$	$\simeq q^2{W}_{{\Sigma }^{\prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$
${𝒪}_9^{TA}$	${\overline{\chi }}_{1}i{\sigma }_{\mu \nu }\frac{q^{\nu }}{m_M}{\chi }_2\overline{𝒩}{\gamma }^{\mu }{\gamma }_5𝒩$+h.c.	$4i\left(\frac{\overrightarrow{q}}{m_M}\times {\overrightarrow{S}}_{\chi }\right)\cdot {\overrightarrow{S}}_{𝒩}$	$4\frac{m_{𝒩}}{m_M}{𝒪}_9$	$q^2{W}_{{\Sigma }^{\prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$
${𝒪}_9^{AT}$	${\overline{\chi }}_1{\gamma }^{\mu }{\gamma }_5{\chi }_2\overline{𝒩}i{\sigma }_{\mu \nu }\frac{q^{\nu }}{m_M}𝒩$+h.c.	$4i{\overrightarrow{S}}_{\chi }\cdot \left(\frac{\overrightarrow{q}}{m_M}\times {\overrightarrow{S}}_{𝒩}\right)$	$-4\frac{m_{𝒩}}{m_M}{𝒪}_9$	$q^2{W}_{{\Sigma }^{\prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$
${𝒪}_{10}^{SP}$	$i{\overline{\chi }}_1{\chi }_2\overline{𝒩}{\gamma }^5𝒩$+h.c.	$i\frac{\overrightarrow{q}}{m_{𝒩}}\cdot {\overrightarrow{S}}_{𝒩}$+h.c.	${𝒪}_{10}$	$q^2{W}_{{\Sigma }^{\prime \prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$
${𝒪}_6^{PP}$	${\overline{\chi }}_1{\gamma }_5{\chi }_2\overline{𝒩}{\gamma }^5𝒩$+h.c.	$-\frac{\overrightarrow{q}}{m_{\chi }}\cdot {\overrightarrow{S}}_{\chi }\frac{\overrightarrow{q}}{m_{𝒩}}\cdot {\overrightarrow{S}}_{𝒩}$	$-\frac{m_{𝒩}}{m_{\chi }}{𝒪}_6$	$q^4{W}_{{\Sigma }^{\prime \prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$
${𝒪}_6^{T^{\prime }{T}^{\prime }}$	${\overline{\chi }}_{1}i{\sigma }^{\mu \alpha }\frac{q_{\alpha }}{m_M}{\gamma }_5{\chi }_2\overline{𝒩}i{\sigma }_{\mu \beta }\frac{q^{\beta }}{m_M}{\gamma }_5𝒩$+h.c.	$4\frac{\overrightarrow{q}}{m_M}\cdot {\overrightarrow{S}}_{\chi }\frac{\overrightarrow{q}}{m_M}\cdot {\overrightarrow{S}}_{𝒩}$	$4\frac{m_{𝒩}^2}{m_M^2}{𝒪}_6$	$q^4{W}_{{\Sigma }^{\prime \prime }}^{\tau {\tau }^{\prime }}\left(q^2\right)$

Table 2. Table 2 : Maximal values in GeV of the effective scale Λ ~ p , e ​ x ​ p m ​ a ​ x subscript superscript ~ Λ 𝑚 𝑎 𝑥 𝑝 𝑒 𝑥 𝑝 \tilde{\Lambda}^{max}_{p,exp} within the parameter space region where the compatibility factors 𝒟 Halo indep. < 1 subscript 𝒟 Halo indep. 1 {\cal D}_{\mbox{Halo indep.}}<1 or 𝒟 M ​ a ​ x ​ w ​ e ​ l ​ l ​ i ​ a ​ n < 1 subscript 𝒟 𝑀 𝑎 𝑥 𝑤 𝑒 𝑙 𝑙 𝑖 𝑎 𝑛 1 {\cal D}_{Maxwellian}<1 for each of the non–relativistic operators of Eq.( 2.4 ).

Non–relativistic operator	Halo–independent	Maxwellian
${𝒪}_4$	163	20.4
${𝒪}_6$	23.0	2.79
${𝒪}_9$	27.5	3.27
${𝒪}_{10}$	3.01	0.32

Table 3. Table 3 : Δ ​ q A 𝒩 Δ subscript superscript 𝑞 𝒩 𝐴 \Delta q^{{\cal N}}_{A} estimations of the axial coupling Γ A superscript Γ 𝐴 \Gamma^{A} = γ μ ​ γ 5 superscript 𝛾 𝜇 subscript 𝛾 5 \gamma^{\mu}\gamma_{5} for q = u , d , s 𝑞 𝑢 𝑑 𝑠 q=u,d,s .

$\Delta q_A^{𝒩}$	$n$	$p$
u	-0.427	0.842
d	0.842	-0.427
s	-0.085	-0.085

Table 4. Table 4 : The same of Table 3 for Δ ​ q 5 𝒩 Δ subscript superscript 𝑞 𝒩 5 \Delta q^{{\cal N}}_{5} .

$\Delta q_5^{𝒩}$	$n$	$p$
u	-108.3	110.55
d	108.60	-107.17
s	-0.57	-3.37

Table 5. Table 5 : Different estimations of the Δ ​ q T p Δ subscript superscript 𝑞 𝑝 𝑇 \Delta q^{p}_{T} = Δ ​ q T ′ p Δ subscript superscript 𝑞 𝑝 superscript 𝑇 ′ \Delta q^{p}_{T^{\prime}} constants for the proton.The last line indicates the energy scale of the measurement/calculation. The corresponding quantities for the neutron are obtainable by assuming isospin symmetry.

$\Delta q_T^p$=$\Delta q_{T^{\prime }}^p$	HERMES[55]	COMPASS[56]	Dyson–Schwinger[57, 58]	Lattice[59]
u	0.57(21)	0.39${}_{-0.12}{}^{}{}_{}{}^{+0.18}$	0.55(8)	0.79(7)
d	-0.18(33)	-0.25${}_{-0.10}{}^{}{}_{}{}^{+0.30}$	-0.11(2)	-0.22(3)
Q2(GeV)	1.0	0.8	2	2

Equations

$$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_{cut}^{lab}<v_{min}^{\ast F},$$

$$\mathcal{H}$$

$$\mathcal{L}=\frac{c_p}{{\Lambda }^2}{\mathcal{O}}_k^{XY}=\frac{c_p}{{\Lambda }^2}\bar{\chi }{\Gamma }^X\chi \bar{N}{\Gamma }^{Y}N,$$

$$\mathcal{L}=\sum _q\frac{c_q}{{\Lambda }^2}\bar{\chi }{\Gamma }^X\chi \bar{q}{\Gamma }^{Y}q,$$

$${\Gamma }_{\gamma \gamma }={\left(\frac{16{\pi }^2{f}_{\pi }^2}{m_{\pi }^2}\right)}^2\frac{{\alpha }_{em}^2{\delta }^9}{512(315{\pi }^9)f_{\pi }^4}\frac{c_4^2}{{\Lambda }^4},$$

$${\Gamma }_{\gamma \gamma }=7.2\times 10^{-56}{\left(\frac{\delta }{\text{mbox}10keV}\right)}^9{\left(\frac{\text{mbox}10GeV}{\tilde{\Lambda }}\right)}^4\text{mbox}GeV,$$

$${\mathcal{L}}_{int,\nu }=\frac{c_{\nu }}{{\Lambda }^2}({\bar{\chi }}_2{\gamma }_{\mu }{\gamma }_5{\chi }_1)(\bar{\nu }{\gamma }^{\mu }{\gamma }_5\nu ),$$

$${\Gamma }_{\nu \nu }=\frac{11{\delta }^5}{768{\pi }^4}\frac{1}{{\tilde{\Lambda }}^4},$$

$${\Gamma }_{\nu \nu }=1.5\times 10^{-33}{\left(\frac{\delta }{\text{mbox}10keV}\right)}^5{\left(\frac{\text{mbox}10GeV}{\tilde{\Lambda }}\right)}^4\text{mbox}GeV,$$

$${\tilde{\eta }}_{0,1}(v_{min})=\frac{{\rho }_{\chi }}{m_{\chi }}{\sigma }_{ref}{\eta }_{0,1}(v_{min})=\frac{\xi {\rho }_{loc}}{m_{\chi }}{\sigma }_{ref}{\eta }_{0,1}(v_{min}),$$

$${\sigma }_{ref}=\frac{(c_k^p{)}^2}{{\Lambda }^4}{\mu }_{\chi \mathcal{N}}^2/\pi ,$$

$${\eta }_0(v_{min},t)$$

$${\eta }_0(v_{min})$$

$${\eta }_1(v_{min})$$

$$S^{[E_1,E_2]}(t)$$

$$S_{0,1}^{[E_1,E_2]}$$

$${\tilde{\eta }}_0(v_{min,2})\le {\tilde{\eta }}_0(v_{min,1})$$

$${\tilde{\Lambda }}_{p,exp}=\frac{\Lambda }{{\xi }^{1/4}|c_p{|}^{1/2}}={\left(\frac{{\mu }_{\chi \mathcal{N}}^2}{\pi (\xi {\sigma }_{ref}{)}_{exp}}\right)}^{\frac{1}{4}},$$

$$<{\overline{\stackrel{~}{\eta }}}_{0,1}>=\frac{{\int }_{v_{min}^{\ast }}^{\infty }d{v}_{min}{\tilde{\eta }}_{0,1}(v_{min}){\mathcal{R}}_{[E_1,E_2]}(v_{min})}{{\int }_{v_{min}^{\ast }}^{\infty }d{v}_{min}{\mathcal{R}}_{[E_1,E_2]}(v_{min})}=\frac{N_{exp}}{{\int }_{v_{min}^{\ast }}^{\infty }d{v}_{min}{\mathcal{R}}_{[E_1,E_2]}(v_{min})},$$

$${\tilde{\eta }}_1^{max}(v_{min})={\tilde{\eta }}_0^{max}(v_{min})=\frac{{\overline{\stackrel{~}{\eta }}}_{1,fit}^{DAMA}}{v_s}\theta (v_s-v_{min}),$$

$${\tilde{\Lambda }}_{p,exp}<{\tilde{\Lambda }}_{p,exp}^{max}={\left(\frac{{\rho }_{loc}{\mu }_{\chi \mathcal{N}}^2}{\pi {\overline{\stackrel{~}{\eta }}}_{fit}^{DAMA}{m}_{\chi }{v}_s}\right)}^{\frac{1}{4}}.$$

$$\widetilde{<\sigma v>}=\frac{1}{{\Lambda }^4}\sum _q{c}_q^2{\widetilde{<\sigma >}}_q,{\widetilde{<\sigma >}}_q=N_q\left[(a{)}_q+\frac{(b{)}_q}{2x_f}\right],$$

$$\xi =\frac{{\widetilde{<\sigma v>}}_0}{\widetilde{<\sigma v>}}=\frac{{\Lambda }^4{\widetilde{<\sigma v>}}_0}{{\sum }_q{c}_q^2{\widetilde{<\sigma v>}}_q}.$$

$${\tilde{\Lambda }}_{p,th}={\left(\frac{{\sum }_q{c}_q^2\widetilde{<\sigma v{>}_q}}{c_p^2\widetilde{<\sigma v{>}_0}}\right)}^{\frac{1}{4}}={\left(\frac{{\sum }_q{r}_q^2\widetilde{<\sigma v{>}_q}}{(c_p/c_u{)}^2\widetilde{<\sigma v{>}_0}}\right)}^{\frac{1}{4}}={\tilde{\Lambda }}_{p,{\Omega }_0}.$$

$${\tilde{\Lambda }}_{p,exp}^{max}>{\tilde{\Lambda }}_{p,{\Omega }_0}>{\tilde{\Lambda }}_{p,{\Omega }_0}^{min},$$

$$1>\xi >{\xi }^{min}=\frac{16m_{\chi }^4{\widetilde{<\sigma v>}}_0}{max{|}_{|c_q|\le 4\pi ,c_n=0}({\sum }_q{c}_q^2{\widetilde{<\sigma v>}}_q)},$$

$${\mathcal{D}}_{\text{mbox}Haloindep.}\equiv \underset{i\in \text{mbox}signal}{\max }\left(\frac{{\overline{\stackrel{~}{\eta }}}_{1i}^{DAMA}+{\sigma }_i}{{\min }_{j\le i}{\overline{\stackrel{~}{\eta }}}_{0j,lim}}\right),$$

$${\mathcal{D}}_{\text{mbox}Haloindep.}\to \max \left({\mathcal{D}}_{\text{mbox}Haloindep.},\frac{N_{k,i}^{expected}}{N_{k,i}^{bound}}\right).$$

$${\mathcal{D}}_{Maxwellian}\equiv \max \left[\underset{k,i}{\max }\left(\frac{{\sigma }_{DAMA}^{min}}{{\sigma }_{ki}^{bound}}\right),\frac{p_{min}}{p}\right].$$

$${\left(\frac{\Lambda }{30\text{mbox}GeV}\right)}^2\frac{4\pi }{|c_u|}to0.0pt\text{raise}0.86108pt<\text{hss}\text{lower}3.87495pt\sim 800.$$

Full text

1 Introduction

Weakly Interacting Massive Particles (WIMPs) are considered the most natural candidates to provide the Dark Matter (DM) in the halo of our Galaxy, and the search for their recoils off nuclear targets represents the most direct way to detect them. In lack of any experimental evidence of physics beyond the Standard Model at the Large Hadron Collider (LHC) the use of non–relativistic Effective Field Theories (EFT)[1, 2] to analyze WIMP direct search data has become popular in the literature[3, 4, 5]. In presence of a putative signal, the question arises of how to shed light on the ultraviolet completion of the EFT, or, at least, of how to use a given EFT to find correlations among different types of WIMP signals. Indeed, this has been discussed in the context of mono–jet signal at the LHC[6, 7], showing that the ensuing constraints at low WIMP masses can be competitive with those from direct searches. Such inferences however imply an extrapolation from the scale of WIMP-nucleus elastic scattering (keV) to the LHC scale (TeV), that becomes questionable if the relevant value of the EFT cut–off scale does not significantly exceed the momentum transfer of the involved physical processes[8, 9, 10]. Similarly, once an effective interaction operator is singled out from direct detection data the question arises if the same operator can yield an acceptable relic abundance through thermal decoupling[3]. Due to the non–relativistic nature of the latter process in this case the involved energy extrapolation is much smaller, since the energy scale relevant to WIMP annihilation at freeze–out is of order $2m_{\chi}$, with $m_{\chi}={\cal O}(GeV)$ the WIMP mass. As a consequence, the use of EFT to calculate the relic abundance appears more reliable, and represents the first necessary step upward in the bottom–up quest to shed light on the ultraviolet completion of the EFT.

In the present paper we wish to address this latter issue, starting from an explicit class of EFT low–energy models that we introduced in a previous paper[11] and suggested by present experimental data . Specifically, our goal is to develop operative bottom–up strategies to test the possibility that the same effective dimension–six operator that explains an experimental excess in WIMP direct detection (at the keV scale) can also yield an acceptable relic abundance through thermal decoupling (at the GeV scale). In particular, we will also address this issue when the WIMP velocity distribution in our Galaxy is not fixed.

Many underground experiments using different materials are strongly pursuing the task of WIMP direct detection, with an impressive improvement of the experimental sensitivities but somewhat confusing and apparently contradictory results. In particular, the DAMA experiment[12] has been measuring for more than 15 years a yearly modulation effect with a sodium iodide target, consistent with that expected from the elastic scattering of WIMPs due to the Earth rotation around the Sun. However, many experimental collaborations using nuclear targets different from $NaI$ and various background–subtraction techniques to look for WIMP–elastic scattering (LUX[13, 14, 15],PANDA[16, 17] XENON100[18], XENON10[19], KIMS[20, 21, 22], CDMS-$Ge$[23], CDMSlite [24], SuperCDMS[25], CDMS II[26], SIMPLE[27], COUPP[28], PICASSO[29, 30], PICO-2L[31, 32], PICO-60[33]) have failed to observe any anomaly so far. This implies severe constraints on the most popular WIMP scenarios used to explain the DAMA excess: a WIMP with scattering cross section off nuclei proportional to the square of the atomic mass number of the target or to the target nuclear spin, and with a galactic velocity distribution given by a Maxwellian. However, when the latter assumptions are relaxed it is possible to show that compatibility between an interpretation of the DAMA effect in terms of a WIMP signal and constraints from other experiment can be recovered [34, 35, 36, 37, 11].

In particular, in the proton–philic Spin–dependent Inelastic Dark Matter (pSIDM) scenario introduced in [11] the DAMA modulation effect is explained by a WIMP which upscatters inelastically[38] to a heavier state and predominantly couples to the spin of protons. In such scenario constraints from xenon and germanium targets are evaded dynamically, due to the suppression of the WIMP coupling to neutrons, while those from fluorine targets are evaded kinematically, 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. The scenario of Ref.[11] is purely phenomenological, and can be mainly considered as a proof of concept of the fact that the parameter space of WIMP direct detection is wider that the one usually assumed. However, it has also some attractive properties: for instance, it predicts in a natural way large modulation fractions compatible to the DAMA experimental data, and provides an explanation of the maximum in the energy spectrum of the modulation amplitude detected by DAMA in terms of WIMPs whose minimal incoming speed matches the kinematic threshold for inelastic upscatters.

In light of this, in the present paper we wish to elaborate more on such scenario, extending the discussion of Ref.[11] in several directions: i) we update experimental constraints including the latest results from LUX[14, 15], PANDA[16, 17] and PICO-2L[32]; ii) we extend the discussion to the most general class of spin–dependent interactions [37] making use of non–relativistic EFT[1, 2]; iii) we discuss the thermal relic density of such particle when its interaction with ordinary matter at the scale of WIMP annihilation in the early Universe is described by the same dimension–six operator responsible for elastic scattering.

The paper is organized as follows: in Section 2 we summarize the main features of the pSIDM scenario, extending it to the most general cases of a spin–dependent coupling (as listed in Table 1); in Section 3 we discuss the DM relic density calculation; in Section 4 we discuss the kind of inferences that an excess in direct detection data can allow on the EFT cut–off scale, showing that the effect of the rescaling of the local density when the DM particle turns out to be subdominant in the Universe limits our access on the latter, but allows to develop a robust halo–independent consistency check on the possibility that the the same effective DM model that explains the excess can fix the thermal density; in Section 5 we show the results of our numerical analysis for pSIDM in the specific class of effective models introduced in Section 2, and in Section 6 we provide our conclusions. In Appendix A we provide the analytic expressions for the thermal average of the WIMP annihilation cross section entering the relic abundance calculation; in Appendix B we summarize the procedure to relate the WIMP–nucleon coupling constants entering direct detection to the WIMP–quark couplings that appear in the calculation of the annihilation cross section; finally, in Appendix C we provide the details of our treatment of the latest experimental constraints from LUX[14, 15], PANDA[16, 17] and PICO-2L[32] (constraints from earlier experiments are implemented as described in Refs. [37, 11]).

2 The pSIDM scenario and DAMA

In this Section we briefly summarize the features of the scenario introduced in Ref.[11] (we refer the reader to such paper for further details).

The most stringent bounds on an interpretation of the DAMA effect in terms of WIMP–nuclei scatterings are obtained by detectors using xenon (LUX[13, 14, 15], PANDA[16, 17]) and germanium (CDMS[23, 24, 25, 26]) whose spin is mostly originated by an unpaired neutron, as well as by the KIMS experiment[20, 21, 22] which uses $CsI$ and thus directly probes the contribution to the DAMA effect from WIMP scatterings off iodine targets. If the WIMP mass is small enough to assume that the DAMA signal is only due to WIMP scatterings off sodium the KIMS constraint can be evaded. Moreover, both sodium and iodine in DAMA have an unpaired proton, so that if the WIMP particle interacts with ordinary matter predominantly via a spin–dependent coupling which is suppressed for neutrons it can explain the DAMA effect in compliance with the bounds from xenon and germanium detectors, whose constraints are strongly relaxed[39, 40]. However this scenario is constrained by droplet detectors (SIMPLE[27], COUPP[28]) and bubble chambers (PICASSO[29], PICO-2L[31, 32],PICO-60[33]) which all use nuclear targets with an unpaired proton (in particular, they all contain ${}^{19}F$, while SIMPLE contains also ${}^{35}Cl$ and ${}^{37}Cl$ and COUPP and PICO-60 use also ${}^{127}I$). As a consequence, this class of experiments rules out a DAMA explanation in terms of WIMPs with a spin–dependent coupling to protons when standard assumptions are made on the WIMP local density and velocity distribution in our Galaxy[40, 31].

In Ref.[11] the alternative approach of IDM was proposed to reconcile DAMA to fluorine detectors. In this class of models 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_{cut}$ the result of the boost in the lab rest frame of some maximal speed value beyond which the WIMP velocity distribution $f(v)$ in the galactic rest frame vanishes (typically $v_{cut}$ is identified with the WIMP escape velocity $v_{esc}$), constraints from droplet detectors and bubble chambers can 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, since in such case WIMP scatterings off fluorine turn kinematically impossible while those off sodium can still serve as an explanation to the DAMA effect. Clearly, this mechanism rests on the trivial observation that the velocity $v_{min}^{*}$ for fluorine is larger than that for sodium.

In Ref. [11] only a standard spin–dependent coupling between the WIMP and the nucleon ${\cal N}=p,n$, ${\cal L}_{int}\ni c^{p}\vec{S}_{\chi}\cdot\vec{S}_{p}+c^{n}\vec{S}_{\chi}\cdot\vec{S}_{n}$ with $c_{n}\ll c_{p}$ was considered (the exact value of $c_{n}/c_{p}$ must be tuned to a different small number depending on the spin–dependent form factor that is used [11]). In the present paper we keep the condition $c_{n}\ll c_{p}$ and extend that analysis to the most general class of spin–dependent interactions[37] by making use of the interaction Hamiltonian which descends from non–relativistic EFT[1, 2]:

[TABLE]

where $t^{0}=1$, $t^{1}=\tau_{3}$ denote the the $2\times 2$ identity and third Pauli matrix in isospin space, respectively, the non–dimensional isoscalar and isovector coupling constants $c^{0}_{k}$ and $c^{1}_{k}$ are related to those to protons and neutrons $c^{p}_{k}$ and $c^{n}_{k}$ by $c^{p}_{k}=(c^{0}_{k}+c^{1}_{k})/2$ and $c^{n}_{k}=(c^{0}_{k}-c^{1}_{k})/2$, $\Lambda$ is a dimension–1 cut–off scale and the operators ${\cal O}_{i}$ are for instance listed in Equations (12) and (13) of [2]: the corresponding interaction operators with an exclusively spin–dependent nuclear response function (i.e. depending only on either the function $W^{\tau\tau^{\prime}}_{\Sigma^{\prime}}$ or $W^{\tau\tau^{\prime}}_{\Sigma^{\prime\prime}}$ in the notation of [11]) are summarized in Table 1, where we follow the naming convention $\bar{\chi}_{1}\Gamma^{a}\chi_{2}\bar{{\cal N}}\Gamma^{b}{\cal N}$+h.c.$\rightarrow{\cal O}_{i}^{ab}$ with $\Gamma^{S}=1$, $\Gamma^{P}=\gamma_{5}$, $\Gamma^{V}=\gamma^{\mu}$, $\Gamma^{A}=\gamma^{\mu}\gamma_{5}$, $\Gamma^{T}$=$\sigma^{\mu\nu}$=$i/2[\gamma^{\mu},\gamma^{\nu}]$, $\Gamma^{T^{\prime}}$=$\sigma^{\mu\nu}\gamma_{5}$. Following this convention the standard spin–dependent coupling corresponds to $\bar{\chi}_{1}\gamma^{\mu}\gamma_{5}\chi_{2}\bar{{\cal N}}\gamma_{\mu}\gamma_{5}{\cal N}$ = $\bar{\chi}_{1}\Gamma^{A}\chi_{2}\bar{{\cal N}}\Gamma^{A}{\cal N}$+h.c.$\rightarrow\vec{S}_{\chi}\cdot\vec{S}_{{\cal N}}\equiv{\cal O}_{4}^{AA}$ with $\Gamma^{A}=\gamma^{\mu}\gamma_{5}$ and ${\cal O}_{4}$=$\vec{S}_{\chi}\cdot\vec{S}_{{\cal N}}$.

3 Relic abundance

A particularly attractive feature of the WIMP scenario is that the interactions that keep the DM and Standard Model particles in thermal equilibrium in the early Universe and that drive the decoupling process fixing the DM relic abundance can be correlated to those which may allow their direct or indirect detection. In particular, in Section 4 we will adopt the “bottom–up” procedure of fixing the effective Hamiltonian of Eq. (2.4) by requiring that it explains the observed DAMA modulation signal, and using it to calculate the relic abundance of particles $\chi_{1}$ and $\chi_{2}$. The total density of the two $\chi_{1}$, $\chi_{2}$ states and the ratio between their abundances depend on the annihilation process $\chi_{1}\chi_{2}\rightarrow\bar{f}f$ of the DM particles to Standard Model fermions $f$ and the decay amplitude for the process $\chi_{2}\rightarrow\chi_{1}$. In order to calculate both processes we make the assumption that the physics of annihilation and decay is dominated at the quark level by the same 6–dimension operator whose non–relativistic limit drives direct detection. As summarized in Table 1 the operators ${\cal O}_{k}$ of the effective Hamiltonian (2.4) can be interpreted as the non–relativistic limits of the interaction terms:

[TABLE]

with XY=AA,VA,TA,AT,SP,PP,T’T’. The above interaction is generated by an analogous WIMP–quark Lagrangian:

[TABLE]

(the relation between the couplings $c_{p}$ and $c_{q}$ is discussed in Appendix B). As far as the decay process is concerned, the Lagrangian of Eq.(3.2) induces the decay $\chi_{2}\rightarrow\chi_{1}\gamma\gamma$ through a quark loop. The scale of the process is of order $\delta~{}{\cal O}(keV)$, and Chiral Perturbation Theory (ChPT) can be used to evaluate it. For the case of the ${\cal O}_{4}^{AA}$ operator it is explicitly given by (see Eq.(3) of Ref.[41]):

[TABLE]

with $m_{\pi}\simeq$ 140 MeV, $f_{\pi}\simeq$ 93 MeV. Substituting numbers one gets:

[TABLE]

with $1/\tilde{\Lambda}^{2}=c_{4}/\Lambda^{2}$. This implies that for the ranges of parameters we will find in Section 5 (\tilde{\Lambda}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 100 GeV, 10 keV$\leq\delta\leq$ 20 keV) the lifetime of the $\chi_{2}$ particle is far larger than the age of the Universe (corresponding in natural units to 1.5$\times$10*-42* GeV*-1*) and for practical purposes $\chi_{2}$ is stable. The same happens for the decay amplitudes calculated for the other models listed in Table 1. In particular, the decay amplitude of model ${\cal O}_{9}^{VA}$ turns out to be the same as that of ${\cal O}_{4}^{AA}$, while models with tensor currents, i.e. ${\cal O}_{9}^{TA}$, ${\cal O}_{9}^{AT}$ and ${\cal O}_{6}^{T^{\prime}T^{\prime}}$ lead to vanishing decay amplitudes at leading order in ChPT; finally, the decay amplitudes for models ${\cal O}_{10}^{SP}$ and ${\cal O}_{6}^{PP}$ turn out to be enhanced by a factor $(f_{\pi}/\delta)^{2}$ and suppressed by a factor $f_{\pi}/m_{\chi}$ compared to Eq. (3.3), respectively, with no practical phenomenological differences with the other models. The bottom line is that for all the models of Table 1, if no other coupling is assumed for the DM particle, our Galaxy is formed by a mixture of $\chi_{1}$ and $\chi_{2}$, which, presumably, are in equal parts and have the same velocity distribution, and direct detection experiments are both sensitive to the upscatters of $\chi_{1}$’s and the downscatters of $\chi_{2}$’s.

This is a problem in the scenario outlined in Section 2, where IDM upscatters off fluorine must be kinematically suppressed in order to comply with the constraints from bubble chambers and droplet detectors, and which requires $v_{min}^{*}>v_{esc}$, because $v_{min}^{*}$ vanishes when $\delta<0$, i.e. for the downscatters of the $\chi_{2}$ particles the condition (2.3) is never achieved. A simple numerical check using the compatibility factors introduced in Section 5 confirms that the downscatters of $\chi_{2}$ rule out such scenario experimentally. A possible way out is adding to the theory a tree–level coupling to neutrinos. In such case, assuming for instance:

[TABLE]

the decay amplitude to standard left–handed neutrinos turns out to be:

[TABLE]

i.e, numerically:

[TABLE]

(where $\tilde{\Lambda}^{2}=\Lambda^{2}/|c_{\nu}|$), so that the $\chi_{2}$ particle can safely decay after decoupling from the plasma and the halo of our Galaxy is only made of $\chi_{1}$ states. Moreover, we notice that Eq.(3.5) allows for a fast decay of $\chi_{2}$ already for $|c_{\nu}|<<|c_{q}|$; in such a case the DM thermal decoupling from the plasma and its relic density are only driven by the couplings to quarks responsible of direct detection, and the correlation between direct detection rates and the relic density is maintained.

In the following we will assume the presence of the coupling of Eq.(3.5) with $|c_{\nu}|<<|c_{q}|$ or an analogous mechanism that depletes the $\chi_{2}$ particles without affecting the relic density. In such case the DM particle relic abundance is given by a standard expression inversely proportional to the quantity $\tilde{<\sigma v>}=a+b/(2x_{f})$, where $<\sigma v>=a+b/x$ is the thermal average of the coannihilation cross section at temperature $T=m_{\chi}/x$, and $x_{f}\simeq 20$ corresponds to the freeze-out temperature (in Appendix A the coefficients $a$ and $b$ for the different models listed in Table 1 are provided). For minimality, we will only assume a coupling of the DM particles to the quarks that drive spin–dependent direct detection, i.e. $f$=$u$, $d$, $s$ for models ${\cal O}^{AA}_{4}$, ${\cal O}^{VA}_{9}$, ${\cal O}^{TA}_{9}$, ${\cal O}^{SP}_{10}$, ${\cal O}^{PP}_{6}$ and $f$=$u$, $d$ for models ${\cal O}^{AT}_{9}$ and ${\cal O}^{TT}_{6}$ (see Appendix B). Combined with the condition $|c_{\nu}|<<|c_{q}|$, this implies that all the quantitative values of the relic abundance discussed in Section 5 will be be upper bounds. Indeed, we will see in Section 5 that this represents the most favorable choice to explain the DAMA effect and get at the same time an acceptable relic abundance through thermal decoupling in a consistent way. In our numerical analysis we will not make any further assumptions on the $c_{q}$ parameters, besides imposing perturbativity ($|c_{q}|<4\pi$) and $c_{n}\ll c_{p}$, and we will neglect the effect of the running of the couplings[42]222In some cases (such as model VA in Table 1[43, 44, 45]) the ultraviolet completion of our scenario requires some level of tuning at high energy, since the running of the couplings generates a non–vanishing SI interaction at low energy.. We conclude this section by pointing out that, obviously, the scenario outlined above does not produce any indirect detection signal.

4 The bottom–up approach

In the present paper the data of DAMA and of the other direct detection experiments will be used to determine estimations or upper bounds of the halo functions $\tilde{\eta}_{0}$, $\tilde{\eta}_{1}$ defined as:

[TABLE]

with $\rho_{loc}$=0.3 GeV/cm2 the observed total DM density in the neighborhood of the Sun, $\rho_{\chi}$=$\xi\rho_{loc}$ the corresponding local density of the $\chi_{1}$ particles (contributing in general a fraction $\xi$ of $\rho_{loc}$), $\sigma_{ref}$ represents a reference WIMP–nucleon cross section that we take as:

[TABLE]

for a given operator ${\cal O}_{k}$, while:

[TABLE]

where the time–dependence of the velocity distribution $f(v,t)$ descends from the rotation of the Earth around the Sun with $\omega=2\pi/(\mbox{365 days})$ and phase $t_{0}\simeq$2 June, and is due to the boost from the galactic to the Earth rest frame, while $v_{min}$ is given by Eq.(2.2). The halo functions of Eqs. (4.4,4.5) are defined in such a way that, for a given experimental setup, the corresponding expected direct detection signal $S^{[E_{1},E_{2}]}(t)$ in the observed energy interval $[E_{1},E_{2}]$ is given by[46, 47, 48]:

[TABLE]

In particular, the response functions ${\cal R}^{[E_{1},E_{2}]}(v_{min})$ are supposed to be known since they depend only on the properties of the detector, while the functions $\tilde{\eta}_{0,1}$ depend both on particle physics (through $\sigma_{ref}$) and on astrophysics (through $\rho_{\chi}$ and $f(v)$) and factorize all our ignorance about the process. Explicit expressions for the response functions ${\cal R}$ for the different experiments considered in this paper are given in [35, 37].

If no assumptions are made about $f(v)$ the two halo functions $\tilde{\eta}_{0}$ and $\tilde{\eta}_{1}$ are subject to the very general conditions:

[TABLE]

The first condition descends from the definition (4.3), that implies that $\tilde{\eta}_{0}(v_{min})$ is a decreasing function of $v_{min}$. The second is a consequence of the fact that $\tilde{\eta}_{1}$ is the modulated part of $\tilde{\eta}_{0}$333The use of representation theorems for distribution functions of the type (4.6) may allows to restrict further the range of the ratio $\tilde{\eta}_{1}/\tilde{\eta}_{0}$ in a statistical way even without specifying the velocity distribution [49].. The last condition with $v_{s}$= $v_{esc}$ reflects the requirement that the WIMPs are gravitationally bound to our Galaxy. We will adopt this choice in our Maxwellian analysis and, in this case, take as the escape velocity of WIMPs in the lab rest frame $v_{esc}$=$v_{esc}^{Galaxy}$+$v_{Sun}$, with $v_{esc}^{Galaxy}$=550 km/sec the escape velocity in the galactic rest frame and $v_{Sun}$=232 km/sec the velocity of the Solar system with respect to the WIMP halo. On the other hand, in our halo–independent analysis $v_{s}$ will represent the maximal value of the $v_{min}$ range corresponding to the DAMA excess. Not that in this case $v_{s}$ will depend on the model parameters $m_{\chi}$ and $\delta$: as explained below it will be adopted in order to get the largest $\eta_{1}(v_{min})$ which does not vanish in the signal range.

Intuitively, the measurement of a signal in a direct detection experiment such as that of DAMA should allow to determine the WIMP–nucleon cross section $\sigma_{ref}$ and in this way the EFT cut–off scale $\Lambda$, shedding light on the scale of the new physics involved in the process. Notice, however, that Eq.(4.1) implies that direct detection allows to get access to the product $\xi\sigma_{ref}$ rather then $\sigma_{ref}$. As a consequence, once an experimental estimate $(\xi\sigma_{ref})_{exp}$ is obtained from the data, using Eq.(4.2) any direct detection experiment is sensitive to the effective scale:

[TABLE]

rather than directly to the physical cut–off scale $\Lambda$.

The firs step to measure $(\xi\sigma_{ref})_{exp}$ is to get experimental estimates of the halo functions $\tilde{\eta}_{0}$ and/or $\tilde{\eta}_{1}$ from the data. In particular, given an experiment with detected count rate or modulation amplitude $N_{exp}$ in the energy interval $E_{1}<E<E_{2}$ the combination[46]:

[TABLE]

can be interpreted as an average of the halo function $\tilde{\eta}_{0,1}(v_{min})$ in an interval $v_{min,1}<v_{min}<v_{min,2}$. The $v_{min}$ interval is defined as the one where the response function ${\cal R}$ is “sizeably” different from zero, including the smearing effect of energy resolution.

Getting $(\xi\sigma_{ref})_{exp}$ from the $<\bar{\tilde{\eta}}_{0,1}>$’s requires to fix the velocity distribution $f(v)$. In this way the function $\eta_{1}(v_{min})$ is known and $\xi\sigma_{ref}$, which is just a normalization factor of the ensuing $\tilde{\eta}_{1}(v_{min})$, can be directly fitted from the data. If, on the other hand, a halo–independent approach is assumed, since $\int f(\vec{v})d^{3}v$=1 and $\tilde{\eta}_{1}\leq\tilde{\eta}_{0}$ the function $\eta_{1}(v_{min})$ can be at least maximized[36] by the choice $\tilde{\eta}_{1}$ =$\tilde{\eta}_{0}$ and $f(\vec{v})=\delta(v_{s}-v_{min})$, with $v_{s}$ the maximal value of the $v_{min}$ range corresponding to the DAMA excess. This corresponds to the largest $\eta_{1}(v_{min})$ which does not vanish in the signal range, and allows to get a lower bound on $(\xi\sigma_{ref})_{exp}$. In this case:

[TABLE]

where the constant value $\bar{\tilde{\eta}}_{1,fit}^{DAMA}$ can be fitted from the DAMA observed modulation amplitudes in a straightforward way444Indeed, due to the large error bars such flat functional form for the $\tilde{\eta}_{1}^{max}(v_{min})$ function is in general not incompatible with the DAMA experimental data.. Plugging Eqs.(4.10), (4.2) and (4.1) in Eq. (4.8) one then gets for the experimental estimate $\tilde{\Lambda}_{p,exp}$ of $\tilde{\Lambda}_{p}$ the upper bound:

[TABLE]

Notice that the halo–independent procedure allows to determine only a lower bound on $\sigma_{ref}$, and so an upper bound on $\tilde{\Lambda}_{p}$.

The quantity above, which is obtained analyzing experimental data, must be compared to the corresponding theoretical expectation for a given model ${\cal O}^{XY}_{k}$ and choice of the physical cut–off scale $\Lambda$ and of the couplings $c_{q}$ (related to $c_{p}$ through Eq.(B.2)). A crucial observation here is that, whenever the DM particle calculated relic density $\Omega$ turns out to be smaller than the observed value $\Omega_{0}$, $\rho_{\chi}$ must be rescaled with respect to the observed local DM density $\rho_{loc}$ by the factor $\xi\equiv\Omega/\Omega_{0}$. To get an expression of the rescaling factor $\xi$ we can express the coannihilation cross section in the form:

[TABLE]

with $(a)_{q}$ and $(b)_{q}$ given in Appendix A. Indicating with $\tilde{<\sigma v>}_{0}=2\times 10^{-9}$ GeV*-2* the value corresponding to the observed relic density, the rescaling coefficient $\xi$ is then given by:

[TABLE]

In this way the theoretical expectation $\tilde{\Lambda}_{p,th}$ of the effective scale turns out to be:

[TABLE]

Notice that in the expression above the dependence on the cut-off scale $\Lambda$ cancels out. Moreover, the second expression in parenthesis shows explicitly that $\tilde{\Lambda}_{p,th}$ depends only on the coupling ratios $r_{q}\equiv c_{q}/c_{u}$ and not on the actual values $c_{q}$. Finally, $\tilde{\Lambda}_{p,th}$ is equal to the quantity $\tilde{\Lambda}_{p,\Omega_{0}}$ that yields $\tilde{<\sigma v>}=\tilde{<\sigma v>}_{0}$, i.e. it corresponds to the observed relic abundance $\Omega_{0}$. This means that a direct detection experiment has in principle no direct access to the physical scale $\Lambda$ and is bound to measure $\tilde{\Lambda}_{p}=\tilde{\Lambda}_{p,\Omega_{0}}$ even if $\Omega\ll\Omega_{0}$555This is at variance with indirect detection where signals scale with $\rho_{loc}^{2}$ and are maximal when $\Omega=\Omega_{0}$[50], irrespective of the actual size of the couplings. Actually, the requirement $\tilde{\Lambda}_{p,exp}$=$\tilde{\Lambda}_{p,th}$ combined with Eqs.(4.11) and (4.14) implies:

[TABLE]

where $\tilde{\Lambda}_{p,\Omega_{0}}^{min}$ represents the minimum of $\tilde{\Lambda}_{p}$ yielding $\Omega_{0}$ when the parameters $r_{q}$ are varied at fixed $c_{n}$=0. In Fig.(1) $\tilde{\Lambda}_{p,\Omega_{0}}^{min}$ is plotted as a function of $m_{\chi}$ for the models of Table 1.

From Eq.(4.15) one can see that the energy scale inferred by the experimental data must yield $\Omega>\Omega_{0}$ (“overclosure”) to be acceptable666This might seem confusing but is simply due to the fact that in a halo–independent approach Eq.(4.10) is used to make inferences on $\tilde{\Lambda}_{p}$; on the other hand the energy scale inferred by the experimental data is supposed to yield $\Omega=\Omega_{0}$ for the correct choice of $f(v)$.. On the other hand, the only direct information on the relic abundance can be obtained by combining the requirement that the EFT of Eq. (3.2) is valid at the scale of the coannihilation process (implying $\Lambda>E_{\chi_{1}}+E_{\chi_{2}}\simeq 2m_{\chi}$) with perturbativity ($|c_{q}|<4\pi$). In this way, taking into account the condition $c_{n}$=0, the following lower bound on $\xi$ can be obtained:

[TABLE]

where in the last expression the denominator is maximized with the conditions $|c_{q}|\leq 4\pi$ and $c_{n}=0$. The quantity above is plotted in Fig.2 as a function of the WIMP mass $m_{\chi}$ for the models of Table 1 for $max(|c_{q}|)=4\pi$ but using instead the same $r_{q}$ coupling ratios employed in Fig.1. From such figure one can see that $\xi\ll$ 1 for all the effective models analyzed in the present paper. This has two consequences: (i) indeed, a positive signal in a direct detection experiment leaves a very large degeneracy on the actual value of the relic abundance $\Omega$; (ii) the constraints from EFT validity and perturbativity from Eq.(4.16) on the allowed parameter space are not significant.

The condition (4.15) is the main result of this Section and will be discussed in the numerical analysis of the pSIDM scenario of Section 5. It should be pointed out that it is not simply a requirement to obtain the correct relic abundance with the same parameters that explain a given experimental excess: instead, it should be interpreted as a more general condition on the possibility that the same effective model that explains such excess can also determine through thermal decoupling the relic abundance, irrespective on its actual value. It should also be pointed out that the degeneracy of direct detection signals on the physical cut–off scale $\Lambda$ discussed in this Section is in principle resolved by LHC signals and by indirect detection (although the specific case of the the scenario discussed in the present paper implies no indirect signals).

5 Results

In this Section we wish to determine the regions in the $m_{\chi}$–$\delta$ parameter space where an interpretation of the DAMA modulation effect in terms of a signal from the WIMPs of the pSIDM scenario outlined in Section 2 is compatible to constraints from other experiments. In order to do so we follow closely the analysis of Ref.[11] (to which we refer for further details) that we extend here to the generalized spin models listed in Table 1. Moreover, in the same parameter space regions we discuss the requirement (4.15), i.e. we show when the same effective model used to interpret direct detection data can be assumed to calculate the relic density through thermal decoupling in a consistent way.

Using Eq.(4.9) it is straightforward, for a given choice of the DM parameters, to obtain estimations $\bar{\tilde{\eta}}_{1i}^{DAMA}$ of the modulated halo function $\tilde{\eta}_{1}(v_{min})$ averaged in appropriately chosen $v_{min,i}$ intervals mapped from the DAMA experimental annual modulation amplitudes, and upper bounds $\bar{\tilde{\eta}}_{0j,lim}$ in the same $v_{min}$ intervals on the functions $\tilde{\eta}_{0}(v_{min})$ from the data of the experiments that have reported null results. Thanks to the condition $\tilde{\eta}_{1}(v_{min})\leq\tilde{\eta}_{0}(v_{min})$ this allows to constrain the DAMA result. In particular, we calculate the response functions ${\cal R}$ as explained in Ref.[11], with the nuclear form factors from [1, 2] and fix $c_{n}/c_{p}$=-0.03[11], which is the value that minimizes the response on xenon using the form factor of Ref. [1, 2] (this is equivalent to taking $c_{n}$=0 as far as Sections 3 and 4 are concerned).

Quantitatively, for a given choice of the WIMP mass $m_{\chi}$ and of the mass difference $\delta$ the compatibility between DAMA and all the other results can be assessed introducing the following compatibility ratio[11]:

[TABLE]

where $\sigma_{i}$ represents the standard deviation on $\bar{\tilde{\eta}}_{1i}^{DAMA}$ as estimated from the data, $i\in\mbox{signal}$ means that the maximum of the ratio in parenthesis is for $v_{min,i}$ corresponding to the DAMA excess, while, due to the fact that the function $\tilde{\eta}_{0}$ is non–decreasing in all velocity bins $v_{min,i}$, the denominator contains the most constraining bound on $\tilde{\eta}_{0}$ for $v_{min,j}\leq v_{min,i}$. The latter minimum includes all available bounds from scintillators, ionizators and calorimeters (see Appendix C of the present paper, Appendix B of Ref.[37] an Appendix A of [11] for a summary of the experimental inputs used in our analysis). Specifically, compatibility between DAMA and the constraints included in the calculation of Eq.(5.1) is ensured if ${\cal D}_{\mbox{Halo indep.}}<1$.

The above procedure cannot be applied to bubble chambers and droplet detectors, which contain different nuclear targets and, being only sensitive to the energy threshold, do not allow to map the corresponding bounds to arbitrary velocity bins. In this case, however, an alternative procedure can be applied, which consists to derive from the measured DAMA modulation amplitudes a piece–wise estimation $\tilde{\eta}^{est}_{0}(v_{min})$ of the smallest non–increasing halo function that can explain the effect (this includes the condition $\tilde{\eta}^{est}_{0}(v_{min}>v_{s})$=0, with $v_{s}$ the maximal value of the $v_{min}$ range corresponding to the DAMA excess), and to use it to calculate for each experiment among $k$=SIMPLE, COUPP, PICO-60, PICASSO and PICO-2L and for each of the corresponding energy thresholds $E_{th,i}$ the expected number of WIMP events $N_{k,i}^{expected}$ and compare it to the corresponding 90% C.L. upper bound $N_{k,i}^{bound}$ (see Appendix B of Ref[37], Appendix A of Ref. [11] and Appendix C of the present paper for further details, and Fig.6 for an explicit example of the $\tilde{\eta}^{est}$ function determination from the data). Then a straightforward generalization of the compatibility factor of Eq.(5.1) is 777In our numerical analysis we use a refined version of Eq.(5.2), given in Eq. (3.9) of [11].:

[TABLE]

In this Section we will also test the compatibility of the DAMA modulation effect when the WIMP velocity distribution is fixed to a standard Maxwellian at rest in the galactic rest frame. Specifically, we adopt the r.m.s. velocity $v_{rms}$=270 km/s, cut it at the escape velocity $v_{esc}^{Galaxy}$=550 km/sec and boost it to the lab rest frame through $v_{sun}$=232 km/s. In this case for a fixed value of the WIMP mass the whole experimental spectrum allows to get directly a single upper bound, or in case of an excess, a single estimation of the combination $\xi\sigma_{ref}$. As far as the DAMA modulation effect is concerned, we estimate an interval $\sigma_{DAMA}^{min}\leq\sigma_{eff}\leq\sigma_{DAMA}^{max}$ for the cross section by minimizing a $\chi$–square of the $\tilde{\eta}_{1}$ function through the corresponding experimental estimations (in this case $\xi\sigma_{eff}$ is just a normalization factor of $\tilde{\eta}_{1}$), and apply a quality check on the corresponding prediction for the modulation amplitudes, requiring that the $p$–value of the minimal $\chi$–square exceeds $p_{min}$=0.05. In the case of null experiments we will use the same energy bins used for the halo–independent approach, and add one last bin containing the whole experimental range analyzed in the experiment. For each experiment $k$ and each energy bin $i$ we will then calculate the corresponding upper bound $\sigma_{ki}^{bound}$. Then for the Maxwellian case we adopt the compatibility factor:

[TABLE]

Also in this case compatibility between DAMA and constraints is ensured if ${\cal D}_{Maxwellian}<1$.

The result of a systematic scan of the compatibility factors ${\cal D}_{\mbox{Halo indep.}}$ and ${\cal D}_{Maxwellian}$ is shown in Fig.3(a), where regions in the $m_{\chi}$–$\delta$ plane enclosed by lines represent configurations with ${\cal D}_{\mbox{Halo indep.}}<1$, and shaded areas those with ${\cal D}_{Maxwellian}<1$ (Fig. 3(b) is a zoom–up of Maxwellian regions). Notice that, as far as the compatibility factor is concerned, only the nuclear scaling law is relevant, and so the relativistic models of Table 1 which correspond to the same non–relativistic effective operator in Eq. (2.4) lead to identical results. Explicitly: ${\cal O}_{4}^{AA}\rightarrow{\cal O}_{4}$, ${\cal O}_{9}^{VA}$, ${\cal O}_{9}^{TA}$, ${\cal O}_{9}^{AT}\rightarrow{\cal O}_{9}$, ${\cal O}_{10}^{SP}\rightarrow{\cal O}_{10}$, ${\cal O}_{6}^{PP}$, ${\cal O}_{6}^{T^{\prime}T^{\prime}}\rightarrow{\cal O}_{6}$. For this reason the relevant discussion is in terms of non–relativistic operators. Indeed, for all of them Fig.3 shows that regions in the parameter space exist where the DAMA effect can be explained in terms of the particle $\chi_{1}$ spin–dependent upscatters off protons in sodium in compliance with all the other constraints. In the same figure the open circle represents a benchmark for which measurements and bounds for the functions $\tilde{\eta}_{0}$, $\tilde{\eta}_{1}$ are shown in detail in Fig.6 when the DAMA effect is interpreted in terms of the non–relativistic operator ${\cal O}_{4}$.

Actually, as far as the standard spin–dependent coupling ${\cal O}_{4}$ is concerned, the region of Fig. 3 represents an update of the analysis of Ref.[11], since we include the latest experimental constraints from PICO-2L[31, 32], LUX[13, 14, 15] and PANDA[16, 17]. Notice that new constraints do not affect the viability of our scenario, thanks to its intrinsic xenon–phobic and fluoro-phobic nature. On the other hand, the regions corresponding to the non–standard interactions ${\cal O}_{9}$, ${\cal O}_{10}$ and ${\cal O}_{6}$ represent an extension of the analysis of Ref. [11]. A first feature that is worth pointing out is that the models ${\cal O}_{9}$ and ${\cal O}_{10}$ have similar behaviors. This is due do the fact that, as can be seen from Table 1, they both scale with the second power of the exchanged momentum $q$, and they depend on the two nuclear response functions $W^{\tau\tau^{\prime}}_{\Sigma^{\prime}}(q^{2})$ and $W^{\tau\tau^{\prime}}_{\Sigma^{\prime\prime}}(q^{2})$, which are related to the two projections of the nuclear spin along the direction perpendicular or parallel to the exchanged momentum, respectively, with the property $W^{\tau\tau^{\prime}}_{\Sigma^{\prime}}(q^{2})\simeq 2W^{\tau\tau^{\prime}}_{\Sigma^{\prime\prime}}(q^{2})$ when $q\rightarrow$ 0, implying a similar scaling law of the cross section off different nuclei. Moreover, the allowed regions shrink when passing from ${\cal O}_{4}$ to ${\cal O}_{9}\simeq{\cal O}_{10}$ and finally to ${\cal O}_{6}$, both for a halo–independent and a Maxwellian analysis. Such behavior is easily understood when one considers that our DM candidate does not interact with xenon, germanium and fluorine, and that for its relevant masses, 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, only iodine in COUPP (with a recoil energy threshold $E_{th}$=7.8 keV which is smaller than that in DAMA and KIMS) is sensitive enough to put constraints on it. We have confirmed this by numerical inspection. Then, scattering off sodium implies smaller exchanged momenta compared to those for scattering off iodine (notice that $q=\sqrt{2M_{N}E_{R}}$), so a larger power of $q$ implies a relative suppression of the expected DAMA signal compared to the same quantity in COUPP, and so a more stringent limit and a smaller allowed parameter space. Indeed for ${\cal O}_{4}$, ${\cal O}_{9}\simeq{\cal O}_{10}$ and ${\cal O}_{6}$ the nuclear response function scales as $q^{n}$, with $n=0,2,4$, respectively.

As explained in Section 4, a discussion of the thermal relic abundance requires the determination from experimental data of the effective scale $\tilde{\Lambda}^{max}_{p,exp}$ defined in Eq. (4.8). As a way of example, contour plots of $\tilde{\Lambda}^{max}_{p,exp}$ are shown in Figs. 4 and 5 for two of the non–relativistic operators, ${\cal O}_{4}$ and ${\cal O}_{9}$, respectively. In both figures the left–handed plot shows the result of the halo–independent procedure described in Eq. (4.11), while the right–handed one corresponds to the same quantity when the WIMP velocity distribution is fixed to a Maxwellian. The situation for all the non–relativistic couplings is summarized in Table 2, where the maximal value reached by $\tilde{\Lambda}^{max}_{p,exp}$ within the corresponding experimentally viable parameter region ($\cal D<$1) is shown for a halo–independent and a Maxwellian analysis. A first observation is that in the halo–independent case the effective scale turns out to be larger than in the Maxwellian one, an obvious consequence of the choice for the modulated halo function $\tilde{\eta}_{1}$ of the expression given in Eq.(4.10): such functional form is numerically much larger than the one in the case of a Maxwellian, and leads to a smaller $\sigma_{ref}$ and a larger $\tilde{\Lambda}^{max}_{p,exp}$. The second feature worth to be pointed out is that, since non–standard interactions with explicit momentum dependence have response functions proportional to powers of factors $(q/m_{\cal N})^{2}<<$1 (with $m_{\cal N}$ the nucleon mass, see Table. (1)), the ensuing suppression in the expected rate off sodium needs to be compensated by smaller $\tilde{\Lambda}^{max}_{p,exp}$ values to explain the observed modulation amplitudes in DAMA. This explains why in Table 2 the effective scale $\tilde{\Lambda}^{max}_{p,exp}$ decreases moving from ${\cal O}_{4}$ to ${\cal O}_{9}\simeq{\cal O}_{10}$ and finally to ${\cal O}_{6}$. Notice that strictly speaking this is a spurious effect since in some of the interactions terms of Table 1 the choice $m_{M}$=$m_{\cal N}$ in the factors $(q/m_{M})$ is arbitrary. However the choice of the scale $m_{M}$ affects direct detection and the relic abundance in the same way, and so does not modify the correlation among the two quantities that is the object of our study.

Moving on to such correlation, the conclusion that can be drawn from Figs.4, 5 and Table 2 is that, when the corresponding ranges of $\tilde{\Lambda}^{max}_{p,exp}$ are compared to the curves for $\tilde{\Lambda}^{min}_{p,\Omega_{0}}$ shown in Fig.1 the inequality $\tilde{\Lambda}^{max}_{p,exp}>\tilde{\Lambda}^{min}_{p,\Omega_{0}}$ (see Eq.(4.15)) can be only verified for the halo–independent analysis of the ${\cal O}_{4}$ model. Indeed, in all other cases $\tilde{\Lambda}^{min}_{p,\Omega_{0}}$ needs to be too large to get $\Omega$=$\Omega_{0}$, an indication that the combination $\sum_{q}N_{q}c_{q}^{2}\tilde{<\sigma>}_{q}$ in the annihilation cross section (see Eq.(4.12)) is also too large. Notice that this implies that our choice to include in the relic abundance calculation only the couplings of the DM particles to the quarks that drive spin–dependent direct detection already includes the most favorable situation and adding more couplings would not change our conclusions. The same argument justifies our choice $|c_{\nu}|\ll|c_{q}|$ at the end of Section 3.

Notice that in Table 2 for the models that do not verify the condition of Eq. (4.15) the validity of EFT at the scale of thermal decoupling, \Lambda\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}2m_{\chi}, needs not to be a requirement, since anyway some alternative mechanism, involving additional degrees of freedom, needs to be advocated for the relic abundance. Such models can explain the DAMA effect provided that the relic abundance is fixed to the observed one by some other mechanism, and that the EFT is valid at the much lower scale of the direct detection process, q\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}100 MeV/c. We find that, when the corresponding values of the effective scale $\tilde{\Lambda}^{max}_{p,exp}$ are converted to the physical cut–off scale $\Lambda$ using Eq. (4.8) such requirement can be met by all the models of Table 2. In particular, for the most favorable choice of parameters in Eq. (4.8) ($\xi=1$ and $|c_{p}|^{1/2}$ maximized in terms of the coupling constants) $\Lambda$ is typically of the order of a few tens of GeV, and always larger that about 800 MeV in the most extreme cases.

As far as model ${\cal O}_{4}$ is concerned, in Fig.4(a) the dashed line represents configuration where $\tilde{\Lambda}^{max}_{p,exp}$=$\tilde{\Lambda}^{min}_{p,\Omega_{0}}$, and the dashed region where $\tilde{\Lambda}^{max}_{p,exp}>\tilde{\Lambda}^{min}_{p,\Omega_{0}}$. Such region represents the only pSIDM parameter space interval, including the generalizations of Table 1, compatible to the assumption that the same effective dimension–six operator that explains the DAMA excess can also yield an acceptable relic abundance through thermal decoupling. To be more specific, in such parameter space $\tilde{\Lambda}^{min}_{p,\Omega_{0}}$ is achieved for $r_{d}\simeq$ 0.47 and $r_{s}\simeq$ -0.3 which, at fixed values of $m_{\chi}$ can be converted in explicit allowed intervals of the couplings $c_{q}$ and the cut–off scale $\Lambda$ using Eq.(4.16) and the relevant curve of Fig. 1. For instance, for $m_{\chi}\simeq$ 15 GeV one has $\xi_{min}\simeq$ 1.5$\times$10*-6*, implying:

[TABLE]

A combination fulfilling the relation above is, for instance, $c_{u}\simeq$ 1.3, $c_{d}\simeq$ 0.66, $c_{s}\simeq$ -0.42 and $\Lambda\simeq$ 270 GeV.

The allowed parameter space of model ${\cal O}_{4}^{AA}$ would shrink if additional couplings besides the minimal ones driving direct detection were included in the relic abundance calculation, and for a halo function $\tilde{\eta}_{1}$ different from Eq.(4.10). On the other hand no compatibility is achieved in the Maxwellian case irrespective of the EFT model.

6 Conclusions

In the present paper we have discussed strategies to make inferences on the thermal relic abundance of a Weakly Interacting Massive Particle (WIMP) when the same effective dimension–six operator that explains an experimental excess in direct detection is assumed to drive decoupling at freeze–out, and applied them to the explicit pSIDM scenario introduced in Ref. [11], a phenomenological set–up containing two DM states $\chi_{1}$ and $\chi_{2}$ with spin–dependent couplings to protons, masses $m_{\chi}$=$m_{\chi_{1}}<m_{\chi_{2}}$ and a small mass splitting $\delta$ that we showed to explain the DAMA effect in compliance with the constraints from other detectors. Compared to Ref. [11] we have updated experimental constraints including the latest bounds from LUX[14, 15], PANDA[16, 17] and PICO-2L[32], and extended the analysis to the most general spin–dependent momentum–dependent interactions allowed by non–relativistic effective theory, listed in Table 1.

To assess the viability of an interpretation of the DAMA modulation effect in terms of pSIDM we have calculated the compatibility factors ${\cal D}_{\mbox{Halo indep.}}$ of Eq.(5.2) and ${\cal D}_{Maxwellian}$ of Eq.(5.3): the resulting allowed regions, defined by the requirements ${\cal D}_{\mbox{Halo indep.}}<1$ and ${\cal D}_{Maxwellian}<1$, are shown in Fig.3. The largest allowed parameter space is obtained for the standard spin–dependent ${\cal O}_{4}$ non–relativistic coupling and for a halo–independent analysis, 10 GeV \mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}m_{\chi}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 27.5 GeV, 12.3 keV \mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}\delta\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 28.2 keV, although for all the models of Table 1 compatible regions exist, also when a Maxwellian velocity distribution is assumed for the WIMPs, albeit in smaller parameter ranges.

The bottom–up problem of calculating the relic abundance (at the GeV scale) fixing the model parameters with direct detection data (at the keV scale) represents the first necessary step upward in the bottom–up quest to shed light on the ultraviolet completion of the EFT. As discussed in Section 5, it is affected by a degeneracy between the WIMP local density $\rho_{\chi}$ and the effective WIMP–nucleon scattering cross section $\sigma_{ref}$ that parameterizes the particle–physics interaction, since $\rho_{\chi}$ must be rescaled with respect to the observed DM density in the neighborhood of the Sun when the calculated relic density $\Omega$ is smaller than the observed one $\Omega_{0}$. As a consequence, a direct detection experiment is not directly sensitive to the physical cut–off scale $\Lambda$ of the effective theory, but on the dimensional combination $\tilde{\Lambda}_{p,exp}$ defined in Eq. (4.8) and predicted by a given effective theory to be equal to the quantity $\tilde{\Lambda}_{p,th}$ defined in Eq.(4.14) and plotted in Fig. 1 for the scenarios of Table 1. Although $\tilde{\Lambda}_{p,th}$ is numerically equal to the interaction strength $\tilde{\Lambda}_{p,\Omega_{0}}$ required to get the observed relic abundance $\Omega_{0}$, the dependence on the physical scale $\Lambda$ cancels out in its calculation, so it is not correlated to the actual $\Omega$. In other words, due to rescaling a direct detection experiment is bound to measure $\tilde{\Lambda}_{p}=\tilde{\Lambda}_{p,\Omega_{0}}$ even if $\Omega\ll\Omega_{0}$. The only direct test on the physical scale is the upper bound on $\Lambda$ corresponding to the condition $\Omega\leq\Omega_{0}$, which should be compared to the lower bound \Lambda\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{>}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}2m_{\chi} required for the Effective Theory to be valid and to $|c_{q}|\leq 4\pi$ from perturbativity; for the specific scenarios discussed in our paper such constraint turns out to be irrelevant (see Fig.2).

The degeneracy discussed above allows to develop the consistency test of Eq. (4.15) for a given experimental excess:

[TABLE]

where $\tilde{\Lambda}_{p,\Omega_{0}}$ is minimized with respect to the couplings ratios $r_{q}=c_{q}/c_{u}$, at fixed WIMP–neutron coupling $c_{n}$=0 (as required by the proton–philic nature of our DM candidate). Such test is not simply a requirement to obtain the correct relic abundance or to avoid overclosure ($\Omega>\Omega_{0}$) when the same parameters that explain a given experimental excess are used in the coannihilation cross section. Instead, it should be interpreted as a more general condition on the possibility that the same effective model that explains such excess can also determine through thermal decoupling the relic abundance $\Omega$ in a consistent way, and irrespective of the actual value of $\Omega$.

In Section 5 we have calculated in a systematic way the scale $\tilde{\Lambda}_{p,exp}^{max}$ from the DAMA data, summarizing the result in Figs. 4 and 5 and in Table 2. In particular, when comparing for each model the maximum value reached by $\tilde{\Lambda}_{p,exp}^{max}$ as provided in Table 2 to the curves for $\tilde{\Lambda}^{min}_{p,\Omega_{0}}$ shown in Fig.1 it is possible to conclude that the inequality $\tilde{\Lambda}^{max}_{p,exp}>\tilde{\Lambda}^{min}_{p,\Omega_{0}}$ can be only verified for the halo–independent analysis of a standard spin–dependent interaction ${\cal O}$=$\bar{\chi}_{1}\gamma^{\mu}\gamma^{5}\chi_{2}\bar{q}\gamma_{\mu}\gamma^{5}q$ +h.c., with no explicit momentum dependence. The ensuing allowed parameter space corresponds to the shades area of Fig.4(a), i.e., approximately, 10 GeV \mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}m_{\chi}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 16 GeV, 17 keV\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}\delta\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}28 keV. From Fig.2 one can see that the actual value of $\Omega$ is basically undetermined in this scenario, 6\times 10^{-7}\Omega_{0}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}\Omega\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}}\Omega_{0}. Such degeneracy on the physical cut–off scale $\Lambda$ could in principle be at least partially resolved by using constraints from accelerator signals. For all the other spin–dependent effective models listed in Table 1 the pSIDM scenario cannot explain the DAMA effect and at the same time provide a thermal relic if the same dimension–six operator drives both direct detection and decoupling at freeze–out, and this conclusion is valid for any choice of the WIMP galactic velocity distribution $f(v)$. Moreover, in no case the same effective coupling can provide a thermal relic and explain DAMA if $f(v)$ is a Maxwellian.

An additional complication of the pSIDM scenario is that for all the models of Table 1, if no other coupling is assumed for the DM particle, the lifetime of the heavy state $\chi_{2}$ is larger than the age of the Universe. In such case our Galaxy would be formed by a mixture of $\chi_{1}$ and $\chi_{2}$. Presumably such states would be present in equal parts and have the same velocity distribution, so that direct detection experiments would be both sensitive to the upscatters of $\chi_{1}$’s and the downscatters of $\chi_{2}$’s. However, as explained in section 2 downscatters of the $\chi_{2}$ particle off fluorine must be suppressed in order to comply with the constraints from bubble chambers and droplet detectors. In our analysis we have assumed the existence of some additional mechanism (such as the coupling of Eq.(3.5) with $|c_{\nu}|\ll|c_{q}|$) that depletes the $\chi_{2}$ particles without affecting the relic density. As a consequence, our model does not produce any indirect detection signal.

Acknowledgments

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. S.S. acknowledges the hospitality of the Instituto de Fisíca Teórica (IFT), Madrid, where part of this work was carried out, and partial support by the Programme SEV-2012-0249 ‘Centro de Excelencia Severo Ochoa’.

Appendix A Thermal coannihilation cross sections

For a Cold Dark Matter particle the thermal average of the coannihilation cross section times velocity at temperature $T$ can be expressed as $<\sigma v>=a+b/x$ with $x\equiv T/m_{\chi}\ll 1$. The relic abundance is inversely proportional to the combination $\tilde{<\sigma v>}=a+b/(2x_{f})$ with $x_{f}\equiv T_{f}/m_{\chi}\simeq 1/20$ and $T_{f}$ the freeze–out temperature. Since $T>>\delta$ in the calculation the mass difference between $\chi_{1}$ and $\chi_{2}$ can be safely neglected. The coefficients $a$ and $b$ for the different models discussed in the present paper are given by $a=1/\Lambda^{4}\sum_{f}N_{f}c_{f}^{2}(a)_{f}$, $b=1/\Lambda^{4}\sum_{f}N_{f}c_{f}^{2}(b)_{f}$, where the sum is over the final fermionic states $f$ with coupling $c_{f}^{2}/\Lambda^{4}\bar{\chi}_{1}\Gamma^{a}\chi_{2}\bar{f}\Gamma^{b}f$ and internal degrees of freedom $N_{f}$, and, using the same notation of Table 1 to indicate the different effective models considered in the present paper, $(a)_{f}$ and $(b)_{f}$ are given by:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The observed DM density corresponds to the reference value $\tilde{<\sigma v>}_{0}=2\times 10^{-9}$ GeV*-2* (in natural units), and the rescaling factor discussed in section 4 is given by $\xi=\tilde{<\sigma v>}_{0}/\tilde{<\sigma v>}$.

Appendix B Nuclear matrix elements

As, for instance, summarized in Section A of [41] it is usually assumed that between the WIMP–nucleon and the WIMP–quark bilinear currents the following relation holds:

[TABLE]

with $\Delta q^{{\cal N}}$ quantities incorporating the complicated non–perturbative physics through which quarks are confined into the nucleons ${\cal N}=n,p$. Usually, such quantities, which most of the times cannot be calculated from first principles but must be estimated from some measured processes depending on the same current, are taken as constants, i.e. independent on the spins and velocities of the nucleon and the quark.

Once the $\Delta q^{{\cal N}}$’s are known, it is straightforward to calculate the effective WIMP–nucleon couplings in terms of the WIMP–quark couplings:

[TABLE]

where the constant $K$ can be directly read–off comparing the second and fourth columns in Table 1.

In particular the axial current $\Gamma^{A}$ is directly related to the nucleon spin (actually, the corresponding $\Delta q^{{\cal N}}_{A}$’s can be directly interpreted as the fractions of the nucleon spin carried by the quark $q$). On the other hand the pseudoscalar current $\Gamma_{5}$ can be related to $\Gamma^{A}$ due to the fact that the latter is not conserved in the nucleon for non–vanishing quark masses [51, 52, 53] so that the corresponding $\Delta q^{{\cal N}}_{5}$ are related to the $\Delta q^{{\cal N}}_{A}$’s[41]. The numerical values of both the $\Delta q^{{\cal N}}_{A}$ and $\Delta q^{{\cal N}}_{5}$’s are for instance given in [54] and, for $q$=$u$, $d$, $s$ are shown in Tables 3 and 4.

As far as the tensor current $\Gamma^{T}$ is concerned, it is related to the transversity distribution of the nucleon [60], which has to do with the contribution of quarks to the projection of the spin perpendicular to the nucleon’s momentum. Only recently the $\Delta q^{{\cal N}}_{T}$’s have been measured by the two experiments Hermes [55] and COMPASS [56], while for a long time their only available estimates had been obtained from from Dyson–Schwinger formalism [57, 58] and lattice QCD [59]. In all cases they are affected by large uncertainties, as shown by the numerical values of the $\Delta q^{p}_{T}$’s for the proton shown in Table 5 (the corresponding ones for the neutron are obtainable by isospin symmetry). Moreover, thanks to the identity $\Gamma^{T^{\prime}}=\sigma^{\mu\nu}\gamma_{5}$=$(i/2)\epsilon^{\mu\nu\alpha\beta}\sigma_{\alpha\beta}$ one has $\Delta q^{{\cal N}}_{T^{\prime}}$=$\Delta q^{{\cal N}}_{T}$. In our quantitative analysis we used the Hermes determination of the $\Delta q_{T}^{\cal N}$ coefficients.

Appendix C Experimental constraints

In the present analysis we include the constraints from LUX[13, 14, 15], PANDA[16, 17], XENON10[19], CDMS-$Ge$[23], CDMSlite [24], SuperCDMS[25], SIMPLE[27], COUPP[28], PICASSO[29, 30], PICO-2L[31, 32] and PICO-60[33]). Compared to Refs.[37] and [11] (to which we refer for the details of the other experiments) in the present paper we include the latest results from LUX[14, 15], PANDA[16, 17] and PICO-2L[32].

As far as the LUX and PANDA analyses are concerned, we use the generalization of the response function ${\cal R}$ as given in Appendix A of Ref.[35].

For LUX we use the events in the lower half of the signal region of Fig. 1 of Ref. [15] with a total exposition of 3.35$\times$105 kg day and the efficiency from Fig.2 of Ref. [15] multiplied by 0.5, and convert recoil energies to photo–electron (PE) numbers using Eq. (3) of [61] with $g_{1}$=0.115 and the $L_{y}$ function taken from Fig. 10 of the same paper.

For PANDA we combine the events in the lower half of the signal regions of Fig. 4 of [16] and Fig. 5 of [17], assuming a total exposition of 3.3$\times$105 kg day and the efficiency from Fig.2 of Ref. [17] multiplied by 0.5. The PANDA collaboration does not provide details on how to convert the recoil energy $E_{R}$ (in keVnr) into the signal $S_{1}$ (in PE). We decide to adopt the same recipe of LUX after checking that it quantitatively reproduces fairly well the correspondence between $S_{1}$ and $E_{R}$ ranges as deduced by the plots in [17].

The latest PICO-2L analysis [32] was performed after an upgrade of the detector that significantly reduced the background compared to [31]. Only the threshold $E_{th}$=3.3 keV was analyzed, with a total exposure of 129.0 kg day and 1 event detected (corresponding to a 95% upper bound of 5.14 events). We use the nucleation probability give in Eq.(B.3) of [37] with $\alpha_{F}$=$\alpha_{C}$=5, and the efficiency of Fig.4 of [31].