Strong New Limits on Light Dark Matter from Neutrino Experiments

TL;DR

This paper demonstrates how neutrino experiments can set new limits on sub-GeV light dark matter by analyzing cosmic-ray interactions and relativistic dark matter flux, expanding the exclusion parameter space.

Contribution

It introduces novel analysis methods and data sources to improve constraints on light dark matter from neutrino experiments, surpassing previous limits.

Findings

Neutrino experiments exclude larger parameter space for sub-GeV dark matter.

Cosmic-ray elastic scattering can produce detectable relativistic dark matter flux.

Strategies to further enhance sensitivity to light dark matter are proposed.

Abstract

The non-detection of GeV-scale WIMPs has led to increased interest in more general candidates, including sub-GeV dark matter. Direct detection experiments, despite their high sensitivity to WIMPs, are largely blind to sub-GeV dark matter. Recent work has shown that cosmic-ray elastic scattering with sub-GeV dark matter would both alter the observed cosmic ray spectra and produce a flux of relativistic dark matter, which would be detectable with traditional dark matter experiments as well as larger, higher-threshold detectors for neutrinos. Using data, detectors, and analysis techniques not previously considered, we substantially increase the regions of parameter space excluded by neutrino experiments for both dark matter-nucleon and dark matter-electron elastic scattering. We also show how to further improve sensitivity to light dark matter.

Equations

$$T_{\chi }=\frac{T_{\mathrm{CR}}^2+2m_{\mathrm{CR}}{T}_{\mathrm{CR}}}{T_{\mathrm{CR}}+(m_{\mathrm{CR}}+m_{\chi }{)}^2/(2m_{\chi })}\left(\frac{1-\cos \theta }{2}\right),$$

$$T_{\chi }^{\max }=\frac{T_{\mathrm{CR}}^2+2m_{\mathrm{CR}}{T}_{\mathrm{CR}}}{T_{\mathrm{CR}}+(m_{\mathrm{CR}}+m_{\chi }{)}^2/(2m_{\chi })}.$$

$$T_{\mathrm{CR}}^{\min }=\left(\frac{T_{\chi }}{2}-m_{\mathrm{CR}}\right)\left(1\pm \sqrt{1+\frac{2T_{\chi }}{m_{\chi }}\frac{(m_{\mathrm{CR}}+m_{\chi }{)}^2}{(2m_{\mathrm{CR}}-T_{\chi }{)}^2}}\right),$$

$$\frac{d{\Phi }_{\chi }}{d{T}_i}=\int \frac{d\Omega }{4\pi }{\int }_{l.o.s.}dl{\sigma }_{\chi i}\frac{{\rho }_{\chi }}{m_{\chi }}\frac{d{\Phi }_i}{d{T}_i}.$$

$$\frac{d{\Phi }_{\chi }}{d{T}_{\chi }}={\int }_0^{\infty }d{T}_i\frac{d{\Phi }_{\chi }}{d{T}_i}\frac{1}{T_{\chi }^{\max }(T_i)}\Theta [T_{\chi }^{\max }(T_i)-T_{\chi }].$$

$$\frac{d\Gamma }{d{T}_k}={\sigma }_{\chi k}{\int }_{T_{\chi }^{\min }}^{\infty }d{T}_{\chi }\frac{1}{T_k^{\max }}\frac{d{\Phi }_{\chi }^z}{d{T}_{\chi }},$$

$$\frac{d{T}_{\chi }}{dx}=-\sum _j{n}_j{\int }_0^{T_r^{\max }}d{T}_r\frac{d{\sigma }_{\chi j}}{d{T}_r}{T}_r,$$

$$\frac{d{T}_{\chi }}{dx}=-\frac{1}{2}\sum _i{n}_j{\sigma }_{\chi j}{T}_r^{\max }.$$

$${\sigma }_{\chi A}={\sigma }_{\chi N}{A}^2{\left[\frac{m_A(m_{\chi }+m_N)}{m_N(m_{\chi }+m_A)}\right]}^2{G}_A^2(Q^2),$$

$$G_A(Q^2)=1/(1+Q^2/{\Lambda }_A^2).$$

$$T_e(T_p)={\int }_0^{T_p}\frac{d{T}_p}{1+k_B⟨d{T}_p/dx⟩},$$

Full text

I Introduction

Though it makes up most of the mass in the universe, dark matter (DM) is only known to interact gravitationally. As a result, its particle mass and scattering cross sections are unknown. Direct-detection experiments, collider searches, and a wide array of cosmological and astrophysical studies have searched for signs of DM interacting with either nucleons or electrons with no clear signals to date dd7 ; dd8 ; dd9 ; dd6 ; Kahrev ; ATLAS17 ; CMS18 ; Lea18 .

Direct detection experiments set strong limits on DM heavier than $\sim 1$ GeV, but rapidly lose sensitivity below this point, because DM that is too light doesn’t carry enough momentum to trigger them. This means different approaches are needed to search for light DM. Cosmology and astrophysics set limits on the DM-nucleon cross section for low masses, but these are of order $\sim 10^{-28}\,$cm2, more than fifteen orders of magnitude weaker than direct detection constraints on GeV-scale DM Nad19 ; Bod18 ; Glu18 ; Sla18 ; Xu18 ; Bho18a ; Bho18b ; Cap19 ; Wad19 . Colliders probe much smaller cross sections, but these limits are more model dependent, and more importantly, there is a ceiling above which DM would interact in the detectors, making traditional missing-energy searches insensitive Bai11 ; Coh15 ; Dac15 . Additional, model-dependent limits may be set based on energy deposition in white dwarfs or neutron stars Bar17 ; Raj17 ; Das19 , Big Bang Nucleosynthesis Krn19 ; Bon19 , DM production in cosmic ray (CR) showers Alv19 , the cooling of stars and supernovae Raf99 ; Fay06 ; Pos08 ; Raf08 ; Guh18 , and Casimir-Polder type forces Fic18 ; Bra18 . There is a large gap between collider and astrophysical/cosmological limits that needs to be probed; see Fig. 1 of Ref. Cap19 .

