Search for dark matter in the form of hidden photons and axion-like particles in the XMASS detector

Hidden photons and axion-like particles are candidates for cold dark matter if they were produced non-thermally in the early universe. We conducted a search for both of these bosons using 800 live-days of data from the XMASS detector with 327 kg of liquid xenon in the fiducial volume. No significant signal was observed, and thus we set constraints on the $\alpha' / \alpha$ parameter related to kinetic mixing of hidden photons and the coupling constant $g_{Ae}$ of axion-like particles in the mass range from 40 to 120 keV/$c^2$, resulting in $\alpha' / \alpha<6 \times 10^{-26}$ and $g_{Ae}<4 \times 10^{-13}$. These limits are the most stringent derived from both direct and indirect searches to date.


Introduction
The existence of Dark Matter (DM) is inferred from many cosmological and astrophysical observations [1]. All indications for its existence so far were based on its gravitational interaction only, E-mail address: xmass .publications11 @km .icrr.u -tokyo .ac .jp. 1 Now at Department of Physics, College of Science and Technology, Nihon University, Kanda, Chiyoda-ku, Tokyo 101-8308, Japan. 2 Now at Center for Axion and Precision Physics Research, Institute for Basic Science, Daejeon 34051, South Korea. and the nature of DM beyond that is still shrouded in mystery. One of the theoretical realizations of DM is Weakly Interacting Massive Particles (WIMPs), which are expected to have masses of roughly 10 GeV/c 2 to a few TeV/c 2 . Although many direct and indirect experiments are searching for WIMPs, no concrete signal has yet been found. Hidden Photons (HPs) and Axion-like Particles (ALPs), which are respectively vector and pseudo-scalar realizations of bosonic super-WIMPs [2], are alternative candidates for DM with expected masses < 1 MeV/c 2  model photon [3]. ALPs arise as pseudo-Nambu-Goldstone bosons that appear generically in all string compactifications [4]. Although a scenario in which they were thermally produced in the early universe was ruled out [5] or disfavored [2], HPs and ALPs can be cold DM candidates that give rise to the observed DM abundance if they were produced non-thermally via the mis-alignment mechanism [4]. HPs and ALPs are experimentally interesting because both of them are absorbed by materials with an interaction analogous to a photoelectric effect [2], transferring their total energy to electrons. Assuming they are cold DM, i.e. they are non-relativistic, the energy they deposit is equivalent to their rest mass. The calculation in Ref. [2] is limited to bosons with keV-scale mass up to roughly 100 keV/c 2 . In this mass region, indirect searches give stringent limits of < 10 −26 on the α /α parameter of HP, which represents the strength of the kinetic mixing, in the mass range from 1 eV/c 2 to 50 keV/c 2 and above 140 keV/c 2 [6]. Around the mass of 90 keV/c 2 , however, the indirect limit is relatively weak as α /α < O (10 −24 ). XMASS already carried out searches around this mass region for HPs and ALPs using data taken in 2010-2012, and had given limits in the mass range of 40-120 keV/c 2 [5]. In this paper, we report an improved result using 800 live-days of data from November 2013 to July 2016 and a fiducial volume containing 327 kg of liquid xenon. The sensitivity of the search is improved by an overall reduction of the background (BG), advances in our understanding of the BG, and a significant increase in available exposure.

