Sensitivity of Solid-Scintillator Detectors to Dark Matter

Abstract

This paper shortly reviews the sensitivities that can be achieved to unambiguously point out the presence of a signal of Galactic origin in dark matter experiments with solid-scintillator detectors. Examples of the experimental sensitivities obtained by exploiting the annual and diurnal modulation signatures are reported with particular regard to the investigations performed in the framework of the DAMA Collaboration. The directionality approach in solid scintillators is also presented and, in particular, the perspectives of the ADAMO project are discussed.

1. Introduction

There is a wide worldwide effort in investigating the dark matter (DM) in the galactic halo through direct detection (for a review see [1]) in terrestrial experiments, in indirect searches (for a review see [2]) and at accelerators (for a review see [3]). In this paper, we focus our attention on the DM direct detection by using model-independent signatures that, as well as other rare searches like the double beta decay processes, the rare nuclear decay processes, etc. allows getting information about either the Standard Model (SM) or beyond the SM of particle physics. In fact, no candidate for DM particle exists in the SM and new Physics in the particle sector must be opened to account for DM. Many DM candidates have been proposed in the literature in theories extending the SM of particle physics (for a review see [4]). On the other side, all the experimental techniques for direct detection of DM particles in terrestrial experiments profit from the large improvements achieved in the low-background technology and of suitable underground laboratories. For a discussion of these techniques see for example [5].

In this paper we describe the achievable sensitivity in time-varying DM signal and we show the results obtained by DAMA/LIBRA in the cases of the DM annual modulation signature and diurnal modulation. Moreover, the case of a novel detector, that exploits the anisotropy of some solid scintillators to the light response of nuclear recoils, is discussed. In particular, the reachable sensitivity for the time-varying DM signal achievable within the ADAMO project is presented.

2. The Achievable Sensitivity in Time-Varying DM Signal

The origin of the DM annual and diurnal modulation signatures and their peculiar features is due to the Earth motion with respect to the DM particles constituting the Galactic dark halo [6,7]. During its motion, a laboratory on the Earth experiences a wind of DM particles, which are trapped in the gravitational field of the Galaxy with a peculiar phase-space distribution function. Many possibilities for the distribution function have been proposed in the literature [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]; many of them profit from new astrophysical observations and new simulations. The DM annual modulation is due to the Earth’s revolution around the Sun, which is moving towards Vega (near the constellation of Hercules) in the Local Standard of Rest (LSR) frame of the Galaxy. In particular, a larger flux of DM particles crosses the Earth around ≃2 June, while a smaller one does around ≃2 December. On the contrary, the DM diurnal modulation is due to the rotation of the Earth around its axis.

The DM annual and diurnal modulation signatures are very distinctive. As regards the DM annual modulation signature, the DM signal must simultaneously satisfy all the following requirements (see for example [1,25]): the rate has to contain a modulation component as a cosine function with a period of one year and a phase that provides the maximum of the signal around ≃2 June. The DM annual modulation is to be found only at low energy, where events induced by DM particles can be present, and can produce events where only one detector module of many actually ”fires” (single-hit events), since the probability that DM particles interact with more than one detector module is negligible. There is also a constraint on the DM annual modulation amplitude: it is $\lesssim$$7\%$ (see later and e.g., [1,25]) of the constant part of the signal in the region of maximal sensitivity for the halo distribution function usually adopted; however, in the literature, some scenarios have been also proposed in which a larger amplitude is expected [26,27,28,29,30]. Thus the DM annual modulation signature does not depend on the nature of the DM particle, has many requirements (as fore-mentioned) and allows a wide range of parameters to be tested in many possible astrophysical, nuclear and particle physics scenarios (see e.g., [1]). Similar considerations can be done for the DM diurnal modulation. In this case, the behavior of the counting rate during the sidereal day must follow a well-defined pattern depending on the location of the detector in the Earth.To explain the DM annual and diurnal modulation signatures, the motion of a detector placed in a laboratory on the Earth surface is described in the Galactic frame. Thus, an observer fixed in it sees the detector moving due both to the rotation of the Earth around its axis and to the revolution of the Earth around the Sun.

The velocity of such a detector can be written as [31]:

$${\overrightarrow{v}}_{\mathrm{lab}}\left(t\right)={\overrightarrow{v}}_{\mathrm{LSR}}+{\overrightarrow{v}}_{\odot }+{\overrightarrow{v}}_{\mathrm{rev}}\left(t\right)+{\overrightarrow{v}}_{\mathrm{rot}}\left(t\right),$$

(1)

where ${\overrightarrow{v}}_{\mathrm{LSR}}$ is the velocity of the LSR due to the rotation of the Galaxy; ${\overrightarrow{v}}_{\odot }$ is the peculiar velocity of the Sun with respect to the LSR; ${\overrightarrow{v}}_{\mathrm{rev}}\left(t\right)$ is the velocity of the Earth due to its revolution around the Sun and ${\overrightarrow{v}}_{\mathrm{rot}}\left(t\right)$ is the velocity of the Earth due to the rotation around its axis at the coordinates (latitude and longitude) of the laboratory. In the Galactic coordinate frame (defined as x axis pointing towards the Galactic center, y axis following the Galaxy rotation and the z axis pointing towards the Galactic North pole), we can write ${\overrightarrow{v}}_{\mathrm{LSR}}=(0,v_0,0)$, where $v_0=(220\pm 50)$ km/s [32,33,34,35] (uncertainty at 90% confidence level (C.L.)), and ${\overrightarrow{v}}_{\odot }=(9,12,7)$ km/s [36]. The Earth velocities due to the revolution and rotation depend on the sidereal time, t.

