First results on dark matter annual modulation from ANAIS-112 experiment

TL;DR

The ANAIS-112 experiment analyzed 1.5 years of data to test the DAMA/LIBRA dark matter modulation claim, finding results consistent with no modulation and supporting future sensitivity goals.

Contribution

First experimental results from ANAIS-112 testing the DAMA/LIBRA modulation claim with null results over 1.5 years.

Findings

Data are consistent with the null hypothesis of no modulation.

Best fits show no significant modulation signal.

Results support the projected 3σ sensitivity in 5 years.

Abstract

ANAIS is a direct detection dark matter experiment aiming at the testing of the DAMA/LIBRA annual modulation result, which standing for about two decades has neither been confirmed nor ruled out by any other experiment in a model independent way. ANAIS-112, consisting of 112.5 kg of sodium iodide crystals, is taking data at the Canfranc Underground Laboratory, Spain, since August 2017. This letter presents the annual modulation analysis of 1.5 years of data, amounting to 157.55 kg$\times$y. We focus on the model independent analysis searching for modulation and the validation of our sensitivity prospects. ANAIS-112 data are consistent with the null hypothesis (p-values of 0.65 and 0.16 for [2-6] and [1-6] keV energy regions, respectively). The best fits for the modulation hypothesis are consistent with the absence of modulation ($S_m$=-0.0044$\pm$0.0058 cpd/kg/keV and -0.0015$\pm$0.0063 cpd/kg/keV, respectively). They are in agreement with our estimated sensitivity for the accumulated exposure, supporting our projected goal of reaching a 3$\sigma$ sensitivity to the DAMA/LIBRA result in 5 years of data taking.

Figures

Equations

$$ϵ(E,d)={ϵ}_{trg}(E,d)\times {ϵ}_{PSA}(E,d)\times {ϵ}_{asy}(E,d),$$

$$R(t)=R_0+R_1\cdot exp(-t/\tau )+S_m\cdot cos(\omega \cdot (t+\varphi )),$$

Full text

There is overwhelming evidence from cosmological and astrophysical observations supporting that a large fraction of the Universe energy-mass budget is not explained in the framework of the standard model of the particle physics, assuming the cosmological standard model 1. The solution to the dark matter (DM) / dark energy (DE) puzzle is probably of a complex nature. In one of the preferred hypothetical scenarios, DM can be understood as a new non-zero-mass particle having a very low interaction probability with baryonic matter. Although proposed candidates span about 45 orders of magnitude in mass, and 60 in cross-section, axions and Weakly Interacting Massive Particles (WIMPs) are among the better motivated 2. Experimental effort devoted to unraveling the nature of the DM particles has been spent, either by direct 3, indirect 4 or accelerator searches 5, which are complementary to each other. Only one experiment, DAMA/LIBRA 6, 7, 8, has provided a long-standing positive result: the observation of a highly statistically significant annual modulation in the detection rate, compatible with that expected for galactic halo dark matter particles. This result has neither been reproduced by any other experiment, nor ruled out in a model independent way. Compatibility among the different experimental results in most conventional WIMP-DM scenarios is actually disfavored 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19. Then, a similar annual modulation search using the same target is mandatory to shed light on the DAMA/LIBRA conundrum, which is the goal of the ANAIS (Annual modulation with NaI Scintillators) experiment. Other efforts sharing ANAIS goal in the international dark matter community are the COSINE-100 experiment, taking data also in dark matter mode at Yang-Yang Underground Laboratory, South Korea 20, 21, 9; and in longer-term, SABRE, aiming at installing twin detectors in Australia and Italy 22, and COSINUS, developing cryogenic detectors based on NaI 23.

An annual modulation in the dark matter interaction rate is expected by the revolution of the Earth around the Sun, which distorts the DM particle velocity distribution function as seen by the detector, typically assumed Maxwellian boosted by the Sun velocity 24, 25. The effect is present unless the DM halo is co-rotating with the Solar System. However, it is strongly dependent on the specific halo model, both in amplitude and in phase. It is natural to assume that the Sun is moving through a locally isotropic DM halo, with the Earth orbiting aside. Consequently, searches are performed for a modulation of DM-like events with a one year period and a well defined phase. On the other hand, preferably an annual modulation analysis should not assume a priori neither the period of the modulation nor the phase, but it should derive them from the data. A full and consistent analysis requires then several years of measurement in very stable conditions. This is the long-term goal of our experiment.

