Search for MeV Electron Recoils from Dark Matter in EXO-200

We present a search for electron-recoil signatures from the charged-current absorption of fermionic dark matter using the EXO-200 detector. We report an average electron recoil background rate of $6.8 \times 10^{-4}\, \mathrm{cts}\,\mathrm{kg}^{-1}\mathrm{yr}^{-1}\mathrm{keV}^{-1}$ above $4\,\mathrm{MeV}$ and find no statistically significant excess over our background projection. Using a total ${}^{136}\mathrm{Xe}$ exposure of $234.1\,\mathrm{kg}\,\mathrm{yr}$ we exclude new parameter space for the charged-current absorption cross-section for dark matter masses between $m_\chi = 2.6\,\mathrm{MeV} - 11.6\,\mathrm{MeV}$ with a minimum of $6\times 10^{-51}\,\mathrm{cm}^2$ at $8.3\,\mathrm{MeV}$ at the $90\%$ confidence level.


I. INTRODUCTION
The significant astrophysical and cosmological evidence for the existence of dark matter combined with the lack of understanding of its particle nature is among the most pressing problems in fundamental particle physics [1]. Dark matter is expected to consist of one or more new, beyond the Standard Model (SM) particles, whose observed interaction with the SM to date is limited to gravitational attraction. Liquid xenon (LXe) time projection chambers (TPCs) are ubiquitous in the field of rare event searches and currently provide the most stringent limits on the scattering cross-section of dark matter in the form of weakly interacting massive particle (WIMPs) [2][3][4]. With the current generation of experiments coming online, they will probe a large fraction of the remaining parameter space for dark matter masses of O(100 GeV) [5][6][7].
However, due to the lack of any confirmed detection to date of well-motivated candidates such as WIMPs or axions [8], alternative dark matter models have received recent theoretical and experimental interest [9]. A new generation of experiments is aiming to explore lower mass WIMP dark matter between 1 MeV and 10 GeV [10]. The need to detect smaller energy deposits becomes an experimental challenge for interactions of dark matter through elastic scattering with nuclei or electrons. In this case, the maximum energy that can be transferred to the detector is the kinetic energy of the dark matter particle. Other interaction mechanisms of dark matter with the SM have been studied, which could produce higher energy deposits, e.g., in the MeV energy range [11][12][13]. Detectors searching for neutrinoless double beta decay (0νββ) are well optimized for this energy range and have demonstrated some of the lowest background rates for any detector technology, along with well-understood background models [14][15][16]. This makes such detectors well suited to perform searches for dark matter interactions depositing energy in the MeV range, in addition to their primary focus of searching for 0νββ.
As a specific example of a dark matter model producing such events, the EXO-200 data is analyzed to search for charged-current absorption of fermionic dark matter [11] by 136 Xe nuclei, for which the energy deposited in the detector depends on the mass of the dark matter particle, which can be much greater than its kinetic energy. This direct detection process results in a unique detector signature with energies of O(MeV), which may have evaded detection in previous direct detection experiments. In this work, we present a search for excess events in the 1 − 8 MeV energy range in the complete EXO-200 data set.