The DM-electron elastic scattering cross section can be probed by direct detection experiments for masses well below a GeV, but these experiments reach a similar kinematic limit around 1 MeV. Cosmological and astrophysical limits on DM-electron scattering exclude cross sections above $\sim 10^{-28}\,$cm2 Ali15 ; Cap19 ; Wad19 (but neutron star heating may provide much stronger constraints Bel19 ). And collider limits (e.g., Refs. Fox11 ; Ess13 ) should have a ceiling analogous to the one in the nucleon case, though as far as we know, it has not been calculated. For the electron case as well, new ideas are needed to close the window between cosmological and collider limits.

Recently, we showed that for allowed values of the DM-proton and DM-electron elastic scattering cross sections, CRs would lose enough energy in collisions with light DM to alter their observed spectra, and used this observation to set strong limits on these cross sections Cap19 . Figure 1 shows a schematic diagram of this scattering process. Following our paper, Refs. Bri18 and Ema18 considered the complementary effect: collisions with CRs can upscatter DM particles to relativistic or near-relativistic speeds, allowing DM particles much lighter than a GeV to produce detectable recoils even in relatively high-threshold detectors. Ref. Bri18 used data from MiniBooNE, Borexino, and Xenon1T to derive large new exclusion regions on sub-GeV DM scattering with nucleons; similarly, Ref. Ema18 used data from MiniBooNE and Super-Kamiokande (Super-K) to set strong new constraints on sub-MeV DM scattering with electrons. (These exclusion regions are complementary to those recently derived based on solar reflection An18 ; Emk17 ; see below.) We emphasize that DM upscattering by CRs is very different from the scenario usually referred to as “boosted” DM, in which energetic DM particles of one species are produced by pair annihilation of a second, heavier DM species Aga14 ; Cha18 ; Ha18 ; McK18 .

Here we expand on the results of Refs. Bri18 and Ema18 by considering detectors and data sets that these papers did not use, and analyzing this additional data more precisely. We use our own numerical code to more carefully model attenuation of the DM flux by nuclei in Earth’s crust. For DM-nucleon scattering, we consider data from Daya Bay, KamLAND, and PROSPECT, and for DM-electron scattering, we consider lower-energy Super-K data than considered in Ref. Ema18 .

II From Cosmic Ray Spectra to Recoil Distributions

In this section, we introduce the basics of CR propagation and compute the recoil spectrum seen in a detector for a given CR spectrum. This overview sets the foundation for our new results in subsequent sections.

II.1 Cosmic Ray Inputs

CRs are energetic charged particles, mostly protons but also including electrons and a range of heavier nuclei, which, at the energies we consider, are thought to largely be accelerated in supernova remnants Bla11 ; Bla13 ; Mor17 ; Byk18 . CRs are energetic enough to upscatter DM particles to relativistic speeds; the observed proton and electron CR spectra at Earth both peak at around 1 GeV, ultra-relativistic energy for electrons and moderately relativistic for protons, and extend many orders of magnitude higher in energy.

CRs are trapped by galactic magnetic fields in a thick, disklike halo around the galactic disk for far longer than it would take to cross the galaxy in a straight line, and their propagation within this halo is described by a diffusion equation. The size of the halo is uncertain, but we adopt the relatively conservative geometry used in Ref. Ema18 , assuming CRs are uniformly distributed in a cylinder with radius $R=10$ kpc and half-height $h=1$ kpc; considering a larger volume, which is probably more realistic (e.g., Ref. Str00 ), would produce stronger limits because the upscattered DM flux would be larger.

The interstellar CR spectrum is different from the spectrum observed at Earth due to solar modulation. For energies below a few GeV, solar modulation suppresses the CR flux at Earth, contributing to the peak mentioned above. Voyager 1 has measured the local interstellar spectra (LIS) of CR electrons and various nuclei down to energies of 1–10 MeV Cu16 , and recent papers have computed the LIS down to the lowest energies of the Voyager electron, proton, and helium spectra Bos17 ; Bos18a ; Bos18b . The LIS has been shown to agree well with the CR spectra elsewhere in the galaxy as inferred by gamma-ray observations, with the CR density increasing somewhat at smaller galactic radii Ack11 ; Yan16 ; Aha18 ; we therefore use for the galactic CR spectra the LIS computed in Refs. Bos17 ; Bos18a . In constraining DM scattering with nuclei, we restrict our attention to the CR proton and helium spectra.

II.2 CR-DM Scattering

Following Refs. Cap19 ; Bri18 ; Ema18 , we assume that both the DM-electron and DM-nucleon cross sections are energy independent, and that DM-electron scattering is isotropic in the center-of-momentum (CM) frame, while DM-nucleus scattering deviates from isotropy only due to nuclear form factors. The assumption of energy independence is a simplification that allows for straightforward comparison between our results and constraints from astrophysics, cosmology, and traditional direct detection. In Sec. VI, we discuss the implications of this choice and possible extensions, such as considering inelastic scattering with nuclei. Following Ref. Bri18 , we expect that the elastic scattering cross section, enhanced by a factor of $A^{2}$, is dominant in the DM energy range we consider.

Typical DM velocities in the galaxy ($\sim 10^{-3}c$) are small compared to the velocities of the CRs we consider, so we treat the DM as being at rest. The kinetic energy transferred to a stationary DM particle of mass $m_{\chi}$ by a CR with mass $m_{\rm CR}$ and kinetic energy $T_{\rm CR}$ is

[TABLE]

where $\theta$ is the CM scattering angle. Consequently, the maximum recoil energy is

[TABLE]

