COSINE-100 and DAMA/LIBRA-phase2 in WIMP effective models

TL;DR

This study analyzes WIMP detection constraints from COSINE-100 and DAMA/LIBRA-phase2 using non-relativistic effective theory, revealing that certain operators allow compatibility between the two experiments' results.

Contribution

It introduces a systematic comparison of WIMP detection results considering effective operators, highlighting the impact of spectral shape and modulation fractions on experimental compatibility.

Findings

Compatibility between DAMA/LIBRA and COSINE-100 depends on particle-physics models.

Certain effective operators predict higher modulation fractions, allowing agreement.

COSINE-100's sensitivity will improve with future analyses.

Abstract

Assuming a standard Maxwellian for the WIMP velocity distribution, we obtain the bounds from null WIMP search results of 59.5 days of COSINE-100 data on the DAMA/LIBRA-phase2 modulation effect within the context of the non-relativistic effective theory of WIMP-nucleus scattering. Here, we systematically assume that one of the effective operators allowed by Galilean invariance dominates in the effective Hamiltonian of a spin-1/2 dark matter (DM) particle. We find that, although DAMA/LIBRA and COSINE-100 use the same sodium-iodide target, the comparison of the two results still depends on the particle-physics model. This is mainly due to two reasons: i) the WIMP signal spectral shape; ii) the expected modulation fractions, when the upper bound on the time-averaged rate in COSINE-100 is converted into a constraint on the annual modulation component in DAMA/LIBRA. We find that the latter effect is the dominant one. For several effective operators the expected modulation fractions are larger than in the standard spin-independent or spin-dependent interaction cases. As a consequence, compatibility between the modulation effect observed in DAMA/LIBRA and the null result from COSINE-100 is still possible for several non-relativistic operators. At low WIMP masses such relatively high values of the modulation fractions arise because COSINE-100 is mainly sensitive to WIMP-sodium scattering events, due to the higher threshold compared to DAMA/LIBRA. A next COSINE analysis is expected to have a full sensitivity for the 5$\sigma$ region of DAMA/LIBRA.

Tables

Table 1. Table 1 : Non-relativistic Galilean invariant operators for dark matter with spin 1 / 2 1 2 1/2 .

${𝒪}_1=1_{\chi }{1}_N$
${𝒪}_2={\left(v^{⟂}\right)}^2$	${𝒪}_9=i{\overrightarrow{S}}_{\chi }\cdot \left({\overrightarrow{S}}_N\times \frac{\overrightarrow{q}}{m_N}\right)$
${𝒪}_3=i{\overrightarrow{S}}_N\cdot \left(\frac{\overrightarrow{q}}{m_N}\times {\overrightarrow{v}}^{⟂}\right)$	${𝒪}_{10}=i{\overrightarrow{S}}_N\cdot \frac{\overrightarrow{q}}{m_N}$
${𝒪}_4={\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{S}}_N$	${𝒪}_{11}=i{\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}$
${𝒪}_5=i{\overrightarrow{S}}_{\chi }\cdot \left(\frac{\overrightarrow{q}}{m_N}\times {\overrightarrow{v}}^{⟂}\right)$	${𝒪}_{12}={\overrightarrow{S}}_{\chi }\cdot \left({\overrightarrow{S}}_N\times {\overrightarrow{v}}^{⟂}\right)$
${𝒪}_6=\left({\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}\right)\left({\overrightarrow{S}}_N\cdot \frac{\overrightarrow{q}}{m_N}\right)$	${𝒪}_{13}=i\left({\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{v}}^{⟂}\right)\left({\overrightarrow{S}}_N\cdot \frac{\overrightarrow{q}}{m_N}\right)$
${𝒪}_7={\overrightarrow{S}}_N\cdot {\overrightarrow{v}}^{⟂}$	${𝒪}_{14}=i\left({\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}\right)\left({\overrightarrow{S}}_N\cdot {\overrightarrow{v}}^{⟂}\right)$
${𝒪}_8={\overrightarrow{S}}_{\chi }\cdot {\overrightarrow{v}}^{⟂}$	${𝒪}_{15}=-\left({\overrightarrow{S}}_{\chi }\cdot \frac{\overrightarrow{q}}{m_N}\right)\left(\left({\overrightarrow{S}}_N\times {\overrightarrow{v}}^{⟂}\right)\cdot \frac{\overrightarrow{q}}{m_N}\right)$

Table 2. Table 2 : Nuclear response functions corresponding to each coupling c j subscript 𝑐 𝑗 c_{j} of the effective Hamiltonian ( 2.6 ), for the velocity–independent and the velocity–dependent components parts of the WIMP response function, decomposed as in Eq.( 2.8 ). In parenthesis the power of q 𝑞 q in the WIMP response function.

