Time-dependent rate of multicomponent dark matter: Reproducing the DAMA/LIBRA phase-2 results

The current paradigm for dark matter direct detection is to assume that the dark sector is solely composed of a single particle species. In this short paper, we make the observation that dark matter comprising both a light and a heavy component that modulate out of phase leads to interesting phenomenology in annual modulation experiments. For an illustrative example, we use the recently released DAMA/LIBRA phase-2 results with a lower energy threshold. Immediately after, it was argued that a one-component spin-independent dark matter explanation of the observed annual modulation is strongly disfavored or excluded unless isospin-violating couplings are invoked. We show that a simple two-component extension can reproduce the observed spectrum without the need to invoke fine-tuned couplings. Using the publicly available DAMA/LIBRA data, we perform a fit of the DAMA/LIBRA energy spectrum of the annual modulation amplitude to a scenario with two dark matter components. We also take into account how gravitational focusing affects the phases of the light and a heavy components differently, which leads to nontrivial effects in the total time-dependent rate. Our results show that there exists a unique solution in agreement with the data in the simplest case of isospin-conserving couplings with equal cross sections. The distinctive features found in this work are crucial for a dark matter interpretation of any observed annual modulation.


Introduction
Dark matter (DM) is one of the greatest mysteries of nature. Until now evidence for its existence stems from its gravitational interactions. There is one potential exception coming from direct detection experiments: the long-standing annual modulation signal observed by the DAMA/LIBRA collaboration [1,2] (denoted by DAMA in the following). The significance for the annual modulation of the complete data set is 12.9 σ. DAMA is unable to distinguish, on an event-by-event basis, electron recoils from nuclear recoils on its sodium iodide (NaI) crystals. However, by measuring the annual modulation (c.f., the total rate) and by performing different checks (only single hit events, phase close to June 2nd, period of 1 year, amplitude, etc.) the DAMA collaboration rejects the possibility that modulated backgrounds (neutrons, muons, solar neutrinos, radon, etc.) are the cause of the signal. Also several independent studies have been performed in the literature without finding any conclusive explanation [3][4][5][6][7]. Until the very latest phase-2 results, under reasonable particle physics and astrophysical assumptions the signal was consistent in both amplitude and phase with that expected from Weakly Interacting Massive Particles (WIMPs). Assuming the Standard Halo Model (SHM) and elastic scattering, the best-fit masses and cross sections are wellknown: a light DM with mass ∼ 10 GeV scattering mainly on sodium (light mass solution), or a heavy DM with mass ∼ 70 GeV scattering mainly on iodine (heavy mass solution) [8][9][10][11][12][13][14][15].
The main issue with a DM interpretation of the DAMA modulation is that no other independent experiments have observed any DM signals, whereas compatibility with DAMA would require such signals to have been seen under currently understood particle physics and astrophysics scenarios. For example, LUX [16], XENON1T [17] and PandaX [18]) appear to exclude a DM origin of the DAMA observations. There is no currently accepted explanation that reconciles DAMA's signal with the absence of results in all other experiments, e.g., see Refs. [12,13,19]. Furthermore, in the last few years halo-independent methods have been designed and applied to DAMA and appear to exclude a DM origin of the results independent of the assumed velocity distribution [20][21][22][23][24][25][26]. This has motivated a large experimental effort to try to reproduce the DAMA experiment with NaI crystals in order to independently either confirm it or reject it [27][28][29][30][31]. The SABRE experiment [29] plans to have a northern and southern hemisphere pair of NaI detectors to search for a seasonal correlation or anti-correlation of any DAMA-like annual modulation signal [32]. Interestingly, the COSINUS experiment [30] aims to also measure the constant rate by developing a cryogenic detector. As studied in Ref. [33], a null result in the latter experiment may rule out model-independently a DM explanation of DAMA.
Very recently, the DAMA collaboration released the long-waited phase-2 results with a lower energy threshold, with three new energy bins below 2 keVee [34]. As first pointed out in version 3 of Ref. [33], and studied also in Ref. [35], the consistency of the DM interpretation of DAMA's signal is under question both for the light and the heavy DM mass solutions mentioned above, even before considering its compatibility with other null results experiments. The issue is that, below 2 keVee the two standard DM solutions behave very differently with decreasing energy: the light DM gives rise to scatterings off iodine, increasing its rate significantly in the case of spinindependent (SI) interactions, while for the heavy DM, that scatters only off iodine, the modulation amplitude decreases giving rise eventually to a phase flip. This was already pointed out in Ref. [36]. With the new data, we also obtain in this work that the one-component spin-independent (SI) DM explanation of the annual modulation is disfavoured or excluded unless isospin-violating (IV) couplings are invoked. In particular, we find that the light DM solution is excluded at 5.3σ, while the heavy one is excluded at 3.0σ, in agreement with Refs. [33,35]. In the case of IV couplings, an adequately tuned destructive interference between the DM couplings to neutrons and protons can still correctly reproduce the spectrum by suppressing the interactions with iodine.
The main observation of this short letter is that the modulation amplitude observed by DAMA at low energies can be reproduced in a natural way by a combination of two DM particles, without the need to invoke fine-tuned IV couplings. As we quantify explicitly below, two DM components give an excellent fit to the data for SI isospin-conserving (IC) couplings, equal to, or better than, that of the one-component scenario with IV couplings. From a theoretical perspective, it is also not difficult to envisage a dark sector, which accounts for 25% of the energy density of the Universe, to consist of more than one light stable state, similar to the way in which the visible sector (which accounts for only 5% of the energy budget) has a few stable particles: protons, electrons, the lightest neutrino, Helium nuclei, etc. Only a few works have studied the case of multi-component DM in direct detection, focusing on constant event rates [37][38][39][40][41][42][43][44]. In the following, we adopt a purely phenomenological approach to try to reproduce the observed modulated signal without going into further details regarding the model building, which would affect the interactions and the abundances of the DM components. This will be studied in a future work [45].
This letter is structured as follows. In Sec. 2 we introduce the relevant notation to describe the DM annual modulation signal in direct detection experiments. We discuss the fitting procedure used in Sec. 3.1. In Sec. 3.2 we show results of the fits of the vanilla one-component scenario with and without IV couplings to the full DAMA data set. We show the main results of this paper for two-component DM in Sec. 3.3. We present best-fit values and two-dimensional confidence level contours for two-component DM. This is done for the simplest case with IC couplings and also for the more general case with IV interactions. We give our conclusions and final remarks in Sec. 4.

