Halo-independent determination of the unmodulated WIMP signal in DAMA: the isotropic case

TL;DR

This paper introduces a halo-independent method to determine the unmodulated dark matter signal in DAMA data, using a profile likelihood approach that accounts for uncertainties in the WIMP velocity distribution, specifically assuming isotropy.

Contribution

It develops a novel, mathematically rigorous framework for extracting the unmodulated WIMP signal from DAMA data without relying on specific halo models, assuming only isotropy in the velocity distribution.

Findings

Estimated unmodulated signal aligns with WIMP modulation expectations.

The modulated amplitude is 5-25% of the unmodulated signal depending on WIMP mass.

Method is applicable to any isotropic velocity distribution in the Galactic frame.

Abstract

We present a halo-independent determination of the unmodulated signal corresponding to the DAMA modulation if interpreted as due to dark matter weakly interacting massive particles (WIMPs). First we show how a modulated signal gives information on the WIMP velocity distribution function in the Galactic rest frame, from which the unmodulated signal descends. Then we perform a mathematically-sound profile likelihood analysis in which we profile the likelihood over a continuum of nuisance parameters (namely, the WIMP velocity distribution). As a first application of the method, which is very general and valid for any class of velocity distributions, we restrict the analysis to velocity distributions that are isotropic in the Galactic frame. In this way we obtain halo-independent maximum-likelihood estimates and confidence intervals for the DAMA unmodulated signal. We find that the estimated unmodulated signal is in line with expectations for a WIMP-induced modulation and is compatible with the DAMA background+signal rate. Specifically, for the isotropic case we find that the modulated amplitude ranges between a few percent and about 25% of the unmodulated amplitude, depending on the WIMP mass.

Tables

Table 1. Table 1 : DAMA modulation amplitudes S m , i subscript 𝑆 𝑚 𝑖 S_{m,i} , estimated unmodulated rates S 0 , i subscript 𝑆 0 𝑖 S_{0,i} , and background+signal rates B i + S i subscript 𝐵 𝑖 subscript 𝑆 𝑖 B_{i}+S_{i} in counts/kg/day/keV. Column 1: electron-equivalent energy bins. Column 2: DAMA modulation amplitudes S m , i subscript 𝑆 𝑚 𝑖 S_{m,i} from [ 2 ] . Columns 3–5: estimated unmodulated spectrum S 0 , i subscript 𝑆 0 𝑖 S_{0,i} for WIMP masses m χ = 5 ​ GeV subscript 𝑚 𝜒 5 GeV m_{\chi}=5~{}{\rm GeV} , m χ = 10 ​ GeV subscript 𝑚 𝜒 10 GeV m_{\chi}=10~{}{\rm GeV} , and m χ = 15 ​ GeV subscript 𝑚 𝜒 15 GeV m_{\chi}=15~{}{\rm GeV} . Column 6: DAMA background+signal rates B i + S i subscript 𝐵 𝑖 subscript 𝑆 𝑖 B_{i}+S_{i} from [ 67 ] (rebinned from 0.25-keVee- to 0.5-keVee-width bins).

$E_i$	$S_{m,i}$	$S_{0,i}$	$S_{0,i}$	$S_{0,i}$	$B_i+S_i$
[keVee]		$m_{\chi }=5$ GeV	$m_{\chi }=10$ GeV	$m_{\chi }=15$ GeV
2.0–2.5	$0.0161\left(39\right)$	${0.513}_{-0.050}^{+0.051}$	${0.333}_{-0.045}^{+0.046}$	${0.227}_{-0.027}^{+0.35}$	$1.029$
2.5–3.0	$0.0260\left(44\right)$	${0.311}_{-0.030}^{+0.031}$	${0.239}_{-0.037}^{+0.039}$	${0.165}_{-0.024}^{+0.33}$	$1.228$
3.0–3.5	$0.0219\left(44\right)$	${0.165}_{-0.016}^{+0.016}$	${0.160}_{-0.030}^{+0.032}$	${0.112}_{-0.020}^{+0.30}$	$1.294$
3.5–4.0	$0.0084\left(40\right)$	${0.0791}_{-0.0078}^{+0.0078}$	${0.103}_{-0.025}^{+0.027}$	${0.073}_{-0.017}^{+0.27}$	$1.140$
4.0–4.5	$0.0107\left(36\right)$	${0.0352}_{-0.0036}^{+0.0035}$	${0.066}_{-0.021}^{+0.022}$	${0.047}_{-0.014}^{+0.24}$	$0.956$
4.5–5.0	$0.0054\left(32\right)$	${0.0148}_{-0.0016}^{+0.0015}$	${0.042}_{-0.016}^{+0.019}$	${0.030}_{-0.011}^{+0.22}$	$0.853$
5.0–5.5	$0.0089\left(32\right)$	${0.00600}_{-0.0065}^{+0.0059}$	${0.026}_{-0.012}^{+0.016}$	${0.018}_{-0.0084}^{+0.20}$	$0.868$
5.5–6.0	$0.0039\left(31\right)$	${0.00236}_{-0.00026}^{+0.00023}$	${0.0156}_{-0.0084}^{+0.013}$	${0.011}_{-0.0059}^{+0.18}$	$0.853$
6.0–6.5	$0.00018\left(308\right)$	${9.04}_{-1.03}^{+0.89}\times {10}^{-4}$	${0.0090}_{-0.0054}^{+0.011}$	${0.0064}_{-0.0038}^{+0.16}$	$0.868$
6.5–7.0	$0.00018\left(281\right)$	${3.41}_{-0.40}^{+0.34}\times {10}^{-4}$	${0.0050}_{-0.0033}^{+0.0094}$	${0.0035}_{-0.0023}^{+0.15}$	$0.860$
7.0–7.5	$0.0015\left(28\right)$	${1.27}_{-0.15}^{+0.12}\times {10}^{-4}$	${0.0026}_{-0.0018}^{+0.0073}$	${0.0019}_{-0.0013}^{+0.13}$	$0.860$
7.5–8.0	$-0.0013\left(29\right)$	${4.67}_{-0.56}^{+0.46}\times {10}^{-5}$	${0.0013}_{-0.0010}^{+0.0054}$	${0.95}_{-0.71}^{+123}\times {10}^{-3}$	$0.890$

Equations

$$S(t)=S_0+S_m\cos [\omega (t-t_0)],$$

$$\tilde{\eta }(v_{\min })=\frac{{\rho }_{\chi }}{m_{\chi }}{\sigma }_{\chi N}{\int }_{|\mathbf{v}|>v_{\min }}\frac{f_{\mathrm{lab}}(\mathbf{v})}{|\mathbf{v}|}{d}^{3}v.$$

$$\int f_{\mathrm{lab}}(\mathbf{v},t)d^{3}v=1.$$

$$S_i(t)=\int {\mathcal{H}}_i(\mathbf{v})f_{\mathrm{lab}}(\mathbf{v},t)d^{3}v,$$

$$\frac{d{R}_{\chi T}}{d{E}_R}=\frac{1}{m_T}\int \frac{d{\sigma }_{\chi T}}{d{E}_R}\frac{{\rho }_{\chi }}{m_{\chi }}v{f}_{\mathrm{lab}}(\mathbf{v},t)d^{3}v.$$

$$\frac{dR}{dE}=\sum _T{C}_T\int {\mathcal{P}}_T(E,E_R)\frac{d{R}_{\chi T}}{d{E}_R}d{E}_R.$$

$${\mathcal{H}}_E(v)=v\frac{{\rho }_{\chi }}{m_{\chi }}\int d{E}_R\sum _T{\mathcal{P}}_T(E,E_R)\frac{C_T}{m_T}\frac{d{\sigma }_{\chi T}}{d{E}_R}.$$