II. THE EXO-200 EXPERIMENT
Between 2011 and 2018, the EXO-200 collaboration searched for the 0νββ of 136 Xe with a total exposure of 234.1 kg yr, setting a lower limit on the 0νββ half-life at T 1/2 > 3.5 × 10 25 yr at 90 % CL [17]. After finishing the first run of the experiment in February 2014 (Phase I) the detector underwent an upgrade before the second run started in May 2016 (Phase II) with lower noise frontend electronics and an increase in the electric field from 380 V cm −1 to 567 V cm −1 . These upgrades led to an improved energy resolution and lower energy threshold in Phase II. This work also uses the same complete EXO-200 dataset, combining Phase I and II.
The detector consisted of a cylindrical time projection chamber (TPC) filled with liquid xenon (LXe) enriched to 80.6 % in 136 Xe, with 19.1 % 134 Xe and the remaining fraction comprised of various other isotopes. The TPC was split into two identical drift regions sharing a common cathode, which along with the surrounding field rings provided an electric field for drifting electrons. Each TPC half had a radius of ∼18 cm and a drift length of ∼20 cm. The TPC was enclosed by a radio-pure thinwalled copper vessel, submerged in a heat transfer cryofluid (HFE-7000 [18]), which provided ∼50 cm of passive shielding, and was maintained inside a vacuum-insulated copper cryostat. Another ∼25 cm of lead around the cryostat provided additional shielding against external radiation. The detector was located inside a clean room at the Waste Isolation Pilot Plant (WIPP) in Carlsbad, New Mexico, under an overburden of 1624 +22 −21 meters water equivalent [16]. An active muon veto system consisted of plastic scintillator panels surrounding the clean room on four sides, which allowed prompt identification of >96 % (>94 %) of the cosmic ray muons passing through the setup in Phase I (II) [19]. Rejecting events coincident with these muons allows suppression of cosmogenic backgrounds [16]. A more detailed description of the detector design and performance can be found in [20,21].
In either drift region, the ionization and scintillation quanta produced by an interaction were collected by a crossed-wire plane at each anode, and by an array of large area avalanche photo-diodes (LAAPDs) behind the wire planes, respectively. The total reconstructed energy of an event was determined by combining the light and charge signals defining a total combined energy variable, significantly improving the energy resolution relative to either signal channel alone [22]. In addition, the combination of light and charge signals provided 3D position reconstruction with an x-y-position resolution of 2.6 mm and a z-position resolution of 0.42 mm [23]. The z-direction is defined along the electron drift direction, while the x/y directions lie in the plane parallel to the anode. More details about the analysis and event reconstruction can be found in [17,24].