The orbital speed of the Earth $|{\overrightarrow{v}}_{\mathrm{rev}}|=V_{\mathrm{Earth}}$ depends weakly on the time because of the ellipticity of the Earth orbital motion around the Sun; in fact, $V_{\mathrm{Earth}}$ varies in the range between 29.3 km/s and 30.3 km/s. Hence, here it is assumed constant and equal to its mean value ≃29.8 km/s.

The rotational speed of the Earth, $|{\overrightarrow{v}}_{\mathrm{rot}}|=V_{\mathrm{r}}$ is given by $V_{\mathrm{r}}=V_{\mathrm{eq}}cos{\varphi }_0$, where ${\varphi }_0$ is the latitude of the considered laboratory and $V_{\mathrm{eq}}$ is the equatorial rotational speed equal to 0.4655 km/s. As an example, $V_{\mathrm{r}}=0.3435$ km/s at Gran Sasso underground Laboratory (LNGS), where ${\varphi }_0={42}^{\circ }$27${}^{\prime }$ N and longitude ${\lambda }_0={13}^{\circ }$34${}^{\prime }$ E.

The expected counting rate in many detectors for DM investigations depends on the detector’s speed in the Galaxy, $v_{\mathrm{lab}}\left(t\right)$. The time-independent contribution is $\left|{\overrightarrow{v}}_{\mathrm{s}}\right|=\left|{\overrightarrow{v}}_{\mathrm{LSR}}+{\overrightarrow{v}}_{\odot }\right|\approx 232\pm 50$ km/s. A Taylor expansion, performed in the smaller time-dependent contributions $\left|{\overrightarrow{v}}_{\mathrm{rev}}\left(t\right)\right|\approx 30$ km/s and $\left|{\overrightarrow{v}}_{\mathrm{rot}}\left(t\right)\right|\approx 0.34$ km/s, allows to determine—at the first order expansion—the detector’s speed in the Galaxy [31]:

$$v_{\mathrm{lab}}\left(t\right)\simeq v_{\mathrm{s}}+V_{\mathrm{Earth}}{A}_{\mathrm{m}}cos\left[\omega (t-t_0)\right]+V_{\mathrm{r}}{A}_{\mathrm{d}}cos\left[{\omega }_{\mathrm{rot}}\left(t-t_{\mathrm{d}}\right)\right],$$

(2)

where $A_{\mathrm{m}}\simeq 0.489$, confirms that the ecliptic is tilted with respect to the Galactic plane of ≃60${}^{\circ }$, and $A_{\mathrm{d}}\simeq 0.671$ (calculated for $v_0=220$ km/s) are “geometrical” factors [31]. The angular velocities are $\omega =2\pi /T$ and ${\omega }_{\mathrm{rot}}=2\pi /T_{\mathrm{rot}}$, where $T=1$ sidereal year and $T_{\mathrm{rot}}=1$ sidereal day. The two phases, $t_0$ and $t_{\mathrm{d}}$, correspond to the time when the annual and diurnal modulations are maximum. The phase of the DM annual modulation is $t_0=t_{\mathrm{equinox}}+73.25\mathrm{solar}\mathrm{days}$, where $t_{\mathrm{equinox}}$ is the spring equinox time (≈21 March). The second term in the previous equation has been calculated for $v_0=220$ km/s and it ranges from 71.76 solar days (for $v_0=170$ km/s) to 74.20 solar days (for $v_0=270$ km/s). For $v_0=220$ km/s it corresponds to $t_0\approx$ 2 June [31]. The phase of the DM diurnal modulation corresponds to $t_{\mathrm{d}}\simeq 14.02$ h local sidereal time and $v_0=220$ km/s; actually this value ranges from 13.94 h ($v_0=170$ km/s) to 14.07 h ($v_0=270$ km/s) [31].

Finally, the expected signal counting rate in a given k-th energy bin can be written, applying a Taylor expansion, as:

$$S_k\left[v_{\mathrm{lab}}\left(t\right)\right]\simeq S_k\left[v_{\mathrm{s}}\right]+{\left[\frac{\partial S_k}{\partial v_{\mathrm{lab}}}\right]}_{v_{\mathrm{s}}}\left[V_{\mathrm{Earth}}{A}_{\mathrm{m}}cos\omega (t-t_0)\right.+\left.V_{\mathrm{r}}{A}_{\mathrm{d}}cos{\omega }_{\mathrm{rot}}\left(t-t_{\mathrm{d}}\right)\right].$$

(3)

For simplicity, we omit higher order terms with no time dependence and higher harmonics of $\omega$ and ${\omega }_{\mathrm{rot}}$. The first term of Equation (3) is the constant part of the signal ($S_0$), while the other two terms are the annual modulation one with amplitude $S_{\mathrm{m}}={\left[\frac{\partial S_k}{\partial v_{\mathrm{lab}}}\right]}_{v_{\mathrm{s}}}{V}_{\mathrm{Earth}}{A}_{\mathrm{m}}$, and the diurnal modulation one with amplitude $S_{\mathrm{d}}={\left[\frac{\partial S_k}{\partial v_{\mathrm{lab}}}\right]}_{v_{\mathrm{s}}}{V}_{\mathrm{r}}{A}_{\mathrm{d}}$.

It is worth to note that the ratio $R_{\mathrm{dy}}$ of the DM diurnal modulation amplitude over the DM annual modulation amplitude is a model-independent constant; it is energy independent being a function of the latitude of the experiment and its order of magnitude is given by the ratio between the rotational speed of the Earth and the orbital speed of the Earth (≃${10}^{-2}$):

