Search for sub-GeV dark matter by annual modulation using XMASS-I detector

A search for dark matter (DM) with mass in the sub-GeV region (0.32-1 GeV) was conducted by looking for an annual modulation signal in XMASS, a single-phase liquid xenon detector. Inelastic nuclear scattering accompanied by bremsstrahlung emission was used to search down to an electron equivalent energy of 1 keV. The data used had a live time of 2.8 years (3.5 years in calendar time), resulting in a total exposure of 2.38 ton-years. No significant modulation signal was observed and 90% confidence level upper limits of $1.6 \times 10^{-33}$ cm$^2$ at 0.5 GeV was set for the DM-nucleon cross section. This is the first experimental result of a search for DM mediated by the bremsstrahlung effect. In addition, a search for DM with mass in the multi-GeV region (4-20 GeV) was conducted with a lower energy threshold than previous analysis of XMASS. Elastic nuclear scattering was used to search down to a nuclear recoil equivalent energy of 2.3 keV, and upper limits of 2.9 $\times$10$^{-42}$ cm$^2$ at 8 GeV was obtained.


Introduction
The nature of dark matter (DM) is a key mystery in cosmology, and detecting it via any force other than gravity is essential for advancing particle physics beyond the standard model. Weakly interacting massive particles (WIMPs) at O(100 GeV) are predicted by theoretical extensions of the standard model, such as the constrained minimal supersymmetric standard model and are strong DM candidates [1]. They have been investigated extensively via nuclear recoil [2,3,4]; however, no significant detections of WIMPs have been confirmed.
Other theories predict a myriad of different DM types, lightmass WIMPs [5], asymmetric DM [6,7,8], or hidden sector DM [9] and many others; the mass of these DM candidates ranges from sub-GeV to a few GeV. Semi-conductor and crystal detectors have searched for these light DM candidates by lowering their nuclear recoil energy thresholds [10,11]. A search via DM-electron scattering by existing detectors have also been performed [12,13]. In addition to these detectors, conventional xenon detectors should also be sensitive to DM with sub-GeV mass [14,15], due to the irreducible contribution of the bremsstrahlung effect accompanying nuclear recoils [14]. The bremsstrahlung effect can occur when DM collides with a nucleus causing it to recoil and accelerate. In the case that a mass of DM particle is 1 GeV, the energy deposited by the bremsstrahlung photon is at most 3 keV. This energy is considerably more than that deposited by elastic nuclear recoil (∼0.1 keV).
In addition to this bremsstrahlung effect, another inelastic effect called the Migdal effect has also been suggested [16]. This effect leads to the emission of an electron from the atomic shell and causes subsequent radiation through the inelastic recoil of DM and nuclei. Although the bremsstrahlung and Migdal effects need both be calibrated experimentally in xenon and cross sections are smaller than that of elastic nuclear recoil (∼10 −6 for Migdal, ∼10 −8 for Bremsstrahlung at 1 GeV), because these inelastic effects lead to larger energy deposition than elastic nuclear recoil, it should be possible to detect sub-GeV DM through these effects.
Moreover, searching for a spin-dependent (SD) interaction utilising these effects is an attractive possibility, since xenon has a larger fraction of odd isotopes than that of other isotopes, such as oxygen [11]. Xenon has two stable odd isotopes, namely 129 Xe and 131 Xe which account for 26.4% and 21.2% of the natural xenon abundance, respectively; oxygen has only 0.04% of odd isotopes. Further theoretical studies are expected to enable the quantitative interpretation of the SD interaction by sub-GeV DM.
This letter reports on the first experimental search for sub-GeV DM (0.32-1.0 GeV) utilizing the bremsstrahlung effect. In the case of xenon, the Migdal effect is accompanied by Mshell electron emission, and the most likely de-excitation energy is 0.66 keV from the 3d orbit. As discussed in Section 4, since our understanding of detector responses is limited to those greater than 1 keV, we focus only on for the signal from the bremsstrahlung effect in this analysis. On the other hand, the search for multi-GeV DM (4-20 GeV) via conventional elastic nuclear recoils [17,18] was performed. For multi-GeV DM search, data with lower energy threshold than in previous studies [17,18] were used to improve sensitivity in the low mass range. These searches were conducted by looking for the annual modulation of the event rate in the XMASS data.