$$\displaystyle f_{\rm lab}({\bf v},t)=f_{\rm gal}\big{(}{\bf u}\big{)},\qquad{\bf u}={\bf v}+{\bf v}_{\odot}+{\bf v}_{\oplus}(t).$$

$$S_i(t)=\int {\mathcal{H}}_i^{\mathrm{gal}}(\mathbf{u},t)f_{\mathrm{gal}}(\mathbf{u})d^{3}u,$$

$$\displaystyle{\cal H}^{{\rm gal}}_{i}({\bf u},t)={\cal H}_{i}\big{(}{\bf u}-{\bf v}_{\odot}-{\bf v}_{\oplus}(t)\big{)}.$$

$$S_{m,i}=\frac{2}{T}{\int }_0^{T}dt\cos [\omega (t-t_0)]S_i(t).$$

$$S_{0,i}=\frac{1}{T}{\int }_0^{T}dt{S}_i(t).$$

$$S_{0,i}$$

$$S_{m,i}$$

$${\mathcal{H}}_{0,i}^{\mathrm{gal}}(\mathbf{u})$$

$${\mathcal{H}}_{m,i}^{\mathrm{gal}}(\mathbf{u})$$

$$-2\ln \mathcal{L}(\{S_{m,j}\})=\sum _{j=1}^N{\left(\frac{S_{m,j}-S_{m,j}^{\exp }}{\Delta S_j^{\exp }}\right)}^2.$$

$$S_{m,j}=\int {\mathcal{H}}_{m,j}^{\mathrm{gal}}(\mathbf{u})f_{\mathrm{gal}}(\mathbf{u})d^{3}u,$$

$$S_{0,i}=\int {\mathcal{H}}_{0,i}^{\mathrm{gal}}(\mathbf{u})f_{\mathrm{gal}}(\mathbf{u})d^{3}u,$$

$${\mathcal{L}}_{\mathrm{p}}(\{S_{0,i}\})=\underset{f_{\mathrm{gal}}\in \mathcal{A}(\{S_{0,i}\})}{\sup }\mathcal{L}(\{S_{m,j}\}).$$

$$-2\Delta \ln {\mathcal{L}}_{\mathrm{p}}(\{S_{0,i}\})=1,$$

$$-2\Delta \ln {\mathcal{L}}_{\mathrm{p}}(\{S_{0,i}\})=-2\ln {\mathcal{L}}_{\mathrm{p}}(\{S_{0,i}\})+2\ln {\mathcal{L}}_{\mathrm{p},\max }.$$

$${\mathcal{L}}_i(S_{0,i})=\underset{\begin{array}{c}{S}_{0,j}\\ j\ne i\end{array}}{\max }{\mathcal{L}}_{\mathrm{p}}(\{S_{0,j}\}),$$

$$-2\Delta \ln {\mathcal{L}}_i(S_{0,i})\le 1,$$

$$-2\Delta \ln {\mathcal{L}}_i(S_{0,i})=-2\ln {\mathcal{L}}_i(S_{0,i})+2\ln {\mathcal{L}}_{\mathrm{p},\max }.$$

$$S_{0,i}^{\inf }(L_0)=\underset{f_{\mathrm{gal}}\in \mathcal{A}(L_0)}{\inf }\int {\mathcal{H}}_{0,i}^{\mathrm{gal}}(\mathbf{u})f_{\mathrm{gal}}(\mathbf{u})d^{3}u$$

$$S_{0,i}^{\sup }(L_0)=\underset{f_{\mathrm{gal}}\in \mathcal{A}(L_0)}{\sup }\int {\mathcal{H}}_{0,i}^{\mathrm{gal}}(\mathbf{u})f_{\mathrm{gal}}(\mathbf{u})d^{3}u,$$

$$\mathcal{L}(\{S_{m,j}\})\ge L_0.$$

$$\text{Extremize }{S}_{0,i}=\int {\mathcal{H}}_{0,i}^{\mathrm{gal}}(\mathbf{u})f_{\mathrm{gal}}(\mathbf{u})d^{3}u$$

$$\int f_{\mathrm{gal}}(\mathbf{u})d^{3}u=1,$$

$$\int {\mathcal{H}}_{m,j}^{\mathrm{gal}}(\mathbf{u})f_{\mathrm{gal}}(\mathbf{u})d^{3}u=S_{m,j}(j=1,\dots ,N),$$

$$\mathcal{L}(\{S_{m,j}\})\ge L_0.$$

$$f_{\mathrm{e}}(\mathbf{u})=\sum _{k=1}^K{\lambda }_k\delta (\mathbf{u}-{\mathbf{u}}_k),$$

Full text

1 Introduction

Discovering the nature of Dark Matter (DM) is one of the most important endeavors in today’s particle physics and cosmology. Weakly Interacting Massive Particles (WIMPs), the most popular and natural DM candidates, provide the correct thermal relic density in the early Universe when their cross section with Standard Model particles is at the level of or not much smaller than weak interaction cross sections. Many experiments are currently trying to exploit this fact to detect the WIMPs supposed to form the dark halo of our Galaxy through their scattering off atomic nuclei in low-background laboratory detectors. An expected feature of halo WIMP scattering is an annual modulation of the scattering rate due to the revolution of the Earth around the Sun [1]. The importance of such yearly modulation for WIMP direct detection rests on the fact that, in absence of a sensitivity to the direction of the incoming particles, it is the only known signature that allows to distinguish a WIMP signal from the background due to radioactive contamination, since the latter may have an energy spectrum indistinguishable from that predicted for WIMPs.

An annual modulation with WIMP characteristics has been claimed for many years by the DAMA experiment [2]. The low–energy event rate in the DAMA sodium iodide scintillators is well represented by a signal of the form

[TABLE]

with $\omega=2\pi/T$, $T=1~{}{\rm yr}$, and $t_{0}\simeq$ 2nd of June, as expected for a nonrotating dark halo of WIMPs. The statistical significance of the DAMA modulation signal exceeds 9 standard deviations and the effect has been recorded through 14 yearly periods. For a typical Maxwellian distribution of WIMP velocities in the Galactic rest frame with rms velocity below 300 km/s, the modulated component $S_{m}$ of the signal is predicted to be less that 10% of the unmodulated component $S_{0}$, in all of the energy bins of the DAMA detected spectrum.