$$R_{\mathrm{dy}}=\frac{S_{\mathrm{d}}}{S_{\mathrm{m}}}=\frac{V_{\mathrm{r}}{A}_{\mathrm{d}}}{V_{\mathrm{Earth}}{A}_{\mathrm{m}}}\simeq 0.0214cos{\varphi }_0;$$

(4)

At the LNGS latitude one has:

$$R_{\mathrm{dy}}=\frac{S_{\mathrm{d}}}{S_{\mathrm{m}}}\simeq 0.016.$$

(5)

The data analysis exploiting the DM annual and diurnal modulation signatures are based on the study of the measured counting rate as a function of time. The initial approach [6,7], that was considered for the DM annual modulation signature, was introduced considering experiments of small masses, high background and poor exposure. Over time the method was extended, in particular, by the DAMA collaboration [37,38,39]. Generally speaking, we can write the time-dependence of the expected signal counting rate in a given k-th energy bin as:

$${\mathcal{S}}_k\left(t\right)={\mathcal{S}}_{0,k}+{\mathcal{B}}_{k}f\left(t\right).$$

(6)

where $f\left(t\right)$ is a periodic function with period T; without loosing generality, we assume that ${\int }_0^{T}f\left(t\right)dt=0$ and ${\int }_0^T{f}^2\left(t\right)dt=0.5$. For the two particular cases we are dealing with here, $f\left(t\right)=cos\omega (t-t_0)$ and $f\left(t\right)=cos{\omega }_{\mathrm{rot}}(t-t_{\mathrm{d}})$, respectively, and in turn ${\mathcal{B}}_k={\mathcal{S}}_{m,k}$ and ${\mathcal{B}}_k={\mathcal{S}}_{d,k}$. In Section 4 the example of the directionality approach using the ZnWO${}_4$ anisotropic scintillators will be addressed in the light of the sensitivity discussed here; in such a case $f\left(t\right)$ is a more complex function with the daily periodicity achievable for example by simulations. Therefore, the ${\mathcal{B}}_k$ are the amplitudes of the studied process.

To get the highest sensitivity, the energy bin $\Delta E$ in the differential energy distribution must be chosen as the best compromise between a sufficient statistics and a good signal-to-background ratio. When the statistics is low, it requires to integrate the events in rather large energy intervals at the expense of the reachable sensitivity. When a large exposure is instead available, the events can be binned on the basis of the arrival time (i index), involved detector module (j index) and energy intervals of width $\Delta E$ (k index): $N_{ijk}$. Each element of $N_{ijk}$ follows a Poisson distribution with expected value:

$${\mu }_{ijk}=\left[b_{jk}+{\mathcal{S}}_{0,k}+{\mathcal{B}}_{k}f\left(t_i\right)\right]M_j\Delta t_i\Delta E{ϵ}_{jk}.$$

(7)

The contribution of background in the k-th energy interval for the j-th detector module is $b_{jk}$, the time interval in which the detector was active on the i-th time bin is $\Delta t_i$, the total efficiency in the k-th energy interval for the j-th detector module is ${ϵ}_{jk}$, and the mass of the j-th detector module is $M_j$. Let us define $w_{ijk}$ = $M_j\Delta t_i\Delta E{ϵ}_{jk}$. One can consider the following variables–the sum run over a single period T–built from the number of counts, $N_{ijk}$ (omitting for simplicity the indexes $j,k$):

$$\begin{array}{ccc}\hfill D& =& \frac{1}{W}\sum _i{N}_i,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill X& =& \frac{2}{W}\sum _i\left[f\left(t_i\right)-\beta \right]N_i,\hfill \end{array}$$

(9)

where:

$$\begin{array}{ccc}\hfill W& =& \sum _i{w}_i,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \beta & =& \langle f\left(t\right)\rangle =\frac{1}{W}\sum _i{w}_{i}f\left(t_i\right),\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill \alpha & =& \langle f^2\left(t\right)\rangle =\frac{1}{W}\sum _i{w}_i{f}^2\left(t_i\right).\hfill \end{array}$$

(12)

In the particular case of “coupled” time bins (i.e., opposite in phase), one obtains $\beta =0$, and in case of a complete period T of measurement $\alpha =0.5$. The D and X variables are stochastic ones since they are function of $N_i$; the latter follows a poissonian distribution with $\langle N_i\rangle =Var\left(N_i\right)={\mu }_i$. These variables can allow interesting studies on the free parameters: the modulation amplitude, $\mathcal{B}$, and the time-independent term, $(b+{\mathcal{S}}_0)$. In fact, in case the variance of the cosine, $(\alpha -{\beta }^2)$, is $\approx 0.5$ for a detector being operational evenly throughout the period T, the expectation values of D and X are $\langle D\rangle =(b+{\mathcal{S}}_0)$ and $\langle X\rangle =\mathcal{B}$. In the general case, they are $\langle D\rangle =(b+{\mathcal{S}}_0)+\beta \mathcal{B}$ and $\langle X\rangle =2(\alpha -{\beta }^2)\mathcal{B}$ (For completeness, another stochastic variable can be used to decouple the two equations: $Z=\frac{1}{W}{\sum }_{i}2\left[\alpha -\beta f\left(t_i\right)\right]N_i$, whose expectation value is directly connected to $(b+{\mathcal{S}}_0)$: $\langle Z\rangle =2(\alpha -{\beta }^2)(b+{\mathcal{S}}_0)$.). Moreover, the variances are: $Var\left(D\right)=\frac{(b+{\mathcal{S}}_0)+\beta \mathcal{B}}{W}$ and $Var\left(X\right)\simeq \frac{4(\alpha -{\beta }^2)(b+{\mathcal{S}}_0)}{W}$ (in the latter we omit the terms negligible when a detector is operational evenly throughout the period T).