Inverting this equation gives the minimum CR energy, $T_{\rm CR}^{\rm min}(T_{\chi})$ required to produce a DM recoil energy $T_{\chi}$:

[TABLE]

where the $+$ applies for $T_{\chi}>2m_{\rm CR}$ and the $-$ applies for $T_{\chi}<2m_{\rm CR}$.

II.3 DM Flux and Spectrum

For the DM density profile, we use a Navarro-Frenk-White (NFW) profile Nav96 with scale radius $r_{s}=20$ kpc and a density at Earth of 0.3 GeV/cm3. This is the same scale radius used in Ref. Ema18 , though they use a local density of 0.42 GeV/cm3. The difference in density has only a mild impact on our results, as discussed below.

The differential flux at Earth of DM upscattered by collisions with CRs of species $i$ (in terms of incident CR energy $T_{i}$) is given by a line-of-sight integral:

[TABLE]

We integrate over the full CR halo with the geometry given above. Note that this flux is in terms of the CR energy, not the DM energy. To convert this into a DM energy spectrum, we integrate over CR energies:

[TABLE]

Figure 2 shows the DM spectra reaching Earth after collisions with either protons (plus helium) or electrons, for several masses. For $m_{\chi}\ll m_{\rm CR}$, the proton-induced spectra show breaks at $T_{\chi}\simeq m_{\chi}$; in the limit where $m_{\rm CR}>m_{\chi}$ and $m_{\rm CR}>T_{\chi}$, the term in the square root in Eq. (3) is $1+T_{\chi}/(2m_{\chi})$, leading to the observed break. For $m_{\chi}\simeq 1$ GeV, the proton and helium form factors begin to matter, causing the break to be slightly lower in energy than $m_{\chi}$. For the electron-induced spectra, because $m_{\rm CR}$ is not necessarily greater than $T_{\chi}$, the corresponding break is roughly compensated by an additional break coming from the factors of $2m_{\rm CR}-T_{\chi}$, which for protons showed up at too high energy to be relevant. The break that does appear is due to the break in the electron LIS, which was included in order to fit the low-energy Voyager data. It is also interesting to note that for light DM, the electron-induced flux is higher at high energy than the proton-induced flux, despite the proton CR flux being higher than that of electrons. This is due to electrons transferring a larger fraction of their energy to light DM than protons because they are closer in mass to the DM.

The incoming DM flux should have significant directional variation: the highest flux should come from the direction of the galactic center, where the line-of-sight integrated DM density is highest. However, scintillator detectors (and typical DM detectors) lack the directional sensitivity to use this information. We discuss below how Super-K’s directional sensitivity could be useful for improving constraints on DM-electron scattering.

II.4 Attenuation of the DM Flux

Direct-detection experiments are blind to DM with sufficiently large cross sections, as it would be stopped from reaching the detectors by scattering in the atmosphere, Earth, and detector shielding. This effect is often neglected for direct detection experiments, as they are designed to probe such low cross sections that Earth is effectively transparent to DM. However, for sub-GeV DM, existing limits are weak enough that for the cross sections we hope to probe, attenuation may be significant.

References Bri18 ; Ema18 ; Dent19 modeled attenuation using a ballistic-trajectory approach that assumes that DM travels in a straight line from the top of the atmosphere to the detector, losing energy as it scatters but not changing direction. This approach is reasonable for DM much heavier than the target particles (as in Ref. Kav17 ), but for light DM ($m_{\chi}\lesssim m_{A}$), it is clearly not well motivated: particles may backscatter into the atmosphere, and even DM that eventually reaches the detector may have a much longer path length than the straight-line distance to the detector Mah17 . On the other hand, Ref. Bri18 neglects the form factor of nuclei in Earth’s crust, which should make their results more conservative. While numerical codes have been used to study attenuation of nonrelativistic DM Mah17 ; Mah17b ; Emk17c ; Emk18 , no such code has previously been used to study CR-upscattered DM.

To model propagation of relativistic or nearly relativistic DM, we wrote our own propagation code, modeled somewhat on the publicly available DMATIS code Mah17 ; DMATIS but including relativistic kinematics and the kinetic energy distribution of CR-upscattered DM. Our purpose here is not to model attenuation in a detector-independent way, but rather to compute the recoil spectrum in a particular detector, which can then be compared to data. In the next section, we describe our code and compare its results with those of the ballistic approximation.

II.5 Target Recoil Distribution

Denoting the DM flux at detector depth $z$ as $d\Phi_{\chi}^{z}/dT_{\chi}$, the differential recoil rate per target particle $k$ (nucleus or electron) is then

[TABLE]

where $T_{k}^{\rm max}$ and $T_{\chi}^{\rm min}$ are obtained from Eqs. (2) and (3) respectively, by replacing $\chi\rightarrow k$ and CR $\rightarrow\chi$.

It is important to note that because $d\Phi_{\chi}^{z}/dT_{\chi}$ contains a factor of the DM-proton or DM-electron cross section, the event rate scales as cross section squared: one factor of the cross section from the DM being struck by a CR, and one factor from the DM interacting in the detector. As discussed below, this makes our final results fairly insensitive to astrophysical uncertainties because they scale only with the square root of astrophysical inputs.

III Propagation of Dark Matter to the Detector

The cross sections we probe are large enough that dark matter may be significantly slowed or deflected by Earth’s crust and atmosphere before reaching a detector. Unlike previous work, which employed an analytic but physically inaccurate ballistic-trajectory method to model attenuation, we develop our own numerical propagation code, using this code to treat attenuation and comparing its results with the ballistic model.

III.1 Ballistic-Trajectory Approximation

As discussed above, studies of DM attenuation in Earth often use a simple, ballistic-trajectory approximation to model the DM energy loss in the crust, atmosphere, and/or detector shielding Sta90 ; Kou14 ; Bri18 ; Ema18 ; Dent19 . In this formalism, the energy loss rate (in units of MeV/cm) is

