1 Introduction↩︎

It has been over two decades since the DAMA/LIBRA NaI experiment claimed the discovery of an annual modulation signal of dark matter (DM) at 6.3\(\sigma\) [@Bernabei:2003za]. In combination with DAMA/LIBRA and DAMA/LIBRA phase2, 2.86 ton-yr of exposure over 22 annual cycles have now yielded a claimed signal over 13\(\sigma\) [@Bernabei:2022ath]. Multiple analyses have shown that the energy dependence of this signal leads to a preference for a particle with a mass around 10 or 50 GeV, depending on whether it scatters primarily off sodium or iodine (e.g. [@Kelso:2013gda; @Baum:2018ekm]). In parallel, however, advances in large direct detection experiments such as SuperCDMS, PICO-60, LZ and XENON1T have left very little room for a DM interpretation, with the current sensitivity being a full six orders of magnitude below DAMA/LIBRA’s region of interest [@PhysRevLett.120.061802; @PhysRevLett.121.051301; @PhysRevLett.121.081307; @PhysRevD.100.022001; @PhysRevLett.131.041002; @PhysRevLett.131.041003]. Nonetheless, the DAMA/LIBRA DM interpretation has persisted on the premise that no experiment to date has formally excluded it by searching for an annual modulation signature with the same technology, namely extremely pure NaI scintillator crystals, leaving the door open to DM particles with some inscrutable affinity towards sodium or iodine. While some exploratory efforts occurred starting from the early 2000s [@Ahmed:2003su; @Kim:2012rza; @DM-Ice:2016snk], none of these experiments were able to completely exclude DAMA/LIBRA in a model-agnostic way.

In recent years, a new generation of NaI experiments has been releasing data: the ANAIS-112 detector [@Amare:2019jul; @Amare:2021yyu; @Coarasa:2024xec; @Amare:2025dfq], located in Canfranc and COSINE-100 [@Adhikari:2018ljm; @COSINE-100:2021xqn; @COSINE-100:2019lgn; @COSINE-100:2021zqh; @Carlin:2024maf], at Yangyang. Neither of these experiments alone has quite reached the 5\(\sigma\) sensitivity to conclusively rule out the NaI-philic DM hypothesis of DAMA/LIBRA. A recent joint analysis [@COSINE-100:2025eyc] combining these datasets into two single-energy bins has excluded a modulation signal compatible with DAMA’s at the level of 4.68\(\sigma\) (3.53\(\sigma\)) in the 1-6 (2-6) keV range.

In addition, the two experiments have separately determined modulation amplitudes for events binned in energy [@Carlin:2024maf; @Amare:2025dfq], assuming the expected phase for DM scattering, a useful result for model discrimination. These can be directly compared to similar results previously reported by DAMA [@Bernabei:2021kdo]. Performing such a comparison will be the goal of this letter. We use the full spectral information reported by the different experiments and determine how well we can reconcile the spectra reported by DAMA/LIBRA and COSINE-100/ANAIS-112 if we assume the modulation is due to DM-nucleon scattering. This comparison is not as trivial as it may seem: given a common nuclear recoil spectrum, the differences in binning and energy resolution can lead to substantially different predictions for the measured modulation spectrum between the various experiments. To remain as general as possible, we adopt two different model-independent approaches: a DM effective field theory in which the nuclear recoil spectrum is calculable and a more general scenario where the nuclear recoil spectrum is parameterized with various functional forms. We undertake fits to the experimental data sets to determine the best-fit parameters of these models, and based on these results we quantify the tension between the experiments.