XMASS detector
The XMASS detector [7] is located at the Kamioka Observatory, which is an underground laboratory in Japan, at a depth of 2,700-m water equivalent. The detector consists of an Inner Detector (ID) and an Outer Detector (OD).
The ID is a liquid xenon scintillator contained in a vacuum insulated copper vessel. The 832 kg of liquid xenon in the sensitive region are surrounded by 642 inward-facing photomultiplier tubes (PMTs). The PMTs are arranged in a copper holder the shape of which is an approximate sphere with a radius of ∼40 cm. The photocathodes of the PMTs cover 62% of that inner surface. Aluminum used as a seal material in the PMT contains some radioactive isotopes (RIs) which were the main BG sources in our previous analysis [5,7]. In 2013 we installed copper covers over these aluminum seals to reduce this BG [8]. Only new data taken after this installation is used in this analysis.
The OD encloses the ID, which is at its center. The OD is a water Cherenkov counter of cylindrical shape, 10 m in diameter and 10.5 m in height, and is read out by 72 20-inch PMTs. The OD is used as an active veto for cosmic-ray muons and serves as a passive radiation shield against environmental neutrons and γ -rays.
The signals from the ID-PMTs are recorded by CAEN V1751 waveform digitizers with a sampling rate of 1 GHz. Data acquisition is triggered if at least 4 PMTs detect signals within 200 ns of each other over 0.2 photoelectron (PE) threshold. The gains of the ID PMTs were monitored with a blue LED embedded in the inner surface of the detector, dimmed to provide single PE signal levels. The detector response to scintillation light was traced by inserting a 57 Co source [9] into the ID every two weeks. The light yield of the detector for 122 keV γ -ray and both the absorption length and the scattering length of the liquid xenon were extracted from this regular calibration. The time evolution of these values was used as input parameters for the Geant4-based XMASS Monte Carlo simulation (MC) [7].

Event selection and vertex reconstruction
Events without associated OD activity (8 or more OD PMTs detect signals over a 0.4 PE threshold) were used for the analysis. We applied our 'standard cuts', which are summarized in [8], to reject events originating from PMT after-pulses and electronic noise. We also required the timing difference to the subsequent event to be > 1 ms, to reject 214 Bi β-ray events from the 214 Bi-214 Po decay sequence in the 222 Rn decay chain.
The event vertex was then reconstructed based on the maximum likelihood evaluation of the observed PE distribution in the ID PMTs [7]. The expected number of PE on each PMT was obtained from MC simulations for various vertex positions r and xenon absorption lengths l abs , and used to calculate a likelihood for each event position. As the regular 57 Co calibration provided the time evolution of l abs , the MC expectations updated accordingly. The position r which maximizes the likelihood was accepted as the event vertex. We required the radius of the reconstructed vertex R rec (≡ | r| with the center of the ID as the origin of our coordinate system) to be < 30 cm. This fiducial volume cut strongly reduces BG events originating from γ -rays or β-rays from RI in/on the detector's inner surfaces and bulk materials; its impact is shown in Fig. 1.
The PE acceptance depends on the vertex position and also on optical parameters such as l abs and light yield in the liquid xenon. We thus use a corrected number of PEs (NPE cor ) as a measure of the energy deposit to avoid such dependencies. NPE cor is defined as: where μ( r, l abs ) is the detection efficiency for the scintillation light generated at r when the absorption length is l abs , and μ( r = 0, l abs ) is the efficiency at the detector center. These efficiencies are calculated from the MC expectations used in the vertex reconstruction.
The ν is the time dependent correction factor for the PE yield at the detector center obtained from the regular 57 Co calibrations.
Changes of ν can be explained by changes in the absorption length of liquid xenon as discussed in [10]. The change of ν over the data taking period are controlled within ±7% from its averaged value.