The dark matter annual modulation signal
In this section, we present the relevant expressions for the direct detection of DM comprising two components with masses m 1,2 (we take m 1 < m 2 such that component 1 is always lighter than component 2), SI cross-sections with protons σ p 1,2 , and local energy densities ρ 1,2 . We impose the constraint ρ 1 +ρ 2 = ρ loc , where ρ loc is the observed DM mass density, which we fix to 0.4 GeV/cm 3 . Following the notation of Ref. [44], it is useful to define the following ratios such that ρ 2 = ρ loc r ρ /(1 + r ρ ). In this work, we focus on the DM annual modulation signal [46,47]. For the case of two-component DM, we can write the amplitude of the annual modulation rate in a detector with target nucleus labelled by (A, Z), for SI interactions, as (see for instance Refs. [15,23,24,48]) Here F A (E R ) is the SI nuclear form factor of element A, for which we use the Helm parametrisation [49,50]. x A is the mass fraction of the target A in the detector. The total modulation amplitude is given by The modulated amplitude for DM consisting of just one component is obtained from Eq. (2) in the limit r ρ = r σ = 0. In the following we use the label α = 1, 2 for the two DM components. In addition to the local energy densities, the astrophysics enters in Eq. (2) through where t J (t J + 0.5) measured in years corresponds to June (December) 2nd, which for minimum velocities v  A (E R ) become negative. In other words, for small enough v (α) m,A , the phase of the modulation flips by 6 months. As we show below, this is precisely what happens for heavy DM components scattering off iodine in DAMA at the lowest energies. Here we defined the halo integral as where v (α) m,A is the minimal velocity of the DM particle α required to produce a recoil of energy E R in element A, and f det ( v, t) describes the distribution of DM particle velocities in the detector rest frame, which we assume to be equal for both DM species. It can be written in terms of the galactic velocity distributions by doing a Galilean boost, is the velocity vector of the Earth in the galaxy rest-frame. We use the SHM, with a Maxwellian velocity distribution with a cut-off at the escape velocity v esc = 550 km s −1 and N esc = erf(z) − 2 √ π ze −z 2 , where z = 3/2 v esc /σ Hα . The velocity dispersion of each DM component, σ Hα , may in principle be different. Indeed, if both components reach equilibrium in the halo, the velocity dispersions are expected to become mass dependent [51]. In this case one can write σ Hα =σ H m/m α , wherem = n α m α / n α , with n α = ρ α /m α the number density of the DM particle α, andσ H ∼ 270 kms −1 is the canonical velocity dispersion. Throughout the paper we generally assume equal velocity dispersions for both components, σ H 1 = σ H 2 =σ H , but we briefly discuss how our results change in the case in which they are given by the mass-dependent relationship provided above.
In order to properly fit a model to the DAMA data we must first take into account the finite energy resolution of the detector. This is done via the convolution of the theoretical modulation amplitude M(E R ) and the differential response function where (E ee ) is the detector efficiency and E ee is the electron equivalent energy, related to the true recoil energy E R through the target-dependent quenching factors E ee = Q A E R . We use Q Na = 0.3 and Q I = 0.09 for the quenching factors of sodium and iodine. We also use the differential response function provided in Ref. [35], so that we can compare our results to the one-component case studied there. The modulation amplitude in bin i, M i , is obtained by averaging the corrected modulation amplitude of Eq. 6 over each bin.
3 Fitting the full DAMA/LIBRA data set In the next section, we briefly discuss details of our numerical scan and the data used. Afterwards, we show results of the fits to the one and the two-component DM scenarios in Secs. 3.2 and 3.3, respectively.