The scattering rate of weakly-interacting DM particles from a nuclear target is given in terms of the differential cross section \(d\sigma/dE_R\) by \[\label{eq:recoil} \frac{dR}{dE_R}=\frac{\rho_0}{m_\chi m_N}\int^{v_{\text{max}}}_{v_{\text{min}}}d^3v v \tilde{f}(\mathbf{v},t)\frac{d\sigma}{dE_R}\tag{1}\] where \(\rho_0\) is the local DM mass density, \(m_\chi\) is the DM particle mass, \(m_N\) is the mass of the recoiling nucleus, \(v_{\max}\) is the escape velocity and \(v_{\min}\) is the minimum velocity needed to cause a nucleus to recoil with energy \(E_R\). \(\tilde{f}(\mathbf{v},t)\) is the DM velocity distribution in the lab frame. The lab-frame distribution is related to the Galactic frame distribution \(f(\mathbf{v})\) by a Galilean boost: \(\tilde{f}(\mathbf{v})=f(\mathbf{v}+\mathbf{v}_{\text{obs}}(t))\), where \(\mathbf{v}_{\text{obs}}(t)=\mathbf{v}_{\odot}+\mathbf{V}_\oplus(t)\), \(\mathbf{v}_{\odot}\) is the velocity of the Sun relative to the DM reference frame and \(\mathbf{V}_{\oplus}\) is the velocity of the Earth about the Sun. In Boreal summer, the addition of the Earth’s orbital velocity with that of the Sun boosts the “wind” of DM in our rest frame, whilst the wind is slower during the Austral summer. This changes the number of particles above the minimum speed required to create an observable nuclear recoil. The resulting annual modulation [@Drukier:1986tm; @Spergel:1987kx; @Freese:1987wu] is parameterised as \[\frac{dR}{dE_R}=A_0+A_1\cos\left[\omega(t-t_0)\right]+... \label{eq:parametrisation}\tag{2}\] Experimental results are reported not in terms of the nuclear recoil energy \(E_R\) but in terms of the measured electron-equivalent energy \(E_{ee}\), which is related to the former by a quenching factor \(Q_T(E_R)\) that depends on both the target \(T\) and the recoil energy. For a NaI detector with two potential nuclear targets \(T=\{\)Na,I\(\}\), the total number of events in a given experimental bin \(i\) is given by \[\begin{align} N_i &= \int_{E_{i,\text{min}}}^{E_{i,\text{max}}} dE_{ee} \epsilon(E_{ee}) \times\\ &\int_0^\infty dE_R \sum_T \xi_T\frac{1}{\sqrt{2\pi}\sigma(E_{ee})}e^{-\frac{(E_{ee}-Q_T(E_R)E_R)^2}{2\sigma(E_{ee})}} \frac{dR_T}{dE_R},\nonumber\label{eq:detectorresp} \end{align}\tag{3}\] where \(\xi_T\) indicates the mass fraction of the target \(T\), \(dR_T/dE_R\) is the predicted recoil energy spectrum for that target, \(E_{i,\text{min}},E_{i,\text{max}}\) define the bin energy range, \(\epsilon\) is the detector efficiency and \(\sigma\) is the detector resolution. As all collaborations report efficiency-corrected results, we will set \(\epsilon=1\) for all detectors.

DAMA/LIBRA phases 1 and 2 [@Bernabei:2021kdo; @Bernabei:2022xgg] have reported results compatible with a DM annual modulation signature in the energy range \(1-6 {\rm \,keV_{\rm ee}}\). However, they have presented results over a larger energy range, with DAMA/LIBRA-phase2 releasing data in the energy range \(0.75-20 {\rm \,keV_{\rm ee}}\) with a total exposure over 8 annual cycles of 1.53 t\(\times\)yr. This is in addition to data from DAMA/NaI and DAMA/LIBRA-phase1, which bring the total exposure in the range \(2-20 {\rm \,keV_{\rm ee}}\) to 2.86 t\(\times\)yr. In this work we fit to all of the available modulation data as presented in [@Bernabei:2021kdo]. The DAMA/LIBRA detector resolution is taken to be \(\sigma(E_{ee}) = a \sqrt{E_{ee}} + b E_{ee}\) with \(a=0.488\sqrt{{\rm \,keV_{\rm ee}}}\) and \(b=0.0091\) [@DAMA:2008bis].

