# Bayesian Analysis of time dependence of DAMA annual modulation amplitude

Srinikitha Bhagvati  and Shantanu Desai 11footnotetext: Corresponding author.

###### Abstract

We implement a test of the variability of the per-cycle annual modulation
amplitude in the different phases of the DAMA/LIBRA experiment using Bayesian
model comparison. Using frequentist methods, a previous study [1] had
demonstrated that the DAMA amplitudes spanning over the DAMA/NaI and the first
phase of the DAMA/LIBRA phases, show a mild preference for time-dependence in
multiple energy bins. With that motivation, we first show using Bayesian
techniques that the aforementioned data analyzed in [1] show a moderate
preference for exponentially varying amplitudes in the 2-5 and 2-6 keV energy
intervals. We then carry out a similar analysis on the latest modulation
amplitudes released by the DAMA collaboration from the first two phases of the
upgraded DAMA/LIBRA experiment. We also analyze the single-hit residual rates
released by the DAMA collaboration to further look for any possible time-
dependency. However, we do not find any evidence for variability of either of
the two datasets by using Bayesian model selection. All our analysis codes and
datasets have been made publicly available.

##  1 Introduction

The identity of cold dark matter is one of the most profound unsolved
mysteries in Cosmology and fundamental Physics. For close to 90 years, we have
known that most of the mass of our galaxy and the rest of the universe is
composed of dark matter [2]. Further observations have shown that it is non-
baryonic, collisionless, and that it decouples while moving at non-
relativistic velocities [3, 4]. Hence, the name “cold dark matter”. The latest
CMB and LSS results indicate that the total density of cold dark matter is
about 25% [5].

A large number of cold dark matter candidates have since been proposed,
spanning many orders of magnitude in mass and interaction strengths [6, 7].
Among those, the Weakly Interacting Massive Particle (WIMP) is one of the most
well motivated and widely studied cold dark matter candidate. This is because
any elementary particle which is a thermal relic from the Big-Bang having
electro-weak scale interactions with ordinary matter satisfies all the
properties needed for a cold dark matter candidate, such as the correct relic
abundance and non-relativistic velocities at the time of decoupling [8, 9,
10].

The detection of WIMPs using coherent elastic scattering with nuclei was
proposed more than thirty years ago [11, 12]. Since then, the experimental
detection of WIMPs and measurements of their couplings with ordinary matter
have been the subject of a large number of experiments in various underground
laboratories [13]. The DAMA experiment is the only one among all these
experiments so far to have claimed to detect dark matter.

The main basis for this claim by the DAMA collaboration is the observation of
an annual cosine-based modulation in the residual single-hit count rates in
the 2-6 keV energy intervals and thereabouts. This observed modulation has the
correct phase (peak in June and minima in December) and amplitude that is in
accord with the theoretical predictions of dark matter halo models [14, 15].
Since its inception in 1995, the statistical significance of the DAMA signal
has been growing steadily with accumulated data over the different phases of
the experiment [16, 17, 18, 19, 20]. The first data release showed only a
$3\sigma$ significance for the annual modulation [21]. This significance
increased to $6.3\sigma$ by the end of the first phase of the DAMA experiment
(DAMA/NaI) [16]. The latest data release in 2018 shows a significance greater
than $8\sigma$ across all the energy intervals, which could contain a WIMP
signal. The net significance in the 2-6 keV energy bin is equal to
$12.9\sigma$, after combining all the accumulated data until the end of phase
2 of the upgraded DAMA experiment (DAMA/LIBRA) [20].

However, other direct detection experiments having different target materials
have failed to confirm this result and their upper limits completely rule out
the allowed regions for WIMP mass and cross-section deduced from the DAMA
annual modulation results under the canonical assumptions for the dark matter
[22, 23, 24, 25]. Experiments having the same target material as DAMA (NaI)
such as COSINE-100 and ANAIS-112, which have recently commenced operations
have also obtained null results [26, 27]. A large part of the DAMA parameter
space has also been ruled out by some indirect dark matter searches [28].
Despite all this, it is important to scrutinize the DAMA dataset to see if it
contains additional signatures which could help distinguish between the dark
matter and an instrumental origin for the observed annual modulation.

In this work, we search for a time dependence of the DAMA best-fit annual
modulation amplitude using Bayesian model comparison. This is a follow-up of a
previous similar study done in [1], who showed using frequentist model
comparison techniques that the annual modulation amplitudes starting from the
first ever phase of the DAMA experiment (starting in 1995) upto 2013 showed
hints of variability. We revisit this question using the latest DAMA data
released in 2018 and use Bayesian model comparison techniques for this
purpose.

The outline of this manuscript is as follows. Section 2 provides an
introduction to Bayesian model comparison techniques that will be used in this
work. In Section 3, we give a brief overview of the analysis and results
obtained for the time-dependence of the modulation amplitude in [1] as well as
our re-analysis of the same data. Section 4 contains the description of the
analysis and results for the model comparison done on the single-hit residual
rates as well as the mean per-cycle modulation amplitude data released by the
DAMA collaboration. Our conclusions are discussed in Section 5.

##  2 Introduction to Model Comparison Techniques

We provide a brief summary of Bayesian model comparison techniques used in
this work. This is the most robust method for model selection, surpassing both
frequentist or information theoretical techniques [29]. More details on
Bayesian model comparison can be found in various recent reviews [30, 31, 29,
32] (and references therein). We have previously applied Bayesian model
comparison techniques to a number of problems in Astrophysics and Cosmology,
including the analysis of data from dark matter experiments such as DAMA,
COSINE-100, and ANAIS-II [33, 34, 35, 36, 37, 38]. For more details, the
reader can refer to [30, 31, 29, 32] (and references therein). Previous
applications of Bayesian model comparison while analyzing the DAMA annual
modulation can be found in [39, 40].