Analysis methods
In order to test the compatibility of a particular DM model with a given set of experimental data, one can employ the method of maximum likelihood. For the case of DAMA we parameterise the likelihood function with a binned Gaussian distribution where N is the number of bins and M i (θ) is the expected modulation amplitude in bin i, defined below Eq. 6, which depends on different DM parameters denoted by θ = (θ 1 , ..., θ n ), where n is the total number of fitted parameters. µ i is the central value of the observed annual modulation  Table 1: Parameters and prior ranges used for the two-component dark matter fits. In the fourdimensional fit, the four parameters are used with κ 1 = κ 2 = r σ = 1.
amplitude in bin i, and σ i its error. Maximising the likelihood in Eq. (7) with respect to the parameters θ is equivalent to minimising -2 times the log-likelihood as follows: This is the familiar chi-square test statistic which, as shown by Pearson [52], in the limit of large µ i follows a χ 2 distribution with N − n degrees of freedom (dof). Since the mean of a χ 2 distribution is the number of degrees of freedom (dof), the goodness of fit can determined by the value of the χ 2 at the minimum, χ 2 min , divided by the number of dof. In our χ 2 fit we use the whole DAMA annual modulation data set, which combines results from DAMA/NaI and DAMA/LIBRA phases 1 and 2. The total exposure is 2.46 ton y. The results are provided in slide 22 of Ref. [34]. We use the data of Tab. I of Ref. [35], which gives the observed modulation amplitude M i in N = 10 bins in the energy range [1,20] keVee.
We use the open source software Minuit [53] to compute the best-fit points, which we give in Tabs. 2 and 3. For completeness, we also compute the corresponding p-values for the fits, as well as the corresponding number of equivalent two-sided Gaussian standard deviations Z. We consider a 'good' fit as one that yields a p-value > 1 × 10 −3 , which corresponds to Z < 1.96. Of course, this is to some extent an arbitrary choice. For the two-component scenarios, we also employ the MultiNest implementation of the nested sampling algorithm [54][55][56]. The latter is well suited for dealing with a highly multimodal target function, which is the case for two-component scenarios. We show results with 10 4 live points and a tolerance of 0.01. For the two-component DM fits we use the priors shown in Tab. 1. We determine the distribution of the profile likelihood ratio (PLR) L/L max throughout the parameter space from the obtained samples, and derive the frequentist 1 and 2σ C.L. contours using pippi [57]. 2