III. EVENT SIGNATURE AND MONTE-CARLO
The charged-current absorption of a fermionic dark matter particle χ by a xenon nucleus will induce a β −decay if the dark matter mass can bridge the gap between the masses of the initial and final nuclei and the mass of the outgoing electron (kinematic threshold). In addition, the daughter nucleus will essentially always be produced in an excited state for the isotopes of interest here due to angular momentum selection rules. The event rate for this process is given by [11,25] summing over all possible excited states k, where m χ is the dark matter mass, ρ χ = 0.3 GeV cm −3 is the local dark matter density, N T,j is the number of targets of a given isotope j, n j is the number of neutrons per target, F k is the Fermi function, p e,k is the momentum of the outgoing β, E e,k the energy of the emitted β, M 2 Aj ,Zj the mass of nucleus A Z X, and |M k | 2 is the spin averaged nucleon level matrix element, all of which are defined in [11]. The details of the derivation of the rate, the Fermi function, and an expression for the matrix element can be found in [11,25]. Because of the higher kinematic thresholds [26] and smaller isotopic fraction, absorption of fermionic dark matter by 134 Xe is negligible relative to absorption on 136 Xe, in EXO-200. The main channel considered is, therefore, χ+ 136 54 Xe → 136 55 Cs * +e − , in which the cesium nucleus is produced in an J P = 1 + excited state, where J is the total angular momentum and P is the parity of the state. In the absence of a detection, such data can constrain the interaction cross-section of this charged-current absorption process, or equivalently, the effective energy scale, Λ, for a beyond the SM operator that could mediate the dark matter interaction where g R is the SU (2) R coupling constant of the theory, M W R = 1 2 g R u is the mass of the mediator (similar to the W-boson as the mediator of charged-current interactions via the weak nuclear force), with u being the vacuum expectation value of the scalar field φ through which the dark matter particle χ interacts with the SM via a Yukawa type coupling.
A major challenge in modeling these signals in EXO-200 comes from the fact that little data exists on the nuclear level structure of 136 Cs, with no sufficiently fast decaying states currently measured between the lowest lying 1 + excited states at 590 keV and its 5 + ground state [27]. Assuming the dark matter absorption preferentially populates the lowest-lying 1 + state at 590 keV, there are only four known states to which the nucleus can decay: three excited states (J P = 9 − , 8 − , 4 + ), and the ground state (J P = 5 + ) [27]. In all cases, the decays are highly forbidden and order-of-magnitude Weisskopf estimates for the 590 keV state give a predicted half-life of 1 s. In this case, the relaxation would appear as a secondary, uncorrelated event in the EXO-200 data, with similar signatures for excitations to higher energy nuclear levels of 136 Cs. However, if unknown intermediate levels of appropriate angular momentum and parity exist, then this half-life could be much shorter, allowing a dark matter absorption event to have a characteristic signature including both the primary β as well as one or more coincident γs. Since neither the decay time nor energy of these γs is currently known, we perform the analysis assuming that the only event signature is the outgoing β. However, as described in Sec. VI, we also consider the case that one or more coincident γs of unknown energy are emitted, and present results that are insensitive to the presence of such γs. We note that the level structure of 136 Cs is under active investigation, with data suggesting the existence of new intermediate levels that may be relevant for such dark matter interactions [28]. In addition, recent shell model calculations predict a decay scheme in which the lowest-lying 1 + state relaxes through levels that include an isomeric state with a lifetime of O(100 µs) [29], which could enable a time-coincidence analysis to separate charged-current absorption events from backgrounds. However, since the existence of these intermediate states has not been conclusively confirmed, we proceed as described above and leave their inclusion in dark matter analyses to future work. The detector response to ionizing radiation is modeled by a detailed Monte Carlo (MC) simulation based on Geant4 [30]. The details of the MC simulation and reconstruction can be found in [17,23]. The absorption event signature of a fermionic dark matter particle is modeled via a single β with a kinetic energy equal to the dark matter mass m χ minus the kinematic threshold energy [11]. The transferred kinetic energy of the dark matter particle results in a O(keV) nuclear recoil of the 136 Xe atom which falls below the energy threshold and is therefore ignored. For a given dark matter mass the injected energy into the detector will be mono-energetic and produce a peak-like signature in the reconstructed energy spectrum. The maximum reconstructed ionization energy spectrum of βs at various energies in Phase II is shown in Figure 1. The more energetic the emitted β, the more bremsstrahlung it emits whose reconstructed energy might fall below the detector threshold, leading to a non-Gaussian low energy tail.

