Limits on Sub-GeV Dark Matter from the PROSPECT Reactor Antineutrino Experiment

If dark matter has mass lower than around 1 GeV, it will not impart enough energy to cause detectable nuclear recoils in many direct-detection experiments. However, if dark matter is upscattered to high energy by collisions with cosmic rays, it may be detectable in both direct-detection experiments and neutrino experiments. We report the results of a dedicated search for boosted dark matter upscattered by cosmic rays using the PROSPECT reactor antineutrino experiment. We show that such a flux of upscattered dark matter would display characteristic diurnal sidereal modulation, and use this to set new experimental constraints on sub-GeV dark matter exhibiting large interaction cross-sections.


I. INTRODUCTION
Despite strong evidence for dark matter's (DM) existence, its particle nature remains unknown, and its identification is one of the most pressing problems in particle physics and astrophysics [1][2][3]. Direct searches for DM, focusing primarily on GeV-scale weakly interacting massive particles colliding with nuclei, are probing ever lower cross sections. However, these searches rapidly lose sensitivity for masses below about 1 GeV: if DM is too light, it does not impart enough momentum to trigger typical detectors. In recent years, a wide range of approaches has been explored in order to probe sub-GeV DM . These include numerous studies of boosted DM, in which light DM is accelerated to high energy through a variety of processes [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39].
One such process that would make light DM detectable is upscattering by cosmic rays (CRs). References [28,29,[34][35][36]39] have explored the experimental signatures of DM particles being struck by CRs, upscattering to high energy, and interacting in direct-detection and/or neutrino experiments. Such analyses have the additional advantage of being insensitive to the cosmogenic DM velocity distribution, a source of uncertainty for traditional direct-detection limits that has received much attention recently [40][41][42][43][44][45][46][47][48]. These analyses have constrained a wide range of parameter space for sub-GeV dark matter, but none of the experiments in question have performed their own analyses aimed at CR-upscattered DM. Using the PROSPECT reactor antineutrino detector, which combines the advantageous features of an on-surface deployment location with excellent particle discrimination capabilities, we have performed the first dedicated experimental analysis constraining sub-GeV DM by considering upscattering by CRs. This analysis, which is also the first to exploit the diurnal sidereal modulation of the boosted DM signal [49,50], addresses regions of parameter space never before probed by terrestrial experiments. While cosmological limits do exist in this parameter space, they are model dependent and thus have faced some debate. Complementary constraints are valuable: while cosmological observables indirectly probe DM scattering in the early Universe, our analysis is based on scattering in the present day, making our analysis more comparable to traditional direct-detection studies than are cosmological limits.
This paper is organized as follows. In Sec. II, we compute the flux of upscattered DM at Earth. In Sec. III, we describe the PROSPECT detector and the experimental search for DM. In Sec. IV, we discuss how we analyze data from PROSPECT to set limits on DM parameter space. In Sec. V, we present our limits.

II. UPSCATTERED DARK MATTER FLUX
If the DM-nucleon scattering cross section σ χN is nonzero, CR nuclei have a chance to collide with DM particles as they propagate in the galaxy. In the DM mass range we consider in this paper, 1 keV < m χ < 1 GeV, CR nuclei carry orders of magnitude more kinetic energy than galactic DM particles, and could upscatter DM particles to high energy. In this section, we compute the spectrum of energetic DM particles reaching Earth after being upscattered by CRs.