[TABLE]

where $T_{r}$ is the energy lost by the DM particle in a collision and the sum is over all relevant target particles (such as various nuclei). For isotropic scattering, $d\sigma_{\chi j}/dT_{r}=\sigma_{\chi j}/T_{r}^{\rm max}$, which gives

[TABLE]

For spin-independent DM-nucleus scattering, the scattering cross section $\sigma_{\chi A}$ is

[TABLE]

where $G_{A}(Q^{2})$ is the nuclear form factor.

This energy loss rate is then used to compute the spectrum at depth $z$, given the spectrum at the top of the atmosphere. However, while the ballistic-trajectory approach is reasonable for sufficiently heavy DM, it is not a good qualitative or quantitative description of the propagation process for light DM. To more accurately model DM propagation through the atmosphere, crust, and detector shielding, we use our own numerical code.

III.2 Description of Numerical Code

Our code models the atmosphere as being composed of 21% oxygen and 79% nitrogen, neglecting additional elements. It models Earth’s crust as being composed of oxygen, silicon, aluminum and iron, with the same abundances used in the DMATIS code DMATIS . Any concrete shielding (relevant for PROSPECT, discussed below) is modeled as having the same composition as the crust.

Our code begins by selecting a kinetic energy from the input DM energy distribution, and initializing a DM particle at the top of the atmosphere, moving toward the detector. In each of the regions defined above, the code computes the mean free path based on the density of material, elemental abundances, and DM-nucleon cross section. It samples from the resulting path length distribution to determine the distance traveled before one interaction, and selects the nucleus encountered in the first scattering with the probability of each nucleus weighted by $A^{2}$ (from the scaling of the cross section) times the abundance (which scales as $1/A$ for a given mass density). It then uses the form factor for that nucleus and the DM energy to determine the scattering angle and energy loss of that collision. For elements in the crust, shielding, and atmosphere, we employ the commonly used Helm form factor Hel56 ; Dud06 .

This process is repeated until the particle either scatters back into the atmosphere, loses too much energy to produce an event in the energy range we consider, or reaches the depth of the detector. If it reaches the depth of the detector, the recoil energy is computed based on the DM energy and the form factor of the target nucleus. This process is repeated for $10^{6}$ particles, and the recoil spectrum plotted. The recoil spectrum is then multiplied after the fact by the probability of a particle actually interacting with a nucleus in the detector.

III.3 Comparison of Numerical Method and Ballistic Approximation

References Mah17b ; Emk18 have compared the results of numerical codes and the ballistic-trajectory approximation, and generally speaking found the ballistic-trajectory approach to be conservative for the detectors and mass ranges they considered. However, these conclusions cannot be naively applied to lower-mass, relativistic DM in higher-threshold detectors. Thus it is necessary to do our own comparison.

As an example, we compute the DM-induced proton-recoil spectrum in Daya Bay, at a mass of 1 MeV and a cross section of $5.0\times 10^{-28}$ cm2, using both the ballistic approximation (neglecting form factors) and our numerical code. We propagate $10^{5}$ DM particles to a depth of 100 m, modeling only particles with kinetic energy up to 1 GeV. We then compare the predicted spectra of DM-induced events in Daya Bay with the observed event spectrum. (In the next section, we provide details of the Daya Bay experiment and describe how we use their data to constrain DM properties.)

Figure 3 shows the spectrum of DM-induced events computed using the ballistic approximation (blue) and our numerical code (black). The ballistic approximation overpredicts the event rate at low energy because it neglects changes in direction during propagation. For low-energy DM particles, say with energy $\sim 200$ MeV when considering MeV DM scattering with silicon, the form factor’s deviation from unity is largely negligible, and particles scatter nearly isotropically. More than half of these particles are backscattered out of the atmosphere, and even the particles that reach the detector’s depth scatter far more times than the ballistic trajectory approach assumes, and lose most of their energy. On the other hand, higher-energy particles (above a couple hundred MeV) are weighted much more heavily toward forward scattering and lower energy loss. Because the ballistic approximation underpredicts the DM-induced event rate where the background rate is lowest, our numerical code finds that the ceiling should be higher than would be computed using the ballistic approximation.

The difference between the ballistic and numerical treatments of attenuation is important not just for computing the ceiling of an exclusion region, but also for determining what data to use to set that ceiling. Ref. Dent19 , using the ballistic approximation, shows that the ceiling of their XENON1T exclusion region would move up noticeably if the threshold were reduced to 1 keV. This is almost certainly due to the high-energy cutoff in the DM spectrum seen in the analytic approximation. However, numerically computing attenuation does not produce this same cutoff, making the threshold energy less important.

IV Constraining the DM-Nucleon Cross Section

In this section, we use the above formalism, along with data from Daya Bay, KamLAND, and PROSPECT, to set new limits on the DM-nucleon cross section $\sigma_{\chi N}$. As above, we consider only elastic scattering on protons in the detectors. We do not consider elastic scattering on carbon nuclei because the effect of quenching is much greater for carbon than for protons (see below), but in Sec. VI we discuss inelastic collisions with various nuclei.

IV.1 Our Improvements Over Previous Work

On top of properly modeling attenuation, we build on previous work by considering additional data and detectors, allowing us to probe additional parameter space. While the exclusion regions of Ref. Bri18 are broad in both mass and cross section, there remains a gap between the MiniBooNE and XENON1T regions, which we almost completely close. Using data from Daya Bay, we almost completely close this gap. KamLAND, the deepest detector we consider, has the lowest ceiling but also the lowest background rate, making it sensitive to cross sections much lower than MiniBooNE and comparable to XENON1T. Finally, we compute projected sensitivities for PROSPECT, whose position at the Earth’s surface will make it sensitive to cross sections higher than MiniBooNE, and JUNO, a future detector with the potential to probe cross sections lower than XENON1T.

