Neutrino experiments probe hadrophilic light dark matter

We use Super-K data to place new strong limits on interactions of sub-GeV Dark Matter (DM) with nuclei, that rely on the DM flux inevitably induced by cosmic-ray upscatterings. We derive analogous sensitivities at Hyper-K and DUNE and compare them with others, e.g. at JUNO. Using simplified models, we find that our proposal tests genuinely new parameter space, allowed both by theoretical consistency and by other direct detection experiments, cosmology, meson decays and our recast of monojet. Our results thus motivate and shape a new physics case for any large volume detector sensitive to nuclear recoils.


Introduction
Dark Matter (DM) is perhaps the strongest evidence for physics beyond the Standard Model (BSM), and thus motivates an intense experimental effort to determine its still elusive properties. Direct detection [1] (DD) experiments constitute one of the leading avenues to probe the non-gravitational interactions of DM particles, and have so far mostly focused on a DM mass range roughly between a GeV and tens of TeV. These searches are now in a mature stage, and no DM signal has been detected so far [2][3][4]. In addition, the LHC has not found clear signs of the BSM physics that, independently of DM, motivates that mass range [5]. Therefore, in the past few years, the community has investigated more in depth other DM mass regimes, and in particular the one of sub-GeV DM (see e.g. the report [6]).
Direct detection experiments rely on the scattering of DM with a target particle on the Earth, whose effects, like the recoil energy of the target, can be detected. The average velocity of DM particles in the solar neighbourhood, v ≈ 10 −3 , implies that for any DM mass there exists an upper limit on the final kinetic energy of the target T , k max = 2v 2 µ 2 DMT /m T × (1 + O(v 2 )), where m T is the target mass and µ DMT is the reduced mass of the DM and the target. Therefore experiments with nuclear targets sensitive to recoil energies keV, like [2][3][4], quickly lose sensitivity for DM masses lighter than about a GeV. In order to test sub-GeV masses, recent efforts have focused on different targets and/or detection strategies, see [6] for a review.
Recently, refs [7,8] realised that another possibility to detect light DM is given by the higher-speed DM component upscattered by cosmic-rays (CRs). The existence of such a component is unavoidable as long as one assumes some DM interactions with the Standard Model (SM), which is nothing more than what is assumed in any direct detection experiment. The larger DM kinetic energy then implies much larger recoil energies at the Earth, so that unexplored regions of sub-GeV DM can be probed both in 'standard' experiments [7] like LUX, PandaX and XENON1T [2][3][4], and with completely different technologies like those available at large neutrino experiments [8]. 1 Focusing now on DM interacting mainly with nucleons, we improve over previous studies [7,[12][13][14] in two respects: i) We derive new limits that are stronger than all existing ones [7,12], by use of Super-K data of protons detected through their Cherenkov light [15,16], and determine analogous sensitivities at Hyper-K and DUNE. The effectiveness of our strategy is due to the volumes and detected nucleon energies of these experiments, both larger than those of any other experiment previously used; ii) We find regions of parameter space where our new limits and sensitivities are stronger not only than all other direct detection techniques, but also than all other constraints we know of (collider, cosmology, etc). We do so in two benchmark 'simplified' DM models, where we add a new particle to mediate the SM-DM interactions below the electroweak scale. Moreover we provide explicit completions of our simplified models above the electroweak scale, that are consistent with other constraints on the new model ingredients, and where our limits and sensitivities still win over all other existing ones. This is a first in the study of detection of CR-upscattered DM, and it allows to define a solid physics case to pursue this novel detection technique at large volume experiments, like neutrino detectors.
To appreciate the novelty of the second item, as well as to better introduce our work, it is useful to summarise here recent related papers. While refs. [7,12] derived new limits on DM interactions with nuclei, they only considered the case of constant cross section. As we stressed already in [8], it is important to include its energy dependence when comparing this detection technique with other ones, because they rely on different energy regimes. Refs [13,14] pursued this direction, but [13] did not compare the constraints with other ones, nor worried about finding an existence proof for their simplified models. Ref [14] did, but found that, in the specific model considered there, the parameter space probed by detection of CR-upscattered DM with Xenon-1T was already entirely excluded by other constraints.
The rest of this paper is organised as follows: in section 2 we present the simplified models of our interest and the resulting DM cross sections, in section 3 we derive new limits on light DM with Super-K data and new related sensitivities at Hyper-K and DUNE, in section 4 we discuss other constraints, in section 5 we prove that our new limits and sensitivities test theoretically-consistent parameter space. We conclude and discuss future directions in Section 6.

Dark Matter Interactions
We consider for definiteness DM as a Dirac fermion χ, singlet under the SM gauge group, with mass m χ . We add to the Standard Model either a new scalar φ with mass m φ , or a new pseudoscalar a with mass m a , to act as mediators of DM-SM interactions via where the index q runs over the SM up, down and strange quarks. For our purpose, we have to convert this interaction with quarks to that with nucleons and eventually with nuclei. In the following, we explain our procedure of the conversion for the scalar and pseudoscalar mediator cases separately. Eqs. (1) and (2) give rise to different regimes of low-energy scattering of DM with nuclei N . At tree level, the leading non-relativistic operators in momentum space, induced by L φ and L a , are respectively (see e.g. [17]) where s N is the nucleus spin, s χ the DM one, q is the momentum exchanged in the χ − N scattering. Therefore 'standard' direct detection techniques test L φ via spin-independent searches, while they are insensitive to L a at tree level. As we will see in Sec. 4, 'standard' DD techniques are sensitive to L a at the one-loop level. We anticipate that, instead, our detection technique is sensitive to both L φ and L a , because CRs have large kinetic energy and the scattering cross section does not suffer from suppression by velocity.