A. Cosmic ray and dark matter inputs
Our analysis depends on the Galaxy's DM density profile, the size of the CR halo, and the spectrum of CRs in the Galaxy. We consider only helium and proton CRs, as heavier nuclei are a small fraction of the flux, and including them would only marginally strengthen our limits.
For the DM density, we assume a Navarro-Frenk-White (NFW) profile [51] with scale radius of r s = 20 kpc (6.17×10 20 m) and a density at Earth of 0.3 GeV/cm 3 . Below, we denote the density ρ χ (r, θ, φ) in spherical coordinates centered at Earth, i.e. r is the distance from Earth and θ and φ denote the direction. As we show below, our results are not strongly affected by the differences between a NFW profile and a more shallow or cored profile near the Galactic Center.
Because we consider CRs scattering with DM throughout the Galaxy, we cannot naively employ the CR spectrum measured at Earth. The solar magnetic field suppresses the flux of low-energy CRs reaching Earth, contributing to the break in the observed spectrum at around 1 GeV/nucleon. Instead, we use the local interstellar spectrum (LIS) computed for protons and helium nuclei in Ref. [52]. We assume that the energy distribution of CRs is independent of position in the CR halo, as the shape of the LIS has been shown by gamma-ray observations to be similar to the CR spectrum in other parts of the Galaxy [53,54].
Galactic magnetic fields prevent CRs from simply streaming out of the Galaxy, binding them in a halo that is often modeled as a cylinder with a half-height of a few kpc. The exact value of the half-height is debated, but in this work we adopt a fairly standard value of 4 kpc (see Ref. [55] and references therein), and assume the CR density is independent of height within this halo and zero outside of it. The radial distribution, meanwhile, can be inferred from gamma-ray observations [56,57]. We use for the radial profile the shallowest curve in Fig. 1 of Ref. [57], as it provides the best fit to relevant astrophysical gamma-ray observations. We neglect anisotropies in the CR flux [58], as they are small compared to statistical uncertainties in the data used (see below). We assume the CR density is zero beyond a cylinder radius of 25 kpc. Thus, our model for the CR density is where g(r cr ) is the radial distribution taken from Ref. [57], ρ 0 is a normalization density fit to the LIS at Earth's position in the Galaxy, and r cr and z are in units of kpc, with the origin at the Galactic Center. In our analysis, we do not consider the change in the Earth's position over the course of the data taking, as the amount the Earth moves in approximately two weeks is small compared to galactic distance scales.