IV. ANALYSIS
The energy region considered here can be divided into two regions based on the relevant backgrounds. In the low energy portion below 3 MeV, the dominant background sources are 2νββ events and γs from the 238 U and 232 Th decay chains. Above 3 MeV the main backgrounds arise from de-excitation γs from cosmogenically produced isotopes. To separate the β arising from dark matter interactions from backgrounds, one can exploit the fact that single MeV-scale βs typically deposit most of their kinetic energy at a single location (i.e., in the most energetic reconstructed charge cluster), in contrast to γ-backgrounds that predominantly Compton-scatter and therefore, on average, have less of their reconstructed energy contained in the most energetic cluster. An additional energy variable, defined as the energy of the most energetic charge deposit in an event, captures this topology difference. Therefore, in combination with the total reconstructed event energy, the energy of the maximum reconstructed charge cluster is used as a second dimension in the likelihood fits discussed in Sec. VI.
The fiducial volume and data quality cuts and the corresponding fiducial mass are the same as in [17], resulting in a total exposure of 234.1 kg yr divided into 117.4 kg yr and 116.7 kg yr in Phase I and Phase II, respectively.
Compared to previous 0νββ analyses from 20] this analysis The energy response of the detector is typically calibrated with external γ-sources placed on the outside of the detector near the cathode, using the full absorption peaks resulting from the decays of 60 Co and daughters of 226 Ra and 228 Th. At the end of the EXO-200 livetime, a composite americium beryllium (AmBe) neutron source was used to produce radiogenically activated 137 Xe, whose β-decay spectral shape was previously studied [31]. In addition, the AmBe source calibration provided a high statistics dataset of various other neutron activated by-products whose de-excitation gammas can be used for an energy calibration up to ∼8 MeV in Phase II only. In addition to the primary γ emitted at 4430 keV ( 9 Be(α, nγ) 12 C), this dataset includes full absorption peaks at 1346 keV ( 64 Cu(e, γ) 64 Ni), 2223 keV ( 1 H(n, γ) 2 H), 4025 keV ( 136 Xe(n, γ) 137 Xe), and 7916 keV ( 63 Cu(n, γ) 64 Cu). The reconstructed energy spectrum resulting from the AmBe source and the corresponding energy calibration results are shown in Figure 2. The lower bound of the energy range considered for this analysis follows previous EXO-200 0νββ analyses [17], whereas the upper bound is limited by the highest available energy calibration peak at 7916 keV.
While the analysis is limited to below 8 MeV to ensure the energy calibration is well-understood, the EXO-200 detector is sensitive to higher energy deposits. The low background data contains 18 (41) events above 8 MeV in Phase I (II) passing all event selection cuts, but which are not accounted for by any of the components in the EXO-200 background model [15,16]. No events above 8 MeV are present in the veto-tagged data containing only events that were collected between 10 µs and 5 ms after any of the 29 muon veto panel triggers. This suggests that these events in the low-background data could arise, e.g., from long-lived cosmogenic activation products, or unaccounted for radiogenic processes not present in the veto-coincident data. The possible effect of such a small, unknown background on the results of this analysis will be discussed in Sec. V.
To maximize the signal efficiency for β events in this energy range, this analysis does not require all ionization signals in an event to be fully reconstructed. Relative to the selection cuts used in 0νββ searches, this recovers events with multiple bremsstrahlung photons, for which there is a substantial probability that at least one lowenergy cluster is detected for which the x and y position can't be fully reconstructed. With all selection criteria included, the Monte Carlo simulation is used to estimate a signal reconstruction efficiency of > 99.7 % for both data-taking phases. Possible systematic errors in this efficiency will be discussed in Sec. V.
MC events from all background and signal components, passing all quality cuts, are used to construct individual two-dimensional probability density functions (PDFs) that depend on the reconstructed total energy and the maximum charge cluster energy of an event. The components entering the background model are the same as in [17], whereas the 0νββ signal is replaced by a monoenergetic β PDF (one for every 50 keV step within the analysis energy window), which is shown in Figure 1 Toy datasets are generated from a background-only fit to the low-background data. These toy datasets are fit with the full background plus signal model, and the upper limit on the number of signal counts for each signal PDF at the 90 % CL is determined from a profile likelihood. The sensitivity is obtained by calculating the median of an ensemble of 1000 upper limits.