Scalar mediator
We write the matrix elements of the quark bilinearqq, between nucleons N = p, n with initial momentum p i and final momentum p f , as where u N is the spinor wavefunction of the nucleon, m N and m q are the nucleon and quark at the leading order in the expansion for small |t| after factorizing the spinor wavefunction, which we extract from [18] as where we have averaged over proton and neutron. Note that the analysis of [18] is based on experimental data and is confirmed by [19,20], while lattice QCD studies in general predict slightly smaller values of f N u and f N d (see, e.g., [21,22]). The scalar nucleon form factor G incorporates the finite size effect of nucleon, which reads [23] G(q 2 ) = 1 We then compute the differential scattering cross sections of DM with nuclei A. For the scalar mediator case, nucleons inside nuclei A coherently contribute to the differential cross section. In our computation, we treat nuclei A as a point particle when we evaluate the kinematical variables, and put the nuclei form factor to incorporate the finite size effect of nuclei A. The differential cross section is then given by 2 where m χ is the DM mass, K f is the kinetic energy of the final state particle of interest (e.g. f = χ when we will compute the DM flux from CR upscatterings, f = p when we will compute the proton recoil spectrum in detectors), K max is its maximal value which will be computed case by case, s is the usual Mandelstam variable, n A is the number of nucleons N in the nucleus A of interest (in this paper we will consider n A = 1, 4 respectively for A = 1 H and A = 4 He). As we mentioned, we evaluate s and t in this equation by the kinematics of nuclei A as a whole. When A is Hydrogen, A = 1 H, this is simply a proton and hence the form factor is given by 2 Another way of including the distribution of nucleons inside nuclei is discussed in [24], which also includes an incoherent contribution. We have checked that the differential cross section computed following [24] gives a similar constraint on the couplings. Since other uncertainties such as astrophysical ones have larger impacts on the final results, eq. (8) is enough for our purpose.
When A is Helium, A = 4 He, the form factor is given by [23] F He (q 2 ) = 1 The coupling g N φ depends on how φ couples to quarks. In the following, we consider the isoscalar coupling g uφ = g dφ = g φ , and assume that couplings to the other quarks vanish. In this case, g N φ is given as

Pseudoscalar mediator
We write the matrix elements of the quark bilinear iqγ 5 q, between nucleons N = p, n with initial momentum p i and final momentum p f , as where h q is defined as at the leading order in the expansion for small |t| after factorizing the spinor wavefunction, that we extract from [25] as for the isovector and isoscalar combinations, respectively. The coupling h s is also computed in [25], but we do not show it here since we take g sa = 0 in our computation below. The pseudoscalar nucleon form factor G q a is evaluated by the partially conserved axial current (PCAC) relation as We infer the axial vector form factor G A from the the isovector results in [25] as 3 This value of Λ a is larger than the world average [26,27], but is motivated by the recent MiniBooNE experiment [28] and the lattice results [25,29]. The induced pseudoscalar form factor G q P is inferred from the pole dominance ansatz as where the parameters are extracted from [25] as for the isovector combination, where m π is the pion mass, and for the isoscalar combination. Finally G G denotes the anomaly form factor which originates from the GG-term in the PCAC relation. Since the isoscalar current is anomalous while the isovector current is not, u+d = 1 and u−d = 0. G G is computed, e.g., in [30]. We have checked that G G is of little significance to our final result, and hence ignore it in the following for simplicity. Since all references we are aware of take the isospin symmetric limit, we assume that the up and down quark masses are equal and are the arithmetic mean of their experimental values when we use the above formulas. Since Helium does not have a net spin, we consider only the scattering with proton for the pseudoscalar mediator case in this paper. The differential scattering cross section with proton is given by In the following, we consider the isoscalar coupling g ua = g da = g a , and assume that couplings to the other quarks vanish. In this case, g N a and F a are given by where we assume the degenerate up and down quark masses as we mentioned above.

New limits and sensitivities
In this section we explain our method to compute the event rates at the terrestrial detectors, and show our new limits on the parameters of our simplified DM models.