B. Cosmic ray-dark matter scattering
As in Refs. [19,28,35], we assume an energyindependent DM-nucleon cross section σ χN . While the feasibility of producing a model with this scaling has been questioned in the literature [37], its use enables useful comparisons to other direct-detection limits. Specific models involving light mediators have also been applied to CR-DM scattering [34,36].
Because CR velocities are so much higher than the velocity of DM bound in the Galaxy, we can treat the DM particles as being at rest. The kinetic energy T χ transferred to a DM particle of mass m χ by a CR of mass m CR and kinetic energy T CR is where θ is the center-of-mass scattering angle [35]. We denote the maximum kinetic energy transferred to a DM particle T max Given the LIS of CR species i, dΦ i /dT i , and the DM density ρ χ , we can compute the spectrum of upscattered DM in terms of CR kinetic energy, The line-of-sight integral extends from Earth to the edge of the CR halo. From here we compute the spectrum in terms of DM kinetic energy T χ : We consider DM boosted to energies between 25 MeV and 1 GeV. The low-energy limit is determined by the analysis threshold (described below), while the highenergy cutoff allows us to neglect quasielastic and inelastic processes [28]. For more discussion and a visualization of the upscattered DM spectrum, see Ref. [35].
In this analysis, we assume the standard cross section scaling for spin-independent DM-nucleon scattering. In the limit of zero momentum transfer, the DM-nucleon cross section σ χN is related to the DM-nucleus cross section σ χA (for a nucleus of mass number A) via the formula where µ χN and µ χA are the DM-nucleon and DM-nucleus reduced masses, respectively. For nonzero momentum transfer, the differential cross section is modified by a form factor F (q 2 ) as For nuclei heavier than hydrogen, we use the standard Helm form factor [59,60], and for protons, we use the dipole form of the hadronic form factor from Ref. [61], both in accordance with Refs. [28,35].
Although this upscattering could have been happening for billions of years, DM boosted to such high energies is no longer gravitationally bound within the galaxy. This prevents a long-term buildup of a high-energy DM in the galaxy, and the flux reaching Earth is relatively small. For DM with a mass of 1 MeV and scattering cross section of 10 −30 cm 2 , the DM flux reaching Earth is approximately 1.26 × 10 −6 cm −2 s −1 sr −1 .
C. Propagation to the detector At the large cross sections we consider, DM may scatter many times while passing through the atmosphere and detector shielding. If the cross section is too large, DM may lose too much energy to be detectable or simply be scattered back out of the atmosphere. The exclusion regions derived in Refs. [28,[34][35][36] all have ceilings, cross sections above which DM would be mostly or entirely blocked from triggering the detector in question. So it is critical to model DM scattering with nuclei as it travels through the atmosphere and detector shielding.
For particles arriving from above the horizon, we simulate propagation through the atmosphere and detector shielding using the same propagation code used in Ref. [35]. This code generates 10 6 -10 7 DM particles at the top of the atmosphere, where each particle is assigned an incoming direction and initial kinetic energy. The kinetic energy is drawn from the spectrum of incoming DM, while the incoming direction is drawn from a distribution of incoming θ and φ, which is in turn based on the direction dependence of the line of sight integral of DM density times CR density. Each particle is then propagated through the atmosphere. The distance to the first collision is drawn from a probability distribution based on the DM's mean free path, the scattering angle is drawn from an isotropic distribution in the center-of-mass frame, and the energy loss and lab-frame scattering angle is computed based on the incoming energy, scattering angle, and target nucleus.
This process is repeated until the particle reaches a sea-level modeled detector location, is scattered out of the atmosphere, or loses too much energy to be detected. The minimal (<1 m water equivalent) concrete overburden provided by the building surrounding the on-surface PROSPECT detector unsurprisingly plays a very minor (percent-level) role in attenuation and down scattering of the DM flux. Reduction in atmospheric overburden due to PROSPECT's ∼ 260 m elevation relative to sea level is similarly negligible.
At the cross sections we consider, any DM particles arriving from below the horizontal direction are completely blocked. This causes the flux at the detector to vary over the course of a day as the Earth's bulk rotates in front of PROSPECT and blocks fluxes from galactic locations with greater or lesser boosted DM density.
The distribution of scattering angles during propagation, and the resulting proton recoil spectra, are determined based on the corresponding nuclear form factors. During propagation, the suppression of the total DMnucleus cross section due to the form factor is neglected, a conservative choice that follows Ref. [28], and has the effect of reducing the flux of very high-energy DM which, as mentioned above, may suffer inelastic collisions. Including such suppression would not change the fact that DM is blocked by the Earth, and would only make our results sensitive to somewhat higher cross sections.
As the flux depends on the location of PROSPECT rel-ative to celestial objects (rather than to the Sun), it is expected that these modulations would exhibit a period of one sidereal day (rather than one solar day). For a given DM mass, cross section, and time, we use the aforementioned propagation code to compute the DM flux at the location of the PROSPECT detector. As detailed below, this allows us to compute the predicted diurnal sidereal modulation of the detected DM signal in PROSPECT and compare it to any time dependence observed in the experiment's data.

III. EXPERIMENTAL SEARCH
The PROSPECT experiment was designed to measure the energy spectrum of antineutrinos emitted from the highly 235 U-enriched High Flux Isotope Reactor (HFIR) reactor at Oak Ridge National Laboratory [62]. Through simultaneous measurements at many <10 m reactordetector distances, PROSPECT is able to probe the existence of short-baseline electron antineutrino disappearance caused by oscillation between active and sterile neutrino states [63,64] with greatest sensitivity in the 1-10 eV 2 scale mass splitting regime [65,66]. By integrating measurements over all baselines, PROSPECT also provides world-leading precision in measuring antineutrino production by 235 U fission products [66,67]. Beyond reactor oscillation and spectrum physics, the PROSPECT detector's on-surface location and powerful particle identification capabilities provide a unique opportunity for DM detection.

A. Experiment and dataset description
The PROSPECT detector consists of an 11×14 array of 1.2 meter long optically isolated segments filled with 6 Li-doped liquid scintillator [68]. Each segment is equipped with a photomultiplier tube (PMT) at each end for collecting scintillation light produced by charged particle interactions. Thin, specularly reflecting segment walls [69] efficiently direct scintillation light towards these PMTs with less than 50% variation in absolute light collection at all points inside a segment. Heavy charged particles producing high ionization density in the PROSPECT scintillator, such as protons, generate a characteristically longer light emission profile than light charged particles, such as electrons, enabling powerful particle identification capabilities via pulse shape discrimination (PSD). The central PSD-capable scintillator detector is surrounded on all sides by tens of cm of passive gamma and neutron shielding. As the primary physics goals of the PROSPECT experiment necessitate it being close to the HIFR core, the detector was deployed on the Earth's surface inside the HFIR building at a location with <1 m water-equivalent overburden. A more detailed description of the PROSPECT experimental layout and detector is available in Ref. [70].
The data analyzed for this paper were collected during ∼14.6 solar days of detector operation from March 16th to March 31st, 2018, and are a subset of the datasets used in Refs. [65,66]. The HFIR reactor was not operating during the entire data-taking period. The triggering of PROSPECT's waveform digitizer (WFD) readout, described in detail in Refs [66,70], occurs when PMTs on one detector segment observe time-aligned waveform features above approximately five photoelectrons in amplitude. To reduce data rates, only sections of digitized waveforms above two photoelectrons in amplitude are recorded to the data stream. With these settings, raw trigger rates were ∼2000 s −1 , with a trigger threshold of ∼75 keV in electron-equivalent energy. The 14bit WFDs, with implemented PMT and electronics gain settings, exhibit linear response below 14 MeV electronequivalent energy, beyond which point electronics saturation results in clipped waveforms. During the data-taking period, one or two PMT channels on 28 of 154 segments had experienced PMT voltage divider instabilities and were turned off. For simplicity, all data from these 28 segments were not considered.

B. Reconstructed physics quantities
PROSPECT low-level data processing, calibration, and physics metric reconstruction is discussed in detail in Ref. [66]. In the present analysis, we take advantage of reconstructed time, segment number, position, energy, and pulse shape variables to select candidate DM interactions. These reconstructed variables are separately assigned to each time-aligned waveform pair in one segment, called a pulse. Time-aligned pulses from different segments are then grouped to form larger data objects, which are referred to as clusters. The reconstructed position of a pulse along a segment (z position) is formed from relative charge integral and arrival time offsets between waveforms; previous analysis has demonstrated a resolution of ∼5 cm or better on reconstructed z positions for pointlike energy depositions.
The reconstructed electron-equivalent energy (MeV) of a pulse is formed using the combined charge information from both PMTs. The absolute scale of reconstructed energy is defined using radioactive gammaray calibration sources and naturally occurring radiogenic and cosmogenic beta-particle and gamma-ray backgrounds [66]. Electron-equivalent energy depositions by heavier charged particles are modeled using Birks quenching parameters [71] fitted to pulse data from triton-alpha products of n-Li capture, radiogenic alpha-particle decays, and collaboration-external measurements of liquid scintillator proton quenching factors [72][73][74]. The accuracy of these external measurements in describing PROSPECT data was further validated using proton recoil spectra generated by deployment of a 241 Am-9 Be fast neutron source inside the central detector.
For pointlike energy depositions in the scintillator, reconstructed pulse energy resolution is dominated by a photo-statistics contribution of approximately (4.5/ √ E)%; reconstructed energy scales are stable to <1% in time and segment number. A reconstructed pulse shape variable is formed from the ratio of late (>44 ns after the leading pulse edge) to total charge (between 12 ns before and 200 ns after the leading pulse edge) in a pulse. For illustration, reconstructed energy versus PSD values for a representative sample of single-pulse clusters are shown in Fig. 1. Three distinct PSD bands are visible, corresponding to electron, proton, and nuclear recoil signatures from low to high PSD value. In the low PSD band, higher event rates at lower energy are clearly visible, including a prominent edge at 2.6 MeV contributed by ambient radiogenic 208 Tl. Low-energy features are also visible in the higher-PSD bands, including a prominent mono-energetic peak from n-Li capture, as well as <1.25 MeV features from 215 Po, 214 Po, 212 Po, and other alpha-particle decays. Higher energy events in all bands arise almost entirely from interactions of cosmic muons and neutrons with the detector. Event rates in the electronlike recoil band far outnumber those in the proton or nuclear recoil bands.

C. Signal selection and background reduction
Using the variables above, we have selected detector signals consistent with DM scattering off of free protons (hydrogen nuclei) in the liquid scintillator. Scattered protons are expected to deposit all their energy in a single segment. For the scattering cross sections probed in this analysis, at most one interaction is expected per incident DM particle, with vertices evenly distributed throughout the detector active volume. As our signal definition, we select reconstructed clusters containing only one pulse, which must have a PSD value consistent with the protonlike recoil band shown in Fig. 1. Allowed PSD ranges were defined by applying three-Gaussian fits to individual energy slices for all single-pulse clusters within 55 cm of the detector z center. Resulting fitted means and 1σ widths are shown in Fig. 2. A representative fit for the 3.0-4.0 MeV energy range is shown in the figure inset. Candidate events are required to have a pulse PSD value within 2σ of the proton recoil band mean, with exact cut values at each energy determined via linear interpolation between fitted values of adjacent energy bins. This cut carries a 5% signal inefficiency, which is propagated through the analysis. PSD parameter distributions are sufficiently time stable as to expect negligible variation in signal efficiency or background contamination over the considered dataset. Due to the expected signal topology described above and the relatively low expected rate of accidentally coincident backgrounds, the single-pulse requirement has negligible associated inefficiency. The single-segment signal topology also ensures that inactive segments do not affect signal selection efficiency or energy response in functioning segments.
As an on-surface detector, PROSPECT is subjected to a high rate of incident cosmogenic neutrons and muons. Both of these particle types are capable of generating single-pulse, high-PSD protonlike signatures in PROSPECT, whether via scattering of primary neutrons from protons in the scintillator, or via scattering of secondary neutrons produced in nearby interactions of those primary particles. To reduce potential backgrounds from these sources, which are expected to dominate after application of all selection cuts, a series of fiducialization and time-coincidence cuts were applied to selected sig-nal candidates. Due to their comparatively higher cross section, incident cosmic or secondary neutrons will preferentially produce proton recoil signatures on the outer edges of the liquid scintillator. Thus, all single-hit clusters reconstructed in segments in the two outermost rows and columns were rejected, as well as those reconstructed further than 20 cm from the z center of the detector. The third-to-bottom row of the detector was also removed from the selection due to higher trigger rates from the detector bottom primarily caused by imperfections in lead shielding coverage.
Cosmogenic time-coincidence veto cuts were optimized to maintain high efficiency (>98%) while still keeping cut lengths substantially longer than the associated physics timescale in question. To reject signals produced by scattering of secondary neutrons, all signals occurring within 5 µs of a preceding muonlike (energy greater than 15 MeV) cluster or 5 µs of a preceding or following proton recoil-like (containing at least one high PSD value pulse) cluster are rejected. All signals occurring less than 500 µs prior to a n-Li capture signal were also rejected. Descriptions of n-Li, neutronlike, and muonlike event class requirements are described in further detail in Ref. [66]. These cosmogenic veto cuts each have an associated dead time of <1.5%, which is corrected for in the analysis. Finally, to reject waveforms truncated by readout window boundaries as well as 212 Po α particle signals in coincidence with preceding 212 Bi γ + β decay signals (0.299 µs half-life), all signals occurring within 2 µs of any other trigger were rejected. This 'pileup' cut has ∼1% associated dead time, which is also corrected for. Negligible (<0.1%) variations in associated veto efficiencies are observed during the data-taking period.
The rare event search performed in this paper rests on the assumption that any observed diurnal sidereal modulation in the rate of detected signal candidates arises from variations in the flux of DM traversing the PROSPECT detector. However, variations in signal-like event rates may also arise from modulations in the flux of incident of cosmic neutrons and muons [66,75,76]: a relative reduction in the flux of incident cosmic neutrons will result in similar reductions in neutron-proton recoils inside of PROSPECT, which are largely indistinguishable from DM-proton recoils. To quantify the level of expected cosmogenic variation in the two week PROSPECT dataset, we use the n-Li capture dataset described above, which occurs at a high rate (∼10 s −1 ) in the detector. We note that PROSPECT has previously demonstrated that differing cosmogenically induced event types, including n-Li captures, show consistent rate fluctuations in response to changing atmospheric pressure [66]. Average n-Li capture rates for each of the 24 sidereal hours in the sidereal day were divided by the total averaged n-Li capture rate to obtain hourly correction factors; these factors were then applied to various signal predictions used to perform the DM exclusion analysis described in the following section. Correction factors for each sidereal hour, depicted in the following section (Fig. 7), have an associated sta-tistical uncertainty of 0.2%, and are all within 1.5% of unity.

D. Final candidate dataset and cross-checks
Efficiency-corrected signal count rates per kg of scintillator obtained after applying all selection cuts are shown in Fig. 3. Event rate reductions from initial protonrecoil-like criteria range from roughly 2 orders of magnitude at high energy to nearly 3 at lower energies, where the PSD requirement largely eliminates previously dominant ambient gamma-ray contributions. Subsequent background cuts contribute almost an additional order of magnitude reduction at most energies. After application of all cuts, DM-like event rates are as low as 5×10 −6 s −1 MeV −1 kg −1 at the highest considered energies in this analysis, or ∼150 d −1 MeV −1 in the 440 kg fiducial detector volume. Below 1.5 MeV, signal rates begin to increase substantially. For this reason, in the DM analysis that follows, only events above 1.5 MeV are considered.
Over the 14.6 solar day dataset, a total of 37522 DMlike signal candidates between 1.5 and 10 MeV (4.8 and 18.5 MeV N R ) are observed. It is expected that this candidate set is dominated by backgrounds consisting of a single proton recoil induced by a scattering cosmic neutron. The clear separation of PSD bands demonstrated in Fig. 2 ensures that electron recoil events, which are both cosmogenic and radiogenic in origin, are subdominant contributors to signal rates. Meanwhile, rates of cosmogenically induced nuclear recoils are lower than that of proton recoils, and are also expected to contribute subdominantly to the final signal dataset.   As shown in Fig. 4, signal events are relatively evenly distributed among the different fiducial detector seg-ments, with per-segment signal rates of ∼0.5 ms −1 between 1.5 and 10 MeV. This figure also illustrates the locations of inactive segments not used in the analysis, as well as the higher rates of signal-like events in nonfiducial segments. To test for unforeseen variations in signal selection efficiency or improperly estimated cosmogenic flux variations, we compared populations of signal-like events occurring only in nonfiducial segments, since these signals, similar to those from the fiducial volume, are expected to be dominated by conventional cosmogenic neutron backgrounds. Two samples of roughly equal live time were formed from events with time stamps between the hours of either 22:00 to 02:00 or 10:00 to 14:00 Greenwich Mean Sidereal Time (GMST), which represent the periods of highest and lowest expected DM fluence through PROSPECT. These nonfiducial datasets contain ∼20 times higher statistics than their counterpart fiducial signal datasets.
The ratio of energy spectra between these two test data samples are plotted versus reconstructed energy in Fig. 5. Good consistency in event rates can be seen across the 1.5-10 MeV energy range of interest within the statistical limitations of the dataset. A flat-line fit to this ratio provides a best fit of 0.987±0.003, in agreement with the 0.988 value expected based on the hourly correction factors calculated above using the cosmic n-Li dataset. This agreement indicates no unexpected diurnal sidereal modulation in background rates within the statistical limitations of this comparatively large dataset. To test for possible modulation in energy spectrum shape, a linear polynomial was fit to the ratio between datasets. The best-fit slope parameter is within 1 standard deviation of zero, as would be expected from an absence of modula- tion. These observations suggest than any diurnal sidereal modulation of conventional origin in the signal DM dataset is negligibly small compared to the statistical uncertainty of the signal dataset. This time-modulation cross-check was also performed for otherwise signal-like events with PSD parameters within 2σ of the the electron recoil band center shown in Fig. 2. This dataset also yielded an event rate ratio, 0.985±0.002, consistent with a lack of unexpected modulation.

IV. DARK MATTER SEARCH RESULTS
PROSPECT records a substantial rate of DM candidate events, most of which are presumably due to Standard Model backgrounds. However, we can still exclude strongly interacting DM by searching for the expected diurnal sidereal modulation of the signal, an effect that has been explored-but not yet used to set limitsin Refs. [49,50,[77][78][79][80]. Because of the large DM cross sections we consider, DM arriving from below the horizon is blocked by the Earth, and of the flux arriving from above the horizon, the further from vertical the initial trajectory is, the more the flux will be attenuated. As the Earth rotates, so does the vertical direction at PROSPECT's location, and the Galactic Center varies from being about 25 degrees above the horizon to being completely blocked by the Earth. As the flux of upscattered DM reaching PROSPECT scales with the line-ofsight integrated DM density times CR density, which is highest in the direction of the Galactic Center, this rotation produces a characteristic modulation over the course of a sidereal day. When a DM particle collides with a proton in the detector, the reconstructed energy in PROSPECT is related to the recoil energy of the proton by an energydependent quenching factor, as discussed in Sec. III. For modeling reconstructed energies of the DM-induced proton recoil signal, we use the "Birk9" fit from Ref. [81]. Figure 6 shows the predicted DM-induced event rate in PROSPECT for a few example DM masses and cross sections representing diverse regions of the parameter space that PROSPECT is sensitive to. Spectra are shown for times integrated from 00:00 to 01:00 GMST, close to when the DM signal is predicted to be strongest. Above a reconstructed energy of 1 MeV, the predicted shape of the DM signal's reconstructed energy spectrum exhibits a gradual downward trend with increasing energy, similar to the observed spectrum. This spectrum is determined just by the incoming DM spectrum (after attenuation) and the kinematics of the collisions, analogous to the computation of the DM spectrum in Eq. 5. While the normalization of the predicted DM signal varies substantially across the DM phase space region of interest, its spectrum shape appears to be relatively consistent across this space. Due to the predicted spectrum's relative flatness and lack of finer-scale features, predicted DM signal event counts in the 1.5-10 MeV reconstructed energy range of interest are largely insensitive to other aspects of detector response, such as PROSPECT's photo-statistics or geometry-dependent energy resolution contributions. As uncertainties in the assumed proton quenching model are correlated across all time bins in the analysis, they also play a negligible role in defining the exclusion limits of this analysis. The latter point was verified by checking for consistency in DM exclusion contours between analyses incorporating each of the different proton quenching models referenced in this paper [72][73][74]81]. As illustrated in Fig. 6, the DM-induced event rate is smaller than the observed rate of signal-like events in PROSPECT in much of the parameter space of interest. This necessitates a time-binned analysis, as mentioned above and in Sec. II. Fig. 7 shows the signal event rate reported by PROSPECT in the 1.5-10 MeV energy range, plotted by time of sidereal day in 24 bins one sidereal hour in width. We note that, due to the shortness of the dataset relative to a solar year and relatively small frequency offset between sidereal and solar time (4 min per solar day), the phase offset between solar and sidereal time in this analysis is roughly consistent. For example, the sidereal time of highest expected DM signal, 23:00 to 00:00 GMST, corresponds to roughly 07:24 to 08:24 and 6:26 to 7:26 Eastern Daylight Time at the beginning and end of the dataset, respectively. Qualitatively, the data show no obvious indications of diurnal sidereal modulation.
To quantitatively determine which regions of DM parameter space produce signal-like event rate modulations inconsistent with the observed data as plotted in Fig. 7, we define the following test statistic: In this test statistic, χ 2 const is defined as a one-parameter flat-line χ 2 fit to the data as binned in Fig. 7; χ 2 DM is defined by a predicted modulating DM contribution specific to each (m χ ,σ χN ) phase space point, added to a fitted flat-line background contribution. Rates for χ 2 const and χ 2 DM predictions are corrected to account for expected percent-level variations in cosmogenically produced signal-like backgrounds, as described in Sec. III; these corrections are illustrated in Fig. 7. This figure also depicts signals for χ 2 const and for χ 2 DM for two test points in dark matter phase space prior to the application of rate correction factors. The black line depicts the best constant fit (minimum χ 2 const /d.o.f. = 35.1/23) with respect to the data, which corresponds to the expected signal from a Standard Model background free from modulating DM effects. The red and blue curves represent χ 2 DM for test points (m χ , σ χN ) = (1 MeV, 3×10 −28 cm 2 ) and (1 MeV, 5×10 −28 cm 2 ), respectively. Minimizing over the remaining parameter gives best-fit background rate contributions of 69.3 ×10 −6 s −1 kg −1 and 65.1 ×10 −6 s −1 kg −1 for these test points, respectively. Comparing this time-independent background rate with the time dependence of the total event rate apparent in Fig. 7, one can see that the DM event rate roughly doubles, from minimum to maximum, over the course of a sidereal day. For example, for a cross section of 5×10 −28 cm 2 , subtracting the background from the total event rate yields a DM event rate that varies from roughly 5.1 ×10 −6 s −1 kg −1 to 10.9 ×10 −6 s −1 kg −1 . Both DMincluding test points provide relatively poor fits to the observed data, with minimum χ 2 DM /d.o.f. of 60.1/23 and 103.1/23 for the red and blue curves, respectively.
By performing similar tests at an array of (m χ , σ χN ) phase space points, we have determined excluded regions of dark matter parameter space using PROSPECT data as shown in Fig. 8. To assign exclusion confidence intervals, we use the Gaussian CL s method [85], which is useful in the context of performing searches for new physics in a continuous parameter space with large sample sizes. The CL s value determined by PROSPECT's dataset x (the data points in Fig. 7) for each phase space point is defined as Here, ∆T (x) is ∆χ 2 (x) (Eq. 8) for PROSPECT's measured dataset, ∆T 0 is χ 2 DM (x H0 ), where x H0 denotes the Asimov dataset following the modulation-free hypothesis, and ∆T 1 is −χ 2 const (x H1 ), where x H1 denotes the Asimov dataset following the dark matter signal for the phase space point in question. Phase space points with CL s values lower than 0.05 are disfavored by the data at the 95% confidence level. The darkly shaded PROSPECT exclusion region in Fig. 8  projection derived by Ref. [35], which assumed both significantly reduced background and improved background modeling. Taking advantage of the daily modulation of the DM signal was crucial to reach this sensitivity. The exclusion's lower limit is defined by the low fraction of incident dark matter flux interacting within the detector, while its upper limit is defined by attenuation of the dark matter flux prior to reaching the active detector region. Due to the relatively similar spectrum shapes between background and signal, negligible additional exclusion power is provided through a finer binning of the statistical analysis in energy. Expanded energy ranges for the analysis also offer limited improvement in exclusion power due to increased background rates at lower energies and low statistics at higher energies.
The strength of the DM exclusion does depend on assumptions about the half-height of the galactic CR halo. We adopted a commonly used value of 4 kpc, but estimates range from roughly 3-7 kpc (see Ref. [55] and references therein). Adjustments of halo half-height within this range result in reduction or expansion of limits at low cross section by less than a factor of 2; exclusions at high cross section are largely unaffected. Similarly, the daily modulation of DM events depends on the assumed DM profile. We tested the robustness of our results by comparing two alternative DM halos: a more concentrated NFW profile with a scale radius of 10 kpc, and an extremely cored model, which follows an NFW profile at large radii but has constant density within about 8 kpc. The more concentrated NFW profile increases the am-plitude of the daily modulation, while the cored profile decreases it, but in both cases only by tens of percent. This produces an O(10%) variation in the strength of our limit at low cross sections.
We have also considered whether the observed hourly signal-like rates in Fig. 7 are consistent with a sinusoidal modulation beyond that allowed by the boosted dark matter signal of central concern in this paper. A best-fit modulation was found with an amplitude 1.46% of the total signal rate and a phase approximately 12 hours behind the DM signal. This fit, which has the phase and amplitude of the sinusoid as free parameters and thus includes 2 fewer degrees of freedom, provided a χ 2 of 30.72, 4.33 below that of the χ 2 const fit described above. A frequentist approach was then employed to determine the strength of this apparent preference towards a modulated signal. Similar ∆χ 2 values were calculated for 10 3 simulated modulation-free PROSPECT datasets with statistical fluctuations matching those of the observed dataset. The ∆χ 2 value of 4.33 derived from the observed data is lower than 22.2% of values from simulated datasets, indicating that the observed data is consistent with a lack of daily modulation.

V. CONCLUSIONS
In summary, we have used a dedicated analysis of 14.6 solar days of the PROSPECT neutrino experiment's reactor-off data to provide new bounds on the nature of dark matter. This search is enabled by PROSPECT's unique experimental configuration, which combines onsurface detector deployment with powerful particle discrimination and event topology reconstruction capabilities. After applying analysis cuts designed to select isolated proton recoil signatures within 440 kg of target liquid scintillator, we have identified 37522 candidate interactions of energetic dark matter upscattered by cosmic rays. As signal detection rates do not exhibit any statistically significant degree of diurnal sidereal modulation, as would be produced by strongly interacting dark matter originating in the galactic halo, we are able to exclude the existence of dark matter over a broad range of sub-GeV parameter space. This new constraint addresses phase space regions beyond those accessible in traditional or low-threshold direct-detection experiments, and reaches cross sections about an order of magnitude larger than those previously probed in other studies of cosmic ray upscattered dark matter. Our limit is complementary to existing constraints from cosmology and structure formation: while these other limits are indirect constraints based on scattering in the early Universe, our result is a direct-detection limit based on scattering in the present day. In the future, longer data collection times and improved background rejection could extend sensitivity substantially at low cross section, but only modestly at high cross section, where useful signatures are sharply cut off by extreme atmospheric attenuation of the incoming DM flux. Mild improvements in sensitivity at high cross section may also be achieved through redeployment at high (>km) elevations.