Expected annual modulation of signal
The annual modulation of the bremsstrahlung signal from the sub-GeV DM is evaluated by following the study in [14]. The differential cross section for such a process is where ω is the bremsstrahlung photon energy, α is the fine structure constant, f (ω) represents atomic scattering factor, µ N is the DM-nucleus reduced mass, v = |v| is the absolute value of the relative velocity between DM and the target v, m N is the nucleus mass, µ N is the DM-nucleus reduced mass, σ S I 0 ≃ A 2 σ n (µ N /µ n ) 2 is the spin-independent DM-nucleus cross section in which σ n is the DM-nucleon elastic cross section, µ n is the DM-nucleon reduced mass, and A is the atomic mass number. The cross section of bremsstrahlung effect is suppressed by the factor of α m 2 N from that of elastic nuclear recoil. The corresponding differential event rate is where N T is the number of target nuclei per unit mass in the detector, ρ χ = 0.3 GeV cm −3 is the local DM mass density [19], m χ is the DM mass, v E is the velocity of the Earth relative to the galactic rest frame. f v (v) is the DM velocity distribution in the galactic frame. It is assumed to be a truncated Maxwellian distribution with escape speed v esc = 544 km/s, most-probable velocity v 0 = 220 km/s and minimum velocity v min = 2ω/µ N [20]. Assuming that the relative velocity between DM and detector varies as {232 + 15 sin 2π(t − φ)/T } km/s [21], in which the phase φ = 152.5 days [19] from January 1st and period T = 365.24 days, we calculated the event rate as a function of bremsstrahlung energy and time. Figure 1 shows the expected bremsstrahlung spectra for 0.5 GeV DM at June and December corresponding to the maximum and minimum v E , respectively, as well as the averaged spectrum. The expected modulation amplitude is about 30% of the average event rate at 1 keV before considering the effect of the detector such as energy nonlinearity or resolution.
The annual modulation in the conventional nuclear recoil signal caused by DM has also been discussed as in [20]. To evaluate the amplitude for this signal, the same calculation in the previous analysis by XMASS was performed [17,18].

XMASS Experiment
The XMASS-I detector is a single-phase liquid xenon (LXe) detector located underground (2,700 meter water equivalent) at the Kamioka Observatory in Japan [22]. The inner detector contains 832 kg of xenon and has a pentakis-dodecahedron structure made of copper that supports 642 Hamamatsu R10789 photomultiplier tubes (PMTs). The quantum efficiency of the R10789 at room temperature is ∼30%. The PMTs cover more than 62% of the inner surface resulting in a large number of photoelectrons per keV detected by the PMTs (PE yield), as it is ∼15 PE/keV for 122 keV γ ray with zero electric field. Here, one PE is defined as the average PE observed at one photon incident to correct for the double PE emission from a PMT in the case of the xenon scintillation [23]. Signals from PMTs are recorded by waveform digitizers (CAEN V1751) with 1 GHz sampling rate. To shield the detector from external neutrons and γ-rays while also providing a muon veto, XMASS-I sits at the centre of a cylindrical water-Cherenkov detector. The Cherenkov detector is 10.5 m in height, 10 m in diameter and has 72 Hamamatsu H3600 PMTs arranged on the inside of its wall. This work used the data collected between November 20, 2013 and June 20, 2017. The xenon was required to maintain a stable operational temperature and pressure. A detailed plot of the LXe temperature and pressure during the first 2.7 years of this dataset are shown in [18], and the values were kept consistently within 0.05 K and 0.2 kPa in the following year. Periods with the problem of data acquisition system or electronics, such as excessive PMT noise, or unstable pedestal levels were removed from the dataset. The dataset has a total live time of 2.8 years, and the exposure is 2.38 ton-years. In addition to this data set, data with a lower energy threshold has also been taken since December 8, 2015. This data, referred to as low threshold data has 0.63 ton-year of exposure, and is used only for multi-GeV analysis. Details are discussed in section 6. In Fig. 2, observed data and simulated signal for bremsstrahlung and nuclear recoils are shown.