${𝐜}_{𝐣}$	$R_{0k}^{\tau {\tau }^{\prime }}$	$R_{1k}^{\tau {\tau }^{\prime }}$	${𝐜}_{𝐣}$	$R_{0k}^{\tau {\tau }^{\prime }}$	$R_{1k}^{\tau {\tau }^{\prime }}$
$c_1$	$M\left(q^0\right)$	-	$c_3$	${\Phi }^{\prime \prime }\left(q^4\right)$	${\Sigma }^{\prime }\left(q^2\right)$
$c_4$	${\Sigma }^{\prime \prime }\left(q^0\right)$,${\Sigma }^{\prime }\left(q^0\right)$	-	$c_5$	$\Delta \left(q^4\right)$	$M\left(q^2\right)$
$c_6$	${\Sigma }^{\prime \prime }\left(q^4\right)$	-	$c_7$	-	${\Sigma }^{\prime }\left(q^0\right)$
$c_8$	$\Delta \left(q^2\right)$	$M\left(q^0\right)$	$c_9$	${\Sigma }^{\prime }\left(q^2\right)$	-
$c_{10}$	${\Sigma }^{\prime \prime }\left(q^2\right)$	-	$c_{11}$	$M\left(q^2\right)$	-
$c_{12}$	${\Phi }^{\prime \prime }\left(q^2\right)$,${\tilde{\Phi }}^{\prime }\left(q^2\right)$	${\Sigma }^{\prime \prime }\left(q^0\right)$,${\Sigma }^{\prime }\left(q^0\right)$	$c_{13}$	${\tilde{\Phi }}^{\prime }\left(q^4\right)$	${\Sigma }^{\prime \prime }\left(q^2\right)$
$c_{14}$	-	${\Sigma }^{\prime }\left(q^2\right)$	$c_{15}$	${\Phi }^{\prime \prime }\left(q^6\right)$	${\Sigma }^{\prime }\left(q^4\right)$

Table 3. Table 3 : Absolute and local minima of the DAMA–phase2 modulation result χ 2 superscript 𝜒 2 \chi^{2} of Eq.( 3.12 ) for each of the couplings c j subscript 𝑐 𝑗 c_{j} of the effective Hamiltonian ( 2.6 ). From Ref.  [ 38 ] .

${𝐜}_{𝐣}$	${𝐦}_{\chi ,\mathrm{𝐦𝐢𝐧}}$ (GeV)	${𝐫}_{\chi ,\mathrm{𝐦𝐢𝐧}}$	$\sigma \left({\text{cm}}^{\mathrm{𝟐}}\right)$	${\chi }_{\mathrm{𝐦𝐢𝐧}}^{\mathrm{𝟐}}$
$c_1$	11.17	-0.76	2.67e-38	11.38
	45.19	-0.66	1.60e-39	13.22
$c_3$	8.10	-3.14	2.27e-31	11.1
	35.68	-1.10	9.27e-35	14.23
$c_4$	11.22	1.71	2.95e-36	11.38
	44.71	-8.34	5.96e-36	27.7
$c_5$	8.34	-0.61	1.62e-29	10.83
	96.13	-5.74	3.63e-34	11.11
$c_6$	8.09	-7.20	5.05e-28	11.11
	32.9	-6.48	5.18e-31	12.74
$c_7$	13.41	-4.32	4.75e-30	13.94
	49.24	-0.65	1.35e-30	38.09
$c_8$	9.27	-0.84	8.67e-33	10.82
	42.33	-0.96	1.30e-34	11.6
$c_9$	9.3	4.36	8.29e-33	10.69
	37.51	-0.94	1.07e-33	15.23
$c_{10}$	9.29	3.25	4.74e-33	10.69
	36.81	0.09	2.25e-34	12.40
$c_{11}$	9.27	-0.67	1.15e-34	10.69
	38.51	-0.66	9.17e-37	13.02
$c_{12}$	9.26	-2.85	3.92e-34	10.69
	35.22	-1.93	2.40e-35	12.47
$c_{13}$	8.65	-0.26	1.21e-26	10.76
	29.42	0.10	5.88e-29	14.28
$c_{14}$	10.28	-0.59	2.61e-26	11.21
	38.88	-1.93	2.19e-27	14.48
$c_{15}$	7.32	-3.58	2.04e-27	12.91
	33.28	4.25	2.05e-33	16.26

Table 4. Table 4 : Minimum value of the modulation fraction ( S m D ​ A ​ M ​ A / S 0 D ​ A ​ M ​ A ) E ′ < 3.5 ​ keVee subscript subscript superscript 𝑆 𝐷 𝐴 𝑀 𝐴 𝑚 subscript superscript 𝑆 𝐷 𝐴 𝑀 𝐴 0 superscript 𝐸 ′ 3.5 keVee (S^{DAMA}_{m}/S^{DAMA}_{0})_{E^{\prime}<3.5~{}\rm keVee} in the three DAMA energy bins for 2 keVee ≤ E ′ ≤ absent superscript 𝐸 ′ absent \leq E^{\prime}\leq 3.5 keVee, where the bulk of the DAMA modulation effect above the COSINE-100 threshold is concentrated

${𝐜}_{𝐣}$			${𝐜}_{𝐣}$
	Low-mass	High-mass		Low-mass	High-mass
	local minimum	local minimum		local minimum	local minimum
$c_1$	0.066	0.054	$c_3$	0.120	0.098
$c_4$	0.065	0.047	$c_5$	0.122	0.059
$c_6$	0.121	0.111	$c_7$	0.097	0.080
$c_8$	0.094	0.072	$c_9$	0.093	0.079
$c_{10}$	0.093	0.085	$c_{11}$	0.094	0.083
$c_{12}$	0.094	0.096	$c_{13}$	0.123	0.139
$c_{14}$	0.126	0.122	$c_{15}$	0.146	0.113

Table 5. Table 5 : Combination of the DAMA/LIBRA–phase1 and the DAMA/LIBRA–phase2 measurements for the modulation amplitudes S m , k subscript 𝑆 𝑚 𝑘 S_{m,k} with statistical errors σ k subscript 𝜎 𝑘 \sigma_{k} used in the present analysis (from Ref.  [ 4 ] ).