In conclusion, the modulation amplitude, $\mathcal{B}$, can be worked out by

$$\mathcal{B}=\frac{\langle X\rangle }{2(\alpha -{\beta }^2)}=\frac{1}{W(\alpha -{\beta }^2)}\sum _i\left[f\left(t_i\right)-\beta \right]\langle N_i\rangle ,$$

(13)

and considering the $N_i$ as the best estimators of $\langle N_i\rangle$ (that is $\langle N_i\rangle =Var\left(N_i\right)\approx N_i$) one has:

$$\begin{array}{ccc}\hfill \mathcal{B}& =& \frac{{\sum }_i\left[f\left(t_i\right)-\beta \right]N_i}{W(\alpha -{\beta }^2)},\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill Var\left(\mathcal{B}\right)& =& \frac{{\sum }_i{\left[f\left(t_i\right)-\beta \right]}^2{N}_i}{W^2{(\alpha -{\beta }^2)}^2}.\hfill \end{array}$$

(15)

To estimate the statistical significance of the presence of modulation one can evaluate the standard deviation associated to the $\mathcal{B}$. As a first approximation without losing generality, we can replace in Equation (14) $N_i$ with $w_i(b+{\mathcal{S}}_0)$:

$$\begin{array}{ccc}\hfill Var\left(\mathcal{B}\right)& \simeq & (b+{\mathcal{S}}_0)\frac{{\sum }_i{\left[f\left(t_i\right)-\beta \right]}^2{w}_i}{W^2{(\alpha -{\beta }^2)}^2}.\hfill \end{array}$$

(16)

After some algebraic calculations using Equations (10)–(12):

$$\sigma \left(\mathcal{B}\right)=\sqrt{Var\left(\mathcal{B}\right)}\simeq \sqrt{\frac{b+{\mathcal{S}}_0}{W(\alpha -{\beta }^2)}}.$$

(17)

Therefore, the obtained sensitivity, defined as $R=\mathcal{B}/\sigma \left(\mathcal{B}\right)$, for observing periodical modulation amplitude is directly proportional to the square root of the generalized exposure (i.e., $W(\alpha -{\beta }^2)$) and inversely proportional to the square root of the mean rate (i.e., $\langle rate\rangle \approx (b+{\mathcal{S}}_0)$). Finally, considering that for a full period T of exposure $\alpha \simeq 0.5$ and $\beta \simeq 0$, one can write:

$$\sigma \left(\mathcal{B}\right)\simeq \sqrt{2\frac{\langle rate\rangle }{W}}.$$

(18)

Of course, when the DM signal contains two or more $f\left(t\right)$ with very different frequencies, as it is the case of annual and diurnal modulations, all these contributions can be treated as independent and all the above considerations still hold separately for each contribution.

In the following, as examples, we consider the cases of the DM annual and diurnal modulation results reported by the DAMA experiments. Only for illustrative purpose, we consider that $\langle rate\rangle \simeq 0.8$ cpd/kg/keV. The exposure, as defined in Equation (10), is $W={\sum }_i{M}_{\mathrm{set}\mathrm{up}}\Delta t_i\Delta E\overline{ϵ}=M_{\mathrm{set}\mathrm{up}}{T}_{\mathrm{tot}}\Delta E\overline{ϵ}$. Typical values of the DAMA experiments are: $\Delta E=4$ and 5 keV, because the DM annual modulation is observed in the (1–6) keV (DAMA/LIBRA–phase2) and (2–6) keV (DAMA/NaI and DAMA/LIBRA–phase1) energy intervals; and the average efficiency is $\overline{ϵ}\approx 0.7$. For the DM annual modulation study the total exposure of DAMA/NaI, DAMA/LIBRA–phase1 and phase2 is $M_{\mathrm{set}\mathrm{up}}{T}_{\mathrm{tot}}=2.46$tons × year [1,40,41], where $M_{\mathrm{set}\mathrm{up}}$ is the total exposed mass and $T_{\mathrm{tot}}$ is the time of data acquisition; thus the expected error on ${\mathcal{S}}_m$ is:

$$\sigma \left({\mathcal{S}}_m\right)\simeq 8\times {10}^{-4}\mathrm{cpd}/\mathrm{kg}/\mathrm{keV}.$$

(19)

Other published results on annual modulation with solid scintillators were recently reported by COSINE-100 [42] and by ANAIS [43]. As deeply argued e.g., in [1] these results have no impact on the observed DAMA annual modulation result, as briefly summarized here. For example, considering the exposures (0.098 and 0.158 tons × year, respectively), the counting rates in the low energy region (3.0 cpd/kg/keV in the 2–6 keV energy region for COSINE-100 and 3.2 (3.6) cpd/kg/keV in the 2–6 (1–6) keV energy region for ANAIS) and the related efficiencies (≃$0.85$ (0.65) in the 2–6 (1–6) keV energy region in either cases), for those experiments, we can estimate by using Equation (18) $\sigma \left({\mathcal{S}}_m\right)\simeq 7\times {10}^{-3}$ cpd/kg/keV for COSINE-100 and $\sigma \left({\mathcal{S}}_m\right)\simeq 5.7\left(6.2\right)\times {10}^{-3}$ cpd/kg/keV in the 2–6 (1–6) keV energy region for ANAIS, respectively. These values are about one order of magnitude larger than those obtained by DAMA, confirming the fact that these experiments are still far from reaching the sensitivity of DAMA. Let us note that the energy intervals in kinetic recoiling energy are not the same, as briefly discussed later.