The starting point of Bayesian Model comparison for a model $M$ with respect
to data $D$ is given by the Bayes Theorem. The model M is usually the
hypothesis which needs to be compared. The Bayes theorem is given by:

$P(M|D)=\frac{P(D|M)P(M)}{P(D)},$ (2.1)

where $P(M|D)$ is the posterior probability and $P(D|M)$ the marginal
likelihood, or sometimes known as the Bayesian Evidence, $P(M)$ is the prior
probability on the model $M$ and $P(D)$ is the probability for the data $D$,
which is usually never used in parameter estimation or model comparison. This
equation can be written as:

$P(D|M)\equiv\int P(D|M,\theta)P(\theta|M)\,d\theta$ (2.2)

where $\theta$ is the parameter vector for the model $M$ and $P(\theta|M)$
denotes the prior on $\theta$ for that model. In order to single out the
preferred model between two models $M_{1}$ and $M_{2}$, we calculate the odds
ratio, defined as the ratio of their posterior probabilities. This odds ratio
for the model $M_{2}$ over the model $M_{1}$ is as follows:

$O_{21}=\frac{P(M_{2}|D)}{P(M_{1}|D)}$ (2.3)

Using Eq. 2.1 and 2.2, this can be simplified as follows:

$O_{21}=\frac{E(M_{2})P(M_{2})}{E(M_{1})P(M_{1})}=B_{21}\frac{P(M_{2})}{P(M_{1})}$
(2.4)

where the term $B_{21}$ is the Bayes Factor, given by the ratio of Bayesian
Evidences for the two models. For equal prior probabilities for both the
models, the Odds Ratio is identical to Bayes Factor. Therefore, we obtain:

$O_{21}=B_{21}=\frac{\int
P(D|M_{2},\theta_{2})P(\theta_{2}|M_{2})\,d\theta_{2}}{\int
P(D|M_{1},\theta_{1})P(\theta_{1}|M_{1})\,d\theta_{1}}$ (2.5)

The Bayes factor is the quantity which is used for Bayesian model comparison.
We note that one key assumption in using the Bayes factor as a proxy for the
odds ratio for model comparison, is that all the models compared have equal
apriori probabilities. The model with the larger value for Bayesian evidence
will be considered as the preferred model. The Jeffrey’s scale is then used to
qualitatively assess its significance [32]. A Bayes Factor $<1$ indicates
negative support for the model for $M_{2}$. A value greater than 10 indicates
moderate evidence for $M_{2}$, whereas a value greater than 100 points to
decisive evidence. In this work, we calculate the Bayesian evidence for our
models using Nestle222http://kylebarbary.com/nestle/, which is a python
package based on the nested sampling algorithm [41].

##  3 Recap of Kelso18 and Re-analysis

In their paper, Kelso et al [1] (K18, hereafter) studied the best-fit
modulation amplitudes of the DAMA/NaI and DAMA/LIBRA Phase-I exposures. For
this purpose, they analyzed the data released by the DAMA collaboration over
the last 14 years starting in 1995 (covering both the DAMA/NaI and DAMA/LIBRA
phases) to explore the time dependence of the modulation amplitude using
frequentist model comparison methods [31]. This would help discriminate
between a dark matter versus an instrumental origin for the DAMA annual
modulation signal. They performed their analyses using the mean-cycle
amplitude data from the 2-4, 2-5, 2-6, 4-5, 4-6, and 5-6 keV energy ranges.

The data used for the 2-4, 2-5, and 2-6 keV energy ranges was obtained from
Figure 3 of Ref. [19], provided by the DAMA collaboration. The mean time for
each cycle is taken relative to January 1, 1995. The mean-cycle data for the
higher energy ranges (4-5, 4-6, and 5-6 keV) were calculated from the
efficiency-weighted average of the mean-cycle amplitudes in the 2-4, 2-5 and
2-6 keV ranges. The main premise behind this is that the data from different
energy ranges modulate with the same period and amplitude, but are
statistically independent. Although this calculation is strictly correct only
when the phase and period are fixed, K18 argued that this approximation is
also valid when the phase and period vary, since the difference in the
modulation amplitude between the fixed versus free period and phase case
differs by less than 0.5$\sigma$.

The best-fit modulation amplitudes for the DAMA/LIBRA Phase-I and DAMA/NaI
exposures are calculated, along with the combined modulation amplitude for
both the experiments. From the results listed in Table-I of K18, one can see
by eye that the amplitudes decrease in every energy range, hinting towards a
decreasing trend in the amplitudes for the two phases of the DAMA experiment.
By using multiple non-parametric statistical tests, K18 concluded that the
null hypothesis is rejected at 2.4$\sigma$ in the 2-6 keV and $p-$values of
21% and 5% in the 2-4 and 2-5 keV energy bins. Therefore, they deduce that one
possible explanation for the results of these tests is that the modulation
amplitude varies with time.

To further quantify this, five parametric models for the time dependence of
the modulation amplitude were analyzed in K18: constant, broken constant with
separate values for DAMA/NaI data and DAMA/LIBRA Phase-I data, linearly time-
dependent, exponentially time-dependent, and a broken exponential model with
independent normalizations for NaI and LIBRA phases. Results from the goodness
of fit test showed that all these models provided reasonable fits for data
from all the energy ranges with the exception of the 5-6 keV range.