Energy (keVee)	$S_{m,k}$	${\sigma }_k$
$1.0-1.5$	0.0242	0.0056
$1.5-2.0$	0.0211	0.0043
$2.0-2.5$	0.0179	0.0023
$2.5-3.0$	0.0197	0.0030
$3.0-3.5$	0.0186	0.0027
$3.5-4.0$	0.0110	0.0026
$4.0-4.5$	0.0109	0.0021
$4.5-5.0$	0.0032	0.0019
$5.0-5.5$	0.0065	0.0019
$5.5-6.0$	0.0059	0.0019
$6.0-6.5$	0.0010	0.0016
$6.5-7.0$	0.0008	0.0017
$7.0-7.5$	0.0009	0.0016
$7.5-8.0$	0.0009	0.0016
$8.0-16.0$	0.0003	0.0004

Equations

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

$$\frac{dR}{d{E}^{\prime }}(t)$$

$$E_{ee}$$

$$\frac{d{R}_{\chi T}}{d{E}_R}(t)=\sum _T{N}_T\frac{{\rho }_{\text{mbox}WIMP}}{m_{\text{mbox}WIMP}}{\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],$$

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

$$\frac{1}{2j_{\chi }+1}\frac{1}{2j_T+1}|\mathcal{M}{|}^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 },c_j^{{\tau }^{\prime }},(v_T^{\perp }{)}^2,\frac{q^2}{m_N^2}\right]W_{Tk}^{\tau {\tau }^{\prime }}(y).$$

$$R_k^{\tau {\tau }^{\prime }}=R_{0k}^{\tau {\tau }^{\prime }}+R_{1k}^{\tau {\tau }^{\prime }}(v_T^{\perp }{)}^2=R_{0k}^{\tau {\tau }^{\prime }}+R_{1k}^{\tau {\tau }^{\prime }}\left(v_T^2-v_{min}^2\right),$$

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

$$N$$

$$|{\vec{v}}_E(t)|=v_{Sun}+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.$$

$${\sigma }_p=(c_j^p{)}^2\frac{{\mu }_{\chi \mathcal{N}}^2}{\pi },$$

$$Asymmetry=(Q_1-Q_2)/(Q_1+Q_2)$$

$$X1=\sum _{100ns}^{600ns}{q}_i/\sum _{0ns}^{600ns}{q}_i$$

$$X2=\sum _{0ns}^{50ns}{q}_i/\sum _{0ns}^{600ns}{q}_i$$

$$X3=\sum _{0ns}^{120ns}{q}_i/\sum _{0ns}^{600ns}{q}_i$$

$$X4=\sum _{100ns}^{150ns}{q}_i/\sum _{0ns}^{600ns}{q}_i$$

$$MT=\sum _{0ns}^{500ns}{q}_i{t}_i/\sum _{0ns}^{500ns}{q}_i$$

$$MTL=\sum _{0ns}^{30ns}{q}_i{t}_i/\sum _{0ns}^{30ns}{q}_i$$

$$MV=\sum _{0ns}^{1000ns}{q}_i{t}_i^2/\sum _{0ns}^{1000ns}{q}_i-(\sum _{0ns}^{1000ns}{q}_i{t}_i/\sum _{0ns}^{1000ns}{q}_i{)}^2$$

$$NC=\text{the number of pulses}$$

$$QC=\sum _{0ns}^{5000ns}{q}_i$$

$$CAT=\text{time of 95\% charge accumulation}(\sum _{0\%}^{95\%}{q}_i),$$

$${\chi }^2(m_{\chi },{\sigma }_p,r)=\sum _{k=1}^{15}\frac{{\left[S_{m,k}(m_{\chi },{\sigma }_p,r)-S_{m,k}^{exp}\right]}^2}{{\sigma }_k^2}$$

Full text

1 Introduction

Since its first presentation, the DAMA/NaI, DAMA/LIBRA–phase1 and DAMA/LIBRA–phase2’s (DAMA for short) observation [1, 2, 3, 4] of an annual modulation signal [5, 6] detected in an array of low–background sodium–iodide crystals has produced continuous speculation about whether or not it is caused by dark matter. Given the substantial statistical significance ($>10\,\sigma$) of the signal, an independent verification of the result with the same sodium–iodide, or NaI(Tl) target crystals is required. Experimental efforts by several groups using the same target medium are on–going and in a few years will confirm the presence of the signal in a model–independent way [7, 8, 9, 10]. On the other hand, a Weakly Interacting Massive Particle (WIMP) [11, 12] dark-matter interpretations for the positive signal have been ruled out by other experiments [13, 14, 15, 16] using time-averaged rates and shapes in the specific context of the Standard Galactic Halo Model [17, 18], assuming a standard spin–independent (SI) or spin–dependent (SD) interaction.

Although preliminary modulation analyses from ANAIS [19] and COSINE-100 [20] have been released, these experiments need several more years of data-taking [21] before they can reach the modulation sensitivity required to probe the DAMA signal. In the meantime, the yearly modulation amplitude and the time–averaged rate can be compared in the context of specific models of WIMP–nucleus interactions. The first such analysis using 59.5 days of COSINE-100 data has excluded the SI WIMP case as a potential model for the DAMA signal interpretation [7].