Typical direct detection experiments use heavy nuclei such as xenon to probe small cross sections, because of the cross section’s commonly assumed scaling with mass number. However, this scaling need not hold in general: DM may couple differently to protons and neutrons, and either way the scaling may break down at large cross sections Dig19 . It is therefore valuable, if possible, to constrain the DM-proton cross section through scattering with individual protons. By focusing on proton recoils in scintillator detectors, we are able to probe the DM-proton cross section directly, making our results complementary to the XENON1T limits.

We compute exclusion regions for these detectors for DM masses from 1 keV to 1 GeV. Masses below about 1 keV are disfavored because of constraints from structure formation Aba05 ; Vie06 ; Boy08 , and fermionic DM cannot be lighter than around 100 eV due to phase space constraints Boy09 ; DiP17 . However, our formalism is perfectly valid for masses far below 1 keV, and could thus constrain models that could evade the cosmological constraints.

IV.2 Proton Recoil Distribution

We consider proton and helium CRs elastically scattering with DM particles, neglecting heavier nuclei, using the LIS computed by Ref. Bos17 for rigidities from 2 MeV to 100 TeV. The cross section is given by Eq. 9, where for protons and helium form factor $G_{A}(Q^{2})$ is

[TABLE]

For a vector current, as considered by Ref. Bri18 , $\Lambda_{p}\simeq$ 770 MeV and $\Lambda_{He}\simeq$ 410 MeV Ang04 . Thus we compute the boosted DM distribution by summing over Eq. (4) hydrogen and helium with the cross section parametrized as above, and inserting the result into Eq. (5). Depending on the energy, the contribution of helium to the upscattered DM flux ranges from negligible compared to the proton contribution, to roughly comparable with it; see Fig. 1 of Ref. Bri18 .

IV.3 New Limits from Daya Bay Data

The Daya Bay Reactor Neutrino Experiment consists of eight antineutrino detectors (ADs) divided among three experimental halls (EHs) inside a mountain, with vertical overburdens ranging from 250 to 860 meters water equivalent (m.w.e.) An12 ; Cao16 . The shallowest of these, EH1, is shallower than XENON1T, Borexino, and KamLAND by about a factor of 10, and deeper than MiniBooNE by around the same factor. It is located near the side of the mountain’s base (see Fig. 1 of Ref. An17 for a diagram of the mountain). We consider only the EH1 ADs, and conservatively assume that all DM arriving from below the horizon is blocked by Earth. Based on the 3-d map, we conservatively assume that the mountain blocks an additional 1/4 of the above-horizon flux, and that the rest is attenuated by 250 m.w.e. of Earth. We make no attempt to reduce the background by a similar angular factor (for Daya Bay or any other experiment), and instead compare this DM flux to the total Daya Bay event spectrum as described below

We use our numerical code described above to model attenuation by Earth’s crust. Because of the detector depth, attenuation in the atmosphere is negligible for Daya Bay and KamLAND, so we set our code to ignore attenuation in the atmosphere.

To compare the DM-induced recoil spectrum to the reported data, we write the recoil spectrum in terms of electron equivalent energy $T_{e}$, given by Refs. Bea02 ; Das11 as

[TABLE]

where $T_{p}$ is the proton recoil energy and $k_{B}$ is the material-dependent Birks’ constant Bir51 . For Daya Bay, we use $k_{B}=0.0096$ cm/MeV for linear alkylbenzene von13 . For KamLAND, we use $k_{B}=0.015$ cm/MeV, as reported by Ref. Yos10 . For PROSPECT, we assume $k_{B}=0.0111$ cm/MeV, following the simpler of the two Birks’ models in Ref. Nor17 .

Ref. An12 shows the event spectrum of the EH1 ADs after muon veto cuts have been applied, taken over three months, in the energy range from 0.3 to 100 MeVee. Here MeVee is short for MeV electron equivalent, meaning the kinetic energy reconstructed from the observed scintillation signal assuming the particle producing it was an electron. We use only single-event data, meaning that no cuts have been applied by looking for the subsequent neutron capture that would be seen in reactor antineutrino events. For Daya Bay (and KamLAND and PROSPECT, discussed below), we treat the signal detection efficiency as unity during the experiment’s specified effective exposure, as the efficiency during livetime is very high. In Fig. 4, we show the Daya Bay spectrum along with the expected spectrum of DM-induced proton recoils for a DM mass of 1 MeV and a range of cross sections.

We consider a particular DM mass and cross section to be ruled out if the predicted DM-induced recoil spectrum is higher than the total measured data at the 90% confidence level at any energy. That is to say, if the predicted DM-induced flux is greater than the measured data at any energy and has less than a 10% probability to fluctuate down to a level equal to or below the measured data, that mass-cross section pair is excluded. Note that this is different than the approach used by Ref. Bri18 . That paper compared the total event rate and DM-induced recoil above 35 MeV proton kinetic energy, rather than the bin-by-bin comparison we make. For cross sections where attenuation is negligible, our method is more sensitive than comparing to the total Daya Bay event rate, because the shapes of the measured spectra and DM-induced spectra differ significantly.

Figure 4 shows Daya Bay data along with predicted DM-induced recoil spectra for a DM mass of 1 MeV and two different cross sections, at the floor and ceiling of our exclusion region. The resulting exclusion region is shown, along with the exclusion region from the KamLAND data, in Fig. 6. This exclusion region almost completely closes the gap between the MiniBooNE and XENON1T limits.

IV.4 New Limits from KamLAND Data