For the DM diurnal modulation study the exposure of DAMA/LIBRA–phase1 has been considered to date [31]: $M_{\mathrm{set}\mathrm{up}}{T}_{\mathrm{tot}}=1.04$ tons × year and, then:

$$\sigma \left({\mathcal{S}}_{\mathrm{d}}\right)\simeq 1.2\times {10}^{-3}\mathrm{cpd}/\mathrm{kg}/\mathrm{keV}.$$

(20)

The values of Equations (19) and (20) will be compared with those obtained by the analyses summarized in the following. Finally, it can be noted that an enhancement of the sensitivity for the DM signal by pursuing the annual and/or diurnal modulation signatures can be obtained not only by increasing the exposure and by decreasing the background rate of the experiments, but also by increasing the average efficiency and by extending the data analysis to lower the energy threshold. In particular, the latter point is important both because it means increasing the $\Delta E$ in the denominator of Equation (18) and because the un-modulated part of the DM signal is expected to have a steep increase at low energy.

3. Practical Examples of Reached Sensitivities: The Case of the DAMA/LIBRA Annual and Diurnal Investigations

The model-independent DM annual and diurnal modulation signatures will be considered in more details in the light of the results of DAMA/LIBRA (see e.g., [1] and references therein) and of the former DAMA/NaI (see e.g., [38,39] and references therein) experiments. These experiments developed highly radio-pure NaI(Tl) target-detectors that ensure sensitivity to a wide range of DM candidates, interaction types and astrophysical scenarios [1].

The DAMA/LIBRA experiment, located deep underground at the Gran Sasso National Laboratory, is composed of about 250 kg highly radio-pure NaI(Tl). At the end of 2010 all the photomultipliers (PMTs) of the set-up were replaced by new Hamamatsu R6233MOD PMTs, having higher quantum efficiency (Q.E.) and lower background with respect to those previously used [44]. The commissioning of this phase of the experiment (DAMA/LIBRA–phase2) was successfully performed in 2011, achieving the software energy threshold at 1 keV (while it was 2 keV for the first phase, DAMA/LIBRA–phase1), and the improvement of some detector’s features [44]. Along with the data taking, within each annual cycle, all the procedures have been the same. The energy calibrations are regularly carried out in the same running condition down to a few keV with X-rays/$\gamma$ sources. In particular, the study of the double coincidences (summing the data over long periods) due to internal X-rays from ${}^{40}$K (at ppt levels in the crystals) allows to calibrate the detectors at 3.2 keV, i.e., rather near to the energy threshold considered for the data analysis. The light response of the detector modules in DAMA/LIBRA-phase2 is from 6 to 10 photoelectrons/keV, depending on the detector module. The overall efficiency as a function of the energy for single-hit events is reported in [44]. In the experiment both single-hit events (where just one of the detector modules fires) either multiple-hit events (where more than one detector module fires) with energy up to a few MeV have been recorded.

In [1,40,41,45] the details of the annual cycles of DAMA/LIBRA–phase2 are reported. The available full annual cycles of DAMA/LIBRA–phase2 are six (exposure of 1.13 tons × year); when also considering the former DAMA/NaI and DAMA/LIBRA–phase1, the cumulative exposure is 2.46 tons × year [1,40,41].

To point out the presence of annual modulation of the single-hit events at low energy, the residual counting rate as a function of the time has been investigated. The time of measurement is discretized in various time intervals, and each experimental point is calculated from the raw rate of the single-hit scintillation events after subtracting the un-modulated part. In Figure 1 the obtained residual rates for single-hit events for DAMA/LIBRA–phase1 and DAMA/LIBRA–phase2 are shown The energy interval is from the software energy threshold of DAMA/LIBRA–phase1 (2 keV), up to 6 keV. The null hypothesis, that is the absence of modulation, is rejected at very high C.L. by ${\chi }^2$ test: ${\chi }^2$/degrees-of-freedom (d.o.f.) = 199.3/102, corresponding to p-value = 2.9 × 10${}^{-8}$.

Other analyses are reported in Refs. [1,40,41,45]. The single-hit residual rates shown in Figure 1 and those of DAMA/NaI have been fitted with the function: $Acos\omega (t-t_0)$, considering a period $T=\frac{2\pi }{\omega }=1$ year and a phase $t_0=152.5$ day (2 June) as expected by the DM annual modulation signature. The obtained ${\chi }^2$/d.o.f. is equal to 113.8/138 and a modulation amplitude $A=(0.0102\pm 0.0008)$ cpd/kg/keV is obtained [1,40,41,45]. When the period and the phase are kept free in the fitting procedure, the achieved C.L. for the full exposure (2.46 tons × year) is 12.9 $\sigma$; the modulation amplitude of the single-hit scintillation events is: $(0.0103\pm 0.0008)$ cpd/kg/keV, the measured phase is $(145\pm 5)$ days and the measured period is $(0.999\pm 0.001)$ year, all these values are well in agreement with those expected for DM particles. The reached confidence level is 12.9 $\sigma$ [1,40,41,45].

Let us note that the errors of the measured modulation amplitudes are well in agreement with those achievable using the theoretical expectation of Equations (17) and (19): $\sigma \left({\mathcal{S}}_m\right)\simeq 8\times {10}^{-4}\mathrm{cpd}/\mathrm{kg}/\mathrm{keV}$.

The absence of any significant background modulation has been verified in energy regions not of interest for DM [1,40,41,45]. These analyses account for any kind of background. No fake events able to reproduce the DM annual modulation signature (that is able to simultaneously satisfy all the peculiarities of the signature and to account for the measured modulation amplitude) are available; for details also see the discussions e.g., in [1] and references therein.