Signal MC
The absorption of HPs and ALPs is analogous to the photoelectric effect. The cross-section of the absorption σ abs for HP can be written in terms of the cross-section of the photoelectric effect σ pe when replacing the photon energy ω by the HP mass m H P , as: where υ is the velocity of the HP, α is the fine structure constant, and α is its analogue for the HP. Assuming a dark matter density of 0.3 GeV/cm 3 , the event rate is expressed as [2]: where A = 131.3 is an atomic mass of xenon of natural composition.
As for the absorption of ALPs, which is also referred to as the axio-electric effect, its cross-section and expected event rate are calculated as [2]: (4) and where m AL P and g Ae are the mass and the axio-electric coupling constant of the ALPs, respectively. The f a is a dimensional coupling constant for ALPs defined as f a ≡ 2m e /g Ae , where m e is the mass of the electron. In this analysis, the expected rate does not depend on DM velocity distribution as indicated in Eqs. (3) and (5). In the following analysis we concentrate on HPs. The result for ALPs can be obtained by replacing the event rate from Eq. (3) by that from Eq. (5). We simulated the photoelectric effect of photons with energies equal to m H P and used this simulation as signal MC. This simulation is identical to the absorption of the HPs, because all the energy of the incident particle, including its rest mass, is transferred to electrons. Photons were generated uniformly inside the ID. We made the momentum directions of the photons simply isotropic, since they do not affect on result of photoelectric effect. The signal MC was made for m H P from 40 to 120 keV/c 2 in steps of 2.5 keV/c 2 . The same selection criteria as described in Sec. 3.1 were applied. The ratio of the number of events passing the selection to the ones generated in the 30-cm-radius fiducial volume is 0.95-0.98 for any of the HP masses. Examples of the expected signal spectrum are shown in Fig. 2.

Background MC
MC samples were prepared for every BG source separately. The BG sources were divided into three groups: RIs dissolved in the liquid xenon, xenon isotopes activated by neutrons, and RIs in/on the detector components.
For the RIs in the liquid xenon, the amount of 214 Pb in the detector, which is a daughter of 222 Rn, was estimated from 214 Bi-214 Po coincidences in the 222 Rn decay chain to be 8.53±0.16 mBq on average throughout the data taking period [11]. The average amount of 85 Kr in the detector was estimated to be 0.25±0.04 mBq using the coincidence of its β emission with Q β = 173 keV followed by a 514-keV γ -ray emission [8]. The amount of the 2νββ decay of 136 Xe ( Q ββ = 2.46 MeV) is estimated from its natural abundance (8.9%) assuming T 1/2 = 2.21 × 10 21 years [12]. A contamination of argon was found from a component analysis of the detector's xenon gas. The β decay of 39 Ar can be a BG source. Also a contamination of 14 C was indicated by its characteristic spectral shape in the observed energy spectrum though its chemical form is not known. The contributions of 39 Ar and 14 C were evaluated from the observed spectrum [11].
The xenon isotopes and their daughters, mainly 131m Xe and 125 I, are thought to be produced in gas phase by thermal neutrons when the gas is outside of the OD water shield, and then returned to the liquid xenon in the detector. The amounts of 131m Xe and 125 I were calculated from the measured thermal neutron flux in the Kamioka Observatory (0.8−1.4) × 10 −5 cm −2 s −1 [13,14].
The amounts of RIs in/on the detector components were obtained from screening with germanium detectors and from the shape of the observed energy spectrum between 30 and 3000 keV without the fiducial volume cut, as summarized in [8].
The uncertainties of the BG amounts were considered as systematic uncertainties in our signal peak search process, as described in Sec. 3.5.

