Direct dark matter search by annual modulation in XMASS-I

A search for dark matter was conducted by looking for an annual modulation signal due to the Earth's rotation around the Sun using XMASS, a single phase liquid xenon detector. The data used for this analysis was 359.2 live days times 832 kg of exposure accumulated between November 2013 and March 2015. When we assume Weakly Interacting Massive Particle (WIMP) dark matter elastically scattering on the target nuclei, the exclusion upper limit of the WIMP-nucleon cross section 4.3$\times$10$^{-41}$cm$^2$ at 8 GeV/c$^2$ was obtained and we exclude almost all the DAMA/LIBRA allowed region in the 6 to 16 GeV/c$^2$ range at $\sim$10$^{-40}$cm$^2$. The result of a simple modulation analysis, without assuming any specific dark matter model but including electron/$\gamma$ events, showed a slight negative amplitude. The $p$-values obtained with two independent analyses are 0.014 and 0.068 for null hypothesis, respectively. we obtained 90\% C.L. upper bounds that can be used to test various models. This is the first extensive annual modulation search probing this region with an exposure comparable to DAMA/LIBRA.


Introduction
There is strong evidence that about 5 times more dark matter exists in the universe than ordinary matter. Despite its prominence, we do not yet know what dark matter is [1]. Among many candidates for dark matter particles, WIMPs are well motivated and have received the most attention to date. However, collider experiments at the LHC do not show any indication for such particles so far [1]. And no experimental indication for a standard WIMP was found in high sensitivity direct search experiments such as LUX [2], XENON100 [3] and SuperCDMS [4] either. On the other hand, that appears to contradict experiments that report signals interpreted as ∼ 10 GeV/c 2 light WIMP dark matter [8,9,10] for many years. In this situation, light mass WIMPs or other dark matter candidates are getting more attention. In fact, XMASS, a high light yield and low background detector, probed this possibility and looked for signals not only from nuclear recoils but also from electrons and gamma rays emanating from interactions of other candidates such as axion-like particles, Super-WIMPs and so on [5,6,7].
The most significant result is that of the DAMA/LIBRA experiment at the Gran Sasso National Laboratory in Italy which indicated an annual modulation signature [11]. The Earth's velocity relative to the dark matter distribution changes as the Earth moves around the Sun and produces such a modulation in the dark matter signal rate. This modulation can be observed with terrestrial detectors [12]. The amplitude of the modulation can be changed from positive (i.e. higher rate in June than in December) to negative at cross-over energy [13] and it is possible to observe this effect if the detector threshold is lower than that energy. For 100 GeV/c 2 WIMP mass and a Xe target, this is about 20 keV nuclear recoil energy and it depends on the WIMP mass and the target materials.
The DAMA/LIBRA experiment reported an observation of event rate annual modulation with a 9σ significance in 1.33 ton·year of data taken over 14 annual cycles with 100 to 250 kg of NaI(Tl) detectors. Their signal may be caused by light WIMPs, or other types of dark matter producing electrons or gamma rays. In such cases, the signal is not observable to direct search experiments if they remove electron events. In this situation, dark matter models, for instance, with interaction via dark matter-electron scattering become well motivated which produce keV energy deposition in the detector because they provide a explanation for the DAMA/LIBRA result while avoiding other direct detection constraints [14,15,16]. Recently, in addition to the WIMP search result [3], an annual modulation search was carried out by the XENON group using only electronic recoil events in their two phase Xe detector with the 34 kg fiducial volume in 224.6 live days data [17]. The result disfavored the interpretation of the DAMA/LIBRA as WIMPelectron scattering through axial-vector coupling. XMASS uses a single phase technology to observe only scintillation light by looking for both types of signals without any electric field. Although XMASS has a modest background rate like that of DAMA/LIBRA, XMASS has a larger mass of 832 kg of liquid xenon and, therefore, is able to reach the DAMA/LIBRA exposure in short time. While the background in this recent modulation study by the XENON experiment is lower, XMASS has a larger target mass and significantly longer exposure time. We will discuss the sensitivity later. Note that XMASS tests this modulation hypothesis with almost half the energy threshold (∼ 1keV) than theirs in a different environment and underground site.