The constant time-independent model was considered as the null hypothesis for
model comparison. It was observed that all the energy ranges, except the 2-4
keV range, show a preference for time-dependence at the 2-3$\sigma$
significance level. The data in the 5-6 keV range show three departures at
$>2.8\sigma$ from the mean value and do not fit any modulation hypothesis.
Therefore, K18 point out that there could be a problem for the data in 5-6 keV
energy range. They conclude that all the time-varying models in the 2-6 keV
range are preferred over a constant amplitude model with a significance
between 2.3$\sigma$ to 2.8$\sigma$. However, a strong preference between the
various time-dependent models could not be discerned using the data available
until then.

###  3.1 Re-analysis of K18 using Bayesian model selection

As described above, K18 used frequentist methods to carry out model comparison
between the time-independent constant model and various time-dependent models
using the per-cycle modulation amplitude data from DMA/NaI and the first phase
of DAMA/LIBRA experiments. Before we move onto model comparison for the
updated data from the second phase of the DAMA/LIBRA experiment, we redo the
analysis using the aforementioned data analyzed in K18 using Bayesian model
comparison for the linear and exponential models. The reason is that these
models showed the largest significance when compared to a constant modulation
amplitude and are representative of the simplest time-varying models. However,
we also report the $\chi^{2}$/DOF for each of the models as a sanity check of
the goodness of fit. A $\chi^{2}$ test is a statistical test used to compare
the observed results of an experiment with the expected results. The
$\chi^{2}$ statistic is defined as

$\chi^{2}=\sum\left(\frac{O_{i}-E_{i}}{\sigma_{i}}\right)^{2}$

where $O_{i}$ and $E_{i}$ are the observed and expected values (for a
particular model) respectively, and $\sigma_{i}$ indicates the error in the
observed values. For all the results in this manuscript, the $\chi^{2}$ values
are obtained from the comparison of the best-fit model to the observed data,
while the $DOF$ corresponds to the number of degrees of freedom, calculated as
the difference in number of data points and the number of free parameters.

We use the data given in Tables 3 and 5 of K18 for the mean per-cycle
modulation amplitudes and mean cycle times. The mean time for each cycle is
taken relative to January 1, 1995. The constant time-independent, linearly
time-dependent and exponentially time-independent models for $S_{m}$ are the
same as those in K18 and are given below:

Constant: $\displaystyle S_{m}$ $\displaystyle=A_{0}$ (3.1) Linear:
$\displaystyle S_{m}$ $\displaystyle=mt+c$ (3.2) Exponential: $\displaystyle
S_{m}$ $\displaystyle=Ae^{-Bt}$ (3.3)

where $t$ denotes time, $A_{0}$, $m$, $c$, $A$ and $B$ are variable parameters
used in the analysis. Equations 3.1, 3.2 and 3.3 describe the three models
used in our analysis. For evaluating the Bayes factors, we choose uniform
priors for the parameters of these models, which are listed in Column-I of
Table 1. The constant model is chosen as the null hypothesis.

| Priors | Priors
---|---|---
Parameters | Per-cycle Modulation Amplitude analysis | Residual Rates Analysis
$A_{0}$ | [-0.05, 0.05] | [-0.05, 0.05]
$m$ | [-0.05, 0.05] | [-0.05, 0.05]
$c$ | [-0.05, 0.05] | [-0.05, 0.05]
$A$ | [-0.05, 0.05] | [-0.05, 0.05]
$B$ | [-0.1 , 0.1] | [-0.01, 0.01]
$\omega$ |  | [0.0015, 0.1281]
$t0$ |  | [0,$2\pi/\omega$]
Table 1: The prior ranges chosen for the parameters used for all the analyses
in this work.

Our results from Bayesian model selection are shown in Table 2. We find that
none of the datasets show a decisive evidence (Bayes factor $>100$) for the
linear or exponential model. The largest Bayes factors which we obtained were
for the exponential model in the 2-5 and 2-6 keV energy ranges of around 11
and 20, respectively. These point to “moderate” evidence for these models
according to Jeffreys’ scale [32]. On the other hand, the Bayes factors for
these models in the 4-5 keV and 4-6 keV energy range are still small, implying
that there is no evidence for time variability in these energy ranges for
these models.

Data | Models | $\chi^{2}/DOF$ | Bayes Factors
---|---|---|---
2-4 keV | Constant | 9.4/13 | -
| Linear | 6.9/12 | 0.05
| Exponential | 7.2/12 | 1.50
2-5 keV | Constant | 13.6/13 | -
| Linear | 6.4/12 | 0.36
| Exponential | 7.5/12 | 11.71
2-6 keV | Constant | 10.7/13 | -
| Linear | 4.3/12 | 0.23
| Exponential | 3.5/12 | 20.44
4-6 keV | Constant | 22.4/13 | -
| Linear | 19.0/12 | 0.05
| Exponential | 15.2/12 | 5.49
4-5 keV | Constant | 19.7/13 | -
| Linear | 13.2/12 | 0.42
| Exponential | 15.0/12 | 4.30
5-6 keV | Constant | 40.0/13 | -
| Linear | 39.8/12 | 0.02
| Exponential | 35.1/12 | 1.28
Table 2: Results for the Bayesian model comparison analysis on the data from
Tables 3 and 5 of K18. The mean time for each cycle is taken relative to
January 1, 1995. The Bayes Factors are calculated for the linear and
exponential models, while considering the constant model as the null
hypothesis, and are listed in Column-IV of this table. We find that 2-5 keV
and 2-6 keV energy intervals show moderate evidence for exponential variation
of the amplitude. All other datasets are consistent with no variability. We
also list the $\chi^{2}$/dof for each of the models/energy ranges as a sanity
check of the goodness of fit.