DM upscattered by Cosmic Rays
Our galaxy is filled with high-energy CRs. Once these CRs have interactions with DM, they can upscatter DM such that DM acquires a large kinetic energy compared to its averaged value inside the galaxy. In the case of the hadrophilic DM, the CR protons and Heliums are the dominant source of such a high-energy component of the DM flux. The DM flux Φ χ upscattered by the CRs is given by where Ω is the solid angle, "l.o.s." stands for the line-of-sight integral, Φ A is the CR flux of the nuclei A in the galaxy, dσ/dK f is the differential cross section between A and χ which we explained in the previous section, and ρ χ is the DM energy density in the galaxy. We include both A = p, He for the scalar mediator case, and include only A = p for the pseudoscalar mediator case. 4 In this paper, we assume for simplicity that the CR flux is homogeneous inside a leaky cylinder centered on the galactic center, and vanishes outside the cylinder. We can then simplify Eq. (22) as where b and l are the galactic latitude and longitude. The J-factor J(b, l) contains the information of the DM-and CR-distributions within the galaxy, and is given by which depends on b and l. In this formula, the integration is bounded by the leaky cylinder, whose radius R and half width h are fixed as R = 10 kpc and h = 1 kpc. We take the DM density profile as the NFW profile [31], where r is the distance from the galactic center. We take r = 8.5 kpc, r c = 20 kpc and, conservatively, ρ χ (r = r ) = 0.42 GeV [32][33][34]. Note that our constraints only mildly depend on the choice of the DM density profile, since the J-factor is linear in ρ χ . We input the CR flux as follows. For K A ≤ 100 MeV, we directly use Voyager's observation with interpolation [35]. The lower thresholds are taken as 22 MeV for proton and 10 MeV for Helium. We assume no CR flux below these energies simply because there is no observation. This assumption does not affect the DM event rates at the neutrino experiments, but may affect other DM direct detection experiments which have lower threshold recoil energy. For 100 MeV < K A < 50 GeV, we use the theoretical computation in [36]. Finally for K A > 50 GeV, we use the fitting formula given in [37]. Note that the fitting formula in [37] underestimates the CR flux compared to Voyager's observation while [36] does not.
In Fig. 1, we show DM fluxes at the Earth's surface for some benchmark values of the parameters. We take m χ = 10 MeV, m φ = 1 GeV and g χφ g uφ = g χφ g dφ = 0.1 for the scalar mediator case, and m χ = 10 MeV, m a = 1 GeV and g χa g ua = g χa g da = 0.1 for the pseudoscalar mediator case, respectively. As a comparison, we also show the standard halo DM flux in the right panel, which is given by where c is the speed of light and v(K χ ) = K χ /2m χ . The velocity distribution function f is given by  where Θ is the Heaviside theta function and Erf is the error function. We took σ = 163 km/sec (corresponding to the peak velocity 230 km/sec) and v esc = 600 km/sec in Fig. 1. As one can see from the figure, although the overall normalization is much smaller compared to the standard halo DM flux, the CR-upscatterred DM flux is extended to far higher energy region. This feature makes it possible to probe this component of the DM by the neutrino experiments that have larger energy thresholds than the standard DM direct detection experiments.

Earth attenuation and events at detectors
The neutrino detectors are located deep inside the Earth in order to reduce backgrounds. As a result, DM can be scattered within the Earth before arriving at the neutrino detectors. We now explain how we implemented this Earth attenuation in our computation. We approximate the energy loss at each scattering as its averaged value for simplicity. We also ignore the change of the direction of the DM velocity at each scattering. This treatment typically gives a slightly weaker bound on the DM cross section than a more refined treatment [38]. With these approximations, the DM kinetic energyK χ at a depth z obeys the following differential equation: where n T is the number density of a target particle T inside the Earth. The initial condition is given byK χ (z = 0) = K χ where the latter is the DM kinetic energy at the Earth's surface.
Since the typical energy scale of the DM events at the neutrino experiments is higher than the nuclear energy scale, we ignore the nuclear structure of the target particles inside the Earth. Then the target particles are protons and neutrons, and we assume their ratio to be 1 : 1. The mass density is taken as ρ p+n = 2.7 g/cm 3 [38], from which the number densities n p and n n are derived. We ignore the form factors of protons and neutrons, which results in a conservative constraint. The DM fluxΦ χ at a depth z is then related to its value at the Earth's surface (computed by Eq. (22)) bȳ which follows from the flux conservation.
As a comparison, we also show the constraints from Xenon1T for the scalar mediator case in the figures below. Since the Xenon1T detector is located deep inside the Earth, we also take into account the Earth attenuation in this case. Contrary to the neutrino experiment case, we treat the nuclei as a whole because of the lower threshold energy of Xenon1T. For simplicity, we assume that all the nuclei inside the Earth are with A = 16, corresponding to oxygen which is the dominant component. We have checked that the results do not change if we instead take A = 12, 20 or 24. We ignored the atomic form factor in this computation to be conservative.
Once we have implemented the Earth attenuation, we are ready to compute the DM event rates at the neutrino detectors. The direction-averaged proton event spectrum at a detector is given by where dN p /dt is the number of the proton events per unit time, and N T is the total number of protons inside the detector. Note that the depth from the Earth's surface z in this formula depends on the solid angle Ω. The depth z is the smallest when DM comes from right above the detector, while it is the largest when DM comes from the opposite side of the Earth. In the case of DUNE, we can also use neutrons as our signal, and hence we use the same formula but replacing p → n to compute the neutron event spectrum.
In Fig. 2, we show the normalized recoil spectra (1/N p )dN p /d log K p without the Earth attenuation. In this figure, we fix the DM and mediator masses as 10 MeV and 1 GeV, respectively. For comparison, the recoil event spectrum for the constant cross section case considered in [7,12] is also shown. As one can see, the spectra for the scalar and pseudoscalar mediator cases are extended to high energy region compared to the constant cross section m φ = 1 GeV, and g χφ g uφ = g χφ g dφ = 0.1. Right: the pseudoscalar mediator case with m χ = 10 MeV, m a = 1 GeV, and g χa g ua = g χa g da = 0.1. Here we take the number of the target particle as that of Super-K, i.e., N T = 7.5 × 10 33 .
case, especially when m χ is smaller. It is because the cross section grows with energy for the scalar and pseudoscalar cases, as long as the energy transfer is smaller than the mass of the mediator. Because of this feature, an experiment with a larger energy threshold such as the neutrino experiments tend to be more sensitive to the cases with the massive mediators.
Comments on the Earth attenuation. In order to visualise the Earth attenuation effects, we plot in Fig. 3 the spectra per unit time d 2 N p /dtdK p at different values of the depth z. The left panel is the scalar mediator case with m χ = 10 MeV, m φ = 1 GeV, and g χφ g uφ = g χφ g dφ = 0.1, and the right panel is the pseudoscalar mediator case with m χ = 10 MeV, m a = 1 GeV, and g χa g ua = g χa g da = 0.1, respectively. It is clear that the Earth attenuation depletes the DM kinetic energy and hence the peak of the recoil spectrum is shifted toward lower energies as z becomes larger. If z, or equivalently the coupling, 5 is too large, the Earth attenuation is so effective that no DM flux reaches the detector. As a result, our constraints on the coupling are bounded both from below and from above. The left-and right-hand plots in Fig. 3 also display a qualitative difference between the scalar and pseudoscalar mediator cases, namely that the peak shifted to lower energies is more pronounced in the latter case than in the former. This feature is due to the fact that the pseudoscalar-mediator cross section, cf. Eq. (20), goes to zero in the limit of no exchanged energy, so that below a certain K p Earth attenuation becomes irrelevant and the DM flux accumulates there. The scalar-mediator cross section, cf. Eq. (8), instead goes to a constant in the limit of no exchanged energy, so that the peak feature at low K p is less pronounced. We have explicitly checked that, as expected from our understanding, the peak feature actually disappears in the scalar-mediator case for a value m χ = 3 GeV (thus different than those in Fig. 3), for which the constant limit in the cross section kicks in for larger values of the exchanged energy.
Another very interesting point, concerning Earth attenuation, is that it can make some constraints stronger, rather than weaker. Let us imagine the situation that a detector is 5 Note that z appears only in the combinations g 2 χφ g 2 uφ z or g 2 χa g 2 ua z.
sensitive to the energy region K min < K p < K max , which satisfies K max < K peak , where K peak is the peak energy of the recoil spectrum without the Earth attenuation. In this case, the Earth attenuation can help to probe the DM by shifting the peak energy to lower such that it falls into the range K min < K peak < K max , before it completely depletes the DM flux. Since the peak energy is O(100 MeV-1 GeV) without the Earth attenuation, this effect makes the constraints stronger for experiments with lower threshold energy such standard DD experiments. On the other hand, this effect does not strengthen the constraints of of experiments with larger threshold energies like Super-K and Hyper-K.