The XMASS experiment
The XMASS detector is located at the Kamioka Observatory (overburden 2700 m.w.e) in Japan. The detailed design and performance are described in [18]. The detector is immersed in a water tank, 10 m in diameter and 10.5 m in height, which is equipped with 72 Hamamatsu H3600 photomultiplier tubes (PMTs), and acts as an active muon veto and a passive radiation shield against neutrons and gamma rays from the surrounding rock. 642 high quantum efficiency (28-40% at 175 nm) Hamamatsu R10789 PMTs are mounted in the liquid xenon detector, an approximate sphere with an average radius of 40 cm. The gain of the PMTs was monitored weekly with a blue LED embedded in the inner surface of the detector. The scintillation light yield response was traced by inserting a 57 Co source [19] into the detector every one or two weeks. The number of events for each source position was about 20,000.
In November 2013, after refurbishing the detector to reduce the radioactive background from the aluminum seal of the PMTs' window that was identified in the commissioning run [18], data taking was resumed with about one order of magnitude improved background by covering these seal parts with plates made of pure copper. The data accumulated between November 2013 and March 2015 were used for this analysis and we selected periods with stable temperature (172.6-173.0 K) and pressure of Xe (0.162 -0.164 MPa absolute). After removing periods of operation with excessive PMT noise or data acquisition problems, the total live time became 359.2 days.
In this paper, two different energy scales were used: 1) keV ee represents an electron equivalent energy incorporating all the gamma-ray calibrations in the energy range between 5.9 keV and 122 keV from 55 Fe, 109 Cd, 241 Am and 57 Co sources by inserting sources into the sensitive volume of the detector. The non-linearity of energy scale was taken into account with those calibrations using a non-linearity model from Doke et al. [20]. Below 5.9 keV, we extrapolated based on this model. We found about 15% energy scale difference from the Noble Element Simulation Technique (NEST) [21] at the threshold energy of 1.1 keV ee (∼8 photoelectrons) in this analysis. 2) keV nr denotes the nuclear recoil energy which is estimated from the light yield at 122 keV by using non-linearity response measurement at zero electric field in [22]. The energy threshold, in this case, corresponds to 4.8 keV nr .

Data Analysis
Events with 4 or more PMT hits in a 200 ns coincidence timing window without a muon veto were initially selected. This resulted in 3.3×10 7 events in the energy region between 1.1 and 15 keV ee . In order to avoid events caused by afterpulses of bright events induced by, for example, high energy gamma-rays or alpha particles, we rejected events occurring within 10 ms from the previous event and having a variance in their hit timings of greater than 100 ns (this selection reduces the number of events to 2.8×10 7 ). A 'Cherenkov cut' removed events which produce light predominantly from Cherenkov emission, in particular from the beta decays of 40 K in the PMT photocathode. Events for which more than 60% of their PMT hits arrive in the first 20 ns were classified as Cherenkov-like events [5] (this selection reduces the number of events to 1.9×10 6 ). Finally, to remove background events that occurred in front of PMT window, we give upper limits on the values of 'Max-photoelectron/Total-photoelectron' where Max-photoelectron and Total-photoelectron are the largest photoelectron counts in one PMT among all PMTs and the total number of photoelectrons in the event, respectively (this selection reduces the number of events to 3.6×10 5 ). These cut values varied as a function of photoelectron from about 0.2 at 8 photoelectrons to about 0.07 at 50 photoelectrons. The count rate for the data after all the cuts is 1.17 (0.028) events/day/kg/keV ee at 1.1 (5.0) keV ee .
The 57 Co calibration data were taken at from z = −40 cm to +40 cm along the center vertical axis of the detector to track photoelectron yield and optical properties of the liquid xenon [18]. A difference of about 10% was observed as the position dependence for this photoelectron yield. The photoelectron yield during the data taking varied about 10%. The absorption and scattering length for the scintillation light as well as the intrinsic light yield of the liquid xenon scintillator are extracted from the 57 Co calibration data the Monte Carlo simulation [18]. With that we found that we can trace the observed photoelectron change in the calibration data as a change as the absorption length, while the scattering length remains stable at 52 cm with a standard deviation of ±0.6%. We then re-evaluate the absorption length and the relative intrinsic light yield to see the stability of the scintillation light response by fixing the scattering length at 52 cm. The absolute absorption length varied from about 4 m to 11 m, but the relative intrinsic light yield (R yield ) stayed within ±0.6% over the entire data taking period. The time dependence of the photoelectron yield affects the efficiency of the cuts. Therefore, we evaluate the absorption length dependence of the relative cut efficiencies through Monte Carlo simulation. If we normalize the overall efficiency at an absorption length of 8 m, this efficiency changes from −4% to +2% over the relevant absorption range. The position dependence of the efficiency was taken into account as a correlated systematic error (∼ ±2.5%). This is the dominant systematic uncertainty in the present analysis. The second largest contribution comes from a gain instability of the waveform digitizers (CAEN V1751) between April 2014 and September 2014 due to a different calibration method of the digitizers used in that period. This effect contributes an uncertainty of 0.3% to the energy scale. Other effects from LED calibration, trigger threshold stability, timing calibration were negligible. The observed count rate after cuts as a function of time in the energy region between 1.1 and 1.6 keV ee is shown in Fig. 2. The systematic errors caused by the relative cut efficiencies are also shown.  To retrieve the annual modulation amplitude from the data, the least squares method for the time-binned data was used. The data set was divided into 40 time-bins (t bins ) with roughly 10 days of live time each. The data in each time-bin were then further divided into energy-bins (E bins ) with a width of 0.5 keV ee . Two fitting methods were performed independently. Both of them fit all energy-and time-bins simultaneously. Method 1 used a 'pull term' α with χ 2 defined as: where R data i, j , R ex i, j , σ(stat) i,j and σ(sys) i,j are data, expected event rate, statistical and systematic error, respectively, of the (i-th energy-and j-th time-) bin. The time is denoted as the number of days from January 1, 2014. K i, j represents the 1σ correlated systematic error on the expected event rate based on the relative cut efficiency in that bin. Method 2 used a covariance matrix to propagate the effects of the systematic error. Its χ 2 was defined as: where N bins (= E bins × t bins ) was the total number of bins and R data(ex) k is the event rate where k = i · t bins + j. The matrix V stat contains the statistical uncertainties of the bins, and V sys is the covariance matrix of the systematic uncertainties as derived from the relative cut efficiency.