Since COSINE-100 uses the same target material, fewer assumptions about how the WIMP elastic scattering rates scale with the target nucleus are needed in order to interpret the result. The dominance of WIMP–sodium or WIMP–iodine scattering events for different WIMP masses can then be determined entirely by kinematics. As a consequence, expected rates can be directly expressed in terms of WIMP-nucleus cross sections. For example, this enabled the Korea Invisible Mass Search (KIMS) experiment, using a CsI(Tl) target, to constrain in a model–independent way WIMP–iodine scattering events in DAMA for WIMP masses above approximately 20 GeV/c2 [13]. On the other hand, for lower WIMP masses, where only WIMP–sodium events are kinematically accessible, the DAMA result could only be probed using different target nuclei. Moreover, the recent DAMA/LIBRA–phase2 result with a lower threshold at 1 keV electron-equivalent (keVee), is more sensitive to WIMP–iodine scattering events for WIMP masses below 20 GeV/c2 and the expected relative numbers of WIMP–iodine and WIMP–sodium scattering events is model dependent. To address these aspects, we use non–relativistic (NR) effective field theory in an interpretation of the DAMA signal with the COSINE–100 data.

COSINE-100 [22, 23, 24] is a joint dark matter search experiment of KIMS [25, 26, 27] and DM-Ice [28, 29] with an array of low radioactive NaI(Tl) crystals at the Yangyang underground laboratory. The experiment is composed of eight encapsulated NaI(Tl) crystals (a total of 106 kg) placed in the middle of a copper box which is filled with 2 tons of liquid scintillator [30] and further shielded by lead and plastic scintillator panels from external radiations. The plastic scintillator panels veto cosmic-ray muons and liquid scintillator actively reduces background radiations that originate from crystals or vicinity of the crystal detectors [23]. The details of the COSINE-100 experimental setup and analysis methods can be found elsewhere [27, 7]. Here, we extend the interpretation of the data collected in Ref. [7] to the general NR effective theory of nuclear scattering for a WIMP of spin 1/2.

The paper is organized by the following structure: Section 2 summarizes the theory of WIMP–nucleus scattering in the NR effective theory of a particle of spin 1/2; Section 3 is devoted to our quantitative analysis; Section 4 contains our conclusions.

2 WIMP–nucleus scattering rates in non–relativistic effective models

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_{exp}$ the live–time exposure 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$ per unit mass in the detector (the sum over $T$ applies in the case of more than one target), while

[TABLE]

where $m_{T}$ is the nuclear target mass and the squared amplitude in parenthesis is given explicitly in Eq.(2.7).

WIMP–nucleus scattering is a non–relativistic process that can be fully described in a non–relativistic Effective Theory approach. In the case of a spin–1/2 DM particle the corresponding Hamiltonian density is given by [31, 32]:

[TABLE]

where, $\mathcal{O}_{j}$’s are non-relativistic Galilean invariant operators which have been collected in Table 1. In the same Table $1_{\chi}$ and $1_{N}$ are identity operators, $\vec{q}$ is the transferred momentum, $\vec{S}_{\chi}$ and $\vec{S}_{N}$ are the WIMP and nucleon spins, respectively, while $\vec{v}^{\perp}=\vec{v}+\frac{\vec{q}}{2\mu_{T}}$ (with $\mu_{T}$ the WIMP–nucleus reduced mass) is the relative transverse velocity operator satisfying $\vec{v}^{\perp}\cdot\vec{q}=0$. In particular, one has $(v^{\perp}_{T})^{2}=v^{2}_{T}-v_{min}^{2}$, where, for WIMP–nucleus elastic scattering, $v_{min}^{2}=\frac{q^{2}}{4\mu_{T}^{2}}=\frac{m_{T}E_{R}}{2\mu_{T}^{2}}$ represents the minimal incoming WIMP speed required to impart the nuclear recoil energy $E_{R}$, while $v_{T}\equiv|\vec{v}_{T}|$ is the WIMP speed in the reference frame of the nuclear center of mass. Moreover $t^{0}=1$, $t^{1}=\tau_{3}$ denote the $2\times 2$ identity and third Pauli matrix in isospin space, respectively, and the isoscalar and isovector (dimension -2) coupling constants $c^{0}_{j}$ and $c^{1}_{j}$, are related to those for 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})$.

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 [31]. So, following Refs.[31, 32], we will not include it in our analysis.

Assuming that the nuclear interaction is the sum of the interactions of the WIMPs with the individual nucleons in the nucleus the WIMP scattering amplitude on the target nucleus $T$ can be written in the compact form:

[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 can be found in Ref. [32]) which depend on the couplings $c^{\tau}_{j}$ as well as the transferred momentum $\vec{q}$ and $(v^{\perp}_{T})^{2}$. In equation (2.7) 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 [31, 32], $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.[32, 33]222Setting $k=M$ and $W^{p,n}_{TM}(q)\equiv(W_{TM}^{00}(q)\pm W_{TM}^{01}(q)\pm W_{TM}^{10}(q)+W_{TM}^{11}(q))/4$, in the case of a standard spin–independent interaction one has $16/(2j_{T}+1)W^{p}_{TM}(q)$=$Z^{2}_{T}F^{2}(q)$ and $16/(2j_{T}+1)W^{n}_{TM}(q)$=$(A_{T}-Z_{T})^{2}F^{2}(q)$, with $Z_{T}$ and $A_{T}-Z_{T}$ the number of protons and neutrons in target $T$, and $F(q)$ the SI nuclear form factor, for which the parameterization in [34] is commonly assumed.. The correspondence between models and nuclear response functions can be directly read off from the WIMP response functions $R^{\tau\tau^{\prime}}_{k}$ [32]. In particular, using the decomposition:

[TABLE]

such correspondence is summarized for convenience in Table 2.

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. The yearly modulation effect is due to the time dependence of the Earth’s speed with respect to the Galactic frame:

[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, $t_{0}$=2 June, $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) while $v_{Sun}$=$v_{0}$+12 km/sec, accounting for a peculiar component of the solar system with respect to the galactic rotation. For the two parameters $v_{0}$ and $u_{esc}$ we take $v_{0}$=220 km/sec [35] and $u_{esc}$=550 km/sec [36]. In the isothermal model the time dependence of Eq. (2.11) induces an expected rate with the functional form $S(t)=S_{0}+S_{m}\cos(2\pi/T-t_{0})$, with $S_{m}>0$ at large values of $v_{min}$ and turning negative when v_{min}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 200 km/s. In such regime of $v_{min}$ and below the phase is modified by the focusing effect of the Sun’s gravitational potential [37], while when $S_{m}\ll S_{0}$ the time dependence differs from a simple cosine due the contribution of higher harmonics [6].

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]

while other experiments put upper bounds on the time average $S_{0}$:

[TABLE]

In the present paper, we will systematically consider the possibility that one of the couplings $c_{j}$ dominates in the effective Hamiltonian of Eq. (2.6). In this case it is possible to factorize a term $|c_{j}^{p}|^{2}$ from the squared amplitude of Eq.(2.7) and express it in terms of the effective WIMP–proton cross section:

[TABLE]

(with $\mu_{\chi{\cal N}}$ the WIMP–nucleon reduced mass) and the ratio $r\equiv c_{j}^{n}/c_{j}^{p}$. It is worth pointing out here that among the generalized nuclear response functions arising from the effective Hamiltonian (2.6) only the ones corresponding to $M$ (SI interaction), $\Sigma^{\prime\prime}$ and $\Sigma^{\prime}$ (both related to the standard spin–dependent interaction) do not vanish for $q\rightarrow$0, and so allow to interpret $\sigma_{p}$ in terms of a long–distance, point–like cross section. In the case of the other interactions $\Phi^{\prime\prime}$, $\tilde{\Phi}^{\prime}$ and $\Delta$ the quantity $\sigma_{p}$ is just a convenient alternative to directly parameterizing the interaction in terms of the $c_{j}^{p}$ coupling. Since we will not consider interferences among different couplings the response functions $W^{\tau\tau^{\prime}}_{Tk}(y)$ for $k$=$\Phi^{\prime\prime}M$, $\Delta\Sigma^{\prime}$ will not play any role in our analysis.

3 Analysis

This analysis uses the COSINE-100 data from October 20, 2016 to December 19, 2016. After application of data quality criteria, the 59.5 live days of good data are used for the results presented here. A total of 11 hours of data did not pass these quality requirement where abrupt high PMT noise triggers and electronic interference triggers were rejected. Total exposure is 6303.9 kg$\cdot$day. During this period, light yield, gain, and other environmental data show stable behavior. The overall crystal PMT gain is changed by less than 1% relative to the beginning of the physics run.

An event is triggered if a photon is observed in each PMT within 200 ns in a crystal. When this happens the data acquisition system reads out the full veto detectors including liquid scintillator and plastic scintillators and other crystal signals simultaneously [24]. The detector stability is checked by using crystal internal gamma calibrations which show consistent results with external source calibrations. The low energy electron signals produced from 60Co calibration of the Compton scattering and tagged by neighboring crystals are used to separate the PMT noise events from data.

First, muon-induced events are rejected by requiring the time difference between muon veto events in the plastic panels and the crystal to be less than 30 ms. This efficiently removes most of muon-induced events that directly pass through the crystals. We, then, require that leading edges of the trigger pulses start later than 2.0 $\mu$s, each waveform contain more than two pulses, and integrated charge below the baseline should be small enough. These reject muon-induced phosphor events and electronic interference events. Next, we demand a single-site condition where neighboring crystals should not have more than four photons and an energy deposit by the surrounding liquid scintillator should be less than 20 keV.

To identify scintillation signals, one must reject two types of backgrounds which are more than the desired signals especially below 20 keVee region. The first class is thin pulses that are originated from PMTs. These noise events are triggered partly from the PMT individually and partly from radioactivities inside the PMT circuitry which make the crystal scintillate. The second class, less often than the first, consists of bell-shaped waveforms that occur sporadically and in a few PMTs only. These bell-shaped pulses are produced due to occasional PMT discharge and the shape of a waveform looks more symmetric than a typical scintillation signal.

The initial rejection algorithms focus on eliminating the thin pulses and other pathological events. We calculate the balance of the deposited charge from two PMTs (Asymmetry : Eq. 3.1 shown in a) of Fig. 1), the charge fraction of 500 ns to 600 ns from the first 600 ns (X1 : Eq. 3.2 shown in b) of Fig. 1), the charge fraction of the first 50 ns to first 600 ns (X2 : Eq. 3.3 shown in c) of Fig. 1), the charge–weighted mean time of pulses within first 500 ns (MT : Eq. 3.6 shown in e) of Fig. 1), the total charge (QC : Eq. 3.10) and the number of pulses (NC : Eq. 3.9). Boosted Decision Trees (BDTs) were trained using aforementioned variables. The electron/gammas signal model is obtained from the energy–weighted 60Co multiple-site distributions and data is used for the noise model. Each crystal is trained for a separate BDT. Six parameters comparing noise-containing data with 60Co multiple signals are shown in Fig. 1.

The definitions of each variable used in the rejection algorithms are following,

[TABLE]

where Q1,2 indicates PMT integrated charges within 5 $\mu$s and qi and ti are waveform amplitudes and times for each 2 ns bin, respectively.