##  4 Analysis and Results

We now present our results on the search for possible time-dependence of the
modulation amplitudes using the latest DAMA data released in 2018. We do this
analysis using two independent datasets, as further described in sub-sections
4.1 and 4.2, below.

###  4.1 Analysis of single-hit residual rates from DAMA/LIBRA Phase-I and
Phase-II

This sub-section entails the analysis of the data for the single-hit residual
rates from the DAMA/LIBRA Phase-I and Phase-II experiments, obtained from [20,
17, 19]. Data from the 2-4 keV and 2-5 keV energy ranges from the Phase-I of
the DAMA/LIBRA experiment, 1-3 keV and 1-6 keV energy ranges in the Phase-II
of the DAMA/LIBRA experiment, and 2-6 keV energy range from both Phase-I and
Phase-II of the DAMA/LIBRA experiment are analyzed. We used the digitized DAMA
data for these energy intervals from our previous work [33]. In the 2-6 keV
energy range, we also combined this data with that from DAMA/NaI experiment
(which was digitized from Ref [17]) and perform the same analysis. The total
exposure for DAMA/NaI, the first phase of the DAMA experiment that ran for
seven years, was $1.08\times 10^{5}$ kg-days (0.29 ton-year). The exposure
then increased over the the first phase of DAMA/LIBRA, the upgraded DAMA
experiment, which registered a total exposure of $3.80\times 10^{5}$ kg-days
(1.04 ton-year) over the course of another seven years. The exposure further
increased in the second phase of the DAMA/LIBRA experiment, registering a
total exposure of $4.11\times 10^{5}$ kg-days (1.13 ton-year). The total
cumulative exposure, spanning DAMA/NaI and DAMA/LIBRA Phase-I and Phase-II, is
2.46 ton-year [20].

The models for the single-hit residual rate term, which we posit for our
analysis are:

Constant $S_{m}$: $\displaystyle f(t)$ $\displaystyle=A_{0}\cos{\omega(t-t0)}$
(4.1) Linear $S_{m}$: $\displaystyle f(t)$
$\displaystyle=(mt+c)\cos{\omega(t-t0)}$ (4.2) Exponential $S_{m}$:
$\displaystyle f(t)$ $\displaystyle=Ae^{-Bt}\cos{\omega(t-t0)}$ (4.3)

where $t$ corresponds to time, $\omega$ and $t_{0}$ are the angular frequency
and phase of the modulation, and $A_{0}$, $m$, $c$, $A$, $B$ represent the
parameters used to model the time-variation of the amplitude $S_{m}$.
Therefore, Eq. 4.2 and Eq. 4.3 contain additional time-dependence on top of
the well known cosine modulation. We chose uniform priors on all these
parameters, which are listed in Column-II of Table 1. The maximum value for
the prior on $\omega$ is taken as $2\pi$ over the median difference in time
data-points. Since the maximum period for which we can detect a potential
modulation is half the total duration of the dataset, we used that to choose
the minimum value for the prior on $\omega$. The ranges for the priors on the
other parameters have been chosen based on the maximum and minimum values for
the residual count rates.

The best-fit plots for the constant, linearly dependent and exponentially
dependent models for $S_{m}$ are shown in Figure 1. The best-fit values used
in the plot are calculated using the median of the samples and the curves are
fitted against the data released by the DAMA collaboration [20, 17, 19].

Data | Models | $\chi^{2}/DOF$ | Bayes Factors
---|---|---|---
DAMA/LIBRA Phase-I: | Constant | $47.0/47$ | -
2-4 keV | Linear | $44.1/46$ | $2.69\times 10^{-4}$
| Exponential | $43.4/46$ | $4.77\times 10^{-2}$
DAMA/LIBRA Phase-I: | Constant | $38.3/47$ | -
2-5 keV | Linear | $34.6/46$ | $2.87\times 10^{-4}$
| Exponential | $34.6/46$ | $9.98\times 10^{-3}$
DAMA/LIBRA Phase-II: | Constant | $52.3/49$ | -
1-3 keV | Linear | $51.7/48$ | $0.68\times 10^{-4}$
| Exponential | $51.6/48$ | $4.39\times 10^{-3}$
DAMA/LIBRA Phase-II: | Constant | $47.7/49$ | -
1-6 keV | Linear | $49.2/48$ | $0.58\times 10^{-4}$
| Exponential | $47.0/48$ | $5.59\times 10^{-3}$
DAMA/LIBRA Phases I and II: | Constant | $69.5/99$ | -
2-6 keV | Linear | $69.5/98$ | $0.19\times 10^{-4}$
| Exponential | $69.4/98$ | $1.92\times 10^{-6}$
DAMA/NaI, DAMA/LIBRA Phases I and II: | Constant | $154.5/136$ | -
2-6 keV | Linear | $147.6/135$ | $1.49\times 10^{-4}$
| Exponential | $149.8/135$ | $9.60\times 10^{-5}$
Table 3: Summary of model comparison results for different models of the
single-hit residual rates as discussed in Section 4.1. The $\chi^{2}$ values
are calculated from the $\chi^{2}$ statistic by comparing the best-fit model
to the data, and $DOF$ is taken as the difference in number of data points and
the number of parameters. The Bayes factors for the linear (Equation 4.2) and
exponentially (Equation 4.3) varying modulation amplitudes have been
calculated by taking the constant amplitude model (Equation 4.1) as the null
hypothesis. Figure 1: Plots of the best-fit curves of the constant, linearly
dependent and exponentially dependent annual modulation amplitudes (on top of
the expected cosine modulaton) from the analysis of residual rates described
in Section 4.1. The best-fits are calculated from the median of the posterior
samples obtained from Nestle used for Bayesian model comparison and the curves
are plotted against the data released by the DAMA collaboration [20, 17, 19].