Results and Discussion
We performed two analyses, one assuming WIMP interactions and the other independent of any specific dark matter model. Hereafter we call the former case the WIMP analysis and the latter a model independent analysis.
In the case of the WIMP analysis, the expected modulation amplitudes become a function of the WIMP mass A i (m χ ) as the WIMP mass m χ determines the recoil energy spectrum. The expected rate in a bin then becomes: where σ χn is the WIMP-nucleon cross section. To obtain the WIMP-nucleon cross section the data was fitted in the energy range of 1.1-15 keV ee . We assume a standard spherical isothermal galactic halo model with the most probable speed of v 0 =220 km/s, the Earth's velocity relative to the dark matter distribution of v E = 232+15 sin2π(t−t 0 )/T km/s, and a galactic escape velocity of v esc = 650 km/s, a local dark matter density of 0.3 GeV/cm 3 , following [13]. In the analysis, the signal efficiencies for each WIMP mass are estimated from Monte Carlo simulation of uniformly distributed nuclear recoil events in the liquid xenon volume. The systematic error of the efficiencies comes from the uncertainty of liquid xenon scintillation decay time of 25±1 ns [5] and is estimated as about 5% in this analysis. The expected count rate for WIMP masses of 7 and 8 GeV/c 2 with a cross section of 2×10 −40 cm 2 for the spin independent case are shown in Fig. 2 as a function of time after all cuts. This demonstrates the high sensitivity of the XMASS detector to modulation. As both methods found no significant signal, the 90% C.L. upper limit by the 'pull term' method on the WIMP-nucleon cross section is shown in Fig. 3. The exclusion upper limit of 4.3×10 −41 cm 2 at 8 GeV/c 2 was obtained. The −1σ scintillation efficiency of [22] was used to obtain a conservative limit. To evaluate the sensitivity of WIMP-nucleon cross section, we carried out a statistical test by applying the same analysis to 10,000 dummy samples with the same statistical and systematic errors as data but without modulation by the following procedure. At first, the time-averaged energy spectrum was obtained from the observed data. Then, we performed a toy Monte Carlo simulation to simulate time variation of event rate of background at each energy bin assuming the same live time as data and including systematic uncertainties. The ±1σ and ±2σ bands in Fig. 3 outline the expected 90% C.L. upper limit band for the no-modulation hypothesis using the dummy samples. The result excludes the DAMA/LIBRA allowed region as interpreted in [8] for the WIMP masses higher than 8 GeV/c 2 . The difference between two fitting methods is less than 10%. The upper limit of 5.4×10 −41 cm 2 is obtained under different astrophysical assumptions of v esc = 544 km/s [24]. The best fit parameters in a mass range between 6 and 1000 GeV/c 2 is a cross section of 3.2×10 −42 cm 2 for a WIMP mass of 140 GeV/c 2 . This yields a statistical significance of 2.7σ, however, in this case, the expected unmodulated event rate exceeds the total observed event rate by a factor of 2, therefore these parameters were deemed unphysical. bands represent the expected 90% exclusion distributions. Limits as well as allowed regions from other searches based on counting method are also shown [2,3,23,8,9,10,5].
For the model independent analysis, the expected event rate was estimated as: where the free parameters C i and A i were the unmodulated event rate and the modulation amplitude, respectively. t 0 and T were the phase and period of the modulation, and t j and ∆t j was the time-bin's center and width, respectively. In the fitting procedure, the 1.1-7.6 keV ee energy range was used and the modulation period T was fixed to one year and the phase t 0 to 152.5 days (∼2nd of June) when the Earth's velocity relative to the dark matter distribution is expected to be maximal. Figure 4 shows the best fit amplitudes as a function of energy for 'pull term' after correcting the efficiency. The efficiency was evaluated from gamma ray Monte Carlo simulation with a flat energy spectrum uniformly distributed in the sensitive volume (Fig. 4 inset). Both methods are in good agreement and find a slight negative amplitude below 4 keV ee . The ±1σ and ±2σ bands in Fig. 4 represent expected amplitude coverage derived from same dummy sample above by the 'pull term' method. This test gave a p-value of 0.014 (2.5σ) for the 'pull term' method and of 0.068 (1.8σ) for the covariance matrix method.
To be able to test any model of dark matter, we evaluated the constraints on the positive and negative amplitude separately in Fig. 4. The upper limits on the amplitudes in each energy bin were calculated by considering only regions of positive or negative amplitude. They were calculated by integrating Gaussian distributions based on the mean and sigma of data (=G(a)) from zero. The positive or negative upper limits are satisfied with 0.9 for , where a and a up are the amplitude and its 90% C.L. upper limit, respectively. The 'pull term' method obtained positive (negative) upper limit of 2.1(−2.1)×10 −2 events/day/kg/keV ee between 1.1 and 1.6 keV ee and the limits become stricter at higher energy. The energy resolution (σ/E) at 1.0 (5.0) keV ee is estimated to be 36% (19%) comparing gamma ray calibrations and its Monte Carlo simulation. As a guideline, we make direct comparisons with other experiments not by considering a specific dark matter model but amplitude count rate. The modulation amplitude of ∼ 2 × 10 −2 events/day/kg/keV ee between 2.0 and 3.5 keV ee was obtained by DAMA/LIBRA [11] and we estimate a 90% C.L. upper limit for XENON100 as 3.7 × 10 −3 events/day/kg/keV ee (2.0-5.8 keV ee ) based on [17] as it was not claimed as a signal. XMASS obtained positive upper limits of (1.7 − 3.7) × 10 −3 events/day/kg/keV ee in same energy region and gives the more stringent constraint. This fact is important when we test the dark matter model.

Conclusions
In conclusion, XMASS with its large exposure and high photoelectron yield (low energy threshold) conducted an annual modulation search. For the WIMP analysis, the exclusion upper limit of 4.3×10 −41 cm 2 at 8 GeV/c 2 was obtained and the result excludes the DAMA/LIBRA allowed region for WIMP masses higher than that. In the case of the model independent case, the analysis was carried out from the energy threshold of 1.1 keV ee which is lower than DAMA/LIBRA and XENON100. The positive (negative) upper limit amplitude of 2.1 (−2.1) × 10 −2 events/day/kg/keV ee between 1.1 and 1.6 keV ee and (1.7 − 3.7)×10 −3 counts/day/kg/keV ee between 2 and 6 keV ee were obtained. As this analysis does not consider only nuclear recoils, a simple electron or gamma ray interpretation of the DAMA/LIBRA signal can also obey this limit.