For the rejection of the bell-shaped pulse events in the Crystal-1 in the later quarter of the data, we trained another BDT (BDTA) using the variance of charge–weighted mean time (MV : Eq. 3.8), the charge ratios of waveform leading edges (X3 : Eq. 3.4 shown in d) of Fig. 1 and X4 : Eq. 3.5), the charge–weighted mean time(MT : Eq. 3.6 and MTL : Eq. 3.7), the charge accumulation time (CAT : Eq. 3.11 shown in f) of Fig. 1) and the energy of the event. These effectively identify the shape distortion compared to the regular scintillation signals. Unlike the previous BDT training, early three quarter of the Crystal-1 data as a good data and the later quarter of the data as a noise-contained data are used for signal and background in the training process and later the same training BDT output is applied to all other crystals.

The single-site event spectrum is obtained after the application of the selection criteria and their efficiencies are measured from the 60Co multiples. On average, event selection efficiency of 70% at 2 keVee is obtained for six crystals except Crystal-5 and Crystal-8 which show higher energy threshold at 4 keVee and 8 keVee, respectively, due to their low light yield. The progression of event selection for Crystal-7 is shown in Fig. 2.

The selection efficiency of these event selections are separately checked with neutron-induced nuclear recoils obtained with a small rectangular (2 cm$\times$2 cm$\times$1.5 cm) NaI(Tl) crystals from the same ingots of the detector crystals in front of 2.42 MeV mono-energetic neutron beams [39]. The efficiency of nuclear recoil events in energies between 2 and 20 keVee is consistent with the 60Co calibration efficiency, which indicates that the selections do not affect the possible WIMP signal region.

To compare the results of COSINE–100 and DAMA/LIBRA–phase2 we start from the best–fit analysis of the DAMA/LIBRA–phase2 modulation effect in terms of NR WIMP effective models in Ref. [38], as summarized in Table 3. In such table each of the NR couplings is assumed to be the only term in the effective Hamiltonian of Eq. (2.6) and the $\chi^{2}$:

[TABLE]

(where we consider 14 energy bins, of 0.5 keVee width, from 1 keVee to 8 keVee, and one high–energy control bin from 8 keVee to 16 keVee) is minimized in terms of the WIMP mass $m_{\chi}$, the neutron–over–proton coupling ratio $r\equiv c^{n}/c^{p}$ and of the effective cross section $\sigma_{p}$ as defined in Eq. (2.14). In Eq. (3.12) $S_{m,k}\equiv S_{m,[E_{k}^{\prime},E_{k+1}^{\prime}]}$ is given by Eq. (2.12), while $S_{m,k}^{exp}$ and $\sigma_{k}$ represent the modulation amplitudes and 1–$\sigma$ uncertainties as measured by DAMA/LIBRA–phase2 [4] and reported in Table 5.

As shown in Table 3 for each NR coupling two local minima of the $\chi^{2}$ are found (low and high– mass) with the low–mass solution corresponding in all cases to the absolute minimum. With the exception of the high–mass minima for the ${\cal O}_{4}$ and ${\cal O}_{7}$ operators, all the $\chi^{2}$ minima of Table 3 are acceptable with $15-3$ degrees of freedom. In particular, the DAMA/LIBRA–phase2 has lowered the energy threshold to 1 keVee, implying that it is sensitive to WIMP–iodine scattering events for WIMP masses below $\simeq$ 20 GeV/c2. In the SI case this requires to highly tune the parameters to suppress the iodine contribution, in order to avoid an otherwise too steeply increasing spectrum at low energy of the modulation amplitudes compared to the data [40]. On the other hand, if the WIMP–nucleus cross section is driven by other operators the fine tuning required to suppress iodine is reduced and/or the hierarchy between the WIMP–iodine and the WIMP–sodium cross section is less pronounced in the first place [38].

The raw and observed WIMP spectra corresponding to the best–fit values of Table 3 for operators $c_{1}$–$c_{15}$ are shown in Fig. 3. The raw spectra are generated by using Eq. (2.4), while the energy resolutions obtained from individual crystal measurements and the DAMA quenching factors (0.3 for Na and 0.09 for I) are applied to the raw signal to create the WIMP signal models.

To test the presence of a WIMP signal in the COSINE–100 data that is consistent to the modulation effect measured by DAMA/LIBRA, we generated WIMP spectra for 5 mass values centered around each of the low–mass and high–mass modulation minima of Table 3.

To extract the WIMP signal from our data, a Bayesian approach with a likelihood function based on the Poisson probability is used. A WIMP is not expected to have multiple scatterings within our detector volume, so our WIMP search window is the 2–20 keVee region in the single–hit spectrum where the average event rate of all crystal is recorded to 3.5 counts/day/kg/keV. All crystals are fitted together with crystal-specific background models and a single WIMP signal at a given mass. The constraints are applied as 1 $\sigma$ Gaussian priors for those obtained from the background understanding. Similarly, systematic parameters that change the shape of the background distributions are added as nuisance parameters.

In this way for each NR coupling and WIMP mass value we produce a posterior probability of the effective cross section $\sigma_{p}$. Examples of the posterior probability versus $\sigma_{p}$ are provided in Fig. 4 for the couplings $c_{1}$, $c_{3}$ and $c_{13}$. For all the couplings $c_{1}$–$c_{15}$ we find no signal, so a 90% confidence level (C.L.) upper limit is obtained by integrating the posterior probability from zero.