The above mentioned COSINE-100 [42] and by ANAIS [43] experiments reported modulation amplitudes in the (2–6) keV energy interval compatible with zero and with the DAMA results: $(0.0083\pm 0.0068)$ cpd/kg/keV [42] and $-(0.0044\pm 0.0058)$ cpd/kg/keV [43], respectively. The errors of these measurements are in agreement with those estimated in Section 2, where this lack of needed sensitivity was already been pointed out. Moreover, let us note that the comparison is done in the same energy interval in keV electron equivalent (keVee). This does not correspond to the same interval in terms of kinetic energy, since the quenching factors of nuclear recoils are not the same. See discussions in [1].

Another approach to study the annual modulation is based on the use of the maximum likelihood method to obtain from the residual rate the modulation amplitudes, ${\mathcal{S}}_m$, as a function of the energy, considering $T=1$ year and $t_0=152.5$ day. The likelihood function of the single-hit events in the k-th energy bin can be written as: ${\mathbf{L}}_{\mathbf{k}}={\Pi }_{ij}{e}^{-{\mu }_{ijk}}\frac{{\mu }_{ijk}^{N_{ijk}}}{N_{ijk}!}$. Hence, by minimizing for each energy bin the logarithm of the likelihood function (the index k is omitted for simplicity): $y=-2ln\left(\mathbf{L}\right)-\mathrm{const}$, the free parameters of the fit ($(b+{\mathcal{S}}_0)$ and ${\mathcal{S}}_m$) can be obtained.

Figure 2 reports the modulation amplitudes for DAMA/NaI, DAMA/LIBRA–phase1 and DAMA/LIBRA–phase2. The total exposure is 2.46 tons × year. The data points below 2 keV are only due to DAMA/LIBRA–phase2 (exposure of 1.13 tons × year). A clear modulation signal is evident in the (1–6) keV energy interval, while the ${\mathcal{S}}_m$ are compatible with zero above 6 keV. This analysis confirms the previous ones. The hypothesis that ${\mathcal{S}}_m$ values in the (6–14) keV energy interval have random fluctuations around zero can be tested by a ${\chi }^2$ analysis: the ${\chi }^2$ is equal to 19.0 for 16 d.o.f. (upper tail probability of 27%). For the case of (6–20) keV energy interval, the ${\chi }^2$/d.o.f. is 42.6/28 (upper tail probability of 4%). The ${\chi }^2$ value is rather large mainly due to two data points, far away from the (1–6) keV energy interval, at 16.75 and 18.25 keV. Excluding only the first and either the points the p-values are 11% and 25%, respectively.

Many other consistency checks have been done; for details see [1] and references therein.

No modulation has been found in any possible source of systematics or side reactions; thus, cautious upper limits on possible contributions to the DAMA/LIBRA measured modulation amplitude have been obtained (see e.g., [1] and references therein). Similar analyses have also been performed for the DAMA/NaI data [38,39]. No systematic effects or side reactions able to account for the whole observed modulation amplitude and to simultaneously satisfy all the requirements of the exploited DM signature have been found. A detailed discussion about all the related arguments can be found in [1] and references therein.

The long-standing annual-modulation evidence measured in DAMA experiments is model-independent, i.e., without any a priori assumption of theoretical interpretations of the identity of DM, of the specifics of its interactions, of dark halo features. It is compatible with a wide set of scenarios regarding the nature of the DM candidate and related astrophysical, nuclear and particle physics [1] (and references therein). In particular, discussions about comparisons with other activities can be found e.g., in [1,25] and references therein.

The DAMA/LIBRA–phase1 data have also been analyzed in terms of the DM diurnal modulation signature due to the Earth’s diurnal rotation around its axis [31]. As shown in Equation (5), the interest in this signature is also because the ratio $R_{\mathrm{dy}}$ of the DM diurnal modulation amplitude over the DM annual modulation amplitude is a model-independent constant at a given latitude whose value at the LNGS latitude is 0.016.

Figure 3 shows the time and energy behaviour of the experimental residual rates of single-hit events both as a function of solar (left) and of sidereal (right) time, in the (2–6) keV energy interval. The used time bin is 1 h (either solar or sidereal, respectively). The null hypothesis (absence of diurnal variation in the residual rate) has been tested by a ${\chi }^2$ test, obtaining the (2–6) keV energy interval ${\chi }^2$/d.o.f. $=25.8/24$ and $21.2/24$ for the solar and sidereal time, respectively [31]. The upper tail probabilities (p-values), calculated by the standard ${\chi }^2$ distribution, are 36% and 63%, respectively [31]. Thus, no diurnal variation with a significance of 95 % C.L. is found. In addition to the ${\chi }^2$ test, an independent statistical test, the runs test [46], has been applied: it verifies the hypothesis that the positive and negative data points are randomly distributed. The lower tail probabilities are equal to: 7% and 78% in the (2–6) keV energy interval for the two cases. Thus, in conclusion, the presence of any significant diurnal variation and of time structures can be excluded at the reached level of sensitivity and upper limits on the diurnal modulation amplitudes can be derived.

The residual rates of the single-hit events as a function of the sidereal time (see Figure 3 right) have been fitted with a cosine function with amplitude $A_{\mathrm{d}}^{exp}$ as free parameter, period fixed at 24 h and phase at 14 h. The obtained diurnal modulation amplitude is compatible with zero: $A_{\mathrm{d}}^{exp}=-(1.0\pm 1.3)\times {10}^{-3}$ cpd/kg/keV (${\chi }^2$/d.o.f. $=20.6/23$, and p = 61%) [31]. Following the Feldman-Cousins procedure [47], the upper limit on the diurnal modulation amplitude can be worked out: $A_{\mathrm{d}}^{exp}\le 1.2\times {10}^{-3}$ cpd/kg/keV (90% C.L.) [31].