KamLAND is located approximately 1 km underground (2700 m.w.e.), comparable to both XENON1T and Borexino. The mountain above KamLAND is wider than it is high, with the main access tunnel running 1.7 km to one side of the mountain and an additional train tunnel running 3 km from KamLAND to another side of the mountain. For this reason, and because KamLAND is surrounded by other mountains, we conservatively consider only DM coming from at least 15$\degree$ above horizontal, and assume that all of the DM we do consider passes through 2 km of rock. As for Daya Bay, we neglect attenuation in the atmosphere.

For KamLAND, we use data from a 8B solar-neutrino search Abe11 , which reports events from 5.5–20 MeVee of visible energy. The signal for a solar-neutrino search is elastic neutrino-electron scattering, unlike a reactor antineutrino search, which looks for a positron followed by a neutron capture. We consider the full spectrum, which includes contributions from other backgrounds.

Figure 5 shows the 8B data along with the DM-induced event spectra for a mass of 1 MeV. We proceed as for the Daya Bay data, considering a DM cross section to be ruled out for a given mass if the DM-induced event rate in the detector is significantly larger than the observed event rate in any energy bin. The resulting KamLAND exclusion region is shown in Fig. 6. Between KamLAND and Daya Bay, we exclude almost the entire XENON1T region through direct scattering on protons.

IV.5 Sensitivity for PROSPECT and JUNO

The PROSPECT reactor neutrino experiment Ash18 is located on Earth’s surface at Oak Ridge National Laboratory. Although its background rate is high, its minimal shielding, less than even MiniBooNE’s overburden, makes it ideal for probing large cross sections. PROSPECT is shielded by a meter of hydrogenous material, half a meter of concrete, and the atmosphere, of which the atmosphere is the dominant contribution. Because the atmospheric column density is at least 10 m.w.e. and because, unlike for neutron scattering, protons are not especially effective shielding compared to other nuclei, we neglect the hydrogenous material. We model the concrete as rock, and because it is a subdominant contribution to the shielding, uncertainties in its composition represent less than a $\sim$10% effect.

We use the reactor-off spectrum of electron-like events in the energy range 0.8–8.8 MeVee from Ref. Ash18 . In the available data, the background rate is too high to set a limit using the same analysis as above. However, the vast majority of the events represented by these data are expected to be electron, not proton, recoils. Applying pulse shape discrimination (PSD) cuts to the data could reduce the event rate by a factor of $\sim$1000 privcomm . We derive a projected sensitivity based on the available data with an assumed factor of 1000 background reduction. Similar to above, we exclude a DM mass and cross section pair if the resulting CR-induced event rate is greater than the total PSD-reduced event rate at the 90% level in any bin. The resulting sensitivity is shown in Fig. 7. This region covers the largest cross sections probed by any direct-detection experiment for $m_{\chi}\ll 1$ GeV.

Additionally, the sensitivity could be further improved with precise background modeling. We have adopted the very conservative approach of comparing the DM-induced event rate to the total event rate in a given bin, but there are known to be standard model backgrounds for PROSPECT (and other detectors). Perfect background modeling would allow us to improve sensitivity by comparing the DM-induced event rate not to the total background, but to the statistical uncertainty on the background. This sensitivity, derived as a 90$\%$ CL and conservatively assuming one week of reactor-off data, is also shown in Fig. 7. The true sensitivity, including systematic as well as statistical uncertainty in background modeling, would lie somewhere in between.

Finally, JUNO is a future reactor neutrino experiment currently under construction, at a depth of about 700 meters and with much larger volume than KamLAND Dju15 . Making the same assumption of perfect background modeling as above, we can compute the sensitivity of JUNO for a one year exposure. We consider only visible energies of 13.5 – 20 MeV, the same energy bin which set our KamLAND limit, and assume that the background rate is the same as that measured by KamLAND in Refs. Abe11 ; Gan11 , which is dominated by atmospheric neutrinos. In fact, the muon-induced neutron background should be much higher, but should only increase the total background by a factor of $\sim 2$. JUNO will also self-shield against neutrons, meaning that a slightly smaller fiducial volume could be used to largely remove this background. The derived 90$\%$ CL sensitivity region is shown in Fig. 7. With only a year of exposure, JUNO can probe the lowest cross sections ever reached by direct detection for a large range of sub-GeV masses.

IV.6 Astrophysical Uncertainties

There is some uncertainty in our results in the local DM density, which is unavoidable for any direct detection experiment. The local DM density is typically taken to be 0.2–0.4 GeV/cm3, and we use a standard value of 0.3 GeV/cm2, so we take the uncertainty on the local DM density to be ${\sim}33\%$. There is also uncertainty in the steepness of the galactic DM density profile, especially close to the galactic center. Models such as the Burkert profile Bur95 have much lower DM density near the galactic center than the NFW profile that we employ. To quantify this uncertainty, we write down an extreme model in which the density profile follows an NFW shape for $r>8$ kpc (the distance from Earth to the galactic center), but is held constant for $r<8$ kpc. We find that the flux of upscattered dark matter from this cored profile is ${\sim}25\%$ lower than for the NFW profile. Adding these uncertainties in quadrature, we find a ${\sim}40\%$ uncertainty on the DM flux due to uncertainties in the DM density profile. Because our limits scale with the square root of the dark matter flux, this would result in only a ${\sim}20\%$ uncertainty on our limits.

For CR inputs, we make strictly conservative choices. As mentioned above, there is evidence that the CR flux is higher closer to the galactic center, which would make our results somewhat stronger. The LIS from Refs. Bos17 ; Bos18a are somewhat lower than the data from Voyager measurements shown in the same papers. And we assume a conservative size for the CR halo, and considering a larger volume would improve our results.

But the uncertainties associated with the DM and CR densities are small compared to the statistical uncertainty on the background, and originate from our choice to compare the DM-induced event rate with the total background. For KamLAND, for example, there is one observed event in the considered energy range, but due to the large statistical uncertainty associated with a single event, we require four DM-induced events. While astrophysical inputs produce a $\sim 20\%$ uncertainty on our limit, the low statistics produce a factor of 2 uncertainty. For Daya Bay, the effect of treating all observed events as potentially due to DM, rather than performing any background subtraction, is similarly large compared to astrophysical uncertainties.