The results for the Bayesian model comparison are summarized in Table 3 for
the different energy ranges. In addition to the Bayes factors, we also list
the $\chi^{2}$/DOF values. We find that both the constant amplitude models as
well all the time-varying amplitudes have $\chi^{2}$/dof close to 1. However,
the Bayes factor for all the time-varying amplitudes are less than or close to
one for all phases of the DAMA experiment. This shows that there is no
evidence for linear or exponential time variation of the modulation amplitude,
when considering the residual count rates.

###  4.2 Analysis of annual mean per-cycle modulation amplitude data

This section discusses the Bayesian analysis of mean per-cycle data for the
modulation amplitude from Fig 14 of [20]. This is an update of the analysis
done in K18 as discussed in Section 3. Although only the data from the 2-4,
2-5, and 2-6 keV energy ranges spanning the DAMA/NaI and DAMA/LIBRA Phase-I
have been released prior to publication of K18, they synthesized the data for
the 4-5, 4-6, and 5-6 keV energy ranges for a fixed period and phase. However,
in the latest data release for the DAMA/LIBRA experiment, the per-cycle
modulation amplitudes for the 4-5 keV and 5-6 keV are provided in Figure 14 of
[20]. We use this data (as given in Table 4) for our analysis.

While K18 analyzed the per-cycle data from DAMA/NaI and the first phase of the
DAMA/LIBRA experiment, the data analyzed here includes the per-cycle
modulation amplitude data from the first and second phase of the DAMA/LIBRA
experiment that began in September 2003. We analyze data from the 1-2, 2-3,
3-4, 4-5 and 5-6 keV energy ranges. Values for the per-cycle amplitudes in the
2-6 and 4-6 keV ranges have not been provided by the DAMA collaboration in
their 2018 data release, and hence have not been analyzed for the DAMA/LIBRA
phase. We also analyze the per-cycle modulation amplitudes for the 1-2 keV
energy range, which have been released only for DAMA/LIBRA Phase-II and are
not available for the DAMA/NaI phase or DAMA/LIBRA Phase-I.

DAMA/LIBRA | $t$ |  $S_{m}$ (1-2 keV) |  $S_{m}$ (2-3 keV) |  $S_{m}$ (3-4 keV) |  $S_{m}$ (4-5 keV) |  $S_{m}$ (5-6 keV)
---|---|---|---|---|---|---
| ($year$) | ($/kg/day/keV$) | ($/kg/day/keV$) | ($/kg/day/keV$) | ($/kg/day/keV$) | ($/kg/day/keV$)
Phase-I |  |  |  |  |  |
Cycle 1 | 1 | - | 0.0184 $\pm$ 0.0088 | 0.0361 $\pm$ 0.0082 | 0.0053 $\pm$ 0.0067 | -0.0049 $\pm$ 0.0067
Cycle 2 | 2 | - | 0.0137 $\pm$ 0.0101 | 0.0180 $\pm$ 0.0096 | 0.0093 $\pm$ 0.0071 | 0.0004 $\pm$ 0.0071
Cycle 3 | 3 | - | 0.0203 $\pm$ 0.0086 | 0.0073 $\pm$ 0.0089 | 0.0191 $\pm$ 0.0069 | -0.0003 $\pm$ 0.0069
Cycle 4 | 4 | - | 0.0284 $\pm$ 0.0079 | 0.0052 $\pm$ 0.0080 | 0.0091 $\pm$ 0.0063 | 0.0019 $\pm$ 0.0063
Cycle 5 | 5 | - | 0.0179 $\pm$ 0.0071 | 0.0033 $\pm$ 0.0071 | 0.0051 $\pm$ 0.0058 | 0.0076 $\pm$ 0.0058
Cycle 6 | 6 | - | 0.0195 $\pm$ 0.0073 | 0.0099 $\pm$ 0.0069 | -0.0009 $\pm$ 0.0055 | 0.0139 $\pm$ 0.0055
Cycle 7 | 7 | - | 0.0224 $\pm$ 0.0067 | 0.0118 $\pm$ 0.0068 | 0.0022 $\pm$ 0.0050 | 0.0035 $\pm$ 0.0050
Phase-II |  |  |  |  |  |
Cycle 1 | 9 | 0.0323 $\pm$ 0.0074 | 0.0202 $\pm$ 0.0064 | 0.0212 $\pm$ 0.0057 | 0.0012 $\pm$ 0.0046 | 0.0076 $\pm$ 0.0046
Cycle 2 | 10 | 0.0147 $\pm$ 0.0081 | 0.0120 $\pm$ 0.0075 | 0.0162 $\pm$ 0.0064 | 0.0075 $\pm$ 0.0050 | 0.0079 $\pm$ 0.0050
Cycle 3 | 11 | 0.0232 $\pm$ 0.0080 | 0.0161 $\pm$ 0.0074 | 0.0076 $\pm$ 0.0062 | 0.0107 $\pm$ 0.0048 | 0.0034 $\pm$ 0.0048
Cycle 4 | 12 | 0.0356 $\pm$ 0.0081 | 0.0199 $\pm$ 0.0073 | 0.0108 $\pm$ 0.0061 | 0.0093 $\pm$ 0.0049 | 0.0026 $\pm$ 0.0049
Cycle 5 | 13 | 0.0095 $\pm$ 0.0084 | 0.0139 $\pm$ 0.0072 | 0.0121 $\pm$ 0.0063 | 0.0033 $\pm$ 0.0047 | 0.0069 $\pm$ 0.0047
Cycle 6 | 14 | 0.0049 $\pm$ 0.0077 | 0.0108 $\pm$ 0.0073 | 0.0079 $\pm$ 0.0061 | 0.0109 $\pm$ 0.0045 | 0.0054 $\pm$ 0.0045