New limits and sensitivities
We now derive our new limits on the DM models by Super-K, and also future sensitivities of Hyper-K and DUNE. We derive the limits and sensitivities in the following way. In deriving them, we use the DM flux averaged over the solid angle. Super-K measured N d = 16 downgoing and N u = 13 upgoing events for 2287.3 days in the energy range 0.485 < K p /GeV < 3.17 (corresponding to the momentum range 1.07 < p p /GeV < 4) [16]. Since our signal is expected to be downgoing in the range of the cross section of our interest, we derive our limits by requiring that where SK = 0.55 is the efficiency [16]. In our computation, we assume that all protons inside the detector contribute to the events because of the high energy threshold, and hence take N SK T = 7.5 × 10 33 , corresponding to an effective volume of 22.5 ktons. The depth of the Super-K detector from the Earth's surface is z SK = 1 km.
Hyper-K is supposed to start taking data in 2027 with its fiducial volume 187 ktons [39]. In order to estimate the upper limit, we assume 5 years of data taking with the same number of background event rate as Super-K. Therefore, we require that where we assume the same signal energy range as Super-K and take HK = 0.55. We take the number of the target particle as N HK T = 6.2 × 10 34 and the depth as z HK = 0.65 km. The DUNE detector will consist of 4 modules, each with 10 ktons of fiducial volume. The first two are expected to be completed by 2024, and be operational by 2026 [40]. Assuming the other two will be installed and operational by 2027, we show sensitivities for 5 years of data taking. Since DUNE relies on the scintillation light, not the Cherenkov light, the threshold energy can be lower than Super-K and Hyper-K. We take it as K p > 50 MeV [41]. We take the upper energy threshold as 10 GeV, although the result is insensitive to this choice as long as it is much larger than 50 MeV. In this energy range, DM dominantly scatters with individual nucleons [42], and hence we take the number of the target particles as N DN T = 2.4 × 10 34 where we summed all protons and neutrons. Concerning background events, we are not aware of any detailed study that is straightforwardly applicable to our case. Therefore, we simply draw the lines by requiring N DN p+n (0.05 < K p /GeV < 10) = 30, The use of XENON1T data have been proposed in [7], the use of KamLAND and future JUNO data in [12]. We have recasted them on the models considered here, with a procedure consistent with that used for our new results, see Sec. 3. Shaded gray areas are, from left to right, limits from BBN, our recast of monojet searches, 'standard' DD, see Sec. 4. The BBN limits are shaded differently than the other ones because, for 1 MeV m χ 10 MeV, one could alleviate them by switching on a coupling to neutrinos, in a way that does not affect any of the other phenomenology presented here, see text for more details. Upper: scalar mediator. Lower: pseudoscalar mediator. In the latter case, limits from meson decay depend further on the coupling to the strange quark. Since this negligibly affects the rest of the phenomenology, we shade only the area that is excluded independently of its value, and display as a dashed lines the ones that are excluded for g sa = g ua (η ) and g sa = −2g ua (η).
in the figures below, where we assume DN = 1, and leave a detailed study on background events as a future work. Finally, the depth of the detector is taken as z DN = 1.5 km [40].
Our results are displayed in Figure 4. There we also show limits on CR-upscattered DM from XENON1T data, computed following [7] 6 , as well as other experimental limits on the models, discussed below in Sec. 4. Their comparison with our new constraints and sensitivities shows that our method can probe a parameter region which is not yet probed by any other experiments. We also show limits from KamLAND data [43] and sensitivities of JUNO [44], computed following [12]. KamLAND observed one event in the bin [13.5 MeVee, 20 MeVee] within 123 kton-day where MeVee is an abbreviation for MeV electron equivalent, and hence we require the DM events smaller than three in this bin. We assume that JUNO has the same background event rate as KamLAND and compute the sensitivity of the bin [13.5 MeVee, 20 MeVee] for the exposure of 5 years with a 20 kton detector. Of course, limits from XENON1T and KamLAND and sensitivities from JUNO have been derived with the same astrophysical inputs used in the rest of this work.
Discussion. The limits and sensitivities that rely on CR-upscattered DM extend to DM masses much smaller than those shown in 4, we do not display that region just because it is anyway excluded by BBN. They also suffer some astrophysical uncertainties, like the precise value of the local DM density and the region of the Milky Way where our input CR spectra are valid. The impact of these uncertainties on limits and sensitivities on CR-upscattered DM, however, is limited: the signal events scale like the square of the cross section (because one hit is necessary to upscatter the DM, and one to detect it), and so as the fourth power of the coupling products on the vertical axis of Fig. 4. The signal events also scale linearly in the astrophysical quantities mentioned above, so that an error as large as a factor of 2 in, e.g., the DM local density, would translate into a shift of less than 20% in the contours shown in Fig. 4. Coming now to comments specific to our new limits and sensitivities, one sees that in the pseudoscalar case DUNE is expected to give marginal improvements with respect to Super-K and Hyper-K, for sub-GeV DM masses. Instead, in the scalar case and over the entire m χ range that we plot, DUNE is expected to improve substantially over Super-K and Hyper-K. This can be explained by noticing that the event recoil spectrum, in the scalar case, stays sizable at lower energies, so that the lower threshold energy of DUNE constitutes an advantage. In the pseudoscalar case instead the recoil spectrum drops significantly at lower energies (see e.g. Fig. 2), explaining why one does not see an analogous improvement at DUNE. Limits from KamLAND and sensitivities from JUNO can be understood in the same way. This property constitutes a concrete example of the necessity to use explicit models, with the resulting relativistic and momentum-dependent cross sections, when comparing the reach of different experiments to CR-upscattered DM. The lower energy threshold of DUNE also explains why its sensitivity improves a lot for DM masses above a GeV, with respect to Super-K and Hyper-K.
Concerning finally XENON1T, it is insensitive to the pseudoscalar mediator case because of its low detection energy. The cross section grows with energy for the massive mediator cases, hence as visible in Fig. 4 neutrino experiments are more sensitive to CR-upscattered DM than standard direct detection experiments, because they have larger detection energies.