The result of our analysis is summarized in Figs. 5 and 6 where, for each NR effective coupling, the DAMA modulation regions at 1–$\sigma$, 3–$\sigma$ and 5–$\sigma$ is compared to the COSINE–100 90% C.L. exclusion limit in the $m_{\chi}$–$\sigma_{p}$ plane. The present analysis is focused on the comparison between DAMA and COSINE-100. However, it is worth reminding here that limits from detectors using target materials different from $NaI$ exclude effective cross sections ranging from one to three orders of magnitude below the DAMA region [38].

In each plot of Fig. 5 the neutron–over–proton coupling ratio $r=c^{n}/c^{p}$ is fixed to the corresponding low–mass best–fit value of Table 3 and the WIMP mass interval is centered accordingly. In Fig. 6 the same is done for the high–mass best fit values of Table 3 (in the latter figure the cases of couplings $c_{4}$ and $c_{7}$ are not included since they do not provide a good fit to the DAMA modulation amplitudes).

Figs. 5-6 show that the exclusion plot on $\sigma_{p}$ from COSINE–100 has a different impact on the DAMA best fit modulation region depending on the specific non–relativistic model. Namely, as far as the DAMA low–mass minima of Fig. 5 are concerned, the tension between DAMA/LIBRA and COSINE–100 is maximal for the couplings $c_{1}$ and $c_{4}$, while for the DAMA high–mass minima of Fig. 6 this happens for the couplings $c_{1}$ and $c_{8}$: in all such cases the 90% C.L. bound from COSINE–100, represented by the (blue) solid line, rules out all the 5–sigma DAMA region shown as the (red) dot–dashed contour. On the other hand, in all other cases the 5–sigma DAMA region is not completely excluded by the corresponding 90% C.L. COSINE–100 upper bound, with two instances ($c_{15}$ at low WIMP mass $c_{13}$ at high WIMP mass) for which all the DAMA modulation region is allowed by the COSINE–100 constraint.

The main motivation of probing the modulation effect claimed by DAMA/LIBRA using a sodium–iodide target with COSINE–100 is to obtain results that depend as little as possible on the unknown properties of the WIMP particle. On the other hand, the model dependence observed in Figs. 5–6 is due to two main reasons: i) the change in the signal spectral shape; ii) the expected modulation fractions $S_{m,[E_{1}^{\prime},E_{2}^{\prime}]}/S_{0,[E_{1}^{\prime},E_{2}^{\prime}]}$ in DAMA.

As far as the spectral shape of the expected WIMP signal is concerned, each of the effective models listed in Table 1 is characterized by a different dependence on the exchanged momentum $q$ (and so on the recoil energy $E_{R}=q^{2}/(2m_{T})$), both through the nuclear response functions $W_{Tk}(q)$ (with $k$=$M$, $\Phi^{\prime\prime}$, $\tilde{\Phi}^{\prime}$, $\Sigma^{\prime\prime}$, $\Sigma^{\prime}$, $\Delta$) and through additional powers of $q$ in the scattering amplitude (as summarized in Table 2). The raw energy spectra in COSINE–100 calculated using Eq. (2.4) are shown with a red solid line in Fig. 3 for each NR operator. Indeed, while in the standard spin–independent and spin–dependent cases (corresponding to $c_{1}$ and $c_{4}$) the expected differential rate is the featureless superposition of two exponentially decaying spectra due to WIMP–iodine and WIMP–sodium scattering events, in the case of other NR operators the WIMP recoil spectrum can show a maximum at low energy that may mimic one of the observed radiation peaks (such as the one due to 40K) potentially affecting the sensitivity to the signal. However, as shown in Fig. 3 with the blue dotted lines, when the quenching factors for Na and I and the crystal resolution are applied to the raw spectra all the expected rates are compressed to lower visible energies, so that this effect is strongly reduced. Moreover the expected rates in COSINE–100 become almost insensitive to WIMP–iodine scattering events, that are driven below the 2 keVee threshold. An example of this effect is provided by the $c_{13}$ coupling for which a slight 0.9 $\sigma$ positive fluctuation from zero is observed in the posterior probability, weakening the bound as show in Fig. 6. This may by partially ascribed to the shape of the signal spectrum, as shown in Fig. 3.

However, a much more important source of model dependence in Figs. 5 and 6 is due to the modulation fractions. In fact, the data used in the present paper (and in the result of Ref. [7]) are sensitive to the time–averaged count rate, so that they are used to put upper bounds on the quantity $S_{0,[E_{1}^{\prime},E_{2}^{\prime}]}$ defined in Eq. (2.13). On the other hand, the DAMA effect WIMP interpretation is in terms of the $S_{m,[E_{1}^{\prime},E_{2}^{\prime}]}$ quantities of Eq. (2.12). Crucially, the ratio of the two quantities depends both on the WIMP velocity distribution (for which we assume here a standard Maxwellian as given in Eq. (2.10)) and on the specific operator assumed to dominate in the Hamiltonian of Eq. (2.6) among those in Table 1. The variation of the modulation fraction with the NR model is shown in Table 4, where for each NR operator we provide the minimum value of the modulation fraction $(S^{DAMA}_{m}/S^{DAMA}_{0})_{E^{\prime}<3.5~{}\rm keVee}$ in the three DAMA energy bins for 2 keVee$\leq E^{\prime}\leq$3.5 keVee, where the bulk of the DAMA modulation effect above the COSINE-100 threshold is concentrated. Such variations are an effect of the same modified spectral features discussed in Fig. 3. In particular, larger modulation fractions appear for WIMP scattering events off sodium targets, which are sensitive to the high–speed tail of the velocity distribution, and when the scattering amplitude is multiplied by powers of the transferred momentum $q$ (as summarized in Table 2) due to the enhanced dependence of the expected rate on the $v_{min}$ parameter.