The single-hit residual rate of DAMA/LIBRA–phase1 as function of the sidereal time has been fitted with the formula $A_{\mathrm{d}}cos\left[{\omega }_{\mathrm{rot}}(t-t_{\mathrm{d}})\right]$; the period is fixed at 24 h and the phase at 14 h, as expected for the DM diurnal effect. The obtained diurnal modulation amplitudes, $A_{\mathrm{d}}$, as functions of energy with 1 keV energy bin are reported in Figure 4.

The $A_{\mathrm{d}}$ values are compatible with zero and show random fluctuations around zero with ${\chi }^2$ equal to 19.5 for 18 degrees of freedom. Moreover, as in the case of DM annual modulation results, the errors of the measured DM diurnal modulation amplitudes are well in agreement with those achievable using the theoretical expectation of Equations (17) and (20): $\sigma \left({\mathcal{S}}_{\mathrm{d}}\right)\simeq 1.2\times {10}^{-3}\mathrm{cpd}/\mathrm{kg}/\mathrm{keV}$.

Finally, taking into account $R_{\mathrm{dy}}$ and the DM annual modulation effect pointed out by DAMA/LIBRA for the single-hit events in the low energy region, the expected value of the DM diurnal modulation amplitude for the (2–6) keV energy interval is $\simeq 1.5\times {10}^{-4}$ cpd/kg/keV. Thus, the effect of DM diurnal modulation, expected because of the Earth diurnal motion on the basis of the DAMA DM annual modulation results, is out the sensitivity of DAMA/LIBRA–phase1 [31]. DAMA/LIBRA–phase2, presently running, with a lower software energy threshold [44] can also offer the possibility to increase sensitivity to such an effect.

4. Example of Reachable Sensitivities for ZnWO${}_4$ Anisotropic Detectors to Measure DM Directionality

The light output and the pulse shape for heavy particles (p, $\alpha$, nuclear recoils) in anisotropic scintillators depend on the direction with respect to the crystal axes. On the contrary, the response to $\gamma /\beta$ radiation is isotropic. This anisotropy property offers the possibility to investigate DM by the directionality approach, which is applicable to those DM candidate particles inducing nuclear recoils. The ZnWO${}_4$ has unique features among the anisotropic scintillators; these make it an excellent candidate for this kind of research. Moreover, plenty of room for improvement of its performance is still present. The possibility of a low background pioneer experiment ADAMO (Anisotropic detectors for DArk Matter Observation) to exploit the directionality approach by using anisotropic ZnWO${}_4$ scintillators deep underground was already discussed in [48].

The directionality approach [49] studies the correlation between the arrival direction of the DM particles, through the induced nuclear recoils, and the Earth motion in the galactic rest frame. In fact, the Earth experiences a wind of DM particles apparently flowing along a direction opposite to that of the solar motion relative to the DM halo. However, because of the Earth’s rotation around its axis, their average direction with respect to an observer fixed on the Earth changes during the sidereal day, as schematically shown in Figure 5. The nuclear recoils produced after a DM scattering are expected to be strongly correlated with the impinging direction of the DM particle. Hence, the study of the nuclear recoils direction offers a powerful tool for investigating the DM in a way largely independent on the assumptions. Unfortunately, the range of recoiling nuclei is very short (of the order of a few mm in low-pressure time projection chambers and typically of the order of m in solid detectors) for practical use in possible experiments aiming at measuring recoil tracks; however, these limitations can be overcome by using anisotropic scintillation detectors.

For an anisotropic detector with the three crystallographic axes fixed in the Earth reference frame, the directions of the possible nuclear recoils induced by the DM particles are expected to be strongly correlated with the sidereal time within the day, while the background events are not. Thanks to the different light responses of ZnWO${}_4$ anisotropic scintillators for nuclear recoils, we can expect that the energy distribution of the DM induced recoils depends on the direction of ${\overrightarrow{v}}_{\mathrm{lab}}\left(t\right)$, the detector’s velocity in the Galaxy, defined in Equation (1). Thus, considering a horizontal coordinate frame located at the North pole, described by the “polar-zenith”, ${\theta }_z\left(t\right)$, and by the “polar–azimuth”, ${\varphi }_a\left(t\right)$, angles, the counting rate is a function of the polar and azimuth angles; an example of this—whose details of the calculations are reported in [48]—is shown in Figure 6.

Therefore, the expected signal counting rate is a function of time, according to a time dependence as that reported in Equation (6), being the time period of $f\left(t\right)$ one sidereal day. The sensitivity of the method can be approached as in the previous Sections.

In particular, the sensitivity reachable by the ADAMO project (200 kg of ZnWO${}_4$), in 5 years of data taking with a duty cycle of 70%, average efficiency of 70%, in 1 keV energy bin and for a typical background level of 0.1 cpd/kg/keV can be calculated by Equation (18). The standard deviation $\sigma \left(\mathcal{B}\right)$ of the amplitude of the variation (here $\mathcal{B}$) due to the effect searched for is $\sigma \left(\mathcal{B}\right)\simeq 1.1\times {10}^{-3}$ cpd/kg/keV. The typical amplitudes of the variation, $\mathcal{B}$, are reported in Figure 6, where it is $\simeq 0.02\frac{\xi {\sigma }_p}{5\times {10}^{-5}\mathrm{pb}}$ for 2 keV energy threshold, and ≃$0.003\frac{\xi {\sigma }_p}{5\times {10}^{-5}\mathrm{pb}}$ for 6 keV energy threshold where ${\sigma }_p$ is the spin-independent elastic scattering cross section of the DM particle on nucleons and $\xi$ is the fraction of the DM local density of the considered candidate. Hence, $\xi {\sigma }_p\simeq 2.8\times {10}^{-6}$ pb and $\xi {\sigma }_p\simeq 1.8\times {10}^{-5}$ pb can be obtained at $1\sigma$ level. These values can be compared with those reported in Figure 7, where curves of sensitivity are reported for different background levels and for a software energy thresholds of 2 keVee. In the same figure the allowed regions obtained in [45] for the DAMA results after the first six annual cycles of DAMA/LIBRA–phase2 are also reported for comparison.