Table 4: The per-cycle modulation amplitudes as determined by the DAMA
collaboration, taken from Figure 36 of Ref [42]. The time axis corresponds to
each annual cycle spanning over 14 years of the DAMA/LIBRA experiment, since
the beginning of the DAMA/LIBRA Phase-I in September, 2003.

The constant time-independent, linearly time-dependent and exponentially time-
independent models for $S_{m}$ used in this analysis are given in Equations
3.1, 3.2 and 3.3. The Bayesian priors used for the constant ($A_{0}$), linear
($m$, $c$), and exponential ($A$, $B$) model parameters remain the same as
before and shown in Table 1. The best-fit plots of this analysis are given in
Figure 2. The best-ft values obtained from the posterior samples provided by
Nestle are shown in Table 5. The $\chi^{2}$/DOF values and the Bayes factors
have been listed in Columns III, and IV of Table 6 respectively for the three
models of $S_{m}$. Unlike K18, we find that $\chi^{2}$/dof estimates is close
to 1 for all the time-varying models. Furthermore, similar to our results for
the residual count rates, we find that the Bayes Factors are close to 1.
Therefore, we do not find any evidence for a linear or exponential dependence
of the (DAMA provided) modulation amplitudes starting from the 2003 phase of
DAMA/LIBRA.

|  | 1-2 keV | 2-3 keV | 3-4 keV | 4-5 keV | 5-6 keV
---|---|---|---|---|---|---
Constant |  $A_{0}$ $[10^{-2}dru]$ | $2.04\pm 0.33$ | $1.80\pm 0.21$ | $1.26\pm 0.19$ | $0.70\pm 0.15$ | $0.61\pm 0.15$
Linear |  $m$ $[10^{-4}dru/year]$ | $15.92\pm 9.08$ | $-5.78\pm 5.38$ | $-5.35\pm 4.87$ | $0.93\pm 3.72$ | $4.38\pm 3.74$
|  $c$ $[10^{-2}dru]$ | $3.79\pm 1.01$ | $2.27\pm 0.47$ | $1.71\pm 0.45$ | $0.63\pm 0.35$ | $0.23\pm 0.35$
Exponential |  $A$ $[10^{-2}dru]$ | $3.62\pm 0.98$ | $2.36\pm 0.52$ | $1.88\pm 0.49$ | $0.76\pm 0.32$ | $0.52\pm 0.23$
|  $B$ $[10^{-2}/year]$ | $5.11\pm 3.01$ | $3.49\pm 2.79$ | $4.92\pm 3.23$ | $0.62\pm 4.79$ | $-1.89\pm 4.47$

Table 5: The best-fit values for the per-cycle modulation amplitudes as given by Equations 3.1, 3.2, and 3.3. The analysis and comparison of these models is described in Section 4.2. $dru$ corresponds to Differential Rate Unit, which is equal to 1 $[/kg/day/keV]$. Data | Models | $\chi^{2}/DOF$ | Bayes Factors
---|---|---|---
DAMA/LIBRA Phase-II: | Constant | $12.5/5$ | -
1-2 keV | Linear | $7.6/4$ | $0.11$
| Exponential | $8.2/4$ | $1.83$
DAMA/LIBRA Phase-I and Phase-II: | Constant | $4.6/12$ | -
2-3 keV | Linear | $3.4/11$ | $0.02$
| Exponential | $3.5/11$ | $0.78$
DAMA/LIBRA Phase-I and Phase-II: | Constant | $15.6/12$ | -
3-4 keV | Linear | $14.4/11$ | $0.02$
| Exponential | $14.2/11$ | $1.34$
DAMA/LIBRA Phase-I and Phase-II: | Constant | $9.7/12$ | -
4-5 keV | Linear | $9.7/11$ | $0.01$
| Exponential | $9.7/11$ | $0.69$
DAMA/LIBRA Phase-I and Phase-II: | Constant | $8.2/12$ | -
5-6 keV | Linear | $6.8/11$ | $0.02$
| Exponential | $7.2/11$ | $0.66$
Table 6: Summary of model comparison results for different models of the per-
cycle modulation amplitudes $S_{m}$ discussed in Section 4.2. The $\chi^{2}$
values are calculated from the $\chi^{2}$ statistic by comparing the best-fit
model to the data, and $DOF$ is taken as the difference in number of data
points and the number of parameters. The Bayes Factors for the linearly and
exponentially varying modulation amplitudes have been calculated using data
released by the DAMA collaboration for the first and second phase of the
DAMA/LIBRA experiment [20]. The constant amplitude model is taken as the null
hypothesis. We do not find any evidence for time modulation of the amplitudes
from this analysis. Figure 2: The best-fit curves for the constant, linear and
exponential time-dependence of the modulation amplitude are plotted over the
per-cycle mean amplitude values. The data for the 1-2 keV energy range has
been released only for DAMA/LIBRA Phase-II. The x-axis represents each annual
cycle spanning over the 14 years of the DAMA/LIBRA experiment, since the
beginning of the DAMA/LIBRA Phase-I in September, 2003 [42].