The DAMA claim has prompted a world–wide effort by other direct detection experiments to confirm or disprove the signal [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The experiments that have by now reached a background level low enough to be sensitive to the DAMA modulation use targets different from sodium iodide. As a consequence, while for standard hypotheses on the WIMP–nucleon cross section and WIMP velocity distribution (i.e., spin–dependent or spin–independent interactions with a truncated Maxwellian distribution) the DAMA signal appears to be in strong tension with constraints from other detectors, when such assumptions are relaxed several WIMP models have been shown to exist for which the yearly modulation effect measured by DAMA can still be reconciled with the non observation of a dark matter signal in other experiments [16, 17, 18, 19]. This shows that there is still a clear need to assess the compatibility of the DAMA excess with other detectors in a model–independent way.

Eliminating the dependence on astrophysics is the underlying goal of the halo–independent approach. Its first formulation [20] was based on the observation that the elastic spin-independent scattering rate of WIMPs in a detector depends on the velocity distribution only through a single velocity integral $\tilde{\eta}(v_{\rm min})$, the same for all experiments,

[TABLE]

Here $m_{\chi}$ is the WIMP mass, $\sigma_{\chi N}$ is the WIMP-nucleon cross section, $\rho_{\chi}$ is the local WIMP mass density and $f_{{\rm lab}}({\bf v})$ is the WIMP velocity distribution in the frame of the detector (the laboratory). The method of [20] has been applied to the comparison of experiments in [21, 22, 20, 23, 24, 25, 26, 27, 28, 29]. It has been generalized to arbitrary WIMP–nucleon interactions and any direct detection experiment (with arbitrary efficiency and energy resolution) in [30] by defining weighted averages of $\tilde{\eta}(v_{\rm min})$ over the range(s) of velocities measured in an experiment. Applications of the latter method to the comparison of experiments for various WIMP–nucleon interactions can be found in [30, 31, 32, 33, 34, 16, 35, 17, 18, 36, 37, 38]. Maximum-likelihood methods to determine the velocity integral and particle physics parameters have been used in [39, 40, 27, 41, 29, 36, 37, 42, 38] and statistical methods to assess the compatibility of experiments have been considered in [41, 43, 37]. Alternative methods to place halo–independent bounds on particle physics parameters have been put forward in [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56].

A weak point of the halo-independent methods above is the way they compare modulated and unmodulated rates. Some authors define two separate velocity integrals $\tilde{\eta}_{0}(v_{\rm min})$ and $\tilde{\eta}_{m}(v_{\rm min})$, one for the unmodulated part $S_{0}$ and one for the modulated part $S_{m}$ in Eq. (1.1), and then proceed to impose either the simple inequality $\tilde{\eta}_{m}(v_{\rm min})<\tilde{\eta}_{0}(v_{\rm min})$ [30, 31, 32, 33, 34, 35, 17, 18, 36, 37] or more sophisticated inequalities valid for smooth distributions [24, 25]. Comparing two separate velocity integrals with proper statistical significance is not straightforward. Other authors replace the modulation amplitude $S_{m}$, which is a coefficient in a Fourier time-series, with half the difference between the maximum and minimum signal during a year, i.e., replace $\tilde{\eta}_{m}(v_{\rm min})$ with $\tilde{\eta}_{1/2}(v_{\rm min})=[\tilde{\eta}_{0}(t_{0})+\tilde{\eta}_{0}(t_{0}+\tfrac{1}{2}T)]/2$ [23, 26, 16]. This replacement is inaccurate in a halo-independent approach, where one must include velocity distributions for which the modulation is not sinusoidal near the threshold region in which the DAMA signal is present (in this case, the theoretical values of $\tilde{\eta}_{1/2}$ and $S_{m}$ may be very different, and without a control on the sinusoidal character of the modulation, it would be inappropriate to compare the theoretical $\tilde{\eta}_{1/2}$ with the measured $S_{m}$).

The main goal of this paper is to show that it is possible to obtain information on the unmodulated signal $S_{0}$ from a measurement of the modulation amplitudes $S_{m}$ without specifying the WIMP velocity distribution (and without assuming two separate velocity integrals or approximating the Fourier coefficient with a difference). For this purpose, (a) we transfer the modulation from being a property of the velocity distribution to being a property of the detector, indeed of the relative motion of the detector with respect to the rest frame of the WIMP population, and (b) we profile the likelihood at fixed $S_{0}$ over all WIMP velocity distributions (a continuum of nuisance parameters) using rigorous mathematical methods based on the theory of linear optimization in spaces of functions that have the dimension of the continuum.

As a first application of the method, which is very general and valid for any elastic or inelastic cross section and any class of velocity distributions, we estimate the DAMA unmodulated signal starting from data on the modulation amplitudes under some simplifying assumptions that have allowed us to explore and understand the difficulties and merits of the method. We restrict the analysis to velocity distributions that are isotropic in the Galactic frame. We apply our method to the DAMA data for the case of WIMP–nucleus elastic scattering and for WIMP masses $m_{\chi}<$ 15 GeV, for which only sodium targets contribute to the expected signal.111At higher WIMP masses the expected rate in DAMA also takes contributions from scattering off iodine. This part of the signal is in principle constrained in a model–independent way by KIMS [6], which uses a CsI detector, and COUPP [11], which uses a CF3I detector, since such detectors employ the same nuclear target (iodine). In this way we will obtain quantitative confidence intervals for the unmodulated components of the WIMP signal in each energy bin, for the first time disentangling the unmodulated signal from the background in a halo–independent way. The $S_{0}$ confidence intervals we find are valid for any WIMP–nucleus interaction in which the 23Na cross section varies negligibly in the 2–4 keVee energy range where the DAMA modulation is present. This includes the standard elastic spin–independent and spin–dependent interactions, with arbitrary ratios of the proton–WIMP and neutron–WIMP coupling constants, but does not include inelastic scattering or effective WIMP–nucleon operators that show an explicit and fast dependence on the WIMP–nucleus relative velocity and/or on the exchanged momentum (like some of those in [57]). Since the properties discussed in the present paper pertain exclusively to sodium targets in DAMA, while other detectors that constrain DAMA use different target nuclei, we will not discuss the latter any further.

The plan of this paper is as follows. In Section 2 we show how to compare the modulated and the unmodulated rates directly in terms of a single velocity distribution function, namely the velocity distribution function in the Galactic rest frame. In Section 3 we present and discuss our method to compute the profile likelihood of the unmodulated signal by using linear optimization theory in the continuum to profile the likelihood over the whole velocity distribution (a continuum of nuisance parameters). Section 4 is devoted to our numerical analysis of the DAMA data for velocity distributions that are isotropic in the Galactic frame. Figs. 8, 9 and 10 and Table 1 of this Section contain our main quantitative results. Finally, the Appendix contains details of the calculation of the modulated and unmodulated Galactic response functions for isotropic velocity distributions.

2 Rates in terms of the Galactic velocity distribution

Let $S_{i}(t)$ be the expected signal in a dark matter detector, where $t$ is time, and the index $i$, which may be continuous, specifies the quantity measured in the experiment, for example detected recoil energy or energy bin or number of photoelectrons. Let $f_{\rm lab}({\bf v},t)$ be the WIMP velocity distribution function in the frame of the detector, normalized to

[TABLE]

The signal $S_{i}(t)$ depends on $f_{\rm lab}({\bf v},t)$ according to the general formula [58]:

[TABLE]

where ${\cal H}_{i}({\bf v})$, called the response function, equals the value the signal would have if all the WIMPs had the same velocity ${\bf v}$. Formula (2.2) can be understood for example by writing the scattering rate per unit target mass $dR_{\chi T}/dE_{R}$ of a WIMP $\chi$ off an isotope $T$ in the target, differential in the nucleus recoil energy $E_{R}$, as a product of the differential cross section $d\sigma_{\chi T}/dE_{R}$ and the WIMP flux $n_{\chi}vf_{\rm lab}({\bf v},t)d^{3}v$ (where $n_{\chi}=\rho_{\chi}/m_{\chi}$ is the $\chi$ number density), the whole quantity divided by the mass $m_{T}$ of the target isotope,

[TABLE]

Here it is understood that the differential cross section $d\sigma_{\chi T}/dE_{R}$ is nonzero only in the kinematically allowed region (e.g., for elastic scattering, only for $E_{R}\leq E_{R}^{\rm max}(v)$ given in equation (4.13) below). Furthermore, if ${\cal P}_{T}(E,E_{R})$ indicates the probability of actually observing an event with observed energy $E$ when a WIMP has scattered off an isotope $T$ in the detector target with recoil energy $E_{R}$, the expected observed event rate per unit target mass $dR/dE$ is given by the convolution

[TABLE]

Here $C_{T}$ is the mass fraction of isotope $T$ in the target. Inserting the scattering rate in equation (2.3) into the latter expression, and exchanging the order of the integrations over $E_{R}$ and ${\bf v}$, leads to the following expression of the response function ${\cal H}_{E}(v)$, where $i=E$, for the differential event rate $dR/dE$,

[TABLE]

The response function depends on the particle physics model for the interaction of the WIMP with the target and includes the probability that a WIMP scattering in the detector is actually observed. The nonzero values of the response function ${\cal H}_{i}({\bf v})$ also indicate the WIMP velocities ${\bf v}$ to which the observed signal $S_{i}(t)$ is sensitive. General expressions of the response functions ${\cal H}_{i}({\bf v})$ for experiments counting number of events in observed energy bins can be found in [58].

One commonly assumes that the response function ${\cal H}_{i}({\bf v})$ is stable, i.e., that it does not depend on time (as already implied in the notation above). One also commonly assumes that the only time dependence in the laboratory velocity distribution $f_{\rm lab}({\bf v},t)$ comes from the motion of the Earth around the Sun or the daily rotation of the Earth. In other words, one assumes that the WIMP velocity distribution in the Galactic frame $f_{\rm gal}({\bf u})$, where ${\bf u}$ is the Galactic WIMP velocity, is stationary (on the time scale of the experiment).

The laboratory velocity distribution $f_{\rm lab}({\bf v},t)$ is related to the Galactic velocity distribution $f_{\rm gal}({\bf u})$ by a Galilean transformation:

[TABLE]

Here ${\bf v}$ is the WIMP velocity relative to the laboratory, ${\bf u}$ is the WIMP velocity relative to the rest frame of the Galaxy, ${\bf v}_{\odot}$ is the velocity of the Sun with respect to the Galactic rest frame, and ${\bf v}_{\oplus}(t)$ is the velocity of the Earth with respect to the Sun.222Since we are interested in the annual modulation, we neglect the daily rotation of the Earth, but our general considerations apply using the velocity of the detector with respect to the Sun in place of ${\bf v}_{\oplus}(t)$. A change of integration variables in Eq. (2.2) gives

[TABLE]

where

[TABLE]

Conceptually, this passage to the Galactic frame means that the time dependence of the signal is a property of the motion of the detector in the Galaxy and not of the (Galactic) distribution function. In particular, the characteristics of a modulated signal are ultimately a property of the detector (composition, energy threshold, motion in the Galaxy) and not of the WIMP velocity distribution in the Galactic frame.

The DAMA modulation amplitude is equal to the coefficient of the $\cos[\omega(t-t_{0})]$ term in the Fourier time-series analysis of the signal. Here $\omega=2\pi/T$ with the period $T=1$ yr, and the time $t=t_{0}$ corresponds to the time of maximum modulation. So:

[TABLE]

The unmodulated signal is the time average of the signal over the course of a year,

[TABLE]

In the Galactic frame, the time dependence is in the response functions, so one can write

[TABLE]

where

[TABLE]

3 Constraining the unmodulated signal

Given $N$ experimental modulation amplitude measurements $S_{m,j}^{\rm exp}$ ($j=1,\ldots,N$) and their 1$\sigma$ errors $\Delta S_{j}^{\rm exp}$, and assuming Gaussian fluctuations, the likelihood function ${\cal L}$, function of the set of parameters $\{S_{m,j}\}$, is given by

[TABLE]

Here, for clarity and simplicity, we only show the dependence of ${\cal L}$ on the expected modulation amplitudes

[TABLE]

which contain all the dependence on the WIMP velocity distribution function $f_{\rm gal}({\bf u})$.

We are interested in constraining the unmodulated signals

[TABLE]

in the given energy bins $i=1,\ldots,N$. For this purpose, we construct a joint profile likelihood ${\cal L}_{\rm p}(\{S_{0,i}\})$ for the $S_{0,i}$ ($i=1,\ldots,N$), treating the velocity distribution $f_{\rm gal}({\bf u})$ as a continuum of nuisance parameters. The joint profile likelihood ${\cal L}_{\rm p}(\{S_{0,i}\})$ is defined as the maximum value of the likelihood function over the set ${\cal A}(\{S_{0,i}\})$ of distribution functions that satisfy Eq. (3.3),

[TABLE]

(Technically, we use the notation sup instead of max because in the infinite-dimensional space of distribution functions it is not automatically guaranteed that there is a distribution that achieves the maximum, although this does happen in our case.)

The maximum-likelihood estimator of the $S_{0,i}$’s then follows as the location of the maximum ${\cal L}_{\rm p,max}$ in the joint profile likelihood ${\cal L}_{\rm p}(\{S_{0,i}\})$, and confidence regions for any combination of the $S_{0,i}$’s can be obtained through the usual profile likelihood procedure. For example, the standard error ellipsoid in the $N$-dimensional $S_{0,i}$ parameter space (or ‘pseudo-ellipsoid’ in the case of non-Gaussian likelihoods) can be obtained by the condition

[TABLE]

where

[TABLE]

Similarly, projections of the standard error ellipsoid onto any pair of two variables $S_{0,i}$ and $S_{0,j}$ can be obtained either by simple geometric construction or by using a profile likelihood further profiled over the other $S_{0,i}$ parameters. In Section 4 we show an example of such 2-dimensional standard error ellipses.

We are interested in particular in finding confidence intervals on each of the $N$ unmodulated signals $S_{0,i}$. For this purpose, we use the profile likelihood function

[TABLE]

which is the joint profile likelihood ${\cal L}_{\rm p}(\{S_{0,i}\})$ further profiled over the $S_{0,j}$ with $j\neq i$, and is thus a function of $S_{0,i}$ only. A 1$\sigma$ confidence interval on a single $S_{0,i}$ can then be obtained through the condition

[TABLE]

where

[TABLE]

Notice that these 1$\sigma$ confidence intervals, sometime called 1$\sigma$ likelihood intervals, have a 68% coverage probability in the limit of large samples when the likelihood is well approximated by a Gaussian, but do not necessarily have a coverage probability of 68% if the likelihood is non-Gaussian.

We examined various ways of computing ${\cal L}_{\rm p}(\{S_{0,i}\})$ and ${\cal L}_{i}(S_{0,i})$, and we have adopted the following procedure. Since both $S_{0,i}$ and $S_{m,i}$ are functionals of the distribution function $f_{\rm gal}$, we can at least conceptually construct a parametric plot of ${\cal L}(\{S_{m,i}\})$ vs. $\{S_{0,i}\}$ by using $f_{\rm gal}$ as the parameter. At each point $\{S_{0,i}\}$ there will be many values of ${\cal L}(\{S_{m,i}\})$, and our goal is to find the maximum of those values, which is ${\cal L}_{\rm p}(\{S_{0,i}\})$. In other words, in this geometrical representation, the joint profile likelihood ${\cal L}_{\rm p}(\{S_{0,i}\})$ is the boundary of all possible values of the likelihood ${\cal L}(\{S_{m,i}\})$ when plotted vs. $\{S_{0,i}\}$. In practice we cannot implement an infinite number of functions $f_{\rm gal}$ (we tried discretizing the distribution function but the maximization procedure did not converge). However, we can think of constructing the boundary of the likelihood values by “rotating the plot by 90 degrees,” i.e., finding the boundary of the $\{S_{0,i}\}$ that have ${\cal L}(\{S_{m,i}\})\geq{\cal L}_{\rm p}(\{S_{0,i}\})$. The latter problem can be written as an extremization problem for $\{S_{0,i}\}$ for which powerful mathematical theorems exist that reduce the infinitely-many functions $f_{\rm gal}$ to a finite number of parameters, making the solution attainable in practice.

For clarity, we illustrate our procedure for ${\cal L}_{i}(S_{0,i})$ only, although we used it for the joint profile likelihood ${\cal L}_{\rm p}(\{S_{0,i}\})$. We write our problem as a linear optimization problem for $\{S_{0,i}\}$ over the distribution functions $f_{\rm gal}({\bf u})$ subject to the constraint that the likelihood function ${\cal L}(\{S_{m,j}\})$ is greater than or equal to a given number $L_{0}$, which we will later vary. Our goal is to find the lower and upper bounds

[TABLE]

over the set ${\cal A}(L_{0})$ of distribution functions that satisfy the constraint

[TABLE]

Varying $L_{0}$, we obtain the lines $S_{0,i}^{\rm inf}(L_{0})$ and $S_{0,i}^{\rm sup}(L_{0})$ in the $S_{0,i}$–${\cal L}$ plane, which we then invert to obtain the profile likelihood ${\cal L}_{i}(S_{0,i})$ as a function of $S_{0,i}$. Figure 1 is a schematic illustration of our procedure.

To compute $S_{0,i}^{\rm inf}(L_{0})$ and $S_{0,i}^{\rm sup}(L_{0})$ we notice that the likelihood function is a function of $f_{\rm gal}({\bf u})$ only through the integrals $S_{m,j}$ in Eq. (3.2). So we rephrase our goal as the following mathematical problem.

[TABLE]

Mathematically, this is an optimization problem of the kind discussed for example in [59, 60]. Specifically, an optimization problem in which the moment set (i.e., the set of distribution functions that obey the moment conditions) is defined by a subset $C(L_{0})\subseteq\mathbb{R}^{N+1}$, namely the subset $C(L_{0})$ in which the moments $(1,S_{m,1},\ldots,S_{m,N})$ satisfy ${\cal L}(\{S_{m,j}\})\geq L_{0}$.

The fundamental theorem in this context [60, 59] states that $S_{0,i}$ achieves its extreme values $S_{0,i}^{\rm inf}(L_{0})$ and $S_{0,i}^{\rm sup}(L_{0})$ on the extreme distributions of the moment set, which are sums with positive coefficients of a finite number $K\leq N+1$ of Dirac delta functions. Here $N+1$ is the number of moment conditions. More precisely, the extreme distributions of the moment set defined by the moment conditions (3.14–3.15) have the form

[TABLE]

where $1\leq K\leq N+1$,

[TABLE]

and the $K$ $N+1$-dimensional vectors $(1,{\cal H}^{{\rm gal}}_{m,1}({\bf u}_{k}),\ldots,{\cal H}^{{\rm gal}}_{m,N}({\bf u}_{k}))$, where the index $k=1,\ldots,K$ specifies the vector, are linearly independent.

Geometrically, the extreme distributions of the moment set are analogous to the vertices of a polyhedron. In the finite-dimensional case, the fundamental theorem states that the maximum and minimum values of a linear function defined over a polyhedron are achieved at one or more vertices of the polyhedron, and thus to find these extrema it suffices to compute the value of the linear function on the vertices. In the continuum case, the fundamental theorem states that the extrema of a linear functional of the distribution (in our case, each $S_{0,i}$) are achieved at one or more “vertices” of the moment set (i.e., at the extreme distributions), and thus to find these extrema it suffices to compute the value of the linear functional on the extreme distributions. The computational advantage is that the fundamental theorem reduces an extremization problem in infinite dimensions (the moment set) into an extremization problem in a finite number of dimensions (the space of extreme distributions, which has dimension at most $(1+d)N^{\prime}$, where $N^{\prime}=N+1$ is the number of moment conditions and $d$ is the dimensionality of the velocity space).

Physically, each delta-function distribution $\delta({\bf u}-{\bf u}_{k})$ in the expression of an extreme distribution, Eq. (3.17), represents a stream of velocity ${\bf u}_{k}$ and zero velocity dispersion. An extreme distribution is a weighted average of streams with weights $\lambda_{k}$. The fundamental theorem allows the computation of $S_{0,i}^{\rm inf}(L_{0})$ and $S_{0,i}^{\rm sup}(L_{0})$ by parametrizing the velocity distribution as a weighted average of no more than $N+1$ streams in velocity space, where $N+1$ is the number of moment conditions, including the normalization condition.333The fundamental theorem is also the rigorous mathematics behind the “interesting” facts that the number of steps in $\tilde{\eta}(v_{\rm min})$ is less than or equal to the number of bins for a binned likelihood [39], or that $N_{O}$ streams suffice to minimize an unbinned likelihood with $N_{O}$ events [40, 17], or that only $1+p+q$ grid points have a nonzero distribution function in the presence of $1+p+q$ moment conditions [56].

The fundamental theorem translates the mathematical problem (3.13–3.16) into the following one.

[TABLE]

In practice this means that at fixed $K$ and given $L_{0}$, the maximal range of the $S_{0,i}$ integral computed using Eq. (3.20) is swept by the $\lambda_{k}$, ${\bf u}_{k}$ parameters that satisfy the constraints (3.22-3.25). The full range of $S_{0,i}$ is then obtained by combining the $N+1$ intervals for $1\leq K\leq N+1$.

It is very important to understand that this method does not in general give the optimal velocity distribution, or the maximum-likelihood velocity distribution. In fact, given a value $L_{0}$ of the likelihood, there are in general many $S_{m,i}$ that have the same likelihood (all those on the likelihood contour level ${\cal L}(\{S_{m,i}\})=L_{0}$). But even if there is only one set of $S_{m,i}$ that corresponds to a given value of the likelihood (and this happens at the point of absolute maximum likelihood for concave likelihood functions), there are in general many velocity distributions with the same moments $S_{m,i}$: some of them are extreme distributions (sums of streams), and some are continuous distributions, or more precisely, continuous linear combinations of sums of streams of the form

[TABLE]

In particular, although the value of the maximum likelihood can be obtained using only sums of streams, there is in general an infinite number of distributions, some discrete and some continuous, that maximize the likelihood. So even if we use extreme distributions to find the extreme values of $S_{0,i}$, it is not correct to think that in general these sums of streams are the only velocity distributions giving those extreme values. The reason we can use the methods of this Section to estimate the unmodulated signal is that we are not interested in finding the optimal velocity distribution but in performing a maximum-profile-likelihood analysis of quantities like $S_{0,i}$ that are integrals of the velocity distribution. For this task, the method described in this Section is adequate and mathematically sound.

4 Analysis: isotropic case

In this Section, we apply the general method described in Section 3 to a specific case: a halo-independent estimate of the unmodulated DAMA signal. Since this is the first implementation of our method, we have made some simplifying assumptions that have allowed us to explore and understand the difficulties and merits of the method itself. First and foremost, to reduce the computing time, we have restricted our analysis to WIMP velocity distributions that are isotropic in the Galactic reference frame, i.e., $f_{\rm gal}({\bf u})=f_{\rm gal}(u)$ with $u=|{\bf u}|$, so that the distribution functions depend on one variable only (the magnitude $u$) instead of three (the components of ${\bf u}$). Under this assumption, the Galactic response functions can be replaced by the angle-averaged Galactic response functions, defined by an average over the directions of the vector ${\bf u}$ as

[TABLE]

Then our extremization problem reads

[TABLE]

Moreover, the isotropic extreme distribution functions are

[TABLE]

where $\overline{f}_{\rm e}(u)=4\pi u^{2}f_{\rm e}(u)$. While mathematically appropriate and well-defined, physically these isotropic extreme distributions do not describe a collection of streams in velocity space but rather some sort of spherical shells in velocity space.

An additional simplifying assumption we make in this first application of our method is to consider only spin-independent scattering off sodium in the DAMA NaI detector. Thus we restrict our analysis to light WIMP masses ($m_{\chi}\leq$ 15 GeV) for which WIMP elastic scattering off iodine is below threshold (for a constant iodine quenching factor $Q_{I}=0.09$ and Galactic escape speeds less than $\sim\!\!580$ km/s). It must be noted that the results of our analysis apply also to cross sections that are not spin-independent but have a mild energy dependence in the 2–4 keVee energy range.

For our analysis we use the $N=12$ DAMA cosine modulation amplitudes in the lowest energy bins in Fig. 8 of Ref. [2]. We list them in Table 1. These measurements were obtained using a total exposure of 1.33 ton yr. The signal is concentrated in the first 6 bins, and the other 6 bins act as a control set with no modulation signal. Data are also available for the sine modulation amplitudes [2] but under our simplifying assumption of isotropic velocity distribution, the sine modulation response functions vanish identically (see Appendix), and thus including them in the likelihood would amount to adding an irrelevant constant. The DAMA Collaboration published also time series of its modulation data [2], with time binnings ranging from   30 days (close to the maxima and minima of the oscillation) and  70 days (close to its equilibrum points). However it has been shown that the corresponding error bars can easily accommodate sizeable distortions from a sinusoidal time dependence of the signal (see for instance the discussion in Section 5 of Ref.[61], relative to the case of a Maxwellian distribution yielding modulation fractions of order unity when the incoming WIMP velocities are very close to the escape velocity), so the ensuing constraint has no impact on our analysis and we neglect it.

For the DAMA response functions off sodium, we take $i$ to be the index of the energy bin in the electron–equivalent energy $E_{ee}$. The latter is related to the recoil energy $E_{R}$ on average by $\overline{E}_{ee}=Q(E_{R})\,E_{R}$, where $Q(E_{R})$ is the quenching factor, with an additional smearing due to a finite energy resolution. We use the DAMA response functions for elastic spin-independent WIMP–nucleus scattering, but our results extend practically unchanged to other elastic interactions in which the velocity and/or energy dependence of the cross section is negligible in the DAMA energy bins (e.g., spin-dependent interactions or other nonrelativistic effective operators that do not show a large variation with recoil energy or WIMP velocity). The DAMA response function for the $i$-th bin with electron–equivalent energies in the range $E_{ee,i}\leq E_{ee}\leq E_{ee,i+1}$ is [58]:

[TABLE]

Here $N_{T}$ is the number of target nuclei in the detector, $M_{\rm det}$ is the mass of the detector, $\Delta E=E_{ee,i+1}-E_{ee,i}$ is the width of the energy bins, $\rho_{\chi}$ is the WIMP density in the neighborhood of the Sun, $\sigma_{\chi T}$ is a reference cross section representing the strength of the WIMP–nucleus interaction,

[TABLE]

is a form factor which depends on the assumed interaction operator, ${\cal G}_{T}(E_{ee},E_{R})$ is the energy resolution smearing function, $\epsilon(E_{ee})$ is an acceptance function including the effect of experimental cuts, and

[TABLE]

is the maximum recoil energy achievable for a WIMP of speed $v$ scattering off a nucleus of mass $m_{T}$ (here $\mu_{\chi T}=m_{\chi}m_{T}/(m_{\chi}+m_{T})$ is the reduced WIMP–nucleus mass).

The reduced response function $\widehat{{\cal H}}_{i}(v)$ has dimensions of velocity. The dimensionless ratio $\widehat{{\cal H}}_{i}(v)/v$ has the immediate physical interpretation as the fraction of an incoming monochromatic flux of speed $v$ that is detected in the $i$-th electron–equivalent energy bin.

We take ${\cal G}_{T}(E_{ee},E_{R})$ to be a Gaussian in $E_{ee}$ centered at $\overline{E}_{ee}=Q(E_{R})\,E_{R}$ and with width $\sigma_{\rm rms}/\mbox{keVee}=0.0091\,(\overline{E}_{ee}/\mbox{keVee})+0.448\,(\overline{E}_{ee}/\mbox{keVee})^{1/2}$. We further assume that ${\cal G}_{T}(E_{ee},E_{R})$ vanishes below the hardware threshold of 1 keVee. We assume a constant quenching factor $Q(E_{R})=0.3$ for sodium.

For the form factor $F(E_{R},v)$ we use the spin-independent form factor of 23Na as given by the Helm form in [62], in correspondence of which $\sigma_{\chi T}$ is the point–like 23Na–WIMP cross section. For the WIMP masses we consider ($m_{\chi}\lesssim 15~{}{\rm GeV}$), this form factor varies negligibly in our analysis: by less than 1% over the 2–4 keVee range where the DAMA modulation is significant, and by $\lesssim 3$% over the whole 2–8 keVee range. By the same token, our analysis applies to all cases in which the variation of the 23Na form factor $F(E_{R},v)$ is negligible. Notice in addition that for such 23Na form factors, any difference in strength between WIMP–proton and WIMP–neutron interactions can be included in the reference cross section $\sigma_{\chi T}$. Thus our analysis applies equally well to 23Na–WIMP elastic scattering that is, for example, any combination of isoscalar and isovector spin-independent interactions, any combination of spin-dependent interactions (for which the 23Na form factors vary by $\lesssim 1$% over the whole 2–8 keVee range), and so on.

We implement the angle-averaged Galactic response functions $\overline{{\cal H}}^{\rm gal}_{0,i}(u)$ and $\overline{{\cal H}}^{\rm gal}_{m,i}(u)$ for elastic spin-independent WIMP–sodium scattering according to Eqs. (A.18–A.19) in the Appendix, with ${\cal H}_{i}(v)$ given by Eq. (4.10) and $v_{\odot}=232$ km/s. Figs. 2 and 3 show these response functions for $m_{\chi}=10~{}{\rm GeV}$. In these figures, we plot the reduced response functions

[TABLE]

The reduced response functions have dimensions of velocity and represent the Galactic response functions normalized in the same way as $\widehat{{\cal H}}_{i}(v)$ in Eq. (4.11). Conceptually, the reduced response functions divided by $u$ give the fraction of WIMPs of Galactic speed $u$ that contribute to the detectable event rate.

As anticipated, we cast the modulation effect as a property of the detector Galactic response function: the only information needed to obtain the reduced Galactic response functions for a given WIMP-nucleus cross section, both modulated and unmodulated, are the experimental properties of the detector and the motion of the detector in the Galaxy.

Once the Galactic response functions are given, the procedure outlined in Section 3, namely solving the extremization problem (3.20–3.25), can be used to estimate the unmodulated signals $S_{0,i}$ from the DAMA data on $S_{m,i}$. The following three considerations are relevant in actually implementing the method of Section 3.

(1) Although no assumptions on $f_{\rm gal}(u)$ are needed in the extremization problem (3.20–3.25), it appears natural to assume that there is a maximum speed for the WIMPs in the Galaxy. Thus we assume that $f_{\rm gal}(u)$ vanishes when $u$ exceeds a maximum velocity $u_{{\rm esc}}$, which we take to be the escape speed from the Galaxy as quoted in [63, 64],

[TABLE]

(This value is within the updated measurements in [65].) We therefore restrict $u\leq u_{{\rm esc}}$.

(2) As seen in Eq. (4.10), and as true in general, the absolute normalization of the response functions contains the unknown factor $\rho_{\chi}\,\sigma_{\chi T}$. However, its value cancels out in the determination of the $S_{0,i}$ from the $S_{m,i}$, essentially because the same factor $\rho_{\chi}\,\sigma_{\chi T}$ appears in both. In principle, one could fix the value of $\rho_{\chi}\sigma_{\chi T}$ that appears in the nonreduced response functions $\overline{{\cal H}}^{{\rm gal}}_{0,i}(u)$ and $\overline{{\cal H}}^{{\rm gal}}_{m,i}(u)$ in the extremization problem (3.20–3.25), solve it as it stands, and then combine the solutions as the value of $\rho_{\chi}\sigma_{\chi T}$ is varied from zero to infinity. Alternatively, and this is the procedure we actually implement, one can use reduced response functions $\widehat{{\cal H}}^{{\rm gal}}_{0,i}(u)$ and $\widehat{{\cal H}}^{{\rm gal}}_{m,i}(u)$, which do not contain the factor $\rho_{\chi}\sigma_{\chi T}$, rescale the coefficients $\lambda_{k}$ to

[TABLE]

so that

[TABLE]

and solve the modified extremization problem in which the normalization condition $\sum_{k=1}^{K}\lambda_{k}=1$ is replaced by the condition

[TABLE]

Since this condition is already contained in the conditions $\widehat{\lambda}_{k}>0$, one of the moment conditions effectively disappears, and the extreme values of the $S_{0,i}$ can be found with sums of up to $N$, instead of $N+1$, streams in velocity space, provided the extreme distributions contain the rescaled coefficients $\widehat{\lambda}_{k}$ instead of $\lambda_{k}$.

(3) Having dropped the normalization condition on the $\lambda_{k}$ as described in the previous paragraph, the $K$ streams ($1\leq K\leq N$) of an extreme distribution must have speeds $u_{k}$ such that the $K$ $N$-dimensional vectors $(\widehat{{\cal H}}_{m,1}(u_{k}),\ldots,\widehat{{\cal H}}_{m,N}(u_{k}))$ are linearly independent. Now the experimental threshold in observed energy (1 keVee in our treatment of DAMA; see our discussion of ${\cal G}_{T}(E_{ee},E_{R})$ after Eq. (4.11)) induces a region below threshold in velocity space, comprised of all speeds $u$ for which the response functions vanish simultaneously, $\widehat{{\cal H}}_{m,i}(u)=0$ ($i=1,\ldots,N$) for $u$ below threshold. For example, in the isotropic case we consider, the threshold speed for the modulated Galactic response functions turns out to be

[TABLE]

Thus, if the velocity of one or more of the $K$ streams is below threshold, the $K$ vectors mentioned above are not linearly independent (one or more of them is the zero vector). On the other hand, the streams with velocity below threshold do not contribute to the $S_{m,i}$ signals at all (indeed, all the response functions are zero for these streams). Thus a sum of $K$ streams in which some streams are below threshold is effectively an extreme distribution with less than $K$ streams. Since we let $K$ vary from 1 to $N$, it is obvious that it is enough for extreme distributions to include only streams above threshold. Therefore we allow only $u>u_{{\rm thr}}$.

The practical implementation of the method described above is conceptually quite simple, although the use of the parametrization (3.17) for the extreme distributions of the moment set requires to explore a parameter space of large dimensionality ($2N=24$ in our $N=12$ case with isotropic Galactic velocity distribution; $4N=48$ if we explored anisotropic Galactic velocity distributions). This kind of task is efficiently performed by using the technique of Markov chains, which makes use of the likelihood function itself to optimize the sampling procedure.444Alternatively one could use grids in velocity space and increase the grid resolution to increase the precision of the computed extreme values, or one could use algorithms to directly maximize and minimize $S_{0,i}$ under the given constraints. We tried some of those methods without success. We chose to generate Monte-Carlo samples of the likelihood to gain confidence on our method. To this aim we use the Markov–Chain Monte Carlo (MCMC) code emcee [66] to generate a large number of sets $\{u,\widehat{\lambda}\}=\{u_{1},\ldots,u_{K},\widehat{\lambda}_{1},\ldots,\widehat{\lambda}_{K}\}$ of Galactic speeds $u_{k}$ and coefficients $\widehat{\lambda}_{k}$ with $1\leq k\leq K$, $1\leq K\leq N$, $u_{{\rm thr}}<u_{k}\leq u_{{\rm esc}}$ and $\widehat{\lambda}_{k}>0$. For each value of $K=1,...,N$, we generate a Markov chain of $5\times 10^{6}$ points using 250 independent walkers and a standard Metropolis-Hastings sampler.

For each MCMC-generated set $\{u,\widehat{\lambda}\}$ we calculate both $\chi^{2}=-2\ln{\cal L}$ and $S_{0,i}$ for $i=1,\ldots,N$ $(N=12)$ and produce $N$ scatter plots of $\chi^{2}$ vs. each of the $S_{0,i}$. The resulting $(\chi^{2},S_{0,i})$ scatter plots are shown in Figs. 4, 5 and 6 for three representative values of the WIMP mass: $m_{\chi}=5~{}{\rm GeV}$, 10 GeV, and 15 GeV, respectively. In each of these figures, the 12 panels correspond to one of the 12 DAMA energy bins between 2 keVee and 8 keVee. In each panel, the boundary of the region covered by points gives the profile likelihood for the $S_{0,i}$ in that panel. More precisely, the boundary in the $i$-th panel is the graph of $\chi_{i}^{2}=-2\ln{\cal L}_{i}(S_{0,i})$, where ${\cal L}_{i}(S_{0,i})$ is the profile likelihood for $S_{0,i}$ with the velocity distribution treated as a continuum of nuisance parameters. We use two colors for the points, red and black. The points in red (black) have $\chi^{2}_{i}\leq\chi^{2}_{i,\rm min}+1$ ($\chi^{2}_{i}>\chi^{2}_{i,\rm min}+1$), where $\chi^{2}_{i,\rm min}$ is the minimum of the $\chi^{2}_{i}$. The value of $S_{0,i}$ where $\chi^{2}_{i}$ is minimum is the maximum-likelihood estimate of $S_{0,i}$. The lowest and highest values of $S_{0,i}$ in the $\chi^{2}_{i}\leq\chi^{2}_{i,\rm min}+1$ region (points in red, which occur at the top boundary of the red region) are the endpoints of the 1$\sigma$ confidence interval for $S_{0,i}$.

To illustrate that our analysis includes the correlations among the $S_{0,i}$’s, in Fig. 7 we show the regions $-2\Delta\ln{\cal L}_{\rm p}(\{S_{0,i}\})\leq 1$ (red points) and $-2\Delta\ln{\cal L}_{\rm p}(\{S_{0,i}\})\leq 3$ (black points) in the plane $S_{0,1}$–$S_{0,2}$ of the unmodulated signals in the first and second energy bins for $m_{\chi}=10$ GeV. The 1$\sigma$ confidence intervals on $S_{0,1}$, for example, are obtained by projecting the red ellipse onto the $S_{0,1}$ axis.

The maximum-likelihood estimates and the 1$\sigma$ confidence intervals of $S_{0,i}$ ($i=1,\ldots,12$), obtained from Figs. 4, 5 and 6 (horizontal blue lines best visible in Fig. 6),555For the higher masses, the sampling becomes sparse at large values of $S_{0,i}$ and the exact value of $S^{\rm sup}_{0,i}(L_{0})$ becomes harder to determine from the figures. This problem is there because the MCMC concentrates points near the maximum-likelihood value of $S_{0,i}$. A few attempted methods of direct minimization of the $\chi^{2}$ at fixed $S_{0,i}$ did not converge. We have performed additional dedicated MCMC runs focused at large values of $S_{0,i}$ and the value of $S^{\rm sup}_{0,i}(L_{0})$ did not change significantly. are listed in Table 1 and plotted versus the electron–equivalent energy $E_{ee}$ in Figs. 8, 9 and 10, for $m_{\chi}=5$, 10, and 15 GeV, respectively. These are the main results of our paper. Table 1 also lists the DAMA modulated amplitudes $S_{m,i}$ (from Fig. 8 of [2]) and the DAMA background+signal rates $B_{i}+S_{i}$ (from Fig. 27 of [67], rebinned from 0.25-keVee- to 0.5-keVee-width bins).

We see that for each $m_{\chi}$, the $S_{0,i}$ decrease with energy, as reasonable for a WIMP signal. We also see that the error bars on the $S_{0.i}$ are small enough to allow a rather good determination of the unmodulated signal, although the error bars become very asymmetric for $m_{\chi}=15~{}{\rm GeV}$. Moreover, an examination of Table 1 leads to the conclusion that the unmodulated signals $S_{0,i}$ are much smaller than the DAMA background+signal measurements. Since it is not trivial to identify what contributes to the DAMA background, we refrain from subtracting an estimated model background like done for instance in [68]. It suffices for us to conclude that the $S_{0,i}$ values we estimate in this paper are reasonable and compatible with the measured DAMA background+signal level.

It is finally interesting to estimate the fraction of the signal that is modulated. In the energy range where the bulk of the DAMA modulation is present, i.e., 2 keVee$<E_{ee}<$4 keVee, we find, for $\chi^{2}_{i}\leq\chi^{2}_{i,\rm min}+1$ (examining the ratio $S_{m,i}/S_{0,i}$ for each point in red in Figs. 4, 5 and 6):

[TABLE]

The modulated amplitude ranges between a few percent and about 25% of the unmodulated amplitude, depending on the WIMP mass. This is in line with expectations for a signal due to dark matter WIMPs.

This conclusion is limited to the case at study, i.e., to velocity distributions that are isotropic in the frame of the galaxy. However it is worth pointing out that solutions with higher modulation fractions are present in the allowed parameter space also in the isotropic case, but are rejected by the DAMA data. This can be seen in Fig. 11, where we plot the ratios between modulated and unmodulated response functions ${\cal H}^{{\rm gal}}_{m,i}/{\cal H}^{{\rm gal}}_{0,i}$ for $m_{\chi}$=10 GeV. In this specific example it is clear that, if $f_{\rm gal}(u)$ is parametrized according to Eq. (3.17) with a set of $u_{j}$ velocities all below $\sim$60 km/s, modulation fractions close to 1 are possible. These configurations, however, have a very small likelihood and do not appear in the 1$\sigma$ confidence intervals. Notice that at even lower velocities, values of the ${\cal H}^{{\rm gal}}_{m,i}/{\cal H}^{{\rm gal}}_{0,i}$ ratios larger than one signal the need of additional harmonics besides the cosine in the description of the modulation time dependence of the signal.

5 Conclusions

We have estimated the unmodulated signal corresponding to the DAMA modulation if interpreted as due to scattering of dark matter WIMPs off 23Na in the DAMA detector. Our analysis, covering WIMPs lighter than $\sim 15$ GeV, is to large extent independent of the dark halo model (we profile the likelihood over velocity distributions that are isotropic in the Galactic rest frame) and of the particle physics model (we cover all 23Na–WIMP interactions in which the 23Na cross section vary little over the 2–4 keVee energy range where the DAMA modulation is significant). The method outlined in Section 3 is however very general and valid for any class of velocity distribution.

We have presented and used an exact and sound mathematical set-up to profile the likelihood over a continuum of nuisance parameters (the whole WIMP velocity distribution). We have used the profile likelihood for each of the unmodulated rates $S_{0,i}$ ($i=1,\ldots,12$ indexing the first 12 DAMA energy bins) to find their maximum-likelihood estimates and their 1$\sigma$ confidence intervals.

Our halo-independent estimates of the unmodulated rates are reasonable and in line with expectations for a signal from WIMP dark matter. The unmodulated rates we estimate give a modulated/unmodulated ratio ranging between a few percent and $\sim 25$%. The unmodulated rates are comfortably below the background+signal level measured by DAMA.

Acknowledgements

This work was completed during Stefano Scopel’s sabbatical semester at the University of Utah in the Spring of 2016, and was presented at various conferences thereafter. We thank the many audience members, and in particular Riccardo Catena, for asking many questions that have helped sharpening the writing of this paper. Paolo Gondolo thanks Sogang University for the financial support and the kind hospitality during his research visits. This work was partially supported by NSF award PHY-1415974 at the University of Utah and by the National Research Foundation of Korea grant number 2016R1D1A1A09917964 at Sogang University.

Appendix A Angle-averaged Galactic response functions

In the application of the optimization method to the modulated/unmodulated DAMA rates in this paper we focus on velocity distributions that are isotropic in the Galactic frame, i.e., for which

[TABLE]

where $u=|{\bf u}|$. In this case, integrals of the form (2.7) become

[TABLE]

where

[TABLE]

is the angle-averaged Galactic response function, and

[TABLE]

is the speed distribution in the Galactic frame normalized to one ($\int_{0}^{\infty}\overline{f}_{{\rm gal}}(u)\,du=1$). A Fourier time-series analysis then gives, in this case of isotropic Galactic velocity distributions,

[TABLE]

where

[TABLE]

For nondirectional dark matter detectors, the laboratory response functions ${\cal H}_{i}({\bf v})$ are isotropic, i.e., ${\cal H}_{i}({\bf v})={\cal H}_{i}(v)$. In this case, Eqs. (2.8), (A.3) and (A.7–A.8) give

[TABLE]

where

[TABLE]

These integrals can be recast as integrals over $v=|{\bf u}-{\bf V}|$ with analytically computed kernels by means of the following manipulations. Let the polar angles $(\theta_{u},\phi_{u})$ of ${\bf u}$ in $d\Omega_{u}=\sin\theta_{u}\,d\theta_{u}\,d\phi_{u}$ be defined with respect to the direction of ${\bf V}$. Then

[TABLE]

Perform the trivial integration over $\phi_{u}$, and change integration variable from $\cos\theta_{u}$ to $v$. Then interchange the order of the $t$ and $v$ integrations. This gives

[TABLE]

Here $\Theta(x)=1$ if $x$ is true and $\Theta(x)=0$ if $x$ is false. The integrals over the time $t$ in Eqs. (A.13) and (A.14) can be performed analytically for a circular Earth orbit. Write

[TABLE]

where $t=0$ when ${\bf v}_{\oplus}$ and ${\bf v}_{\odot}$ form their smallest angle $\beta$, i.e., ${\bf v}_{\odot}\cdot{\bf v}_{\oplus 1}=v_{\odot}v_{\oplus}\cos\beta$ and ${\bf v}_{\odot}\cdot{\bf v}_{\oplus 2}=0$. Then ${\bf v}_{\odot}\cdot{\bf v}_{\oplus}(t)=v_{\odot}v_{\oplus}\cos\beta\cos(\omega t)$ and

[TABLE]

where

[TABLE]

In terms of these variables,

[TABLE]

where we define

[TABLE]

The functions $K_{0}(x,y)$ and $K_{m}(x,y)$ are symmetric in $x$ and $y$, and the integrals appearing in them can be expressed in terms of the incomplete elliptic functions $E(\phi|m)$ and $F(\phi|m)$.

[TABLE]

where

[TABLE]

and

[TABLE]

We used Eqs. (A.18) and (A.19) to compute the DAMA angle-averaged Galactic response functions.

We notice lastly that the angle–averaged Galactic response functions for the $\sin(\omega t)$ Fourier mode vanish identically. This follows from the observation that the formula for the angle–averaged sine modulation response functions is analogous to equation (A.14) with the replacement of $\cos(\omega t)$ with $\sin(\omega t)$ in the numerator but the same dependence on $V$. Since $V$ depends on $\cos(\omega t)$ only, the integrand is thus a product of an odd and an even function of $\omega t$ and vanishes when integrated over a whole period.