The anisotropic scintillators were suggested for the first time as detectors to study the directionality signature in the DM investigation in [51]. At a later time, the anthracene scintillators were considered as a good possibility [50]. However, only in 2013 the ZnWO${}_4$ anisotropic scintillators were proposed [48] (Some preliminary activities were also carried out by various authors; for details see [48,50,51] and references therein). In fact, the ZnWO${}_4$ is a very promising crystal for the DM directionality approach offering suitable features [48,52,53,54,55,56,57,58,59,60,61].

After those papers, thanks to the promising features of ZnWO${}_4$ anisotropic scintillators, the collaboration between the DAMA group of INFN and the INR-Kyiv group [53,54,62,63] has been developing several ZnWO${}_4$ detectors in the last years. The first crystals were grown and produced by the Institute for Scintillation Materials (ISMA, Kharkiv, Ukraine); later on, the DAMA/INR-Kyiv groups have been profiting from a collaboration with the Nikolaev Institute of Inorganic Chemistry (Novosibirsk, Russia) to carry on an R&D for producing ultra-radiopure ZnWO${}_4$. This R&D is based on the use of the low-thermal gradient Czochralski technique in a platinum crucible; it is still ongoing [64]. The crystals produced within these collaborations and R&D were measured in the DAMA/R&D underground facility at LNGS [53,54,62,63]. All these measurements and R&D show that the ZnWO${}_4$ scintillators are competitive detectors for experiments based on the approach of the DM directionality. In particular, it has been confirmed that the light output and the time profile of the scintillation pulses for heavy particles (p, $\alpha$, nuclear recoils) are dependent on the direction of such particles relative to the crystal axes. On the contrary, such dependence was not observed for $\gamma /\beta$ radiation [52,61]. In addition, the ZnWO${}_4$ anisotropic scintillators have shown a pulse shape discrimination capability between $\gamma /\beta$ radiation and heavy particles ($\alpha$). The R&D works have also shown the possibility to produce single ZnWO${}_4$ crystals with masses of some kg [64]. ZnWO${}_4$ have a high atomic weight and, due to the presence of three target nuclei with very different masses (Zn, W and O), these scintillators—as also the NaI(Tl)—are sensitive to both small and large mass DM candidates.

The ZnWO${}_4$ scintillators, developed in the R&D works, have shown a very good level of radio-purity; the measured upper limits are: 20 Bq/kg for the ${}^{40}$K, in the range (0.17–1.3) Bq/kg for the ${}^{228}$Th, and 2 Bq/kg for the ${}^{226}$Ra [53,56]; further radio-purification of ZnWO${}_4$ crystal scintillators is still feasible and one of the issues of the present R&D works. The measurements at LNGS confirmed that the light outputs of the crystals are relatively high at room temperature: ≃20% of the Na(Tl) scintillator and it can be further improved at low-temperature [65]. To study this feature, a small cryostat is currently under test at LNGS; it will allow a stable working temperature around −50 ${}^{\circ }$C.

The anisotropy of the ZnWO${}_4$ to nuclear recoils is a crucial issue for DM investigation with the directionality approach. For this purpose, specific measurements were performed by using a 7.99 g mass ZnWO${}_4$ crystal scintillator of $(10\times 10\times 10.4)$ mm${}^3$, in the framework of the ADAMO project [48,61]. The size of the crystal was chosen so small to avoid multiple scatterings of neutrons. The crystal was obtained by a second crystallization procedure using low-thermal gradient Czochralski technique from zinc tungstate crystals made from tungsten oxide additionally purified by double sublimation of tungsten chlorides [61]. The crystallographic axes, identified by the producer, were experimentally verified. The crystal was irradiated with monochromatic neutrons at ENEA-CASACCIA and with $\alpha$ particles at LNGS [61]. In the measurements with neutrons, the scattered neutrons were tagged by suitable neutron detectors; in such a way, the quenching factors of nuclear recoils were determined for three different neutron scattering angles (i.e., nuclear recoils energies) along the different crystallographic axes. Hence, the anisotropy of the light response for nuclear recoils in the ZnWO${}_4$ crystal scintillator was determined at 5.4 standard deviations [61].

5. Conclusions

In this paper, the annual modulation, the diurnal modulation and the directionality technique, as main signatures that can be exploited in order to investigate the presence of dark matter particles in the Galactic halo, have been discussed and presently achieved sensitivities have been summarized. Achieved experimental results obtained by DAMA, exploiting such signatures have been also reported. Further improvements in sensitivity to dark matter annual and diurnal modulation signatures are expected with the presently running DAMA/LIBRA–phase2 experiment; moreover lowering the software energy threshold below 1 keV can also offer the possibility to increase sensitivity to such effects. Efforts of the DAMA collaboration towards this direction are in progress. The perspectives of anisotropic solid scintillators, studied in the framework of the ADAMO project, have also been outlined.