MC treatment and its uncertainty
Several corrections, which were mainly based on calibration data, were applied to the MC as summarized in Table 1. For reference purpose we label our five corrections C1 through C5. The correction C1 follows the standard MC treatment in XMASS [11]. In addition, the corrections C2-C5 were introduced to include our knowledge of the various BG components in the energy range used in this analysis. The uncertainty in each of these corrections was used as a systematic uncertainty in the fitting process described in Sec. 3.5.
The correction C1 is for the non-linear scintillation efficiency of xenon, which is important for the signal search as a change of non-linearity moves the signal peak position and deforms the BG spectrum shapes. The non-linearity was taken into account in our MC based on the model described in [15]. As with the XMASS standard treatment [11], it was then calibrated using several γ -ray sources; 55 Fe (5.9 keV), 241 Am (59.5 keV γ -ray, 17.8 keV Table 1 Corrections to the MC and their errors. For the non-linearity correction C1, the results of our five calibration points were connected using f model as described in the text. As for the energy and position resolution corrections C2 and C3, the results at 59.3 and 122.1 keV were interpolated with a function of energy A ⊕ B/ √ E where the operator ⊕ represents quadratic sum, and the parameters A and B were determined by requiring the function to connect the results at the two calibration points.

ID kind of correction
(energy of γ -rays or events used to derive the correction) correction factor and its error (≡ C m ± δC m used in Eq. (8)  Np's X-ray, and ∼30 keV escape peak of Xe's X-ray), and 57 Co (122.1 keV γ -ray and 59.3 keV X-ray from the tungsten housing of the source). The results were expressed relative to the 122.1 keV 57 Co γ -rays. These calibration points were linearly interpolated to model the energy dependence of the non-linearity as shown in Fig. 3, and the observed energy in the MC was scaled according to the model function. The measurement error of each calibration point was considered as a systematic uncertainty. The uncertainty of the energy dependence model was also taken into account. In this analysis, in addition to the linear interpolation, we modeled the energy dependence with several interpolation methods (spline, and polynomial interpolation), and also fitted the calibration points with polynomial functions of various order, as well as using combinations of linear interpolation and fitting. Examples of such model functions ( f model ) are shown in Fig. 3. As described in Sec. 3.5, the significance of signals was evaluated including such model uncertainties.
The resolutions for energy (NPE cor ) and position (R rec ) were evaluated from the 59.3 keV and 122.1 keV peaks in the 57 Co calibration runs. The energy resolution in the data was evaluated to be 8.3 ± 0.8% at 59.3 keV and 3.9 ± 0.2% at 122.1 keV, which were worse than in the MC. The quadratic subtraction of the resolution in the MC from the one in the data was 3.8% at 59.3 keV and 1.1% at 122.1 keV. The energy spectrum of the MC was additionally smeared to compensate for this difference (correction C2). The po-sition resolution in the data was evaluated to be 6.6 ± 0.5 mm for 59.3 keV X-rays and 4.5 ± 0.1 mm for 122.1 keV γ -rays at the fiducial volume radius (= 30 cm). The resolution in the MC was slightly better, the quadratic subtraction of which from the one in the data was 2.6 mm (1.3 mm) for 59.3 keV X-rays (122.1 keV γ -rays). The efficiency for the fiducial volume cut depends on the position resolution. The efficiency in the MC was thus corrected according to this difference (correction C3), but compared to the other corrections its impact on the analysis is not significant.
Depending on the period of data taking, eight to ten PMTs among all 642 PMTs were not operational because of their high noise rates or electrical problems. Our MC simulation shows that events generated on the detector surface near these dead PMTs are often mis-reconstructed within the fiducial volume. Such events come from RIs in the PMTs or 210 Pb and its decay products in the copper cover for the aluminum seal, and contribute in particular to the low energy region [8]. By artificially masking normal PMTs in the data, we found that the reconstructed vertices ( r) of such events concentrate around the axes connecting the detector center with the dead PMTs, and that reconstruction moves them towards the center of the detector. The probabilities of mis-reconstruction in the data and the MC were estimated using this feature [8] for two energy regions separately. According to the difference between the probabilities, the event rate in the MC was increased by 7% for the 441 < NPE cor < 515 region and 19% for the 515 < NPE cor < 588 region (correction C4 and C5).

Signal peak search
We can extract a potential signal from the data by comparing the observed energy spectrum with the combined predictions of signal and BG MC including their respective uncertainties.
The energy spectrum after applying all the selections is shown in Fig. 4. The peak around NPE cor = 2400 came from residual 131m Xe after calibration with 252 Cf, which was useful as a ref- To search for the HP signal, the histogram of the observed energy spectrum was fitted with the sum of a putative signal and the BG MC spectra. The range of the histogram used for fitting was 440-2650 NPE cor (corresponding to ∼ 30-180 keV γ -ray energy), and divided equally into 150 NPE cor bins. The chi-square was defined as: where the summation is taken for all the RIs of the three BG categories described in Sec. 3.3. The R j-th BG is the expected event rate from the j-th RI and p j is its scale parameter whose initial value is unity. The χ 2 sys is a penalty term to handle systematic uncertainties. It is defined as: where the δ p j are the 1σ uncertainties of the BG estimates described in Sec. 3.3. The first summation on the right side constrains the parameters p j around ±1σ from unity during the fit. Because the amounts of 14 C and 39 Ar in xenon are obtained directly from the fit to the data, these RIs are not included in this summation while they are included in the summation in Eq. (7). The second summation is a penalty term related to the five special MC corrections C1-C5 described in Sec. 3.4. The m-th correction factor C m listed in Table 1 R H P = R H P (m H P , g Ae ; E , C 1 , ..., C 5 ; f model ) . (10) The χ 2 f it was minimized separately for every 2.5 keV/c 2 step in HP mass between 40 and 120 keV/c 2 and for each step of α /α in line shows the 90%CL constraint presented in this paper. The black line shows our previous result [5]. The blue, magenta, green, and orange lines are limits reported by the XENON100 [18], the Majorana Demonstrator [19], the LUX [20], and the PandaX-II [21]. The dotted, dashed, and dash-dotted lines in light blue color are constraints from indirect searches derived from red giant stars (RG), diffuse γ -ray flux, and horizontal branch stars (HB), respectively [6].
the α /α > 0 region, by fitting the p j , C m , E , and f model . The optimization of these parameters except for f model was done using ROOT TMinuit [16]. As for f model , each f model was tested separately, and the one giving the smallest χ 2 f it was chosen for each mass and each α /α. This method corresponds to handling the model function shape as one of the fitting parameters [17]. The sensitivity obtained with this method is conservative compared to sticking with one model. We thus obtained the χ 2 f it profile as a function of α /α for each mass. The minimum of the profile is the most probable α /α parameter for that mass.

Result
The best fit result for HP with m H P = 85 keV/c 2 is shown in Fig. 4, where the minimum χ 2 f it /NDF = 131/122 with α /α = 1.1 × 10 −26 . No significant signal was found at any HP mass. The difference between the minimum χ 2 f it and the χ 2 f it with no signal was at most 1.62. We thus set the 90% confidence level (CL) constraint on α /α from the relation: where a and a 90 denote α /α and its constraint, respectively. The constraint for each mass is shown in Fig. 5. Compared to our previous work, the constraints improved by a factor of 10-50. The constraints from other direct and indirect searches are also shown in the figure. Our result gives the most stringent limit in the mass range from 40 to 120 keV/c 2 . The indirect limits for HP around 90 keV/c 2 are α /α < O (10 −24 ), relatively weak compared to the higher and lower mass regions (α /α < O (10 −27 ) for m H P ≥ 200 or ≤ 40 keV/c 2 ). Our limit, α /α ≤ 2 × 10 −27 -6 × 10 −26 , bridges that region of relative weakness for the indirect limits. The same result converted to a constraint on g Ae for ALPs using Eqs. (3) and (5) is also shown in the figure. LUX and PandaX-II present limits for ALPs below 20 keV/c 2 . Our result covers the higher mass region. While XENON100 and the Majorana Demonstrator report limits over a wider mass region, our constraint g Ae < a few 10 −13 is the best limit for m AL P > 40 keV/c 2 .

Conclusion
The searches for HPs and ALPs, which are candidates for cold DM, were conducted using 800 live-days XMASS data with 327 kg xenon in a 30-cm-radius fiducial volume. We searched for signal peaks from HPs or ALPs interactions analogous to the photoelectric effect in the energy spectrum around tens of keV, where the event rate in XMASS is ∼ 5 × 10 −4 day −1 kg −1 keV −1 . With the absence of any significant signal, we set the most stringent upper limits for the parameter for kinetic mixing α /α of the HP and for the coupling constant g Ae of the ALPs in the mass range from 40 to 120 keV/c 2 .