##  5 Conclusions

Around three years ago, K18 analyzed the DAMA annual modulation data collected
over 14 cycles (starting from the DAMA/NaI phase in 1995) in order to discern
a possible time dependence in the modulation amplitude. Their results indicate
that all time varying amplitudes are preferred over a constant amplitude by
upto 2-3$\sigma$ in almost all the energy ranges they analyzed. If the DAMA
annual modulation is caused by astrophysical dark matter, one would not expect
the modulation amplitude to change.

Motivated by these considerations, we revisit this issue using the latest DAMA
data. We first re-analyze the modulation amplitudes used in K18 by carrying
out Bayesian model comparison of the linear and exponential time-varying
amplitudes versus the constant amplitude. Our results for this re-analysis can
be found in Table 2. We find that although none of the models show a decisive
evidence for any linear or exponential variation, the data in 2-5 keV and 2-6
keV energy intervals show Bayes factors of around 11 and 20, respectively,
pointing to moderate evidence for a time-variation.

We then update this analysis using both the cumulative residual event rates as
well as the updated per cycle modulation amplitude provided by the DAMA
collaboration in their latest paper [20] starting from the DAMA/LIBRA phase in
2003. The analysis of residual event rate was done by an augmented fit,
consisting of a varying modulation amplitude in addition to the cosine term.

Our results for the Bayesian model comparison using the single-hit residual
rates and the per-cycle modulation amplitudes are summarized in Table 3 and
Table 6, respectively. Using both the datasets, we find that the Bayes factors
for all the time-varying amplitudes with respect to the constant amplitude are
less than or close to one. Therefore, we do not find any evidence for a
decreasing modulation amplitude as hinted from the K18 results for the
DAMA/LIBRA phase of the DAMA experiment.

To promote transparency in data analysis and enable the reader to reproduce
these results, we have uploaded all our codes datasets used in this analysis
at https://github.com/srinixbhagvati/Time-dependence-of-DAMA-LIBRA-modulation-
amplitude