$\rm ANAIS{\hbox{-}}112$, consisting of 112.5 kg of NaI(Tl) detectors, was installed in 2017 at the Canfranc Underground Laboratory, LSC, in Spain. The $\rm ANAIS{\hbox{-}}112$ set-up undergoes a different residual cosmic ray flux and environmental conditions than DAMA/LIBRA (800 m versus 1400 m rock overburden, for instance). Consequently, the potential confirmation of a modulation with same phase and amplitude would be very difficultly explained as an effect of backgrounds or systematics. $\rm ANAIS{\hbox{-}}112$ experimental details have been recently reported, as well as the performance of the first year’s operation 26, analysis of backgrounds 27, and sensitivity prospects for a five-years operation 28; here we just briefly summarize the most relevant features of the experimental apparatus.

$\rm ANAIS{\hbox{-}}112$ uses nine NaI(Tl) modules produced by Alpha Spectra Inc. in Colorado (US)‡‡‡http://www.alphaspectra.com/. These modules have been manufactured from 2012 to 2017, and shipped to Spain avoiding air travel in order to prevent cosmogenic activation of the module materials. Each crystal is cylindrical (4.75” diameter and 11.75” length), with a mass of 12.5 kg, and it is housed in OFE (Oxygen Free Electronic) copper. This encapsulation has a Mylar window allowing low energy calibration using external gamma sources. It incorporates two quartz optical windows to couple the photomultipliers (PMTs). All PMT units, as well as all relevant materials used in the building of the detectors, have been screened for radiopurity using HPGe detectors in the low background facilities at LSC. Their contribution to the experiment background has been estimated 29, 27 and included in our background model (see below). We would like to emphasize the outstanding light collection measured for the nine modules, at the level of 15 photoelectrons (phe) per keV 26, 30. $\rm ANAIS{\hbox{-}}112$ is calibrated every two weeks using external 109Cd sources: all the nine modules are simultaneously calibrated using a multi-source system which minimizes down time periods. Background events from the decay of 40K and 22Na in the crystal bulk, associated to 3.2 and 0.9 keV energy depositions, and selected by coincidence with an energy deposition in a second module around 1461 and 1275 keV, respectively, are also used to improve the accuracy of the calibration down to the energy threshold 26.

The $\rm ANAIS{\hbox{-}}112$ shielding consists of 10 cm of archaeological lead, 20 cm of low activity lead, an anti-radon box (continuously flushed with radon-free nitrogen gas), an active muon veto system made up of 16 plastic scintillators designed to cover top and sides of the whole ANAIS set-up and 40 cm of neutron moderator (a combination of water tanks and polyethylene blocks). In the design of the muon veto system we follow a tagging strategy instead of a hardware vetoing. The goal is twofold: on the one hand, to discard events in the NaI(Tl) crystals coincident with muon veto triggers. On the other hand, to analyse eventual correlations between muon hits in the plastic scintillators and events in the NaI(Tl) crystals, specially in the region of interest (ROI), $\rm[1{\hbox{-}}6]$ keV§§§Energy will be shown in electron equivalent units throughout this letter.

The $\rm ANAIS{\hbox{-}}112$ electronic chain and data acquisition system (DAQ) have been described in Refs. 31, 26. Each PMT charge signal is independently processed and divided into: (1) a trigger signal; (2) a low energy (LE) signal that goes to the digitizers which sample the waveforms at 2 Gs/s with high resolution (14 bits); and (3) a high energy (HE) signal, conveniently attenuated. The trigger of each PMT signal is done at phe level, while the single module trigger is done by the coincidence (logical AND) of the two PMT triggers in a 200 ns window. The global trigger is the logical OR of the nine modules trigger signals. Trigger efficiency is close to 100% down to the analysis threshold established at 1 keV 26.

$\rm ANAIS{\hbox{-}}112$ started taking data in the DM mode on August, 3rd, 2017. It has accumulated almost 1.5 years of data-taking time in quite stable conditions by February, 12th, 2019. Total live time available for the annual modulation analysis is 527.08 days: 341.72 (first year) + 185.36 (half second year). This implies a live time of 94.5% (95.2%), dead time of 2.9% (2.2%), and down time of 2.6% (2.6%) for the first (second) year of data taking, respectively. The down time is mainly due to the periodical calibration runs carried out using low energy gamma sources. We remove events arriving less than one second after the last muon veto trigger, correcting the total live time by subtracting one second per muon veto trigger, so the live time used for the annual modulation analysis which follows is 511.16 days.