IV.7 Future Ways to Improve Sensitivity

It is apparent from Fig. 6 that detectors at different depths cover different regions of parameter space, but share a similar shape: the ceiling of each region is set by attenuation in the atmosphere and Earth, the floor is set by detector backgrounds, and the high-mass end of each region is set by how they join. Because detectors at different depths are sensitive to different ranges of cross section, considering additional detectors may provide useful new coverage. We have considered only two of the Daya Bay ADs in setting our exclusion regions, and considering the other (deeper) detectors could push Daya Bay’s sensitivity to lower cross sections. Probing cross sections above the PROSPECT sensitivity region would require detectors with minimal atmospheric shielding, such as balloon-, rocket-, or satellite-based experiments.

V Constraining the DM-Electron Cross Section

In this section, we apply the general treatment presented in Sec. II to DM-electron scattering. For our new results, we use lower-energy data than Ref. Ema18 , and show that doing so yields a larger signal to background ratio at a given DM mass and cross section, thus producing stronger limits.

V.1 Electron Recoil Distribution

The case of DM-electron scattering is simpler than DM-nucleon scattering because we consider only one CR species, and because there is no form-factor suppression. Additionally, we do not have to convert between nuclear recoil energy and electron recoil energy in the detector. So we obtain the DM flux directly from Eqs. (4) and (5), with the DM-electron cross section $\sigma_{\chi e}$ being the relevant cross section.

In computing the DM flux at Super-K, we consider only DM arriving from at least 15$\degree$ above horizontal, and neglect attenuation for this flux. Over the entire DM mass range we consider, the lowest cross section probed by Ref. Ema18 is orders of magnitude below the ceiling caused by attenuation in Earth. We restrict our attention to cross sections too small to be excluded by Ref. Ema18 , so attenuation is negligible. As a result, we can insert the computed DM spectrum directly into Eq. (6) to get the recoil spectrum per target electron, and multiply this by $7.5\times 10^{33}$ electrons Ema18 to get the recoil spectrum in Super-K.

We note, however, that the ceiling of our analysis could be lower if $\sigma_{\chi N}$ or the DM-proton cross section $\sigma_{\chi p}$ is closely related to $\sigma_{\chi e}$, such that increasing $\sigma_{\chi e}$ increases attenuation by nuclei as well (see, e.g., Ref. Emk17b ).

V.2 Our Improvements Over Previous Work

As seen in our Fig. 8, the DM-induced recoil spectrum is steeply falling, meaning that for a given cross section, a higher signal rate is observed at lower energy. Below, we show that considering data below 100 MeV produces significantly stronger limits than previously derived. We also derive projected limits from Hyper-Kamiokande that are stronger than the previously derived projections for Super-K and DUNE.

V.3 New Limits from Super-K Data

In the previous analysis, Ref. Ema18 considered Super-K data above 100 MeV in order to stay above the atmospheric neutrino background. We instead consider energy as low as 10 MeV, a range where the atmospheric neutrino background falls with decreasing energy while the predicted DM spectrum rises sharply, in order to obtain a higher signal-to-background ratio. We consider the spectrum of electron-recoil events from a diffuse supernova neutrino background (DSNB) search performed using Super-K data from the SK-I exposure, from April 1996 to July 2001 Bay11 . During this time period, 239 electron-recoil events are reported in the energy range from 16–88 MeV visible energy (related to our work, these data have also been used to set constraints on dark radiation Cui17 ). Note that because we consider electron-recoil events, there is no need to include quenching, as was done in the previous section to convert nuclear recoil energy to electron-equivalent energy. We follow the same approach as for DM-nucleon scattering, considering a DM mass and cross section ruled out if the predicted DM-induced event rate is higher than the measured event rate at the 90% level at any energy.

Some fraction of DM-induced events may be lost to analysis cuts. We obtain the signal detection efficiency as a function of recoil energy by interpolating Fig. 10 of Ref. Bay11 , and multiply our computed event spectrum by the efficiency in order to get the spectrum of DM-induced events passing analysis cuts. As the data only extend down to 16 MeV, we take the signal efficiency to be zero below this energy.

Figure 8 shows the Super-K data along with the (efficiency-corrected) spectrum of CR-induced DM events for a mass of 1 MeV and several values of the cross section. The increasingly large cross sections, as labeled, correspond to our limit, the halfway point between our limit and the previous limit from Ref. Ema18 , and the limit from Ref. Ema18 itself. The resulting limit is shown in Fig. 9, along with the limits from Ref. Ema18 , as well as astrophysical, cosmological, direct detection and solar reflection limits.

Our limit is tighter than the limit set by Ref. Ema18 by a factor of at least $\sim 4$, depending on the mass. It is comparable to the Super-K sensitivity curve derived in Ref Ema18 , and close to their DUNE sensitivity curve. This means that we have already excluded most of the parameter space which was previously thought to be reachable only with careful background subtraction, restriction to events coming from the galactic center, and possibly several years of DUNE data.

V.4 Sensitivity for Hyper-Kamiokande

Hyper-Kamiokande is a proposed neutrino experiment that will consist of two large water Cherenkov detectors, with a total fiducial volume 16.8 times that of Super-Kamiokande Abe16 . Assuming that the observed background will be similar to that of Super-Kamiokande, we can compute the 90$\%$ CL sensitivity in the same way as we did for PROSPECT and JUNO. We assume an exposure 16.8 times larger than the Super-Kamiokande exposure used to set our limit. This projection, shown in Fig. 9, is stronger than previous projections for Super-K and DUNE and at least an order of magnitude stronger than previous constraints.