Limits and sensitivities on cross sections
In Figure 5 we translate our results, for the scalar mediator case, to limits and sensitivities on the spin-independent non-relativistic cross section for DM scattering with protons where µ χp = m χ m p /(m χ + m p ), see Sec. 4.3 for more details. First, this is a quantity that could be more familiar to some readers. Second, this now allows us to comment on why the use of simplified models is, for our new constraints, more appropriate than the one of an EFT.
Our new strategy probes parameter space allowed by other limits, like monojet at the LHC and standard DD, only for mediator masses around a GeV, while only DUNE will probe new parameter space for mediator masses of 3 GeV. When the mediator mass is around a GeV, the energy dependence in the denominators of the cross sections (see e.g. Eq. (8)) cannot be neglected as it is in an EFT approach. Indeed, if it could be neglected, then our new limits in the left-and right-hand plots of Fig. 5 would coincide, while they are different by a factor as large as ∼ 2 in the case of Super-K and Hyper-K. They are still different, but slightly less, in the case of DUNE, as can be understood by the fact that the energies relevant for DUNE are smaller than those relevant for Super-K and Hyper-K. The even lower energies relevant for KamLAND, JUNO and especially XENON1T are, then, the reason why those lines almost does not change by changing the mediator mass from 1 to 3 GeV. To summarise, our new limits and sensitivities really test new parameter space in a region where the use of an EFT is not fully justified, so that they should be accompanied by the value of the mass of the mediator of the interactions between DM and the SM. We conclude this section with a comment on the use of our limits and sensitivities, and more in general on all those coming from CR-upscattered DM. For detection techniques that involve only energies much lower than the DM and nucleon masses, one typically finds several models that map to a given assumption on the energy-dependence of the cross section. For example, in 'standard' DD, the model we studied of fermion DM with scalar mediators realises a cross section which is constant in the energies relevant for those experiments. Explicitly, at those experiments the Mandelstam variable t is negligible with respect to m χ , m A and m φ in Eq. (8). This is not the case for CR-upscattered DM, so that the cross-section depends on the relevant energies. Therefore, the use of the latter limits and sensitivities, shown in colour in Figure 5, is perfectly consistent with the assumption of cross section constant at smaller energies. More than that, our results proves that they are actually the stronger existing constraint over a significant range of DM masses. On the contrary, the assumption of a cross section constant in the energies relevant for CR-upscattered DM, as done in a relatively wide portion of the literature on the subject, at present does not stand on the same solid footing. While that regime in energy could be approximately realised in the case of scalar DM with scalar mediated interactions, one should assess if the resulting limits on CR-upscattered DM (see e.g. [12]) test parameter space that is allowed by monojet searches, whose shape in the σ NR − m χ plane would be the same as in Fig. 5.