## References

  * [1] C. Kelso, C. Savage, P. Sandick, K. Freese and P. Gondolo, _Examining the time dependence of DAMA’s modulation amplitude_ , _European Physical Journal C_ 78 (2018) 223 [1710.03770].
  * [2] G. Bertone and D. Hooper, _History of dark matter_ , _Reviews of Modern Physics_ 90 (2018) 045002 [1605.04909].
  * [3] G. Jungman, M. Kamionkowski and K. Griest, _Supersymmetric dark matter_ , _Physics Reports_ 267 (1996) 195 [hep-ph/9506380].
  * [4] G. Bertone, D. Hooper and J. Silk, _Particle dark matter: evidence, candidates and constraints_ , _Physics Reports_ 405 (2005) 279 [hep-ph/0404175].
  * [5] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. & Astrophys._ 641 (2020) A6 [1807.06209].
  * [6] L. Roszkowski, _Particle dark matter - A theorist’s perspective_ , _Pramana_ 62 (2004) 389 [hep-ph/0404052].
  * [7] J. L. Feng, _Dark Matter Candidates from Particle Physics and Methods of Detection_ , _Ann. Rev. Astron. Astrophys._ 48 (2010) 495 [1003.0904].
  * [8] B. W. Lee and S. Weinberg, _Cosmological lower bound on heavy-neutrino masses_ , _Physical Review Letters_ 39 (1977) 165.
  * [9] H. Goldberg, _Constraint on the Photino Mass from Cosmology_ , _Physical Review Letters_ 50 (1983) 1419.
  * [10] G. Steigman, B. Dasgupta and J. F. Beacom, _Precise relic WIMP abundance and its impact on searches for dark matter annihilation_ , _Physical Review D._ 86 (2012) 023506 [1204.3622].
  * [11] A. Drukier and L. Stodolsky, _Principles and applications of a neutral-current detector for neutrino physics and astronomy_ , _Physical Review D._ 30 (1984) 2295.
  * [12] M. W. Goodman and E. Witten, _Detectability of certain dark-matter candidates_ , _Physical Review D._ 31 (1985) 3059.
  * [13] D. Bauer, J. Buckley, M. Cahill-Rowley, R. Cotta, A. Drlica-Wagner, J. L. Feng et al., _Dark matter in the coming decade: Complementary paths to discovery and beyond_ , _Physics of the Dark Universe_ 7 (2015) 16 [1305.1605].
  * [14] A. K. Drukier, K. Freese and D. N. Spergel, _Detecting cold dark-matter candidates_ , _Physical Review D._ 33 (1986) 3495.
  * [15] K. Freese, J. Frieman and A. Gould, _Signal modulation in cold-dark-matter detection_ , _Physical Review D._ 37 (1988) 3388.
  * [16] R. Bernabei, P. Belli, F. Cappella, R. Cerulli, F. Montechia, F. Nozzoli et al., _Dark Matter search_ , _Nuovo Cimento Rivista Serie_ 26 (2003) 1 [astro-ph/0307403].
  * [17] R. Bernabei, P. Belli, F. Cappella, R. Cerulli, C. J. Dai, A. D’Angelo et al., _First results from DAMA/LIBRA and the combined results with DAMA/NaI_ , _European Physical Journal C_ 56 (2008) 333 [0804.2741].
  * [18] R. Bernabei, P. Belli, F. Cappella, R. Cerulli, C. J. Dai, A. D’Angelo et al., _New results from DAMA/LIBRA_ , _European Physical Journal C_ 67 (2010) 39 [1002.1028].
  * [19] R. Bernabei et al., _Final model independent result of DAMA/LIBRA-phase1_ , _Eur. Phys. J._ C73 (2013) 2648 [1308.5109].
  * [20] R. Bernabei, P. Belli, A. Bussolotti, F. Cappella, V. Caracciolo, R. Cerulli et al., _First Model Independent Results from DAMA/LIBRA-Phase2_ , _Nuclear Physics and Atomic Energy_ 19 (2018) 307 [1805.10486].
  * [21] R. Bernabei, P. Belli, R. Cerulli, F. Montecchia, M. Amato, G. Ignesti et al., _Dark Matter Search with the DAMA experiments_ , _Nuclear Physics B Proceedings Supplements_ 87 (2000) 67.
  * [22] LUX collaboration, _Results from a search for dark matter in the complete LUX exposure_ , _Phys. Rev. Lett._ 118 (2017) 021303 [1608.07648].
  * [23] PandaX-II collaboration, _Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment_ , _Phys. Rev. Lett._ 119 (2017) 181302 [1708.06917].
  * [24] XENON collaboration, _Dark Matter Search Results from a One Ton-Year Exposure of XENON1T_ , _Phys. Rev. Lett._ 121 (2018) 111302 [1805.12562].
  * [25] CDEX collaboration, _Limits on Light Weakly Interacting Massive Particles from the First 102.8 kg ${\times}$ day Data of the CDEX-10 Experiment_, _Phys. Rev. Lett._ 120 (2018) 241301 [1802.09016].
  * [26] G. Adhikari, E. Barbosa de Souza, N. Carlin, J. J. Choi, S. Choi, M. Djamal et al., _Strong constraints from COSINE-100 on the DAMA dark matter results using the same sodium iodide target_ , _arXiv e-prints_ (2021) arXiv:2104.03537 [2104.03537].
  * [27] J. Amare, S. Cebrian, D. Cintas, I. Coarasa, E. Garcia, M. Martinez et al., _Annual Modulation Results from Three Years Exposure of ANAIS-112_ , _arXiv e-prints_ (2021) arXiv:2103.01175 [2103.01175].
  * [28] Super-Kamiokande collaboration, _Search for dark matter WIMPs using upward through-going muons in Super-Kamiokande_ , _Phys. Rev. D_ 70 (2004) 083523 [hep-ex/0404025].
  * [29] S. Sharma, _Markov Chain Monte Carlo Methods for Bayesian Data Analysis in Astronomy_ , _Ann. Rev. Astron. Astrophys._ 55 (2017) 213 [1706.01629].
  * [30] Ž. Ivezić, A. Connolly, J. Vanderplas and A. Gray, _Statistics, Data Mining and Machine Learning in Astronomy_. Princeton University Press, 2014.
  * [31] M. Kerscher and J. Weller, _On model selection in cosmology_ , _SciPost Physics Lecture Notes_ 9 (2019)  [1901.07726].
  * [32] R. Trotta, _Bayesian Methods in Cosmology_ , _ArXiv e-prints_ (2017) [1701.01467].
  * [33] A. Krishak, A. Dantuluri and S. Desai, _Robust model comparison tests of DAMA/LIBRA annual modulation_ , _JCAP_ 2020 (2020) 007 [1906.05726].
  * [34] A. Krishak and S. Desai, _An independent assessment of significance of annual modulation in COSINE-100 data_ , _Open J. Astrophys._ (2019)  [1907.07199].
  * [35] A. Krishak and S. Desai, _An independent search for annual modulation and its significance in ANAIS-112 data_ , _Progress of Theoretical and Experimental Physics_ 2020 (2020) 093F01 [1910.05096].
  * [36] A. Krishak and S. Desai, _Model comparison tests of modified gravity from the Eöt-Wash experiment_ , _JCAP_ 2020 (2020) 006 [2003.10127].
  * [37] H. Singirikonda and S. Desai, _Model comparison of $\Lambda$ CDM vs Rh=c t using cosmic chronometers_, _European Physical Journal C_ 80 (2020) 694 [2003.00494].
  * [38] R. Agrawal, H. Singirikonda and S. Desai, _Search for Lorentz Invariance Violation from stacked Gamma-Ray Burst spectral lag data_ , _JCAP_ 2021 (2021) 029 [2102.11248].
  * [39] A. Messina, M. Nardecchia and S. Piacentini, _Annual modulations from secular variations: not relaxing DAMA?_ , _JCAP_ 2020 (2020) 037 [2003.03340].
  * [40] G. Grilli di Cortona, A. Messina and S. Piacentini, _Migdal effect and photon Bremsstrahlung: improving the sensitivity to light dark matter of liquid argon experiments_ , _JHEP_ 11 (2020) 034 [2006.02453].
  * [41] J. Buchner, _Nested Sampling Methods_ , _arXiv e-prints_ (2021) arXiv:2101.09675 [2101.09675].
  * [42] R. Bernabei, P. Belli, A. Bussolotti, F. Cappella, V. Caracciolo, R. Cerulli et al., _The DAMA project: Achievements, implications and perspectives_ , _Progress in Particle and Nuclear Physics_ 114 (2020) 103810.