COSINE-100 recently released results with \(6.4 \mathrm{yr}\) of exposure [@Carlin:2024maf], where they find no evidence of an annual modulation signal and a greater than \(3\sigma\) tension with DAMA/LIBRA in the energy range \(1-6 {\rm \,keV_{\rm ee}}\). In our fits, we include their reported modulations in the range \(0.75-20 {\rm \,keV_{\rm ee}}\) and make use of the detector resolution presented in [@Kang:2019uuj]: \(\sigma(E_{ee}) = \sqrt{a E_{ee} + b E_{ee}^2}\) with \(a = 0.081483 {\rm \,keV_{\rm ee}}\), \(b = 0.001885\).

ANAIS-112 recently presented binned modulation amplitude data from their 6-year exposure dataset [@Amare:2025dfq]. This data is incompatible with the DAMA/LIBRA modulation signal in the \(1-6 {\rm \,keV_{\rm ee}}\) range at nearly \(4 \sigma\) confidence level. We fit to their reported modulation amplitudes over the range \(1-20 {\rm \,keV_{\rm ee}}\), and take the detector resolution to be \(\sigma(E_{ee}) = a + b \sqrt{E_{ee}}\) with \(a = -0.08 {\rm \,keV_{\rm ee}}\), \(b = 0.378 \sqrt{{\rm \,keV_{\rm ee}}}\) [@Amare:2018sxx].

Equally important is the nuclear recoil quenching factor \(Q_T(E_R)\), defined as the ratio of scintillation light yield produced by nuclear recoil to that of electron recoil at the same energy. We take recently-measured energy-dependent quenching factors for sodium from Ref. [@Carlin:2024maf] and iodine from Ref. [@Lee:2024unz], for all experiments. We will, however, examine the impact of alternative quenching factors.

2 Quantifying the tension between the datasets↩︎

We follow the approach from [@Maltoni:2003cu] to quantify the tension between the datasets of different experiments.

Consider two experiments, \(A\) and \(B\), that measure a number of observables \(m_A\) and \(m_B\), respectively, and a model that describes the signal using \(n\) parameters \(c_i\). We define \(\chi_X^2(c_i)\) as the \(\chi^2\) obtained from the data of experiment \(X=A,B\). We denote the best-fit point obtained by fitting the data of experiment \(X\) alone as \(c_i^X\), so that \(\chi_X^2(c_i)\) has a minimum at \(c_i=c_i^X\). Similarly, \(c_i^{A+B}\) is the best-fit point obtained by fitting the data of both experiments simultaneously. In this case, the function \[\chi_{A+B}^2(c_i) = \chi_A^2(c_i) + \chi_B^2(c_i)\] has a minimum at \(c_i = c_i^{A+B}\), and the quantity \[\delta\chi^2 = \chi_{A+B}^2(c_i^{A+B}) - \chi_A^2(c_i^{A}) - \chi_B^2(c_i^{B}) \label{eq:tensionstat}\tag{4}\] follows a \(\chi^2\) distribution with \(n\) degrees of freedom.

This can be used to determine the level of statistical inconsistency of the two datasets, under the assumption that the chosen model describes the signals with the right level of precision to reproduce the data of the experiments. We will use this statistical test to assess the compatibility between DAMA/LIBRA and the combination of COSINE-100 and ANAIS-112. We will employ two qualitatively different approaches to derive the tension between the DAMA/LIBRA and combined COSINE-100 and ANAIS-112 datasets.

3 Detailed particle astrophysics test↩︎