A blind analysis strategy in three levels has been followed: first, we calculate pulse parameters, the time since the last muon veto, and we apply the peak-finding algorithm to identify individual phe in low energy pulses; second, we calibrate the energy response of every detector at LE and HE 26; third, we optimize the pulse shape cuts and calculate their efficiency, being the LE variable hidden for events corresponding to single hits (M1 events). Only 10% of the data were unblinded for general background assessment and fine tuning of procedures 26; those data were randomly chosen (34 days amounting to 32.9 days live time) and time evolution was kept hidden until the data unblinding, presented in this letter.

Calibration procedures and efficiency corrections applied in the following have been derived as done in Ref. 26 for the first year, continuing with same procedure in the half second year added in this analysis. Events in the ROI are selected, after energy calibration, by imposing the following criteria: single hit events (M1); a pulse shape cut combining the fraction of the pulse area in [100-600] ns after the event trigger, defined following 32, and the logarithm of the mean time of the distribution of the individual phe arrival times in the digitized window 20; events having an asymmetric light sharing between the two PMT signals are removed by imposing a cut on the number of peaks identified in each PMT (asymmetry cut). The total detection efficiency, $\epsilon(E,d)$, calculated independently for every detector $d$ as a function of the energy, $E$, can be written 26 as

[TABLE]

where the trigger efficiency $\epsilon_{trg}(E,d)$ is calculated from Monte Carlo (MC) simulations, and the efficiencies of the pulse shape, $\epsilon_{PSA}(E,d)$, and asymmetry, $\epsilon_{asy}(E,d)$, cuts are evaluated from the 3.2 and 0.9 keV events selected by the coincidence with the high energy gammas following 40K and 22Na decays, and 109Cd calibration events, respectively, accumulated for all the analysed exposure. Total detection efficiency ranges from 0.15 to 0.35 at 1 keV, depending on the detector, increases up to 0.8 at 2 keV and is nearly 1 at 4 keV for all the modules. Statistical errors in the total efficiency vary from 3 to 8% at 1 keV down to 1% at 6 keV. Comparing different methods for the efficiency calculation we have also estimated a systematic uncertainty that amounts up to 20% at 1-1.2 keV and is negligible above 1.5 keV. More details can be found in 26.

The background model for all the nine detectors used in the $\rm ANAIS{\hbox{-}}112$ set-up has been developed. It is based on MC simulations using the measured activity in external components and in crystals, including cosmogenic products, quantified in dedicated, independent measurements using different analysis techniques. It provides a good overall description of measured data at all energy ranges above 2 keV and at different analysis conditions (coincidence or anticoincidence) 27. In the ROI the background is dominated by the emissions from the crystals themselves, in particular, 210Pb (32.5%) and 3H (26.5%) continua, and 40K (12%) and 22Na (2.0%) peaks are the most significant contributions. Short-lived isotopes activated cosmogenically are still present in the bulk of the last received crystals 33, 34, contributing as background in the ROI, specially in the [3-5] keV region. However, from 1 to 2 keV there is a large fraction of our background lacking from explanation 27. It could have as origin non-bulk scintillation leaking through our event selection criteria.

The time evolution of the rate of those events surviving all the cuts during the first year and a half of $\rm ANAIS{\hbox{-}}112$ operation is shown for different populations in panels a)-f) of Figure 1. Data from all the modules have been added together and corrected by the corresponding efficiencies. The two lower panels correspond to two different energy windows in the ROI: $\rm[1{\hbox{-}}6]$ keV, a), and [3-5] keV, b), while the upper panels show the evolution of control populations for which no modulation is expected: c) [6-20] keV, d) double-coincidence events (M2) in energy region $\rm[1{\hbox{-}}6]$ keV, and coincident events attributed to e) 40K, and f) 22Na decays in the crystals¶¶¶Low energy M2 events having an energy deposition in a second detector in a window around the corresponding high energy gamma. The uppermost panel, g), presents the evolution of muon-related low-energy events before the cuts (M1 events in $\rm[1{\hbox{-}}6]$ keV region arriving less than 1 second after a muon veto trigger). As pointed out in Ref. 26, many low energy events are registered after a muon passage through the detector. They likely originate in the long tail of the scintillation pulse produced by a large muon energy deposition, which is able to trigger many times the DAQ system and produce fake low-energy events. This population is clearly fluctuating in time and will be studied in a future work.