V.5 Future Ways to Improve Sensitivity

Unlike in the case of DM-nucleon scattering, we do not consider detectors at different depths for DM-electron scattering, nor would doing so improve on existing limits. Super-K is already very deep, and higher cross sections are already covered: though it is not shown in Fig. 9, Ref. Ema18 also derived a constraint based on MiniBooNE data, which covers cross sections above the ceiling they computed for the Super-K region. Still larger cross sections are constrained by considering CR downscattering Cap19 and cosmological observations Ali15 . We refer the reader to Ref. Ema18 for a plot of the higher cross section constraints.

Additionally, in a detector that has directional sensitivity, as Super-K does, it is possible to search for events coming only from the direction of the galactic center. Ref. Car15 has studied the morphology of a hypothetical signal coming from CR-DM interactions, and the signal distribution they compute could be useful for such a search. These possibilities have both been discussed in Ref. Ema18 .

A coming improvement to the Super-K detector is the addition of water-soluble gadolinium salt, which will allow tagging of antineutrino events, greatly reducing backgrounds Bea03 ; Lab18 . Another is the further reduction of spallation backgrounds, which are the dominant background in Super-K from 6–18 MeV, based on the cuts proposed by Refs. Li14 ; Li15a ; Li15b . Background reduction at low energy, where the predicted DM-induced event rate is largest, will improve sensitivity.

VI Directions for Future Work

Here we mention new directions for testing low-mass DM, beyond the discussions above about improving the sensitivity of this paper’s results.

We first note that the energy dependence of the cross sections may be non-constant, contrary to the energy-independent case we assume, which would change the exclusion regions. If the cross sections increase with energy, as is likely, this would improve the sensitivity, although it may also lower the ceiling of the corresponding region. See Ref. Dent19 ; Bon19 for initial explorations. More generally, it will be interesting to develop concrete models for these light DM scenarios; that could also provide connections to constraints from fixed-target experiments Dia13 ; Bat14 ; Iza14 ; Kah14 ; Sop14 ; Col15 ; deN16 ; Ban17 ; Agu18 , as done in Ref. Bon19 .

We have only considered elastic scattering on protons and electrons, but other target nuclei and additional signals could also be considered. Carbon recoils in liquid scintillator would be heavily quenched, but would benefit from coherent enhancement of the cross section, possibly allowing neutrino detectors to probe lower cross sections. Carbon nuclei in a scintillator-based detector could also be excited by a collision with high-energy DM, and the de-excitation would produce a characteristic 15.11 MeV gamma-ray signal. Additional signals such as bremsstrahlung photons, considered by Refs. Kou16 ; McC17 ; Ake18 ; Bel19 ; Liu19 , and the Migdal effect, as discussed by Refs. Ibe17 ; Dol17 ; Ake18 ; Arm19 ; Bel19 ; Liu19 , could be seen in low-threshold detectors.

A DM particle with sufficient energy could also perhaps expel a neutron from a carbon nucleus. The resulting neutron capture and subsequent decay of the unstable 11C nucleus would be a distinctive delayed-coincidence signal Gal04 . An analogous delayed coincidence signal could be seen in Super-K, with a neutron expelled from an oxygen nucleus, leading to both a neutron capture signal and the subsequent gamma-ray deexcitation of an excited 15O nucleus Lan95 . For even higher-energy DM, collisions with nuclei could produce pions, and the subsequent pion decays could be additional signals to search for. Future work could also consider inelastic and quasi-elastic contributions to the attenuation, although we expect these effects to be subdominant in the energy range that drives our limits.

VII Conclusions

We have expanded and improved upon previous studies of cosmic ray-upscattered dark matter, in which dark matter that is nominally too light to be detected is accelerated by cosmic rays to detectable energies. We have developed the first numerical propagation code to model attenuation of light, relativistic dark matter, resulting in more accurate ceilings which allow us to probe higher cross sections than the previously used ballistic approximation. We have excluded parameter space for dark matter-proton scattering never before probed by direct detection, and identified detectors that can expand this coverage even further. We have also substantially increased the parameter space probed directly by dark matter-proton scattering, making our results complementary to xenon-based limits. Finally, we have improved the sensitivity of Super-K to dark matter-electron scattering, for some masses superseding even the future projections derived in previous work, and derived the sensitivity of the future Hyper-Kamiokande detectors.

Though our limits are based on the astrophysics of cosmic rays and the galaxy’s dark matter halo, our results are fairly insensitive to astrophysical uncertainties. Unlike traditional direct-detection searches, the observed event rate scales with the square of the cross section, so that our limits depend only on the square root of the cosmic ray flux or dark matter density. Uncertainties in the interstellar cosmic ray spectra, cosmic ray halo size, and dark matter profile are suppressed by this square root and thus largely negligible.

Future work could lead to significant improvement on our results, extending sensitivity to both higher and lower cross sections. Significant background reduction and precise background modeling would lead to strengthened constraints on both dark matter-proton and dark matter-electron scattering. With such developments, as well as consideration of specific dark matter models, cosmic rays will only become more powerful as a probe of sub-GeV dark matter, constraining the parameter space between collider and cosmological limits.

Acknowledgments

We are grateful to the anonymous referee for insightful comments and suggestions, which substantially improved the paper. We are grateful for helpful discussions with Nathaniel Bowden, Jun Cao, Eric Carlson, Timon Emken, Rouven Essig, Glennys Farrar, Karsten Heeger, Thomas Langford, Bryce Littlejohn, M. Shafi Mahdawi, Pieter Mumm, Kenny Ng, Annika Peter, Pranava Teja Surukuchi, Lindley Winslow, and Xingchen Xu, and especially Torsten Bringmann, Yohei Ema, Maxim Pospelov and Filippo Sala.

CVC and JFB are supported by NSF grant PHY-1714479.