Calibration
The gain of each PMT was monitored by measuring single PE using a blue LED attached to the inner surface of the detector. This LED is flashed once per second, and gain of each PMT was calculated based on the weekly averaged LED data. The PE yield was tracked by inserting a 57 Co source into the detector every one or two weeks. These calibration processes are described in detail in [17,18,22,24]. The PE yield, absorption and scattering length for the scintillation light as well as the number of generated LXe scintillation photons per keV (light yield), are evaluated from the 57 Co calibration data with the help of a Monte Carlo (MC) simulation. In the simulation, two PE emissions are also taken into account. The variation in PE yield can be explained by changes of the absorption length in the LXe [18]. To reduce this change of PE yield, xenon gas has been purified continuously by circulating through hot metal getters since March 2015. The standard deviation of the PE yield was ±2.4% and ±0.5% before and after the circulation has been started, respectively. In this letter, two different energy scales, "keV ee " and "keV nr ", are used to indicate the electron-equivalent energy and nuclear recoil energy, respectively. These are different from those used in the previous analysis [17,18] below 5.9 keV ee and 3 keV nr as new calibrations were performed in this low energy region as explained as follow.
For the electron-equivalent energy, the non-linearity of the light yield (scintillation efficiency) along energy was taken into account using the model from Doke et al. [25] with corrections based on the result of calibration. The scintillation efficiency below 5.9 keV was calibrated using the L-shell X-ray escape peaks measured during calibration with an 55 Fe source. These escape peaks distribute energy in 1.2-2 keV, and the weighted mean energy of these escape peaks was 1.65 keV. Figure 3 shows the distribution of the number of PEs for the escape peak. The scintillation efficiency at 1.65 keV was evaluated by comparing these escape peaks in the data (solid blue line) and total MC (solid red line) considering systematic uncertainties such as the source assembly with its shadowing and reflection effects, trigger efficiency, the choice of fitting functions. The dashed red histogram represents the PE distribution only for escape peaks, whereas the green line represents the PE distribution for tail component from 5.9 keV peak, which was caused by the shadowing effect of the calibration source. The tail component was also modelled with parameters and simultaneously fitted because of the uncertainty. Total MC distribution was then calculated as the summation of these two components. Considering all the systematic and statistical uncertainties, the scintillation efficiency at 1.65 keV was estimated to be 39 +4 −4 % of that of 122 keV. As the result of this calibration, the energy scale at 1.65 keV became 20% lower than the previous scale used in [17,18]. The energy threshold for sub-GeV DM analysis via bremsstrahlung was set to 1.0 keV ee , since the uncertainty below that energy considerably increases. The scintillation efficiency at 1 keV ee was estimated to be 31 +7 −4 % of that of 122 keV.
In addition to the scintillation efficiency, detector resolution was also calibrated using these peaks. The resolution of the detector at 1.65 keV was estimated from the calibration measurement to be 40% , and Gaussian smearing was applied to MC to reproduce the data. This extra smearing was (17±10)%. The 10% uncertainty was mainly due to the surface roughness and reflection of the source.
The nonlinear response for nuclear recoil with energy over 3 keV nr was estimated using the scintillation efficiency at zero electric field in [26]. The LUX group conducted a nuclear recoil calibration [27] using neutrons from a deuterium-deuterium beam at 180 V/cm, the resultant scintillation efficiency for nuclear recoil is used to estimate the response for nuclear recoil energy below 3 keV nr . The existence of an electric field in [27] reduces the light yield. The amount of the reduction due to electric field was considered to be level of 10% [28,29]. Although the XMASS detector is operated under zero electric field, we used the unaltered results with 10% uncertainty, a typical reduction amount. The energy threshold for multi-GeV DM analysis via nuclear recoil is set to 2.3 keV nr such that we could suppress an impact of the systematic error caused by the flasher events explained in Section 6, to be smaller than other errors. At this energy, a 50% trigger efficiency of the signal simulation (8 GeV) was obtained; this threshold corresponds to 2.3 PEs. The scintillation efficiency at this energy was changed to 8.5% from 6.5%.