One can notice that in the same energy intervals of Table 4 the time–averaged spectrum in DAMA is above $\simeq$ 0.8 events/day/kg/keVee (see Fig. 1 in [4]). An analysis of the DAMA unmodulated data similar to the one that we perform in the present paper for COSINE–100 is not available and beyond our capability (since it would require a detailed understanding of the background, possibly from Monte Carlo simulations, and access to the systematics of the experiment). If, instead, a vanishing background is conservatively assumed in DAMA, no constraint on the modulation effect can be obtained. In fact the minimal modulation fraction of Table 4 ($\simeq$ 0.065 events/kg/day/keVee) implies an upper bound S_{m}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{<}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt\sim}} 0.05 events/kg/day/keVee, which is above the corresponding measured values of the modulated amplitudes (see Table 5).

As far as COSINE–100 is concerned, a given upper bound on the time–averaged rate $S_{0}$ is converted into a bound on the yearly modulated component $S_{m}$ in DAMA whose strength is inversely proportional to the expected modulated fraction $S_{m}/S_{0}$. In Table 4 one can see that the smallest values for $(S^{DAMA}_{m}/S^{DAMA}_{0})_{E^{\prime}<3.5~{}\rm keVee}$ (and so the strongest limits on the modulated amplitudes) correspond to the standard SI and SD couplings ($c_{1}$ and $c_{4}$) at the level of about 7%. Indeed, as shown in Figs. 5 and 6 such values are small enough for the corresponding COSINE–100 exclusion plots to exclude all the DAMA 5–$\sigma$ region. However, for several models ($c_{3}$, $c_{5}$, $c_{6}$, $c_{13}$ and $c_{14}$ at low WIMP mass and $c_{14}$ at high WIMP mass) the modulation fraction $(S^{DAMA}_{m}/S^{DAMA}_{0})_{E^{\prime}<3.5~{}\rm keVee}$ is as high as $\simeq$ 12%, and about half of the corresponding 5 sigma C.L. DAMA regions are allowed. For most other models the modulation fraction is at an intermediate value $\simeq$0.10 so that most of the DAMA 5-sigma region is excluded. Notice that the two highest numerical values in Table 4 ($\simeq$ 0.14 for the high–mass DAMA best–fit of $c_{13}$ and $\simeq$ 0.15 for the low–mass DAMA best–fit of $c_{15}$) correspond to the two less–constraining cases in Figs. 5 and 6 for which all the DAMA modulation region is allowed by the COSINE–100 bound. Therefore, the lower energy threshold of COSINE-100 would improve the bound at low WIMP masses because WIMP–iodine scattering events in the energy range 1 keVee$\leq E^{\prime}\leq$2 keVee drive $S_{m}/S_{0}$ to lower values. The case of $c_{5}$ in Fig. 6 is a peculiar one: as discussed in [38] the corresponding effective operator leads to a velocity–dependent cross section for which the $\chi$–square can saturate to a constant, acceptable value at large WIMP masses. This explains the peculiar elongated shape for $c_{5}$ in Fig. 6.

4 Conclusions

Assuming a standard Maxwellian for the WIMP velocity distribution, in the present paper we have discussed the bounds from the null WIMP search result of the COSINE-100 experiment on the DAMA/LIBRA–phase2 modulation effect within the context of the non–relativistic effective theory of WIMP–nucleus scattering. To this aim we have systematically assumed that one of the effective operators allowed by Galilean invariance dominates in the effective Hamiltonian of a spin–1/2 DM particle.

We find that, although DAMA/LIBRA and COSINE–100 use the same sodium–iodide target, the comparison of the two results still depends on the particle–physics model. This is mainly due to two reasons: i) the WIMP signal spectral shape; ii) the expected modulation fractions, when the upper bound on the time–averaged rate in COSINE–100 is converted into a constraint on the yearly modulated component in DAMA/LIBRA. We find that the latter effect is the dominant one. In particular, for several effective operators we find that the expected modulation fractions are larger than in the standard spin–independent or spin–dependent interaction cases. As a consequence, for such operators compatibility between the modulation effect observed in DAMA/LIBRA and the null result from COSINE–100 is still possible. COSINE-100 has been taking stable data for more than 2.5 years and 1 keVee threshold analysis is forthcoming. This would improve the bound at low WIMP masses because WIMP–iodine scattering events in the energy range 1 keVee$\leq E^{\prime}\leq$2 keVee drive $S_{m}/S_{0}$ to lower values.

Acknowledgments

The Sogang group’s 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. Jong–Hyun Yoon acknowledges support from the Magnus Ehrnrooth Foundation. We thank the Korea Hydro and Nuclear Power (KHNP) Company for providing underground laboratory space at Yangyang. The work of COSINE-100 is supported by: the Institute for Basic Science (IBS) under project code IBS-R016-A1, Republic of Korea; UIUC campus research board, the Alfred P. Sloan Foundation Fellowship, NSF Grants No. PHY-1151795, PHY-1457995, DGE-1122492 and DGE-1256259, WIPAC, the Wisconsin Alumni Research Foundation, Yale University and DOE/NNSA Grant No. DE-FC52-08NA28752, United States; STFC Grant ST/N000277/1 and ST/K001337/1, United Kingdom; and Grant No. 2017/02952-0 FAPESP and CAPES Finance Code 001, Brazil.