V. SYSTEMATICS
The assessment of systematic errors follows the same general approach as in [17] and includes an 1. uncertainty in the signal-specific detection efficiency caused by discrepancies in the shape of data and MC PDFs 2. uncertainty in the activity of radon in the LXe as determined in stand-alone studies via measurement of time-correlated decays 3. uncertainty in the overall event detection efficiency due to possible errors in the estimated event reconstruction and selection efficiencies The first systematic error arises from shape disagreements between data and MC PDFs and is propagated to the fit as a Gaussian constraint on the normalization of the signal PDF, as was done in [17]. We performed a fit to the veto-tagged datasets in Phase I and II and treated the bin-by-bin ratio between the data and the best-fit model in each fit dimension as a possible error, reweighting the shape when building the PDFs for the likelihood fit to the low-background data to correct for the spectral distortions. Toy datasets were generated from these modified background model PDFs with the signal PDF present and were fit against the original unweighted PDFs. The bias in the number of fitted signal counts as a function of the injected number of signal counts is used to construct a conditional Gaussian constraint for the signal-only PDF normalization that has the following functional form It consists of a component proportional to the number of signal counts N with a proportionality factor a and a signal-independent component b. The signal-independent component quantifies shape errors in the background model that lead to a fixed bias in the number of signal counts, while the signal-dependent term quantifies a relative error that would be expected to scale with the signal size. The energy dependence of the parameters a and b were found to be approximately described by linear and exponential functions, respectively, and is shown in Figure 3. For Phase II an average between the bias from the spectral shape error in the muon veto-coincident data and the AmBe calibration data was used. The second uncertainty, quantifying possible error on the Rn activity in the LXe, was estimated in [16] to be 10 % and is only applied to the Rn-related backgrounds. For the last item, this analysis assumes the same systematic error for the overall signal efficiency of 3.0 % (2.9 %) for Phase I (II) as was used in [17]. This systematic error may provide a slightly conservative estimate for the reconstruction efficiencies since β-like events above 3 MeV are easier to identify with smaller reconstruction errors due to their higher signal-to-noise. As mentioned, the observed events above 8 MeV suggest that a small background component not included in the existing EXO-200 background model could be present. If such backgrounds arise from e.g., Compton scattering of high energy γs, this background component may also extend into the energy region below 8 MeV that is considered in this analysis. We estimated the effect that such an unknown background might have on our results by including an additional freely floating flat background PDF into the fit model and reevaluating the sensitivity. The relative difference in the sensitivity with and without this additional background is taken as a systematic error due to the possible incompleteness of the background model and was calculated to be ∼0.1 %. While a small energy dependence is observed, the largest systematic error over all energies is used as a conservative estimate.
The signature of charged-current absorption events is similar to a solar neutrino interaction with the detector, where in the case of the charged-current interaction (ν + 136 Xe → 136 Cs * + e − ) 136 Cs might also be produced in an excited state. In contrast, an elastic scattering process will only result in a single β in the final state. The solar neutrino background from both processes is estimated to produce 1 event throughout the livetime of the experiment within the entire energy region for this MeV. The systematic error from spectral shape disagreement is energy-dependent and is shown over the entire energy range in Fig. 3.
analysis [29] and can therefore be safely ignored.
Lastly, the systematic error on the energy scale is evaluated by taking the maximum spread in the calibrated energy for various possible calibration functions that interpolate between the measured photo-peaks of known energy (e.g. interpolating functions consisting of 1st or 2nd order polynomials, with or without a constant offset). The uncertainty in our energy calibration was estimated to be ∼2 % (∼0.5 %) in Phase I (II). The larger error in Phase I arises from the lack of AmBe calibration data, which provides known energy peaks constraining the calibration functions to the 8 MeV upper energy threshold. This systematic is propagated as a Gaussian constraint on the γ-scale, which is a freely floating parameter in the likelihood fit that scales the energy of all PDFs representing interactions of γ events via a common multiplicative factor. This floating γ-scale can compensate for possible differences between calibrations and low-background data, but the best fit value was found to be consistent with unity within less than 0.2 %. An additional β-scale is introduced into the fit that allows the energy scale of PDFs arising from β-like events to float independently from γ-like PDFs, which allows the fit to correct possible differences in the energy response between βs and γs in the detector. However, this parameter was similarly found to be consistent with unity to within less than 1 %, in good agreement with previous studies [32]. Table I shows a summary of the systematic errors considered for this analysis.