In panels a) and b) of Figure 1 we observe a relevant decrease in the rate at the ROI, which amounts up to 8% in the analysed period. It is caused by cosmogenically activated isotopes and is well reproduced qualitatively in all energy windows by our background model 27, shown in Figure 1 in green dashed lines. It is normalized by a factor, $f$, to be more easily compared with the measured rates. Agreement in [3-5] keV region is also quantitatively very good, $f$=1.04, as the background time evolution in this region is dominated by short-lived cosmogenic isotopes. However, our background model underestimates the rate in $\rm[1{\hbox{-}}6]$ keV region ($f$=1.28), because from 1 to 2 keV our background model does not reproduce the measurement 27, as commented above.

We perform a least-squares fit on panels a) to d) to a constant term plus an exponential function, to account for the mentioned background reduction, following our background model 27. The $\chi^{2}/NDF$ and the values of the fitted parameters are shown in the figure. Good fits are obtained, with p-values larger than 0.20 in all cases. The 40K (3.2 keV) and 22Na (0.9 keV) M2 populations, on the other hand, are modelled by a constant and an exponential decay, respectively. We derive consistent results (p-value = 0.22) in the first case, and a reasonable agreement with the 22Na half-life ($T_{1/2}=3.01\pm 0.36$ y, p-value = 0.07) in the second.

In this letter we present our modulation results in two energy regions: $\rm[2{\hbox{-}}6]$ keV and $\rm[1{\hbox{-}}6]$ keV, to allow direct comparison with the DAMA/LIBRA results. The values of the modulation amplitude observed by DAMA/LIBRA, $S_{m}^{DAMA}$, are $0.0102\pm 0.0008$ and $0.0105\pm 0.0011$ cpd/kg/keV in the full exposure for $\rm[2{\hbox{-}}6]$ keV and using only phase-2 data for $\rm[1{\hbox{-}}6]$ keV energy region, respectively 8. We expect results derived from $\rm[2{\hbox{-}}6]$ keV to be more robust because our data selection efficiencies strongly go down below 2 keV, increasing the risk to be affected by unknown systematics.

We evaluate the statistical significance of a possible modulation in our data by a least square method in the time-binned data. The efficiency-corrected rate of events surviving the cuts in $\rm[1{\hbox{-}}6]$ and $\rm[2{\hbox{-}}6]$ keV energy regions is modelled as

[TABLE]

where $R_{0}$ and $R_{1}$ are free parameters and $\tau$ is fixed to the value obtained from our background model in the corresponding energy range. We also fix the period ($\omega=2\pi/365$ d $=0.01721$ rad d*-1*) and the phase ($\phi=-62.2$ d, corresponding the cosine maximum to June, 2 when taking as time origin August 3), while Sm is fixed to 0 for the null hypothesis and left unconstrained (positive or negative) for the modulation hypothesis. This allows a direct comparison with the results from the DAMA/LIBRA analysis with 1 free parameter 8. We present the best fit for both hypothesis for ${\rm 10{\hbox{-}}day}$ time binning in Figure 2. In order to highlight the presence or absence of modulation, we plot the data with the constant and exponential terms subtracted. For the sake of comparison, in the plot we show the modulation measured by DAMA/LIBRA (green lines).

In both energy regions the null hypothesis is well supported by the $\chi^{2}$ test, with $\chi^{2}/NDF=48.0/53$ for the $\rm[2{\hbox{-}}6]$ keV (p-value = 0.67) and $\chi^{2}/NDF=62.0/53$ for the $\rm[1{\hbox{-}}6]$ keV regions (p-value = 0.18). The best fits for the modulation hypothesis are $S_{m}=-0.0044\pm 0.0058$ cpd/kg/keV and $-0.0015\pm 0.0063$ cpd/kg/keV for $\rm[2{\hbox{-}}6]$ keV and $\rm[1{\hbox{-}}6]$ keV, respectively. In both cases, p-values are slightly lower than those of the null hypothesis (0.65 and 0.16, respectively). The best fits are incompatible at 2.5$\sigma$ (1.9$\sigma$) with the DAMA/LIBRA signal.