Other limits
We describe here other constraints on DM with interactions with nuclei, and constraints specific to the models of Eqs. (1) and (2). They are shown in Figures 4 and 5.

Cosmology
If the dark sector is light, of order MeV, it can modify the thermal history in the early universe. Recent studies of the BBN and CMB constraints on MeV-scale dark matter thermally coupled to the SM plasma in the early universe include [45][46][47][48][49][50][51][52][53].
BBN. In our benchmark points, we assume that m φ , m a ∼ O(1) GeV, and that these mediators decay well before BBN, thus not affecting BBN nor CMB. However, the dark matter can affect BBN (and CMB, see later) and hence can be constrained if m χ O(10) MeV. Although the precise bound depends on the specific couplings to the SM particles [53], to be conservative we simply assume that BBN requires m χ > 10 MeV. By switching on a coupling of the mediators to neutrinos, with appropriate size, one could relax the BBN limit to m χ O(1) MeV [51,53], and this is why we extend the range of some Figures to m χ = 1 MeV and we shade BBN differently than other constraints.
Dark Matter abundance. By using the only couplings of Eqs. (1), (2) and for the coupling values of our interest, DM is in chemical equilibrium with the early-universe bath as long as it contains quarks. After the QCD phase transition takes place, DM pairs could still annihilate into two (scalar-mediator) or three (pseudoscalar-mediator) pions. Therefore, achieving the correct relic abundance via thermal freeze-out would require m χ m π or m χ 3m π /2. One would then conclude that smaller masses overclose the universe, and are therefore excluded. We note however that a small coupling of the mediator to neutrinos (perhaps even to electrons in the scalar case) would allow to achieve the correct relic abundance via thermal freeze-out even for DM masses smaller than O(m π ), without worsening the BBN constraints discussed so far, and without introducing other ones. We then refrain from displaying lines corresponding to the correct relic abundance in Figs. 4 and 5, because they would depend on free parameters unrelated to the phenomenology presented there. In addition, if the DM annihilation cross section is larger than the thermal one (e.g. for m χ enough larger than m π or if we add a small coupling to neutrinos), then the symmetric DM component annihilates away leaving only a possible initial asymmetric component, dependent on other more UV physics.
CMB. Since CMB constraints are absent when DM is asymmetric enough [54,55], we refrain from displaying them in Figs. 4 and 5, because as we have seen they can be dealt with by physics that is unrelated to the phenomenology that is the focus of this work. If one was instead interested in the case where the DM abundance is set by the only couplings of Eqs. (1) and (2), so that as we have argued m χ > O(m π ), then i) in the scalar mediator case, CMB would not exclude the currently allowed regions of parameter space, because the DM annihilations into mesons and quarks (and mediators, if the masses allow) are all p-wave suppressed (see e.g. [56]); ii) in the pseudoscalar mediator case, CMB would exclude the model, because DM annihilations are instead in s-wave. Exceptions exist, e.g. if the thermal abundance is set by resonant annihilations [57], but we do not explore them further here.

