First Probe of Sub-GeV Dark Matter Beyond the Cosmological Expectation with the COHERENT CsI Detector at the SNS

The COHERENT collaboration searched for scalar dark matter particles produced at the Spallation Neutron Source with masses between 1 and 220~MeV/c$^2$ using a CsI[Na] scintillation detector sensitive to nuclear recoils above 9~keV$_\text{nr}$. No evidence for dark matter is found and we thus place limits on allowed parameter space. With this low-threshold detector, we are sensitive to coherent elastic scattering between dark matter and nuclei. The cross section for this process is orders of magnitude higher than for other processes historically used for accelerator-based direct-detection searches so that our small, 14.6~kg detector significantly improves on past constraints. At peak sensitivity, we reject the flux consistent with the cosmologically observed dark-matter concentration for all coupling constants $\alpha_D<0.64$, assuming a scalar dark-matter particle. We also calculate the sensitivity of future COHERENT detectors to dark-matter signals which will ambitiously test multiple dark-matter spin scenarios.

The COHERENT collaboration searched for scalar dark matter particles produced at the Spallation Neutron Source with masses between 1 and 220 MeV/c 2 using a CsI[Na] scintillation detector sensitive to nuclear recoils above 9 keVnr. No evidence for dark matter is found and we thus place limits on allowed parameter space. With this low-threshold detector, we are sensitive to coherent elastic scattering between dark matter and nuclei. The cross section for this process is orders of magnitude higher than for other processes historically used for accelerator-based direct-detection searches so that our small, 14.6 kg detector significantly improves on past constraints. At peak sensitivity, we reject the flux consistent with the cosmologically observed dark-matter concentration for all coupling constants αD < 0.64, assuming a scalar dark-matter particle. We also calculate the sensitivity of future COHERENT detectors to dark-matter signals which will ambitiously test multiple dark-matter spin scenarios. Introduction: Standard model (SM) fermions only account for ≈ 20% of cosmologically observed mat-galaxies [2] and galaxies within clusters [3,4]. Despite continuing improvement in our understanding of the gravitational effects of DM, its particle nature has not been determined. Searches for traditional weakly interacting massive particle (WIMP) DM have not yet found a positive signature [5][6][7][8][9]. Further, experimental sensitivity is rapidly approaching a "neutrino fog" [10] of background from coherent elastic neutrino-nucleus scattering (CEvNS) [11] events from astrophysical neutrino sources which will hinder progress.
In response, the interest in sub-GeV DM particles, too light to be observed in many conventional WIMP detectors, has increased recently. Cosmological observations suggest that such DM could not interact with SM matter through the weak force [12]. However, sub-GeV hidden-sector DM particles could interact with SM fermions mediated by a "portal" particle [13][14][15][16]. These proposed hidden sector particles are viable DM candidates.
If sub-GeV DM exists, these particles would be produced at accelerators. Beam-dump experiments have already begun to survey the possible parameter space [17][18][19][20][21][22][23] with more experiments planned [24,25]. Searches for accelerator-produced DM are of particular interest as the DM particles are relativistic so that the scattering cross section is relatively spin-independent [26].
Experiments capable of seeing low-energy nuclear recoils associated with CEvNS can search for an analogous coherent elastic DM-nucleus scattering process [27,28]. As the cross section scales according to the square of the proton number Z 2 , such an experiment can achieve competitive sensitivity with relatively low mass. Further, as sub-GeV interactions with the SM are mediated by a new force particle, the hadronic and leptonic DM couplings may be radically different, thus calling for experimental efforts to constrain both. While the majority of previous constraints test the coupling between DM and leptons, CEvNS experiments are sensitivity to coherent DM-nucleus scattering, thus probing the quark coupling.
In this paper, we present the first search by CO-HERENT for accelerator-produced DM particles. This uses data collected by our CsI detector which measured CEvNS [29,30] at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory [31]. We focus on a single benchmark model of scalar DM particle, χ, mediated by a vector portal particle, V , [16] with masses m χ and m V , respectively. In this model, V kinetically mixes with the SM photon with a coupling ε. DM particles can then be produced through V → χχ decay with a coupling α D . The DM scatterign cross section at thermal freeze-out, and thus relic DM abundance in the modern universe, depends on the single parameter, Y [32], defined as We thus adopt this parameter when presenting our results for convenient comparison to other measurements. Our analysis is restricted to scalar DM in the mass range 1 < m χ < 220 MeV/c 2 . A spin ½ particle is also viable as a DM candidate, though these scenarios would require lower couplings to match the cosmologically observed concentration. With improved sensitivity in the future, we will explore constraints on Majorana and pseudo-Dirac fermion DM, but currently focus on scalar DM. The COHERENT CsI detector at the SNS: The SNS operates a 1.4 MW proton beam incident on a mercury target running at 60 Hz. For our detector operations, the SNS maintained a beam-pulse width of 378 ns FWHM and an average proton energy of 0.984 GeV. With its high beam power, the SNS could produce an enormous flux of DM particles through proton bremsstrahlung and hiddensector decays of π 0 and η 0 mesons produced in the target. Neutrinos from the accelerator, produced by the decay of π + particles which formed as protons stop in the target, induce CEvNS in our detectors, one of our principal backgrounds to dark-matter detection. This neutrino flux includes a prompt component from π + → µ + ν µ and a delayed component, τ = 2.2µs, from µ + → e + ν eνµ decay.
We operate several detectors in "Neutrino Alley" at the SNS, a basement hallway with sufficiently low backgrounds to allow for neutrino measurements. One of our detectors was a 14.6 kg CsI[Na] scintillating crystal [33,34], commissioned in 2015, which made the first observation of CEvNS [29]. This detector was decommissioned in 2019. During its run, this detector collected 13.99 GWhr, 3.20 × 10 23 protons on target, of beam data. The detector was situated 19.3 m from the beam target, 90 • off-axis from the beam direction. The light was collected by a single Hamamatsu R877-100 photomultiplier (PMT) sampled at a rate of 500 MS/s. We assembled shielding with multiple materials to moderate both environmental γ and neutron activity.
The light yield of the detector was determined to be 13.35 photo-electrons (PE) / keV-electronequivalent (keV ee ). Two sources were used for calibration: a 59 keV ee γ line from 60 Co decay and a 57.5 keV ee peak from 127 I(n,γ) which includes a small, quenched nuclear recoil. The light yield was shown to be uniform across the crystal by taking calibration data at multiple locations along the detector length.
Dark matter events in the CsI detector: We use the BdNMC [35] simulation package to predict the DM flux in Neutrino Alley along with the scattering rate and kinematics within our detectors. BdNMC is versatile, calculating DM production and detection through several channels. Coherent elastic DMnucleus scattering has been implemented specifically for CEvNS experiments.
The dominant production channels for portal particles at the SNS are π 0 → γ + V decay, η 0 → γ + V decay, and p + N → p + N + V bremsstrahlung. Production from π 0 decay, η 0 decay, and proton bremsstrahlung dominate for DM masses below 40 MeV/c 2 , between 40 and 130 MeV/c 2 , and above 130 MeV/c 2 , respectively. We do not have sensitivity for m χ > 220 MeV/c 2 , beyond which bremsstrahlung is kinematically forbidden. With a GEANT4 [36] simulation, we predict 0.107 ±10% π 0 produced per proton incident on the target [37]. For our most sensitive masses, π 0 decay dominates the sensitivity. The production of portal particles from π 0 decay is given by the branching ratio for m V < m π [38] which is proportional to the expected DM flux. Though the beam energy at the SNS, T p ≈ 0.98 GeV, is slightly lower than the production threshold for p + p → p + p + η 0 production, there are η mesons produced in the target due to the Fermi momentum of mercury [39]. A calculation of this sub-threshold production [40] suggests that about 0.002 ±30% η 0 are produced per π 0 at the SNS. BdNMC predicts the timing of scattering events which typically scatter within a few ns of the speed-of-light-delayed DM production in the target. As this is a small delay, we assume all DM we study travels at the speed of light.
To lowest order, the differential cross section in recoil energy, E r , is where α is the electromagnetic fine structure constant, m N is the nuclear mass, and p χ and E χ are the incident DM momentum and energy. The event rate is given by the flux × the cross section and thus depends on the couplings as ∝ ε 4 α D ∝ Y 2 /α D . The scattering model used for our sensitivity estimates presented in [41] had a calculation error with the definition of Q 2 = 2m N E r , described in [27], that has now been fixed. We have confirmed event rates predicted by BdNMC using this new model with an independent, cross-check calculation from COHER-ENT.  Figure 1. The expected spectrum of coherent DM scattering signal in CsI from both galactic and SNS-produced DM for a mass of 25 MeV/c 2 , near our optimal sensitivity. Though the rate for galactic DM is higher, the recoil energies are far below threshold while we select 26% of DM produced at the SNS.
For m χ = 25 MeV/c 2 , the expected average recoil energy is 9 keV, just at our analysis threshold. The spectra of interacting DM in our CsI detector are shown in Fig. 1 for both galactic and SNS-produced DM assuming a DM mass at our peak sensitivity, m χ = 25 MeV/c 2 . The prediction of the galactic recoil spectrum assumes a local DM density of 0.3 GeV/cm 3 near Earth with a Boltzmann speed distribution with v/c ≈ 0.001 [42]. Though fewer interactions are expected for SNS-produced DM, the typical recoil energy is higher than for galactic DM by a factor of 10 6 allowing for detection of 26% of the SNS-produced DM.
The detector response for DM recoil events is assumed to be the same as for CEvNS, described in [30], apart from quenching at high recoil energies. All data used to fit our quenching model were taken at E rec < 70 keV nr . This is sufficient to cover all CEvNS recoils; however a small percentage of DMinduced recoils lie beyond this point. For recoil energies above 70 keV nr , we assume a constant quenching factor, (9.8 ± 1.8)%, which is the quenching and uncertainty implied by our fit at 70 keV nr .
Data analysis: We performed a search for light DM particles in our CsI data collected during SNS operations. The analysis was blinded, defining all selection cuts, uncertainties, and fitting methods before determining the observed data spectrum. The DM scattering model, however, was updated after box-opening to correct an error discovered in the coherent cross section. The corrected version is given by Eqn. 3.
We used the same event reconstruction used to determine the CEvNS cross section [30]. We also applied the same event selection, except that the highest recoil energy analyzed was increased from 60 PE to 250 PE to capture the most energetic recoils expected from high-mass DM interactions. The recoil energy binning was also re-optimized for ideal DM separation from the backgrounds. The analysis binning was not reoptimized after box opening and updating the cross section model. Steady-state accidentals (SSBkg) and CEvNS interactions are the dominant backgrounds. A small number of beamrelated neutron (BRN) and neutrino-induced neutron (NIN) events were also accounted for. Neutron rates and uncertainties were determined from EJ-301 [43] liquid scintillator data collected before detector commissioning [30].
The CEvNS cross section was fixed to the SM prediction and allowed to float within the form-factor uncertainty, 3.4%. We also included systematic uncertainties from neutrino flux, background normalization, threshold, selection efficiency, and quenching that are calculated as described in [30,44] and propagated to the DM signal prediction when appropriate. The form-factor, efficiency, and quenching uncertainties included both normalization and shape uncertainties.
Our DM prediction, parameterized by DM mass, m χ , and coupling, Y , was added to our expected SSBkg, BRN, NIN, and CEvNS backgrounds. We tested DM masses between 1 and 220 MeV/c 2 which covers the range where COHERENT has competitive sensitivity.
The timing of observed events, t rec , is critical for this result. We measure t rec relative to a trigger signal from the accelerator indicating protons on target, and the neutrino pulse arrives in our detector roughly 150 ns later. The DM region of interest (ROI) is defined as 0.25 ≤ t rec < 0.75 µs, determined by selecting the optimal timing region for separating the DM signal from backgrounds with a s/ √ b figure of merit. Over 92% of DM recoils but only 25% of CEvNS events are expected in this interval. Most neutrino events are delayed relative to the DM events by τ µ = 2.2 µs. These delayed events can be used to constrain systematic uncertainties in situ to improve the precision of background estimates within the DM ROI.
The data was binned in two dimensions: recoil energy and recoil time. For each DM mass and coupling, we performed a binned log-likelihood fit profiling over nuisance parameters relating to systematic uncertainties. As Y depends on two couplings, α D and ε, we conservatively fix α D = 0.5 in our constraint [27]. Lower values reduce both the expected event rate and Y . Since the DM scattering cross section at freeze-out is ∝ Y , when assuming a fixed DM relic abundance, a decrease in α D must correspond to an increase in ε such that the total event rate increases. Thus, lower values of α D give tighter bounds on the dark matter model. For α D > 0.5, the scattering cross section, and thus event rate, increases due to higher order diagrams, thus improving constraints. Y also depends on two masses, m χ and m V . As the production and scattering both only depend on m V , increasing the value of m V /m χ with a fixed m V does not affect our event rate but does give lower values of Y and tighter bounds. We thus assume m V /m χ = 3 according to convention near the threshold for on-shell V decay though there is viable parameter space at lower values. For a given value of m χ , we calculate the ∆χ 2 (Y ) curve relative to the best-fit DM coupling. Allowed values of Y are determined according to the Feldman-Cousins prescription [45] with 90% confidence.
Results: We selected 5142 events in the analysis region of 0 ≤ E rec < 250 PE and 0 ≤ t rec < 6 µs. For each DM mass tested, the best-fit was identical, preferring no DM events in each case with a fit χ 2 /dof = 103.0/120. Our observed best-fit, SSBkgsubtracted spectra, both in the DM timing ROI and the CEvNS background timing region, are shown in Fig. 2 Table I. A summary of prior prediction and best-fit event rates for each background and the 90% limit for 25 MeV/c 2 DM.
The fit prefers slightly fewer CEvNS events than predicted after profiling over nuisance parameters. This is consistent with our CEvNS measurement using the same dataset [30]. Our critical ∆χ 2 values used to construct 90% confidence intervals on N DM are 2.1 − 2.3 depending on m χ . These are slightly lower than those expected from Gaussian statistics due to the proximity to the boundary N DM ≥ 0. At our peak sensitivity, m χ = 25 MeV/c 2 , we determined there are < 15.8 DM events in our sample to 90% CL, though the constraint on the number of DM scatters in our dataset depends on the mass assumption.
Our constraint on DM parameters for α D = 0.5 is shown in Fig. 3 along with our projected sensitivity and other current constraints. The parameters that yield the relic abundance for scalar DM [46] is also shown. With a small 14.6 kg detector, we improve constraints on Y for 11 < m χ < 165 MeV/c 2 by  up to 5× suggesting that future, large-scale CEvNS detectors will be successful in ambitiously limiting light DM models. With the current dataset, we can reject coupling parameters consistent with cosmological DM for masses between 20 and 33 MeV/c 2 assuming α D = 0.5. The constraint is strongest at m χ = 25 MeV/c 2 where we can eliminate the scalar target for all α D < 0.64.
Additionally, as there are few accelerator-based searches for DM that test the DM-quark coupling, we also compare our constraint to both astroparticle and accelerator-based searches of light DM sensitive to the quark coupling. Comparisons to astroparticle results are made by averaging the coherent DM-nucleus cross section given in Eqn. 3 with couplings determined by our constraint over the velocity distribution expected for the DM halo near Earth [42]. This result is also shown in Fig. 3. Our constraint improves on all constraints of the DMquark coupling for masses below 166 MeV/c 2 where COHERENT data probes more than an order of magnitude of previously untested parameter space. At higher masses, astroparticle experiments exploiting the Migdal effect [47][48][49] dominate [50,51] with an additional constraint from CRESST-III [52].
As our constraint depends on our particular choice of α D , we can explore this parameter by constraining the values of α D for which we reject the relic abundance at 90%, as shown in Fig. 4. For a given DM mass, the relic abundance is given by a fixed value of Y . Decreasing α D while holding Y fixed at the relic abundance increases ε ∝ 1/ √ α D such that the overall signal rate expected in COHERENT, which scales like ε 4 α D ∝ 1/α D , increases. Thus, we show the lower bounds for allowed α D . For scalar DM, we constrain the cosmological abundance with very conservative choices of α D . However, if DM is a Majorana or a pseudo-Dirac fermion, significant parameter space remains. In the future, with larger detectors sensitive to lower nuclear recoils, CEvNS  . The values of αD for which we can reject the cosmologically observed DM concentration as a function of DM mass compared to constraints shown in Fig. 3. All parameter space below the contour is excluded. We include scenarios where DM is scalar, a Majorana fermion, or a pseudo-Dirac fermion. A dashed line at αD = 0.5, the value assumed in Fig. 3, is also drawn. data can probe fermion DM models at α D ≈ 0.5, which favor Y values up to 20× lower.
Conclusion: We have exploited CEvNS data collected by our decommissioned CsI[Na] detector at the SNS to search for hidden-sector DM particles produced in the beam. This dataset, in addition to making the most precise measurement of CEvNS to date, has considerably improved on current constraints for accelerator-produced DM particles with masses between 11 and 165 MeV/c 2 . Additionally, this result improves on constraints of the DM-quark coupling for masses below 166 MeV/c 2 . In particular, this is the first result to test scalar DM for even the conservative choice of α D = 0.5 in the studied mass range. The data also constrains Majorana and pseudo-Dirac DM, but constraints on these scenarios will not be as exhaustive until future data is collected. We have developed powerful methods for understanding background rates by exploiting timing information. In the future, these techniques will constrain systematic uncertainties in-situ allowing much more stringent searches. In particular, future argon and CsI detectors placed in the STS beam have significant potential to discover an excess of DM scatters in currently unexplored parameter space independent of DM mass and spin phenomenology.