One-component dark matter
We perform fits of one-component DM (1DM) to the DAMA data. We do this for SI interactions with both IC and IV couplings. The former scenario involves just 2 free parameters, while the 2 We have checked the consistency of the best-fit points obtained with Multinest and with Minuit by adequately choosing the search ranges in the latter. The reason is that Multinest is well suited to finding different local minima in multi-dimensional parameter spaces, but, unless a very low tolerance and a high efficiency are used, which would make the scans very computationally expensive, the best-fit points may not be as accurate as those given by Minuit. Therefore in the tables we show those obtained with the latter.  Table 2: Results of the one-component dark matter fit to DAMA data. The DM mass m is in GeV and the scattering cross section with protons σ p 1 is in units of 10 −40 cm 2 . The different fits are for both the light and the heavy dark matter mass solutions, which correspond to scatterings mainly in Na and I, respectively. We show isospin-conserving (IC) and isospin-violating (IV) couplings. The dashes refer to the parameters that are fixed to 1 and therefore are not varied in the fit. The last two columns provide the p-value and the number of Gaussian-equivalent standard deviations Z (two-sided).
latter has 3. We show the results for the light and the heavy DM mass solutions. In Tab. 2 we show the best-fit points obtained, as well as the values of χ 2 min /dof, the p-values and Z. As can be observed, unless IV couplings are invoked, there is no good fit to DAMA.
It is quite easy to visualise how a combination of both light and heavy DM could give a spectrum that reproduces the observed signal in the lowest energy bins. This is the main idea of this paper, which we study in the next section.

Two-component dark matter
In this section we perform a fit of the annual modulation amplitude generated by two DM particles (2DM) to the DAMA data. We first take IC couplings and fix r σ = 1, leaving r ρ free. These models have 4 free parameters. The energy density and the cross section enter identically in the rate, M (α) ∝ ρ α σ α . Therefore, apart from the overall normalisation, one can fix r σ and consider r ρ as a free parameter. We have also checked that the other scenario (fixing r ρ = 1 and leaving r σ free) gives a similar fit with a cross section that is roughly twice that of the case of r σ = 1 because r ρ 3, confirming the expectations. We then do a general fit with free r ρ and r σ and IV couplings, which involves 7 free parameters.
We show in Tab. 3 the results of the fits of two-component DM (2DM) to the annual modulation signal observed by DAMA. We show the best-fit points of the different fits performed, as well as the values of χ 2 min /dof, the p-values and the Gaussian-equivalent number of standard deviations Z. In the first two rows we show the two IC best-fit solutions, while in the last one we show the more general case of IV ones. For IC couplings two solutions exist. One can observe that solution B, with the larger mass splitting between the two DM particles, is a better fit than solution A. Although the general IV model (last row) produces a smaller χ 2 min than both IC cases, its larger number of free parameters yields a smaller p-value and hence a worse fit. Therefore, in the following we only focus on the four-parameter IC cases.
Let us also briefly mention how our results change in the case of mass-dependent velocity dispersions. To investigate their impact, we performed a dedicated Minuit fit around the param-    eters given by the IC solution A in Tab. 3, in which the velocity dispersions were allowed to be mass-dependent. We obtain χ 2 min /dof = 10.69, with best-fit values m 1 = 19.66 GeV, m 2 = 75.8, σ p 1 = 0.13 and r ρ = 3.02. Therefore, the results are very similar to the mass-independent velocitydispersion case.
In Fig. 1 we show with a solid black line the binned modulation amplitude for the best-fit points for two-component DM in the case of IC couplings for solution A (left panel) and solution B (right panel). We also show the individual contributions in dotted green (component 1) and dashed blue (component 2) lines, and the DAMA experimental results as red points. We have checked that solution A corresponds to scattering of both DM components dominantly off iodine. The fact that the lighter component scatters dominantly off iodine with a negligible contribution from sodium is due to two factors: first, the smaller quenching factor in iodine compensates its larger mass, translating into smaller v For energies below such phase flips, the modulation amplitude of component 2 becomes negative and therefore there is a partial cancellation in the combined modulation amplitude between the individual DM contributions. As can be seen, such two-component DM combinations reproduce very well the modulated spectrum observed by DAMA.
In Fig. 2 we show the profile likelihood ratio (PLR) density L/L max overlaid with 1 and 2σ C.L. contours. We only show four illustrative slices of the full four dimensional parameter space. One immediately notices two distinct regions of high PLR density which correspond to the two solutions A and B, which we indicate in the plot. There is a slight preference for solution B (see Tab. 3). One can observe that the 1 and 2σ regions are quite extended, with a significant range of parameters being able to reproduce the DAMA results.
In the m 1 -m 2 plane (top left panel ) one can observe that, for each of the solutions, one of the DM components has a narrow range of masses (for example component 1 for solution B, which always gives the largest contribution at all energies) while the other one (component 2 for solution B) spans a large mass range, with a contribution that is subdominant. There is a sizable range of mass splittings that provide a good fit. In the log 10 (r ρ )-log 10 (σ p 1 ) plane (top right panel ) one can observe that the two solutions are aligned in a band with negative correlation, i.e., the larger the cross section the smaller the ratio of energy densities. 3 In the log 10 (r ρ )-m 2 plane (bottom left panel ) we see that it is solution B (with lighter m 1 and heavier m 2 ) that has the larger cross section. In the log 10 (r ρ )-m 2 plane (bottom right panel ) the large degeneracy in the case of solution B between m 2 and r ρ is clearly visible. It is due to the fact that, for m 2 m I , the minimum velocity of the second component in iodine, v (2) m,I , approaches a constant value, so that its contribution to the amplitude only depends on the ratio r ρ /m 2 .