Collider
Monojet at LHC. Dark Matter pairs can be produced at the LHC together with SM radiation, in processes like qq → gφ(a) → gχχ where, depending on its mass, the mediator φ (a) can either be off-shell, or be produced on-shell and then decay to DM. Thus the parameter space can be constrained by searches of monojet plus missing energy like [58,59], that we recast here to our model. We use FeynRules 2.3.32 [60] to generate the model files, and MadGraph 5 2.7.3 [61,62] to generate events pp → χχ+jets, up to three jets. We then shower with Pythia 8.2 [63] switching on MLM merging [64], and use Delphes 3.4.1 [65] to simulate the ATLAS detector. We impose weak cuts at generation level, the most important ones being htjmin = 180 and transverse momentum of each DM particle larger than 60 GeV. We then use Madanalysis5 [66] to reject events after detector simulation that do not satisfy the cuts described in [58]: we demand that the leading jet transverse momentum, p T (j 1 ), be larger than 150 GeV, that there is not a fifth jet with p T (j 5 ) larger than 30 GeV, and that the missing transverse energy (MET) be larger than 200 GeV. This defines the inclusive region called IM0 in [58] (we find that additional cuts on the azimuthal separation between the leading jet and MET have no further impact in reducing our events). Regions from IM1 to IM12 are then defined by the cut MET/GeV > 250, 300, 350, 400, 500, 600,... 1200. By using the merged cross sections given by Pythia 7 and the efficiencies of the various cuts given by MadAnalysis 5, we finally obtain the signal cross sections after cuts for each of the 13 regions above.
We then compare these cross sections with the model-independent limits, on the visible (i.e. multiplied by acceptance and efficiency) cross sections, given in Table 7 of [58] for each of the regions IM0,...IM12 above. We translate the strongest of these limits into an upper limit on g uφ = g dφ and g ua = g da . By doing so we are implicitly assuming that our simulation accounts for the entire acceptance and efficiency: in this sense our limits can be seen as aggressive ones, which is a conservative choice from the point of view of determining the parameter space open for the new searches at neutrino experiments proposed in this paper. We repeat the procedure for mediator masses of 1 and 3 GeV, and DM masses from 0.01 to 10 GeV, finding the limits displayed in Figs. 4 and 5, which for m χ < m a,φ /2 read g ua,uφ = g da,dφ 0.06. In the pseudoscalar case and for a selected subset of parameters, we repeat the procedure also for g uφ = −g dφ and g ua = −g da , and find comparable limits. This was expected because that relative sign matters only in interference diagrams, which are present only for multiple jet emissions and thus are less important. It was also expected that our limits are the same in the scalar and pseudoscalar case and that, for m χ < m a,φ /2, they are independent of the DM mass, because DM is produced from mediator decays and so what matters is the mass of the latter. As a check of our results, we found comparable limits also when we generated events with stronger generation level cuts (htjmin = 380 and transverse momentum of each DM particle larger than 120 GeV), so to have more statistics in the signal regions with more aggressive cuts. Finally, we checked our limits with another simulation carried out independently.
Meson invisible decay. The hadrophilic mediators and dark matter particles can be produced, and hence can be probed, by meson decays. Since the mediators couple to light quarks, light mesons such as π, K and η are relevant in our case. We are interested in the relatively heavy mediator mass, m φ , m a 1 GeV, and hence these mediators cannot be directly produced from the decay of π, K and η. Still the mesons can decay into dark matter particles and one can constrain the parameter region from their invisible decay. This is particularly relevant for the pseudoscalar mediator case since the pseudoscalar mediator has the same quantum numbers as neutral mesons and can mix with them.
Because of the coupling with light quarks (2), the pseudoscalar mediator a can mix with the neutral mesons. The chiral Lagrangian with the mediator a is given by M = diag(m u + ig ua a, m d + ig da a, m s + ig sa a).
Here f π = 93 MeV,m = (m u + m d )/2 = 3.5 MeV [67] and the eta prime decay constant is taken to be equal to the pion decay constant. The last term in Eq. (35) gives an anomalous mass term to η . By expanding this Lagrangian, one finds the following mass mixing terms between a and the mesons In the limit of small g's, the mixing angles are given as The total decay widths of neutral mesons are Γ π 0 = 7.7 × 10 −3 keV (τ π 0 = 8.5 × 10 −17 s), Γ η = 1.3 keV and Γ η = 1.9 × 10 2 keV [67]. The constraints on their invisible decay are Br(π 0 → invisible) < 2.7 × 10 −7 , Br(η → invisible) < 1.0 × 10 −4 and Br(η → invisible) < 6 × 10 −4 [67]. The invisible decay rates are given by We then obtain the following constraints on the couplings: Thus, the meson invisible decays in general put a severe constraint on the pseudoscalar mediator model. However, we note that some of these constraints can be evaded once we switch on the couplings to the strange quark. For instance, one can consider the SU(3) singlet case, g ua = g da = g sa . In this case, one can evade the mixing with π 0 and η, and only the η invisible decay applies. Alternatively one can take g sa = −2g ua = −2g da , which makes only the η invisible decay relevant. Since our constraint is insensitive to g sa , we shade the parameter region covered by both the η and η invisible decays in Figure. 4. The individual constraints are indicated as the dashed lines, implying that it depends on g sa which line is actually relevant. For these lines, we assume that g sa = g ua = g da for the η invisible decay and g sa = −2g ua = −2g da for the η invisible decay, respectively. Yet another option would be the isovector coupling g ua = −g da and g sa = 0. In this case, the π 0 invisible decay puts a severe constraint for m χ < m π 0 /2, but there is no meson decay constraints for m χ > m π 0 /2. However, our limits and sensitivities at neutrino experiments also get weaker due to the exact cancellation of the form factor in the large |t| limit (see Eqs. (17) and (18)). Since we could not derive a meaningful constraint on the couplings within our treatment, we do not seek this possibility any further. 8 If the mediators are lighter, m φ , m a 1 GeV, they can be directly produced by the decay of mesons. In this case, the dark sector can be probed directly by meson decays with invisible final-state particles, and also by beam-dump experiments where the mediators are produced on-shell by the decay of mesons, see e.g. [68] for the scalar mediator case. The models then suffer from severe constraints and we do not explore that mass region in this paper. 8 Within our treatment, only DUNE can put a constraint on the couplings in the isovector case. This is partly because our treatment of the Earth attenuation is too conservative since we ignore the form factor in the computation of the Earth attenuation. However, even if we include the form factor, the derived constraints are not strong, so we do not discuss this possibility in this paper.
Constraints on other couplings. As shown in Eqs. (1,2), in this paper, we assume the mediators couple only to light quarks. A rich phenomenology would of course arise if one switched on couplings of the mediators to other SM particles, see e.g. [69][70][71][72]. Here we just briefly mention that a coupling to the top quark would induce, at one-loop, a bsφ/a vertex constrained by B → Kφ/a(→ invisible) [73,74], and that bottom and charm couplings would be constrained by Υ → γφ/a(→ invisible) and J/ψ → γφ/a(→ invisible) [75].

