First results on dark matter annual modulation from ANAIS-112 experiment
J. Amaré, S. Cebrián, I. Coarasa, C. Cuesta, E. García, M. Martínez, M.A. Oliván, Y. Ortigoza, A. Ortiz de Solórzano, J. Puimedón, A. Salinas, M.L. Sarsa, P. Villar, J.A. Villar
FFirst results on dark matter annual modulation from ANAIS-112 experiment
J. Amar´e,
1, 2
S. Cebri´an,
1, 2
I. Coarasa,
1, 2
C. Cuesta,
1, 3
E. Garc´ıa,
1, 2
M. Mart´ınez,
1, 2, 4
M.A. Oliv´an,
1, 5
Y. Ortigoza,
1, 2
A. Ortiz de Sol´orzano,
1, 2
J. Puimed´on,
1, 2
A. Salinas,
1, 2
M.L. Sarsa ∗ ,
1, 2
P. Villar,
1, 2 and J.A. Villar †1, 2 Laboratorio de F´ısica Nuclear y Astropart´ıculas, Universidad de Zaragoza, C/ Pedro Cerbuna 12, 50009 Zaragoza, Spain Laboratorio Subterr´aneo de Canfranc, Paseo de los Ayerbe s.n., 22880 Canfranc Estaci´on, Huesca, Spain Present Address: Centro de Investigaciones Energ´eticas,Medioambientales y Tecnol´ogicas, CIEMAT, 28040, Madrid, Spain Fundaci´on ARAID, Av. de Ranillas 1D, 50018 Zaragoza, Spain Present Address: Fundaci´on CIRCE, 50018, Zaragoza, Spain
ANAIS is a direct detection dark matter experiment aiming at the testing of the DAMA/LIBRA annual modu-lation result, which standing for about two decades has neither been confirmed nor ruled out by any other exper-iment in a model independent way. ANAIS-112, consisting of 112.5 kg of sodium iodide crystals, is taking dataat the Canfranc Underground Laboratory, Spain, since August 2017. This letter presents the annual modulationanalysis of 1.5 years of data, amounting to 157.55 kg × y. We focus on the model independent analysis searchingfor modulation and the validation of our sensitivity prospects. ANAIS-112 data are consistent with the nullhypothesis (p-values of 0.65 and 0.16 for [ ] and [ ] keV energy regions, respectively). The best fits forthe modulation hypothesis are consistent with the absence of modulation ( S m = − . ± . − . ± . σ sensitivity to the DAMA/LIBRA resultin 5 years of data taking. There is overwhelming evidence from cosmological and as-trophysical observations supporting that a large fraction of theUniverse energy-mass budget is not explained in the frame-work of the standard model of the particle physics, assumingthe cosmological standard model [1]. The solution to the darkmatter (DM) / dark energy (DE) puzzle is probably of a com-plex nature. In one of the preferred hypothetical scenarios,DM can be understood as a new non-zero-mass particle hav-ing a very low interaction probability with baryonic matter.Although proposed candidates span about 45 orders of mag-nitude in mass, and 60 in cross-section, axions and WeaklyInteracting Massive Particles (WIMPs) are among the bettermotivated [2]. Experimental effort devoted to unraveling thenature of the DM particles has been spent, either by direct [3],indirect [4] or accelerator searches [5], which are complemen-tary to each other. Only one experiment, DAMA/LIBRA [6–8], has provided a long-standing positive result: the observa-tion of a highly statistically significant annual modulation inthe detection rate, compatible with that expected for galactichalo dark matter particles. This result has neither been re-produced by any other experiment, nor ruled out in a modelindependent way. Compatibility among the different exper-imental results in most conventional WIMP-DM scenarios isactually disfavored [9–19]. Then, a similar annual modulationsearch using the same target is mandatory to shed light on theDAMA/LIBRA conundrum, which is the goal of the ANAIS(Annual modulation with NaI Scintillators) experiment. Otherefforts sharing ANAIS goal in the international dark mattercommunity are the COSINE-100 experiment, taking data alsoin dark matter mode at Yang-Yang Underground Laboratory, ∗ Corresponding author: [email protected] † Deceased
South Korea [9, 20, 21]; and in longer-term, SABRE, aim-ing at installing twin detectors in Australia and Italy [22], andCOSINUS, developing cryogenic detectors based on NaI [23].An annual modulation in the dark matter interaction rate isexpected by the revolution of the Earth around the Sun, whichdistorts the DM particle velocity distribution function as seenby the detector, typically assumed Maxwellian boosted by theSun velocity [24, 25]. The effect is present unless the DM halois co-rotating with the Solar System. However, it is stronglydependent on the specific halo model, both in amplitude and inphase. It is natural to assume that the Sun is moving through alocally isotropic DM halo, with the Earth orbiting aside. Con-sequently, 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 analysisshould not assume a priori neither the period of the modula-tion nor the phase, but it should derive them from the data.A full and consistent analysis requires then several years ofmeasurement in very stable conditions. This is the long-termgoal of our experiment.ANAIS-112, consisting of 112.5 kg of NaI(Tl) detectors,was installed in 2017 at the Canfranc Underground Labora-tory, LSC, in Spain. The ANAIS-112 set-up undergoes a dif-ferent residual cosmic ray flux and environmental conditionsthan DAMA/LIBRA (800 m versus 1400 m rock overburden,for instance). Consequently, the potential confirmation of amodulation with same phase and amplitude would be very dif-ficultly explained as an effect of backgrounds or systematics.ANAIS-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 afive-years operation [28]; here we just briefly summarize themost relevant features of the experimental apparatus.ANAIS-112 uses nine NaI(Tl) modules produced by Alpha a r X i v : . [ a s t r o - ph . I M ] M a r Spectra Inc. in Colorado (US) ∗ . These modules have beenmanufactured from 2012 to 2017, and shipped to Spain avoid-ing air travel in order to prevent cosmogenic activation of themodule materials. Each crystal is cylindrical (4.75” diameterand 11.75” length), with a mass of 12.5 kg, and it is housedin OFE (Oxygen Free Electronic) copper. This encapsulationhas a Mylar window allowing low energy calibration using ex-ternal gamma sources. It incorporates two quartz optical win-dows to couple the photomultipliers (PMTs). All PMT units,as well as all relevant materials used in the building of the de-tectors, have been screened for radiopurity using HPGe detec-tors in the low background facilities at LSC. Their contribu-tion to the experiment background has been estimated [27, 29]and included in our background model (see below). We wouldlike to emphasize the outstanding light collection measuredfor the nine modules, at the level of 15 photoelectrons (phe)per keV [26, 30]. ANAIS-112 is calibrated every two weeksusing external Cd sources: all the nine modules are simul-taneously calibrated using a multi-source system which min-imizes down time periods. Background events from the de-cay of K and Na in the crystal bulk, associated to 3.2 and0.9 keV energy depositions, and selected by coincidence withan energy deposition in a second module around 1461 and1275 keV, respectively, are also used to improve the accuracyof the calibration down to the energy threshold [26].The ANAIS-112 shielding consists of 10 cm of archaeo-logical lead, 20 cm of low activity lead, an anti-radon box(continuously flushed with radon-free nitrogen gas), an ac-tive muon veto system made up of 16 plastic scintillators de-signed to cover top and sides of the whole ANAIS set-up and40 cm of neutron moderator (a combination of water tanksand polyethylene blocks). In the design of the muon veto sys-tem we follow a tagging strategy instead of a hardware veto-ing. The goal is twofold: on the one hand, to discard eventsin the NaI(Tl) crystals coincident with muon veto triggers. Onthe other hand, to analyse eventual correlations between muonhits in the plastic scintillators and events in the NaI(Tl) crys-tals, specially in the region of interest (ROI), [ ] keV † .The ANAIS-112 electronic chain and data acquisition sys-tem (DAQ) have been described in Refs. [26, 31]. Each PMTcharge signal is independently processed and divided into: (1)a trigger signal; (2) a low energy (LE) signal that goes to thedigitizers which sample the waveforms at 2 Gs/s with highresolution (14 bits); and (3) a high energy (HE) signal, conve-niently attenuated. The trigger of each PMT signal is done atphe level, while the single module trigger is done by the co-incidence (logical AND) of the two PMT triggers in a 200 nswindow. The global trigger is the logical OR of the nine mod-ules trigger signals. Trigger efficiency is close to 100% downto the analysis threshold established at 1 keV [26].ANAIS-112 started taking data in the DM mode on August,3 rd , 2017. It has accumulated almost 1.5 years of data-takingtime in quite stable conditions by February, 12 th , 2019. To- ∗ † Energy will be shown in electron equivalent units throughout this letter tal live time available for the annual modulation analysis is527.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 dueto the periodical calibration runs carried out using low energygamma sources. We remove events arriving less than one sec-ond after the last muon veto trigger, correcting the total livetime by subtracting one second per muon veto trigger, so thelive time used for the annual modulation analysis which fol-lows is 511.16 days.A blind analysis strategy in three levels has been followed:first, we calculate pulse parameters, the time since the lastmuon veto, and we apply the peak-finding algorithm to iden-tify individual phe in low energy pulses; second, we calibratethe energy response of every detector at LE and HE [26];third, we optimize the pulse shape cuts and calculate their ef-ficiency, being the LE variable hidden for events correspond-ing to single hits (M1 events). Only 10% of the data wereunblinded for general background assessment and fine tuningof procedures [26]; those data were randomly chosen (34 daysamounting to 32.9 days live time) and time evolution was kepthidden until the data unblinding, presented in this letter.Calibration procedures and efficiency corrections applied inthe following have been derived as done in Ref. [26] for thefirst year, continuing with same procedure in the half secondyear 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 frac-tion of the pulse area in [100-600] ns after the event trigger,defined following [32], and the logarithm of the mean time ofthe distribution of the individual phe arrival times in the digi-tized window [20]; events having an asymmetric light sharingbetween the two PMT signals are removed by imposing a cuton the number of peaks identified in each PMT (asymmetrycut). The total detection efficiency, ε ( E , d ) , calculated inde-pendently for every detector d as a function of the energy, E ,can be written [26] as ε ( E , d ) = ε trg ( E , d ) × ε PSA ( E , d ) × ε asy ( E , d ) , (1)where the trigger efficiency ε trg ( E , d ) is calculated fromMonte Carlo (MC) simulations, and the efficiencies of thepulse shape, ε PSA ( E , d ) , and asymmetry, ε asy ( E , d ) , cuts areevaluated from the 3.2 and 0.9 keV events selected by the co-incidence with the high energy gammas following K and Na decays, and
Cd calibration events, respectively, ac-cumulated for all the analysed exposure. Total detection ef-ficiency ranges from 0.15 to 0.35 at 1 keV, depending on thedetector, increases up to 0.8 at 2 keV and is nearly 1 at 4 keVfor all the modules. Statistical errors in the total efficiencyvary from 3 to 8% at 1 keV down to 1% at 6 keV. Comparingdifferent methods for the efficiency calculation we have alsoestimated a systematic uncertainty that amounts up to 20% at1-1.2 keV and is negligible above 1.5 keV. More details canbe found in [26].The background model for all the nine detectors used in theANAIS-112 set-up has been developed. It is based on MCsimulations using the measured activity in external compo-nents and in crystals, including cosmogenic products, quanti-fied in dedicated, independent measurements using differentanalysis techniques. It provides a good overall description ofmeasured data at all energy ranges above 2 keV and at differ-ent analysis conditions (coincidence or anticoincidence) [27].In the ROI the background is dominated by the emissions fromthe crystals themselves, in particular,
Pb (32.5%) and H(26.5%) continua, and K (12%) and Na (2.0%) peaks arethe most significant contributions. Short-lived isotopes acti-vated cosmogenically are still present in the bulk of the lastreceived crystals [33, 34], contributing as background in theROI, specially in the [3-5] keV region. However, from 1 to2 keV there is a large fraction of our background lacking fromexplanation [27]. It could have as origin non-bulk scintillationleaking through our event selection criteria.The time evolution of the rate of those events surviving allthe cuts during the first year and a half of ANAIS-112 op-eration is shown for different populations in panels a)-f) ofFigure 1. Data from all the modules have been added to-gether and corrected by the corresponding efficiencies. Thetwo lower panels correspond to two different energy windowsin the ROI: [ ] keV, a), and [3-5] keV, b), while the upperpanels show the evolution of control populations for which nomodulation is expected: c) [6-20] keV, d) double-coincidenceevents (M2) in energy region [ ] keV, and coincident eventsattributed to e) K, and f) Na decays in the crystals ‡ . Theuppermost panel, g), presents the evolution of muon-relatedlow-energy events before the cuts (M1 events in [ ] keV re-gion arriving less than 1 second after a muon veto trigger). Aspointed out in Ref. [26], many low energy events are regis-tered after a muon passage through the detector. They likelyoriginate in the long tail of the scintillation pulse produced bya large muon energy deposition, which is able to trigger manytimes the DAQ system and produce fake low-energy events.This population is clearly fluctuating in time and will be stud-ied in a future work.In panels a) and b) of Figure 1 we observe a relevant de-crease in the rate at the ROI, which amounts up to 8% in theanalysed period. It is caused by cosmogenically activated iso-topes and is well reproduced qualitatively in all energy win-dows by our background model [27], shown in Figure 1 ingreen dashed lines. It is normalized by a factor, f , to be moreeasily compared with the measured rates. Agreement in [3-5] keV region is also quantitatively very good, f =1.04, asthe background time evolution in this region is dominated byshort-lived cosmogenic isotopes. However, our backgroundmodel underestimates the rate in [ ] keV region ( f =1.28),because from 1 to 2 keV our background model does not re-produce the measurement [27], as commented above.We perform a least-squares fit on panels a) to d) to a con-stant term plus an exponential function, to account for thementioned background reduction, following our backgroundmodel [27]. The χ / NDF and the values of the fitted param- ‡ Low energy M2 events having an energy deposition in a second detector ina window around the corresponding high energy gamma R a t e ( c pd / k g / k e V ) days after August 3, 20170 200 400 600 a) M1 [1-6] keV rate = 2.94 + 0.87exp(-t/1217.8) (cpd/kg/keV)/NDF = 59.4 / 52 [pval=0.22] c f = 1.28 b) M1 [3-5] keV rate = 3.02 + 0.58exp(-t/420.9) (cpd/kg/keV)/NDF = 50.0 / 52 [pval=0.55] c f = 1.04 c) M1 [6-20] keV rate = 1.73 + 0.13exp(-t/313.7) (cpd/kg/keV)/NDF = 58.0 / 52 [pval=0.26] c d) M2 [1-6] keV rate = 0.25 + 0.05exp(-t/326.4) (cpd/kg/keV)/NDF = 48.6 / 52 [pval=0.61] c K (M2) [2-5] keV e) 0.001) (cpd/kg/keV) – rate = (0.091 /NDF = 59.4 / 54 [pval=0.29] c Na (M2) [0-2] keV f) rate = 0.12exp(-t/1583.6) (cpd/kg/keV)/NDF = 68.6 / 53 [pval=0.07] c -related events [1-6] keV m g) M1 FIG. 1: Rate of events corresponding to different populations alongthe first year and a half of ANAIS-112 operation, calculated in10-day binning. Events surviving all the filters are shown in pan-els a)-f): in energy region [ ] keV, a); [ ] keV, b); [ ] keV,c); M2 events in energy region [ ] keV, d); M2 events in energyregion [ ] keV in coincidence with a high energy gamma in a sec-ond detector in the K window, e); and M2 events in energy region [ ] keV in coincidence with a high energy gamma in a second de-tector in the Na window, f). Events arriving less than one secondafter a muon veto trigger before the cuts are shown in panel g). Fitsare shown as red solid lines, and corresponding fit parameters andchi-squared and p-values are also given in the plot. Green dashedlines correspond to the background model, normalized according toa factor, f , which is also given in the plots. eters are shown in the figure. Good fits are obtained, withp-values larger than 0.20 in all cases. The K (3.2 keV)and Na (0.9 keV) M2 populations, on the other hand, aremodelled by a constant and an exponential decay, respec-tively. We derive consistent results (p-value = 0.22) in thefirst case, and a reasonable agreement with the Na half-life( T / = . ± .
36 y, p-value = 0.07) in the second.In this letter we present our modulation results in two en-ergy regions: [ ] keV and [ ] keV, to allow direct com-parison with the DAMA/LIBRA results. The values of themodulation amplitude observed by DAMA/LIBRA, S DAMAm ,are 0 . ± . . ± . [ ] keV and using only phase-2 data for [ ] keV energy region, respectively [8]. We expect resultsderived from [ ] keV to be more robust because our data se-lection efficiencies strongly go down below 2 keV, increasingthe risk to be affected by unknown systematics.We evaluate the statistical significance of a possible modu-lation in our data by a least square method in the time-binneddata. The efficiency-corrected rate of events surviving the cutsin [ ] and [ ] keV energy regions is modelled as R ( t ) = R + R · exp ( − t / τ ) + S m · cos ( ω · ( t + φ )) , (2)where R and R are free parameters and τ is fixed to thevalue obtained from our background model in the correspond-ing energy range. We also fix the period ( ω = π /
365 d = . − ) and the phase ( φ = − . m is fixed to 0 for the null hypothesis andleft unconstrained (positive or negative) for the modulationhypothesis. This allows a direct comparison with the resultsfrom the DAMA/LIBRA analysis with 1 free parameter [8].We present the best fit for both hypothesis for 10-day timebinning in Figure 2. In order to highlight the presence or ab-sence of modulation, we plot the data with the constant andexponential terms subtracted. For the sake of comparison, inthe plot we show the modulation measured by DAMA/LIBRA(green lines). Days after August 3, 2017 /NDF = 62.0 / 52 [pval=0.16] c fi – mod hyp: Sm = (-0.0015 /NDF = 62.0 / 53 [pval=0.18] c fi null hyp DAMA mod hyp: Sm = 0.0105 (cpd/kg/keV) [1-6] keV Days after August 3, 2017 -0.100.10.2 /NDF = 47.4 / 52 [pval=0.65] c fi – mod hyp: Sm = (-0.0044 /NDF = 48.0 / 53 [pval=0.67] c fi null hyp DAMA mod hyp: Sm = 0.0102 (cpd/kg/keV) [2-6] keV c pd / k g / k e V FIG. 2: ANAIS-112 data in the energy windows [ ] keV (bottompanel) and [ ] keV (top panel) surviving all the cuts and efficiencycorrected [26]. Data is displayed after subtracting the constant andexponential functions fitted to Equation 2. Fits are also shown in thesame way, both in the modulation (3 free parameters) and the nullhypothesis (2 free parameters). χ and p-values displayed allow thecomparison of both hypothesis, and DAMA/LIBRA results on mod-ulation amplitude in both energy windows are shown in green [8]. In both energy regions the null hypothesis is well sup- ported by the χ test, with χ / NDF = . /
53 for the [ ] keV (p-value = 0.67) and χ / NDF = . /
53 for the [ ] keV regions (p-value = 0.18). The best fits for the modu-lation hypothesis are S m = − . ± . − . ± . [ ] keV and [ ] keV,respectively. In both cases, p-values are slightly lowerthan those of the null hypothesis (0.65 and 0.16, respec-tively). The best fits are incompatible at 2.5 σ (1.9 σ ) withthe DAMA/LIBRA signal.The statistical significance of our result is determined bythe standard deviation of the modulation amplitude distribu-tion, σ ( S m ) , which would be obtained in a large number of ex-periments like ANAIS-112 with the present exposure. Then,we quote our sensitivity to DAMA/LIBRA result as the ra-tio S DAMAm / σ ( S m ) , which directly gives in σ units the C.L.at which we can test the DAMA/LIBRA signal. At present,our result σ ( S m ) = . ( . ) cpd/kg/keV for [ ] keV( [ ] keV) corresponds to a sensitivity of 1.75 σ (1.66 σ ) tothe DAMA/LIBRA signal. In Ref. [28] we found an analyti-cal expresion to calculate σ ( S m ) at a given exposure from themeasured background and detection efficiency. Figure 3 (darkblue lines) displays our sensitivity projection calculated fol-lowing Ref. [28] for the two studied energy ranges, whereasthe blue bands represent the 68% uncertainty in S DAMAm as re-ported in Ref. [8]. In the calculation we take into account theANAIS-112 live time distribution, the background reductionexpected due to decaying isotopes and the statistical error inthe detection efficiency. The black dots are the sensitivitiesderived in this work, including a systematic error estimatedby changing the time-binning from 1 to 20 days, and consid-ering the systematics in the efficiency [26]. The results per-fectly agree with our estimates, confirming the ANAIS-112projected sensitivity to the DAMA/LIBRA result. A 3 σ sen-sitivity should be at reach in 4-5 years of data-taking. real time (y) real time (y) ) s C . L . ( FIG. 3: ANAIS-112 sensitivity to the DAMA/LIBRA signal in σ C.L. units (see text) as a function of real time in the [ ] keV (upperpanel) and [ ] keV (lower panel) energy regions. The black dotsare the sensitivities derived in this work, σ ( S m ) . The blue bandsrepresent the 68% C.L. DAMA/LIBRA uncertainty [8]. Finally, Figure 4 presents the best fit amplitudes, S m , cal-culated per 1 keV energy bins, from 1 to 20 keV (bottompanel, black dots), together with the DAMA-phase-2 modu-lation amplitudes extracted from Ref. [8] (blue triangles). Thetop panel shows the p-values for the null (open squares) andmodulation hypothesis (closed circles) for every energy bin.All the modulation amplitudes in the ROI are compatible with0 and, in general, p-values for the null hypothesis are slightlylarger than those of the modulation hypothesis. The 1 σ and2 σ bands shown in the figure are obtained following Ref. [28]for the present ANAIS-112 exposure.In summary, to test the DAMA/LIBRA annual modulationresult in a model independent way, an analysis of the first yearand half of data from ANAIS-112 experiment has been pre-sented. It is compatible with the sensitivity estimates donein Ref. [28], and confirms our goal of reaching sensitivity at3 σ level in five years (see Figure 3) to DAMA/LIBRA re-sult. We derive best fits for the modulation amplitude S m = − . ± . − . ± . [ ] and [ ] keV energy regions, respectively, compatiblewith the absence of modulation.This work has been financially supported by the Span-ish Ministerio de Econom´ıa y Competitividad and the Eu-ropean Regional Development Fund (MINECO-FEDER) un-der grants No. FPA2014-55986-P and FPA2017-83133-P,the Consolider-Ingenio 2010 Programme under grants Mul- energy (keV) S m ( c pd / k g / k e V ) -0.04-0.0200.020.04 s s p - v a l ue null hypmodulation hyp FIG. 4: Modulation amplitude per 1 keV energy bins combiningdata from all the modules. We show for reference the correspond-ing DAMA/LIBRA result [8], and the 1 σ and 2 σ ANAIS-112 bandsfollowing the analysis done in Ref. [28]. Top panel compares the p-values of the fits to Equation 2 with those corresponding to the nullhypothesis for every energy bin. tiDark CSD2009-00064 and CPAN CSD2007-00042 and theGobierno de Arag´on and the European Social Fund (Groupin Nuclear and Astroparticle Physics and I. Coarasa predoc-toral grant). We thank the support of the Spanish Red Con-solider MultiDark FPA2017-90566-REDC. We acknowledgethe technical support from LSC and GIFNA staff. ProfessorJ.A. Villar passed away in August, 2017. Deeply in sorrow,we all thank him for his dedicated work and kindness. [1] G. Bertone and D. Hooper, Rev. Mod. Phys. , 045002 (2018),arXiv:1605.04909.[2] L. Baudis, J. Phys. G43 , 044001 (2016).[3] J. Liu, X. Chen, and X. Ji, Nature Phys. , 212 (2017),arXiv:1709.00688.[4] J. Conrad and O. Reimer, Nature Phys. , 224 (2017),arXiv:1705.11165.[5] O. Buchmueller, C. Doglioni, and L. T. Wang, Nature Phys. ,217 (2017).[6] R. Bernabei et al. (DAMA), Eur. Phys. J. C56 , 333 (2008),arXiv:0804.2741.[7] R. Bernabei et al., Eur. Phys. J.
C73 , 2648 (2013),arXiv:1308.5109.[8] R. Bernabei et al., Universe , 116 (2018), [At. En-erg.19,307(2018)], arXiv:1805.10486.[9] G. Adhikari et al., Nature , 83 (2018).[10] M. Kobayashi et al. (XMASS) (2018), arXiv:1808.06177.[11] D. S. Akerib et al. (LUX), Phys. Rev. D98 , 062005 (2018),arXiv:1807.07113.[12] K. Abe et al. (XMASS), Phys. Rev.
D97 , 102006 (2018),arXiv:1801.10096.[13] E. Aprile et al. (XENON), Phys. Rev. Lett. , 101101 (2017),arXiv:1701.00769.[14] C. Savage, K. Freese, P. Gondolo, and D. Spolyar, JCAP ,036 (2009), arXiv:0901.2713.[15] E. Aprile et al. (XENON100), Science , 851 (2015),arXiv:1507.07747.[16] J. Herrero-Garcia, JCAP , 012 (2015), arXiv:1506.03503.[17] S. Baum, K. Freese, and C. Kelso, Phys. Lett.
B789 , 262 (2019),arXiv:1804.01231.[18] S. Kang, S. Scopel, G. Tomar, and J.-H. Yoon, JCAP , 016(2018), arXiv:1804.07528.[19] J. Herrero-Garcia, A. Scaffidi, M. White, and A. G. Williams,Phys. Rev.
D98 , 123007 (2018), arXiv:1804.08437.[20] K. W. Kim et al. (KIMS), arXiv:1806.06499 (2018).[21] G. Adhikari et al., Eur. Phys. J.
C78 , 107 (2018),arXiv:1710.05299.[22] C. Tomei (SABRE), Nucl. Instrum. Meth.
A845 , 418 (2017).[23] A. G¨utlein et al., Nucl. Instrum. Meth.
A845 , 359 (2017).[24] K. Freese, J. A. Frieman, and A. Gould, Phys. Rev.
D37 , 3388(1988).[25] K. Freese, M. Lisanti, and C. Savage, Rev. Mod. Phys. , 1561(2013), arXiv:1209.3339.[26] J. Amar´e et al., accepted for publication in Eur. Phys. J. C(2018), arXiv:1812.01472.[27] J. Amar´e et al., submitted to Eur. Phys. J. C (2018),arXiv:1812.01377.[28] I. Coarasa et al., accepted for publication in Eur. Phys. J. C(2018), arXiv:1812.02000. [29] J. Amare et al., Eur. Phys. J. C76 , 429 (2016),arXiv:1604.05587.[30] M. A. Oliv´an et al., Astropart. Phys. , 86 (2017),arXiv:1703.01262.[31] M. Oliv´an, Ph.D. thesis, Universidad de Zaragoza (2016). [32] R. Bernabei et al., Nuclear Instruments and Methods in PhysicsResearch A , 297 (2008), arXiv:0804.2738.[33] P. Villar et al., Int. J. Mod. Phys. A33 , 1843006 (2018).[34] J. Amar´e et al., JCAP1502