Analysis and results for sub-GeV DM
Event selection was applied in two stages that we referred to as standard and likelihood cuts [18]. The standard cut eliminates events that are indicative of electric noise, afterpulses, or Cherenkov emissions inside the quartz window of PMTs rather than physical interactions in the detector. Following the standard cut, we applied the likelihood cut on the basis of PE hit patterns, which removes background events occurring in front of a PMT window or near the detector wall.
The treatment of systematic uncertainties was the same as in [18]. The dominant systematic uncertainty in this analysis was associated with the variation in the PE yield during exposure. As discussed in section 4, the variation in the LXe absorption length causes a variation in the PE yield. This variation both distorts the spectrum and changes the cut efficiency. These effects were corrected based on the calculation of the relative change in the spectrum using MC simulations. To correct for each time/energy bin of measured data, MC simulations with corresponding absorption lengths derived from 57 Co calibration in each period were generated. Using these simulation results, the correction factors for the corresponding time/energy bins were calculated. These correction factors for each bins are referred to as the relative efficiency.
As it is explained in [30], the main source of background in these energy regions is 238 U and 210 Pb contained in the sealing material between the quartz window and metal body of each PMT. Since the relative efficiency depends on the spectrum shape of the expected background, the uncertainties were evaluated by comparing reasonable background models. This uncertainty of the background contributed the most to the systematic error in the relative efficiency, 1.2% and 2.5% at 1 and 5 keV ee , respectively. Note that these errors of the count rate have a correlation between each energy and time bin. The next-leading contribution came from the gain instability in the waveform digitizers between April 2014 and September 2014. During that period, a different calibration method was used for the digitizers. This variation contributed an extra uncertainty of 0.3% to the energy scale. Other contributions from the uncertainty in the PMT gain calibration using a LED, trigger-threshold stability and timing calibration were negligible.
The dataset was divided into 86 time bins (t bins ) with roughly 15 live days in each bin. The data in each time bin was further divided into energy-bins (E bins ) with bin width of 0.5 keV ee . For the DM search through the bremsstrahlung effect, the data was fitted in the energy range from 1.0 to 20 keV ee .
Minimum-χ 2 fitting was performed in the annual modulation analysis. In this analysis, the 'pull method' [31], one of the two different methods in previous analyses [17], was used to fit all energy and time bins simultaneously and to treat the correlated errors. The χ 2 function is defined as follows: where R data i, j , R ex i, j , are the data and expected number of events for the i-th energy and j-th time bins after considering the efficiency of all event selections, respectively. σ(stat) i, j and σ(sys) i, j are the statistical and systematic uncertainty of the expected number of events, respectively. The 'pull terms', α and β k represent the size of the systematic uncertainties that have correlations in energy bins or time bins. α is overall size of the relative efficiency errors common for all energy bins. Therefore, the error size of each bin changes simultaneously during the fit procedure. α = 1 (−1) corresponds to the 1 σ (−1 σ) correlated systematic error on the expected event rate. β k is the k-th systematic uncertainty of the signal simulation caused by the properties of LXe.
The uncertainties for scintillation time constants and the scintillation efficiency for the electron-recoil signal were considered. These uncertainties correlatively alter the signal spectrum between energy bins. For time constants, two components referred to as fast and slow component were used on the basis of the γ-ray calibration of the XMASS-I detector [32]. These were 2.2 and 27.8 +1.5 −1.0 ns, respectively, with the fast component fraction of 0.145 +0.022 −0.020 . For the scintillation efficiency, the uncertainty described in section 4 was used. We assumed that the signal efficiency below 1.0 keV ee is zero because of the uncertainty in the scintillation efficiency. The effect of the uncertainty of the energy resolution is much smaller than that of scintillation efficiency and is negligible. The expected number of events R ex i, j (α, β) is then expressed as follows: where t j and ∆t j are the center and width of the j-th time bin, respectively; σ χn is the DM-nucleon cross section; ǫ b i, j (α) and ǫ s i, j (α) are the relative efficiencies for the background and signal, respectively. To account for the changing background rates from long-lived isotopes, we added a simple linear function with slope B b i and constant C b i in the i-th bin. The source of the decay was considered as 210 Pb, which has a half-life of 22.3 years. A s i (β) represents the amplitude, and C s i (β) represents the unmodulated component of the signal in the i-th energy bin. In this analysis, the signal efficiencies for each DM mass were estimated using the MC simulations of uniformly injected photons from the bremsstrahlung effect in the LXe volume. The unmodulated component and amplitude of the signal spectrum were calculated for a particular cross section and mass of DM. The sub-GeV DM analysis was conducted for DM masses between 0.32 and 1.00 GeV. Figure 4 shows the observed event rate with the best fit and expected time valuation for 0.5 GeV at 1.0-1.5 and 1.5-2.0 keV ee . The search for DM mass more than 1 GeV via this bremsstrahlung effect has not been performed because the assumptions for the signal calculation in [14], such as that for form factor were not proper. The deviation was ∼0.3% and ∼3% at a maximum momentum transfer of 1 and 3 GeV DM, respectively. The best fitted cross section from the data was -1.4 +1.3 −1.6 ×10 −33 cm 2 at 0.5 GeV. The best fitted χ 2 /NDF was 3333.8/3188, and pull parameter α was 0.6 for 0.5 GeV.
Considering that we found no significant signal, the 90% confidence level (CL) upper limit on the DM-nucleon cross section σ up was calculated by the Bayesian approach [19]: where P is the probability function defined as follows: The result of the DM search via the bremsstrahlung effect is shown in the sub-GeV region of Fig. 5. The expected sensitivity for the null-amplitude case is calculated by using the statistical samples. They were generated based on the event rate obtained from a fitted result of data with only background components decreasing linearly in time, as described in [17,18]. When generating these statistical samples, data for each period and each energy bin was fitted without the signal amplitude in the first step. Thereafter, the expected number of events in each period was calculated while considering systematic errors such as relative efficiency. Finally, the Poisson fluctuation of the number of events was calculated for each energy bin, on the basis of the livetime of each period. One thousand sets of statistical samples were generated, and the 90% CL upper limit sensitivity was calculated for each sample. The 90% CL sensitivity for DM at 0.5 GeV was 2.4 +1.2 −0.8 ×10 −33 cm 2 (the range containing 68% of statistical samples) and our upper limit was 1.63 ×10 −33 cm 2 (p-value: 0.27).