The first is to use concrete assumptions for the astrophysics and particle physics defined in Eq. 1 . The DM distribution \(f(\mathbf{v})\), is well-described by a Maxwellian velocity distribution in the halo frame, with a peak velocity \(v_{0} = 240\, \pm \,8\) km s\(^{-1}\) [@Reid:2014boa], cut off at the escape velocity \(v_{\rm{\text{esc}}} = 528 \pm 25\) km s\(^{-1}\), based on Gaia data [@Deason:2019kgj]. The halo parameters are fixed to their central values, as varying these parameters is expected to affect all experiments in the same way. The potential impact of varying the velocity distribution will be implicitly captured by the more-model-independent parametrization discussed later. On the particle physics side, we consider elastic scattering of DM on nuclei. However, instead of simply considering spin-independent and spin-dependent scattering, we adopt a more model-agnostic approach by considering an effective field theory (EFT) allowing for different types of interactions between DM and Standard Model particles. This approach avoids the need to specify the detailed microphysics of DM interactions, as long as the relevant energy scale is below the cutoff scale \({\rm \Lambda}\). For elastic scattering, we can simply set \({\rm \Lambda}\) equal to the hadronic scale, i.e.\(\Lambda=2\) GeV without loss of generality. (A non-relativistic EFT approach was taken with earlier data [@Kang:2019uuj; @COSINE-100:2019wga; @COSINE-100:2021xqn]. Relativistic operators were also fit to the earlier data in [@Kang:2019dbr], but there only one operator was considered at a time, in contrast to our global fit.)

To construct the EFT, we further assume that the DM particle is a Dirac fermion and a singlet under the Standard Model gauge group. At the non-relativistic level, the interactions of all DM candidates map onto a finite set of Galilean-invariant operators. In this work we restrict ourselves to spin-1/2 dark matter for concreteness. While the interactions of scalar and vector dark matter arise from different relativistic operator bases, and hence lead to different correlations among the non-relativistic operator coefficients, the fermionic operators we consider already generate a broad set of relevant non-relativistic interaction structures. In particular, they include the dominant spin-independent and spin-dependent responses, as well as momentum- and velocity-suppressed interactions. We therefore expect that extending the analysis to spin-0 or spin-1 dark matter would not qualitatively change our conclusions, although it could shift the precise best-fit regions. Sec. 4 covers more generic parametrisations of the signal. Following the notation of Refs. [@Bishara:2017nnn; @Brod:2017bsw], we write the interaction Lagrangian for the theory as \[\mathcal{L}_{\rm{int}} = \sum_{a,d} \dfrac{ \if\relax\detokenize{d}\relax \mathcal{C}_{a} \else \mathcal{C}_{a}^{(d)} \fi }{{\rm \Lambda}^{d-4}} \if\relax\detokenize{d}\relax \mathcal{Q}_{a} \else \mathcal{Q}_{a}^{(d)} \fi \,,\] where \(\if\relax\detokenize{d}\relax \mathcal{Q}_{a} \else \mathcal{Q}_{a}^{(d)} \fi\) is a particular effective operator involving DM and Standard Model fields, \(d\geq 5\) is the mass dimension of the operator and \(\if\relax\detokenize{d}\relax \mathcal{C}_{a} \else \mathcal{C}_{a}^{(d)} \fi\) is the dimensionless Wilson coefficient associated to \(\if\relax\detokenize{d}\relax \mathcal{Q}_{a} \else \mathcal{Q}_{a}^{(d)} \fi\). The full Lagrangian density for the theory is then \(\mathcal{L} = \mathcal{L}_{\rm{SM}} + \mathcal{L}_{\rm{int}} + \overline{\chi}\left(i\cancel{\partial}-m_\chi\right)\chi\,,\) such that the free parameters of the theory are the DM mass \(m_\chi\), and the set of dimensionless Wilson coefficients \(\{ \if\relax\detokenize{d}\relax \mathcal{C}_{a} \else \mathcal{C}_{a}^{(d)} \fi \}\).

The phenomenology of DM in this model will generically be dominated by the lowest-dimension operators, unless their Wilson coefficients are accidentally suppressed or vanish, in which case higher-dimensional operators may become relevant. Under this assumption, we limit ourselves to operators with \(d\le 6\). At dimension 5, there are the two dipole operators \[\begin{align} \if\relax\detokenize{5}\relax \mathcal{Q}_{1} \else \mathcal{Q}_{1}^{(5)} \fi &= \frac{e}{8\pi^2} (\overline{\chi} \sigma_{\mu\nu} \chi) F^{\mu\nu} \,, \\ \if\relax\detokenize{5}\relax \mathcal{Q}_{2} \else \mathcal{Q}_{2}^{(5)} \fi &= \frac{e}{8\pi^2} (\overline{\chi} i \sigma_{\mu\nu} \gamma_5 \chi) F^{\mu\nu} \, , \end{align}\] where \(F_{\mu\nu}\) is the electromagnetic field strength tensor and \(e\) is the electromagnetic charge. These operators give rise to long-range interactions, i.e.steeply-falling recoil spectra.