VI. RESULTS
Following [17], the likelihood fits are initially performed independently in each data-taking phase, minimizing a binned NLL function and then profiling over the signal PDF at a given energy. Hence, the combined NLL at a given energy is simply the sum of both NLLs, and the efficiency and livetime corrected profile likelihood curves, which are a function of the dark matter interaction crosssection, can be added to obtain the combined upper limit at 90 % CL. The determination of the 90 % percentile of FIG. 5. Limit on the absorption cross-section of fermionic dark matter by 136 Xe nuclei at the 90 % CL. We show our main 2D fit result in addition to a de-excitation γ agnostic 1D fit result (see text for more details). The grey region is excluded by direct constraint from searches at the LHC [33].
the test statistic is assuming that Wilks' approximation holds [34]. A comparison with the obtained upper limit using MC-based test statistic distribution at an example dark matter mass of 5 MeV shows that the results of this work produce a ∼10 % weaker limit, and thus the assumption of Wilks' approximation does not have a significant effect on the limit (and may be slightly conservative).
The combined median 90 % CL sensitivity is similarly determined by combining the profile likelihoods of 1000 toy datasets from each phase. The combined 90 % upper limit and sensitivity are shown in Fig. 4, together with the 68 % and 95 % percentiles around the median sensitivity. We calculate the local p-value as the fraction of toy datasets under the no-signal hypothesis whose ∆NLL is at least as large as the one obtained from the data fit. The largest local significance was found to be ∼ 2σ at 2.7 MeV. We, therefore, find no statistically significant evidence for dark matter in our data. Overall, the observed limit lies within the ±2σ distribution around the median sensitivity over most of the mass range. In the region between 5 MeV and 7 MeV, the limit is slightly stronger than the expected ±2σ region. This could arise from either a statistical under-fluctuation in the backgrounds at these energies or, possibly, a small overestimate of the background in this region relative to that assumed in the toy Monte Carlo studies. Figure 5 shows the 90 % CL exclusion limit on the interaction cross-section as a function dark matter mass m χ . In addition to the primary result using a two-dimensional likelihood fit consisting of the total combined energy and the maximum charge energy in an event, we also performed the fits using only the latter variable. These 1D fits have reduced sensitivity to a signal but provide an analysis that does not strongly depend on the presence of additional de-excitation γs, since the primary β cluster will contain the largest energy for nearly all events over the mass range considered, resulting in only marginal differences in the signal PDF shape. We also include limits from constraints on these models from collider experiments at the Large Hadron Collider (LHC) at CERN [33], demonstrating that relative to collider searches low-background detectors can perform competitive searches for certain dark matter models. Additional constraints arise from indirect detection searches due to the decay of such dark matter particles that become more stringent at higher dark matter masses [35].
In Fig. 6 we provide the total combined energy spectrum of the low-background data from 1 MeV to 8 MeV, along with the best fit to the data with a background-only model, and the bin-wise residual counts for Phase I and Phase II. The excellent spectral shape agreement and linear detector response to both β-like and γ-like events are a result of more than ten years of work to build, run, and understand the EXO-200 detector and its response. These spectra, published here for the first time, may also constrain other physics beyond the SM that is outside the scope of the present work.
This analysis represents the first search for the absorption of MeV fermionic dark matter in a liquid xenon detector via a charged-current interaction, with a total 136 Xe exposure of 234.1 kg yr. As no statistically significant evidence was observed, we exclude new parameter space for these models at the 90 % CL for dark matter masses between 2.6 MeV and 11.6 MeV. These results are complementary to searches for neutral-current absorp- FIG. 6. Total combined energy spectrum of the low-background data and the best fit (blue) to the data with a background-only model in Phase I (a) and Phase II (c). The resulting residual number of counts between data and the best from Phase I and Phase II is shown in panels (b) and (d), respectively. No significant disagreement is seen between the data and the background-only model. The average electron-recoil background rate above 4 MeV is 4.0 × 10 −4 cts kg −1 yr −1 keV −1 in Phase I and 6.8 × 10 −4 cts kg −1 yr −1 keV −1 in Phase II with a total exposure of 117.4 kg yr and 116.7 kg yr, respectively. tion of fermionic dark matter reported in [36,37] and charged-current absorption of light fermionic dark matter in [38]. Furthermore, a detailed understanding of the detector response is demonstrated. The low backgrounds achievable with the liquid xenon technology can enable future rare event searches such as 0νββ, direct detection of WIMP dark matter, or other novel dark matter interaction mechanisms [39][40][41].