Analysis and results for multi-GeV DM
An additional search for multi-GeV DM signals from elastic nuclear recoil was conducted. The analysis was mostly identical to that of the sub-GeV DM search, but data on energy less then 1.0 keV ee were analysed using nuclear recoils as low as 2.3 PE (∼ 2.3 keV nr , ∼ 0.5 keV ee ). This type of data, the low threshold data, has been recorded since December 8, 2015 with three PMT hit trigger. The total exposure of the data was 0.63 ton-years. The signal efficiency after all the data selection was improved from 5% and 10% to 10% and 15% at the lowest energy bin (2.3 -4.8 kev nr ) for 4 GeV and 8 GeV DM, respectively. This improvement of the trigger condition occurred due to the decrease of dark hits of each PMT. Average dark hits for each PMT were approximately 15 Hz at earlier periods and decreased to approximately 5 Hz during the operation. After the several data-taking tests, we were able to record stable data with the three-PMT hit triggers.
The primary uncertainty in the low-threshold data came from a weak light emission of the PMTs with a one PE. From the measurement for several PMTs in room temperature, the probability of the emission per a one PE was ∼0.3 -1.0%. Given that the light emission occurs even after dark hits, changes in the dark hits for each PMT directly change the event rate around the threshold. Thus, an additional condition for the run selection was applied to suppress this uncertainty; periods where the dark-hit rates for individual PMTs as well as the total dark-hit  [11] and CRESST-II [33], which are searching for the elastic nuclear recoil signals, are shown in each colour. The black line shows the result of the nuclear recoil search at 4-20 GeV. For comparison, results from CDMS-Si [34], CDMSLite [10], SuperCDMS [35], LUX [3], XENON1T [2], PandaX-II [4], DAMA/LIBRA [36,37], and XMASS-I [18], DarkSide-50 [38], and the liquid scintillator experiment by Collar [39] are shown for each colour. The green and yellow bands for each result show the ± 1 σ and ± 2 σ expected sensitivity of 90% CL upper limits for the null-amplitude case, respectively. rate among all the PMTs changed more than 500 Hz from the nominal values were removed from the analysis. Furthermore, the event with this light emission has characteristic timing and angular distributions of hit PMTs; the time difference between the PMT emitting the light and other PMTs receiving the light after emission distributed more than 35 ns and the latter PMTs were located within 50 degrees from the former PMT. Therefore, if any pair of hits in the events agrees with these conditions, the event was eliminated from the analysis. This event selection, referred to as a flasher cut, was applied only for three PMT hit events, and the uncertainty due to the weak flash effect after this cut is 0.4% at maximum.
The χ 2 and expected event rate functions for the time variation fitting are the same as those in the sub-GeV DM analysis except for the energy range. Most of the uncertainty for elastic nuclear recoil signal is discussed in [18], only the uncertainty of the xenon scintillation efficiency for nuclear recoil is different. As discussed above in section 4, the measurements for energy below 3 keV nr in [27] are considered.
From the multi-GeV DM analysis, we obtained the best-fit cross section between 4 and 20 GeV DM mass. The best-fit cross section is -3.8 +2.0 −4.5 × 10 −42 cm 2 at 8 GeV, and no significant signal was found in this analysis including other mass. Because of this, a 90% CL upper limit on the DM-nucleon cross section was determined. The 90% CL sensitivity at 8 GeV was 5.4 +2.7 −1.7 × 10 −42 cm 2 , and the upper limit was 2.9 × 10 −42 cm 2 (p-value: 0.11). The result of the DM search via the nuclear recoil signal is plotted in the multi-GeV region of Fig. 5. The upper limits and allowed regions determined by other experiments are also shown.
Compared with the result from the previous analysis of XMASS data [18], the result of the present analysis is approximately 6.7 times better at 8 GeV. Because both the lowthreshold data and the new scintillation efficiency below 3 keV nr in [27] improve the sensitivity. The search for DM mass below 3 GeV was not performed via nuclear recoil. This is because the maximum recoil energy is below 1 keV nr , which is the lowest calibrated energy in [27].

Conclusion
We carried out the annual modulation analysis for XMASS-I data to search for the sub-GeV and multi-GeV DM via the bremsstrahlung effect and elastic nuclear recoil, respectively. The former search limits the parameter space of DM with a mass of 0.5 GeV to below 1.6 × 10 −33 cm 2 at 90% CL. This is the first experimental result for a sub-GeV DM search focused on annual modulation and bremsstrahlung photons emitted by inelastic nuclear recoils. The additional search for the multi-GeV DM with the lower threshold data obtained a limit for the parameter space of DM with a mass of 8 GeV to below 2.9 × 10 −42 cm 2 at 90% CL.