At dimension six, we consider the operators \[\begin{align} \tag{5} \if\relax\detokenize{6}\relax \mathcal{Q}_{1,q} \else \mathcal{Q}_{1,q}^{(6)} \fi &= (\overline{\chi} \gamma_\mu \chi)(\overline{q} \gamma^\mu q)\,, \\ \if\relax\detokenize{6}\relax \mathcal{Q}_{2,q} \else \mathcal{Q}_{2,q}^{(6)} \fi &= (\overline{\chi} \gamma_\mu \gamma_5 \chi)(\overline{q} \gamma^\mu q)\,, \\ \if\relax\detokenize{6}\relax \mathcal{Q}_{3,q} \else \mathcal{Q}_{3,q}^{(6)} \fi &= (\overline{\chi} \gamma_\mu \chi)(\overline{q} \gamma^\mu \gamma_5 q)\,, \\ \if\relax\detokenize{6}\relax \mathcal{Q}_{4,q} \else \mathcal{Q}_{4,q}^{(6)} \fi &= (\overline{\chi} \gamma_\mu \gamma_5 \chi)(\overline{q} \gamma^\mu \gamma_5 q)\,. \tag{6} \end{align}\] The first two operators give rise to spin-independent interactions, while the last two give rise to spin-dependent interactions. Moreover, \(\if\relax\detokenize{6}\relax \mathcal{Q}_{1,q} \else \mathcal{Q}_{1,q}^{(6)} \fi\) and \(\if\relax\detokenize{6}\relax \mathcal{Q}_{4,q} \else \mathcal{Q}_{4,q}^{(6)} \fi\) are independent of the momentum transfer and the DM velocity, while \(\if\relax\detokenize{6}\relax \mathcal{Q}_{2,q} \else \mathcal{Q}_{2,q}^{(6)} \fi\) and \(\if\relax\detokenize{6}\relax \mathcal{Q}_{3,q} \else \mathcal{Q}_{3,q}^{(6)} \fi\) are suppressed in the low-velocity limit. Together these operators therefore capture a wide range of different possibilities for elastic scattering.

The effective operators are defined at the scale \({\rm \Lambda}=2\) GeV, where the Higgs, \(W\) and \(Z\) bosons as well as the top, bottom and charm quarks have been integrated out. We do not consider interactions with leptons, which do not give rise to nuclear scattering at tree-level. Following the assumption of Minimal Flavour Violation, we take the Wilson coefficients for operators involving the down and strange quarks to be equal, but we allow the Wilson coefficients for operators involving up quarks to differ. We parameterize the relative couplings between the \(u\) and \(d\)-type quarks by angles \(\theta^{(6)}_a\), in the following form: \[\begin{align} C^{(6)}_{a,d} &= C^{(6)}_a \sin \theta^{(6)}_a = C^{(6)}_{a,s}\,, \\ C^{(6)}_{a,u} &= C^{(6)}_a \cos \theta^{(6)}_a \; . \end{align}\] In the EFT setup, we therefore have 11 model parameters to consider.

The fits are performed with GAMBIT[@gambit; @gambit_addendum] and its DarkBit [@GAMBITDarkMatterWorkgroup:2017fax] module, using Diver 1.3¹ [@ScannerBit] to explore the parameter space, DirectDM [@Bishara:2016hek; @Bishara:2017pfq; @Bishara:2017nnn] to match the effective operators introduced above onto the non-relativistic operators relevant for nuclear scattering, and DDCalc² [@GAMBIT:2018eea] to evaluate nuclear form factors, calculate the differential event rate, and evaluate the experimental likelihoods.