Direct detection
The parameter space is constrained by conventional DD experiments. In the scalar mediator model, the dark matter can scatter with nucleons by a tree-level φ exchanging diagram. The spin-independent non-relativistic cross section is given as [68] where µ = (1/m χ + 1/m N ) −1 is the reduced mass for dark matter-nucleon scattering, A and Z are the mass and atomic numbers of the target particle, and y φpp and y φnn are effective φ-nucleon couplings defined as We take m u = 2.16 MeV, m d = 4.67 MeV, f N u = 0.0199 and f N d = 0.0431 for both N = p, n. In the pseudoscalar mediator model, a tree-level contribution is suppressed in the nonrelativistic limit both by the spin and by the velocity as shown in Eq. (3). The dark matter can still scatter with nucleons by one-loop a exchanging box diagrams [76][77][78][79][80]. Here we follow [80] to evaluate this contribution. The spin-independent non-relativistic cross section from the box diagrams is given as Here C q is the Wilson coefficient for a twist-0 operator, and C (1) q and C (2) q are the Wilson coefficient for twist-2 operators. They are given as where G, X 001 , and X 111 are loop functions defined as The quantities q N (2) andq N (2) are the second moments for quark distribution functions for nucleons whose values we extract from [80]. We combine the results from NEWS-G [81], DarkSide-50 [82], CDMSlite [83], CRESST-I II [84], XENON1T (S2 only) [85], and XENON1T (S1+S2) [4], and show the DD bound in Figures. 4 and 5. In the case of our interest, CRESST-III, DarkSide-50 and XENON1T give the strongest constraints in the relatively small, middle and large dark matter mass ranges, respectively.

Theoretical Consistency
The quark interactions in Eqs. (1) and (2) are effective descriptions in the phase where the electroweak symmetry is broken. A possibility to UV-complete the interactions with couplings g uφ and g ua is by adding to the model two left-handed Weyl fermions U L and U † R , with U L,R in the same SM gauge representations of the SM field u R . We then add to the SM Lagrangian (in Weyl notation, with H the SM Higgs doublet) and, depending on the mediator under consideration, L UV a = g U a iaU L u † R + h.c. .
Upon the Higgs taking VEV, Eq. (57) induces a mass mixing between u L and U L which, via Eqs. (58), (59), gives where v 246 GeV and M U = y 2 U v 2 /2 + M 2 U is the physical mass of the Dirac fermion (U L , U R ) (mass eigenstate up to relative orders y u v/M U , with y u the SM Yukawa coupling). LHC limits on new vector-like quarks constrain M U 1.3 TeV [86,87], so that g uφ,ua 0.14 g U φ,Ua y U .
As our new limits and sensitivities probe values of g uφ,ua smaller than 0.1, they test a region compatible with a fully perturbative UV completion (i.e. with g U φ and y U smaller than one). The UV completion of the couplings g dφ , g da , g sφ , g sa is entirely analogous.

Conclusions and Outlook
The quest for new experimental tests of sub-GeV Dark Matter has been a booming field in recent years. One such test [7,8] consists in relying on the DM sub-component that is upscattered by cosmic rays, and that thus possesses a larger kinetic energy than DM in the virialised halo. The existence of this component, which is unavoidable as long as DM possesses any interaction with the SM, opens the possibility to test sub-GeV DM with 'standard' DD experiments like XENON1T [7], and with large neutrino detectors built for other purposes, like Super-K [8].
Within this context and focusing on DM interactions with nucleons, in this paper we obtained the following new results: by the use of Super-K data on protons detected via their Cherenkov emission [15,16], we derived new limits on DM-nucleon cross sections that are stronger than all previous ones relying on CR-upscattered DM [7,12]. We also estimated related sensitivities at Hyper-K and DUNE; we casted our limits and sensitivities on the parameter space of two simplified models, where DM is a fermion and its interactions with quarks are mediated by a new scalar or pseudoscalar particle, cf. Eqs. (1) and (2). This allowed us to put on solid grounds our proposal to test light DM at large neutrino experiments, because i) we took into account the energy dependence of the DM scattering cross sections with nuclei, which is necessary whenever the energies relevant for different detectors are different (as it is the case with Super-K vs DUNE vs XENON1T); ii) we compared our new limits and sensitivities with other tests of these models, like 'standard' DD, BBN, monojet and meson decays, see Sec. 4. Note that the last two tests require to specify the mediator of DM interactions and its mass. As a byproduct of this comparison, we recasted monojet searches to DM models with scalar and pseudoscalar mediators, and we assessed the status of limits on sub-GeV pseudoscalar coupling to quarks.
iii) we assessed the theoretical consistency of the parameter space we test, by explicitly building a UV completion that is free from other constraints, see Sec 5.
Our main results are shown in Fig. 4 in the parameter space of the simplified models. For the reader's convenience, in Fig. 5 we translated them to the usual plane of DM mass vs spinindependent non-relativistic cross section. We stress that the use of our limits and sensitivities in that plane, together with others that assume a cross section constant in the energies relevant for other detection techniques, is not only consistent, but stands on very solid footing, see Sec. 3.4 for more details.
Our results open several future directions of exploration. First and foremost, they motivate the realisation of our proposed searches at Super-K, Hyper-K and DUNE. They also pave the way to study the impact of these searches, and their interplay with other experimental probes, in other simplified models of hadrophilic DM. Finally, it would be interesting to study the case where DM couples to both leptons and hadrons. We plan to come back to some of these directions in future work.