The statistical significance of our result is determined by the standard deviation of the modulation amplitude distribution, $\sigma(S_{m})$, which would be obtained in a large number of experiments like $\rm ANAIS{\hbox{-}}112$ with the present exposure. Then, we quote our sensitivity to DAMA/LIBRA result as the ratio $S_{m}^{DAMA}/\sigma(S_{m})$, which directly gives in $\sigma$ units the C.L. at which we can test the DAMA/LIBRA signal. At present, our result $\sigma(S_{m})=0.0058~{}(0.0063)$ cpd/kg/keV for $\rm[2{\hbox{-}}6]$ keV ($\rm[1{\hbox{-}}6]$ keV) corresponds to a sensitivity of 1.75$\sigma$ (1.66$\sigma$) to the DAMA/LIBRA signal. In Ref. 28 we found an analytical expresion to calculate $\sigma(S_{m})$ at a given exposure from the measured background and detection efficiency. Figure 3 (dark blue lines) displays our sensitivity projection calculated following Ref. 28 for the two studied energy ranges, whereas the blue bands represent the 68% uncertainty in $S_{m}^{DAMA}$ as reported in Ref. 8. In the calculation we take into account the $\rm ANAIS{\hbox{-}}112$ live time distribution, the background reduction expected due to decaying isotopes and the statistical error in the detection efficiency. The black dots are the sensitivities derived in this work, including a systematic error estimated by changing the time-binning from 1 to 20 days, and considering the systematics in the efficiency 26. The results perfectly agree with our estimates, confirming the $\rm ANAIS{\hbox{-}}112$ projected sensitivity to the DAMA/LIBRA result. A 3$\sigma$ sensitivity should be at reach in 4-5 years of data-taking.

Finally, Figure 4 presents the best fit amplitudes, $S_{m}$, calculated per 1 keV energy bins, from 1 to 20 keV (bottom panel, black dots), together with the DAMA-phase-2 modulation amplitudes extracted from Ref. 8 (blue triangles). The top panel shows the p-values for the null (open squares) and modulation hypothesis (closed circles) for every energy bin. All the modulation amplitudes in the ROI are compatible with 0 and, in general, p-values for the null hypothesis are slightly larger than those of the modulation hypothesis. The 1$\sigma$ and 2$\sigma$ bands shown in the figure are obtained following Ref. 28 for the present $\rm ANAIS{\hbox{-}}112$ exposure.

In summary, to test the DAMA/LIBRA annual modulation result in a model independent way, an analysis of the first year and half of data from $\rm ANAIS{\hbox{-}}112$ experiment has been presented. It is compatible with the sensitivity estimates done in Ref. 28, and confirms our goal of reaching sensitivity at 3$\sigma$ level in five years (see Figure 3) to DAMA/LIBRA result. We derive best fits for the modulation amplitude $S_{m}=-0.0044\pm 0.0058$ and $-0.0015\pm 0.0063$ cpd/kg/keV, in the $\rm[2{\hbox{-}}6]$ and $\rm[1{\hbox{-}}6]$ keV energy regions, respectively, compatible with the absence of modulation.

This work has been financially supported by the Spanish Ministerio de Economía y Competitividad and the European Regional Development Fund (MINECO-FEDER) under grants No. FPA2014-55986-P and FPA2017-83133-P, the Consolider-Ingenio 2010 Programme under grants MultiDark CSD2009-00064 and CPAN CSD2007-00042 and the Gobierno de Aragón and the European Social Fund (Group in Nuclear and Astroparticle Physics and I. Coarasa predoctoral grant). We thank the support of the Spanish Red Consolider MultiDark FPA2017-90566-REDC. We acknowledge the technical support from LSC and GIFNA staff. Professor J.A. Villar passed away in August, 2017. Deeply in sorrow, we all thank him for his dedicated work and kindness.

• Bertone and Hooper [2018]

G. Bertone and D. Hooper, Rev. Mod. Phys. 90, 045002 (2018), eprint arXiv:1605.04909.

• Baudis [2016]

L. Baudis, J. Phys. G43, 044001 (2016).

• Liu et al. [2017]

J. Liu, X. Chen, and X. Ji, Nature Phys. 13, 212 (2017), eprint arXiv:1709.00688.

• Conrad and Reimer [2017]

J. Conrad and O. Reimer, Nature Phys. 13, 224 (2017), eprint arXiv:1705.11165.

• Buchmueller et al. [2017]