The best-fit values for all parameters for each scan to the modulation spectrum from DAMA/LIBRA, both ANAIS-112 and COSINE-100, or all three experiments combined, are shown in Table 1. The corresponding \(\chi^2\) values and resulting tension is shown in Table 2. The combination of ANAIS-112 and COSINE-100 is found to rule out the DM nuclear recoil interpretation of the DAMA signal with a significance of \(5.1\sigma\).

Our best fit to the DAMA data alone prefers interactions via the \(C^{(6)}_{3}\) operator, which induces spin- and momentum-dependent interactions, with a DM mass of 184 GeV. We find an almost equally good fit at a mass of approximately 40 GeV. Fitting to ANAIS+COSINE, as well as to all three experiments, yields small couplings to all operators, i.e. no DM signal.

Restricting the nuclear recoil energy range in DAMA below 7 keV\(_{\mathrm ee}\) provides a best fit mass of 36 GeV. However, this restriction does not change the qualitative conclusions regarding the tension between DAMA and other experiments, which remains above 5\(\sigma\), hence we choose to present the fit over the full range for completeness.

It has been speculated that quenching factors may be specific to individual crystals, rather than to the material itself. Since DAMA has never measured the energy dependence of its crystal’s quenching factor, they assume a constant value of \(0.3\) for sodium and \(0.09\) for iodine [@Bernabei:1996vj]. Measurements from other groups have consistently shown a strong energy-dependence, with sodium quenching factors never being higher than \(\sim 0.2\), and going as low as \(0.1\) below 5 keV [@Carlin:2024maf]. Nonetheless, we have repeated all scans using the DAMA quenching factors for the DAMA response modeling, and still find that this only modestly reduces the tension to 4.63\(\sigma\). The final column of table 2 shows the resulting \(\chi^2\) values.

Fig. 2 shows the expected spectrum seen at DAMA/LIBRA (left), ANAIS-112 (centre) and COSINE-100 (right), if the signal is fit to that experiment only (orange), or to all three simultaneously (red). The large discrepancy between the data points (black) and the red curves for all three experiments illustrates the tension—in this EFT model, the strong preference of DAMA for modulation at low energies leads to similar predictions for the other experiments, a scenario which the data do not support.

In addition to the tension present in the modulation data, there is also tension in the predicted total rate. The best-fit point in the EFT setup predicts a total event rate in COSINE corresponding to around 8000 signal events in the bin [1, 1.5] keV, which should be compared to an observed number of 32294 events and a background expectation of \(31572 \pm 1074\) events based on an exposure of 97.7 kg years [@COSINE-100:2021xqn]. Hence, while the total rate predicted by the best-fit point cannot be excluded considering the total rate alone, it is in strong tension with data once the background model is taken into account. However, since these data points are not statistically independent from the modulation data, they cannot be easily combined with the tension presented above.

4 How generic is the tension?↩︎

The use of specific particle physics and astrophysics assumptions means that any derived tension between DAMA/LIBRA and other experiments is necessarily model-dependent. We therefore also adopt a second approach which simply replaces the right-hand side of Equation 1 by a generic functional form that is not constrained by any known particle, nuclear or astrophysics. We consider two sets of model-independent parameterisations: 1) A modulated nuclear recoil signal injected in \(2N\) uniform \(E_R\) bins between 1 and 80 keV³ allowing independent bin amplitudes between recoils with sodium and iodine, and 2) a polynomial \(\times\) exponential parameterisation, inspired by the expected DM signal, but leaving coefficients completely free: \[\begin{align} A_1 &=& \sum_{T = \mathrm{Na,I}} \left(c^T_0 + c^T_1 E_R + c^T_2E_R^2 + ...\right)e^{-d_T E_R} , \label{eq:polyexp} \end{align}\tag{7}\] where \(d_T\) and the \(c_i^T\) are free parameters. We use Diver 1.3 embedded in GAMBIT-light⁴ to maximise the likelihoods of the resulting spectra with respect to the DAMA/LIBRA data, combined ANAIS-112 and COSINE-100 datasets, and the combination of all three experiments, and compute the tension statistic as in Eq. 4 . Demanding that the spectra be produced only by scattering with Na or I further worsens the tension in all cases. Resulting best fit \(\chi^2\) values, p-values, and tension statistics are presented in the appendix.