Conclusions
The observed annual modulation in the low energy bins below 2 keVee in the Phase-2 run DAMA/ LIBRA results is no longer compatible with a one-component dark matter interpretation in the simplest spin-independent isospin-conserving scenario. The reason is easy to understand: a light DM particle that scatters predominantly off sodium produces a spectrum that steeply increases at low recoil energies due to iodine scatterings becoming allowed; conversely, a heavy DM mass produces a spectrum that falls off at the lowest energy bins. The issue is that the observed values in the lowest energy bins happen to lie, roughly speaking, between the expected modulated spectrum of a light and that of a heavy DM state. It then seems natural that a two-component model comprising both a light and heavy component can revive the DM interpretation of the low threshold data. Indeed, below a certain energy value, the modulation amplitude of the heavy DM state becomes negative, so that in the combined modulation amplitude the individual DM contributions counterbalance each other to some extent. Of course, their masses and their relative local energy densities (or similarly, the relative strength of their interactions) play an important role regarding the compatibility of the two-component DM scenario as a solution to the observed spectrum by DAMA.
In this short letter we performed a fit of two-component DM to DAMA data. Contrary to the case of one-component DM where now only IV couplings provide a good fit, we find very good agreement to the data also for IC couplings. Two different solutions with qualitatively different spectral behaviour are found, shown in Fig. 1. On one of them (A), scatterings are predominantly on iodine, with a crossing between the spectrum of the two individual DM components at roughly 2 keVee. On the other one (B), the lighter component scatters off both sodium and iodine at the lowest energies and the individual spectra do not cross. The results of the fits are summarised in Tab. 3 and Fig. 2, which involve reasonable values of the relative energy densities and the cross section. Although model-dependent assumptions regarding the velocity dispersions or the relic abundance of the different DM components may modify the results, the main conclusions are expected to hold, that is, that the DAMA spectrum can be well reproduced by two DM particles. 4 3 Therefore the figures in the log 10 (rρ)-m1 and the m2-log 10 (σ p 1 ) planes can be easily obtained from the ones in Fig. 2 and we do not show them. 4 In this work we focused on the amplitude of the annual modulation. Let us also mention that the phase of the Two-component DM therefore looks like a natural solution to the first part of the DAMA puzzle: the compatibility of the spectrum with that expected from DM under standard astrophysical and particle physics assumptions. The second part of the DAMA puzzle, that is, the compatibility of DAMA data with other null results is not solved. It is interesting to speculate, however, how it might be somewhat improved in the scenario considered here.