O. Buchmueller, C. Doglioni, and L. T. Wang, Nature Phys. 13, 217 (2017).

• Bernabei et al. [2008]

R. Bernabei et al. (DAMA), Eur. Phys. J. C56, 333 (2008), eprint arXiv:0804.2741.

• Bernabei et al. [2013]

R. Bernabei et al., Eur. Phys. J. C73, 2648 (2013), eprint arXiv:1308.5109.

• Bernabei et al. [2018]

R. Bernabei et al., Universe 4, 116 (2018), [At. Energ.19,307(2018)], eprint arXiv:1805.10486.

• Adhikari et al. [2018a]

G. Adhikari et al., Nature 564, 83 (2018a).

• Kobayashi et al. [2018]

M. Kobayashi et al. (XMASS) (2018), eprint arXiv:1808.06177.

• Akerib et al. [2018]

D. S. Akerib et al. (LUX), Phys. Rev. D98, 062005 (2018), eprint arXiv:1807.07113.

• Abe et al. [2018]

K. Abe et al. (XMASS), Phys. Rev. D97, 102006 (2018), eprint arXiv:1801.10096.

• Aprile et al. [2017]

E. Aprile et al. (XENON), Phys. Rev. Lett. 118, 101101 (2017), eprint arXiv:1701.00769.

• Savage et al. [2009]

C. Savage, K. Freese, P. Gondolo, and D. Spolyar, JCAP 0909, 036 (2009), eprint arXiv:0901.2713.

• Aprile et al. [2015]

E. Aprile et al. (XENON100), Science 349, 851 (2015), eprint arXiv:1507.07747.

• Herrero-Garcia [2015]

J. Herrero-Garcia, JCAP 1509, 012 (2015), eprint arXiv:1506.03503.

• Baum et al. [2019]

S. Baum, K. Freese, and C. Kelso, Phys. Lett. B789, 262 (2019), eprint arXiv:1804.01231.

• Kang et al. [2018]

S. Kang, S. Scopel, G. Tomar, and J.-H. Yoon, JCAP 1807, 016 (2018), eprint arXiv:1804.07528.

• Herrero-Garcia et al. [2018]

J. Herrero-Garcia, A. Scaffidi, M. White, and A. G. Williams, Phys. Rev. D98, 123007 (2018), eprint arXiv:1804.08437.

• Kim et al. [2018]

K. W. Kim et al. (KIMS), arXiv:1806.06499 (2018).

• Adhikari et al. [2018b]

G. Adhikari et al., Eur. Phys. J. C78, 107 (2018b), eprint arXiv:1710.05299.

• Tomei [2017]

C. Tomei (SABRE), Nucl. Instrum. Meth. A845, 418 (2017).

• Gütlein et al. [2017]

A. Gütlein et al., Nucl. Instrum. Meth. A845, 359 (2017).

• Freese et al. [1988]

K. Freese, J. A. Frieman, and A. Gould, Phys. Rev. D37, 3388 (1988).

• Freese et al. [2013]

K. Freese, M. Lisanti, and C. Savage, Rev. Mod. Phys. 85, 1561 (2013), eprint arXiv:1209.3339.

• Amaré et al. [2018a]

J. Amaré et al., accepted for publication in Eur. Phys. J. C (2018a), eprint arXiv:1812.01472.

• Amaré et al. [2018b]

J. Amaré et al., submitted to Eur. Phys. J. C (2018b), eprint arXiv:1812.01377.

• Coarasa et al. [2018]

I. Coarasa et al., accepted for publication in Eur. Phys. J. C (2018), eprint arXiv:1812.02000.

• Amare et al. [2016]

J. Amare et al., Eur. Phys. J. C76, 429 (2016), eprint arXiv:1604.05587.

• Oliván et al. [2017]

M. A. Oliván et al., Astropart. Phys. 93, 86 (2017), eprint arXiv:1703.01262.

• Oliván [2016]

M. Oliván, Ph.D. thesis, Universidad de Zaragoza (2016).

• Bernabei et al. [2008]

R. Bernabei et al., Nuclear Instruments and Methods in Physics Research A 592, 297 (2008), eprint arXiv:0804.2738.

• Villar et al. [2018]

P. Villar et al., Int. J. Mod. Phys. A33, 1843006 (2018).

• Amaré et al. [2015]

J. Amaré et al., JCAP 1502, 046 (2015), eprint 1411.0106.