Using the tension statistic defined above, more than 12 bins (6 Na, 6 I) are required for the tension between DAMA/LIBRA and ANAIS-112 and COSINE-100 to fall below 5\(\sigma\). With 20 bins, the tension reduces to 2\(\sigma\). If we only allow for a signal in Na or I, the tension does not fall below 5\(\sigma\) for 10 or fewer bins. Indeed, we find that the tension is only reduced when the contributions from Na and I in each bin have opposite signs (such that the one contribution peaks in summer, while the other peaks in winter), allowing for delicate cancellations to simultaneously fit the overall signal shape in all experiments thanks to their different binning, and therefore different mapping from nuclear to electronic recoil.

The DM-inspired parameterisation from Eq. 7 also leads to a reduced tension as the number of parameters is increased. When including only 4 free parameters (i.e. \(c^{T}_0 e^{-d_T E_R}\) for both sodium and iodine), the tension is above the 5\(\sigma\) level, but it drops modestly to 4.5, 4.1 and 3.2\(\sigma\) as linear, quadratic and cubic terms in energy are added, respectively.

As in the binned case, the best-fit points to the combined data tend to be pathological, preferring a signal with a positive modulation amplitude for Na, and negative for I. In other words, while recoils off sodium peak in the summer, nuclear recoils on iodine would need to peak in winter, or vice-versa. Fig. 2 shows the resulting signals in all three experiments in the case of the 10-parameter fit. The NR signal that leads to these results is in Fig. 3.

It is worth noting that this solution implies that the individual modulation amplitudes \(A_1^T\) in both sodium and iodine must be much larger than the actually observed modulation amplitude \(A_1 = \sum_T A_1^T\). However, the positivity of the time-dependent event rate in Eq. 2 implies that the average rate \(A^T_0\) in each target must satisfy \(A^T_0 \gtrsim |A^T_1|\). If \(A_1^T\) has opposite sign for sodium and iodine, it follows that the total rate \(A_0\) must be much larger than the observed modulation amplitude \(|A_1|\). As discussed above, the total rate can be constrained using independent measurements from COSINE-100 [@COSINE-100:2021xqn]. This self-consistency requirement would place these scenarios under further strain.

5 Conclusion↩︎

We have revisited the tension between the DAMA/LIBRA, ANAIS-112 and COSINE-100 DM annual modulation datasets. For a Dirac fermion interacting with Standard Model particles through dimension-five and dimension-six operators, DAMA/LIBRA is in a 5.10 \(\sigma\) tension with the other two experiments. This conclusion is now independent of the fact that direct search DM experiments with other target nuclei would exclude the DAMA/LIBRA excess. If instead a generic parameterisation of the modulation amplitude vs nuclear recoil energy is used, the tension is still greater than 5\(\sigma\) unless there is a significant degree of fine-tuning, resulting in a cancellation of the contributions from sodium and iodine. It remains interesting to consider a Southern hemisphere NaI experiment, such as the forthcoming SABRE South experiment [@SABRE:2022twu], which may shed further light on the mystery, in particular because it is expected to place a world-leading limit on the total rate of DM scattering in NaI detectors.

Appendix↩︎

We present additional details on the results of our fits in Table 3.

---

1. https://github.com/diveropt/Diver/releases/tag/v1.3.0↩︎

2. The rate calculations and likelihoods for the modulation experiments included in this paper will be part of the the upcoming DDCalc 3.0 public release.↩︎

3. This range in nuclear recoil energy ensures that for quenching factors ranging from 5% to 30%, we cover the signal range in electron-equivalent energy.↩︎

4. GAMBIT-light is a lightweight version of GAMBIT, available at github.com/GambitBSM/gambit_light_1.0.↩︎