No room to hide: implications of cosmic-ray upscattering for GeV-scale dark matter

The irreducible upscattering of cold dark matter by cosmic rays opens up the intriguing possibility of detecting even light dark matter in conventional direct detection experiments or underground neutrino detectors. The mechanism also significantly enhances sensitivity to models with very large nuclear scattering rates, where the atmosphere and rock overburden efficiently stop standard non-relativistic dark matter particles before they could reach the detector. In this article, we demonstrate that cosmic-ray upscattering essentially closes the window for strongly interacting dark matter in the (sub-)GeV mass range. Arriving at this conclusion crucially requires a detailed treatment of both nuclear form factors and inelastic dark matter-nucleus scattering, as well as including the full momentum-transfer dependence of scattering amplitudes. We illustrate the latter point by considering three generic situations where such a momentum-dependence is particularly relevant, namely for interactions dominated by the exchange of light vector or scalar mediators, respectively, and for dark matter particles of finite size. As a final concrete example, we apply our analysis to a putative hexaquark state, which has been suggested as a viable baryonic dark matter candidate. Once again, we find that the updated constraints derived in this work close a significant part of otherwise unconstrained parameter space.


Introduction
The strategies to search for a dark matter (DM) component in the Universe are nowadays extremely varied, targeting many possible gravitational and non-gravitational properties such as the DM mass or standard model (SM) couplings [1]. In astrophysical, cosmological, and laboratory settings, this broadband approach has yet to conclusively reveal any nongravitational signatures. However, via both indirect and direct searches, the very wide DM model space has been significantly restricted. The focus of this article concerns the reach of the generic class of experiments aiming to directly detect DM through a possible DM-nucleon coupling [2], known as direct detection facilities. Currently, world-leading examples of this setup include e.g. LUX-ZEPLIN (LZ) [3], PandaX-4T [4], and Xenon-1T [5], which set the strongest limits in the DM mass m χ vs. spin-independent nuclear coupling σ SI parameter space.
The sensitivity of a given direct detection experiment is controlled by a number of factors. Firstly, the event rate Γ N scales with the number of DM particles that have a sufficiently large kinetic energy. Specifically, the DM energy must be large enough to induce a nuclear recoil that can trigger a signal above the detector threshold. Secondly, the rate also scales linearly with the DM-nucleon cross section dσ χN /dT N , at least in the above examples, where T N is the nuclear recoil energy. Thirdly, as in any count-based experiment, this signal rate should be compared to some background event rate to derive a statistically significant detection threshold. Notably, in direct detection facilities, the background rates are typically extremely low as necessitated by the small expected signal rates, although there are some important exceptions, such as a dedicated CRESST surface run [6].
The standard target for these experiments is the DM in the Galactic halo, which has characteristic velocities of the order v χ ∼ 10 −3 c and in any case cannot exceed the Galactic escape velocity v esc ∼ 540 km/s [7,8]. For a given DM mass m χ , there is hence unavoidably a maximum DM kinetic energy available to excite nuclear recoil signals of the order T N ∼ m 2 χ v 2 esc /m N . For some DM mass m min χ this recoil energy must fall below the detectable threshold, and the experimental sensitivity drops to zero. For experiments such as Xenon, PandaX and LZ, it is well-known that this cut-off lies around the GeV-scale, corresponding to a detectable threshold in the keV range. As such, even though these detectors have impressive reach -currently down to the level of spin-independent cross sections of σ SI ∼ 10 −47 cm 2 [3][4][5], and even approaching the neutrino floor [9,10] with ongoing searches -there is ample motivation (and hence, in fact, both experimental and theoretical activity) for methods to probe the sub-GeV mass range [11,12]. This describes the first "window" in which DM can hide -it could just be that DM has a small mass out of the reach of direct detection experiments. There is yet another window at large values of the cross section σ SI , however, which will be a key focus of this article. This arises due to the fact that if DM interacts too strongly, then it can actually be the case that DM particles are unable to reach the detectors due to the attenuation of the flux in the atmosphere or the rock overburden [13][14][15]. This typically becomes the main prohibitive factor for cross sections at the level of σ SI 10 −28 cm 2 [16].
There have been a number of promising experimental proposals to probe these two open windows. Attempts to extend the sensitivity to DM-nucleus interactions into the sub-GeV realm include searches for Migdal electrons [17,18] or bremsstrahlung photons [19], accompanied by an intense low-threshold direct detection program in the development of novel detector concepts (for a recent review, see Ref. [12]). Cross sections sufficiently large for DM to scatter inside the Earth before reaching underground detectors, on the other hand, can be probed by surface runs of conventional direct detection experiments (like the one performed by the CRESST collaboration [6]), or by targeting the expected diurnal modulation in the signal in this case [20,21]. As far as this work is concerned, however, we will be interested in the role played by the irreducible astrophysical flux of highly boosted DM that originates from cosmic ray collisions with DM particles in the Galactic halo (CRDM). This was pointed out only relatively recently [22,23], and subverts the issue of a loss in sensitivity by noting that a sub-dominant component of DM with velocities well above those in the Galactic halo can produce a detectable signal even if it is very light, i.e. for DM masses (well) below 1 GeV. The sub-dominant nature of the flux naturally introduces a trade-off with the interaction rates that can be probed, quantitatively resulting in limits at the level of σ SI ∼ 10 −31 cm 2 [22]. Interestingly, CRDM does not only probe previously open parameter space at small DM masses but also results in bounds extending into the relevant regime of the second open window described above. After this initial work pointed out the advantages of considering such a boosting mechanism, a large number of further analyses have addressed various aspects of the production , attenuation [46,47], and detection [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62] of astrophysically boosted DM. For a recent comprehensive (re-)analysis of all of these aspects see, e.g. Xia et al. [63], who stressed in particular that form-factor suppressed attenuation in the overburden seemingly allows us to exclude cross sections much larger than σ SI ∼ 10 −28 cm 2 .
This article builds on this literature in three important ways: firstly, we point out that when DM acquires such large energies, inelastic scattering in the rock overburden above detectors such as Xenon-1T will at some point become the dominant attenuation mechanism. As such, to avoid being over-optimistic in terms of how much parameter space is excluded, we show how to include this physical effect in a self-consistent manner and derive the resulting bounds. Secondly, we broaden the applicability of these limits to models that are more realistic for DM with sub-GeV masses, moving beyond simplified contact interactions to interactions mediated by vector or scalar mediators, or DM that has some internal structure. Finally, we argue that with these improvements, and when taking into account fully complementary constraints from cosmology, there is generically no remaining open parameter space left unconstrained for nuclear cross sections exceeding 10 −30 cm 2 , for DM masses in the entire MeV to GeV range. We demonstrate that possible loopholes to this statement -still allowing an open window at larger cross sections -require a combination of (i) questioning the principal ability of CRESST to probe DM masses down to the published limit of m χ = 140 MeV [6] and (ii) choosing a rather narrow range of mediator masses m φ ∼ 30 MeV (or finite DM extent r χ ∼ 10 fm). For our numerical analysis throughout the article, we use the package DarkSUSY [64]. The improved CRDM treatment reported in this work, including also updated cosmic ray fluxes and a more sophisticated use of form factors in the attenuation part, will be included in the next public release of the code.
The rest of the article is organized as follows: we start in section 2 by briefly reviewing the production of CRDM and the attenuation of the subsequent flux on its way to the detector, establishing our notation and setting up the basic formalism that our analysis relies on. In the next two sections, we discuss in more detail how to model nuclear form factors (section 3) and the impact of inelastic scattering (section 4) on the attenuation of the flux. In section 5, we consider a number of generic options for the Q 2 -and s-dependence of the scattering amplitude that are more realistic than assuming a constant cross section. We complement this in section 6 with the analysis of a specific example, namely a baryonic DM candidate that has been argued to evade traditional direct detection bounds despite its relatively strong interactions with nuclei. We conclude and summarise our results in section 7.

Cosmic-ray upscattering of dark matter
We describe here, in turn, how initially non-relativistic DM particles in the Galactic halo are up-scattered by cosmic rays (CRs), how the flux of these relativistic CRDM particles is attenuated before reaching detectors at Earth, and how to compute the resulting elastic scattering rate in direct detection experiments.
Production: The basic mechanism that we consider is the elastic scattering of CR nuclei N , with a flux of dΦ N /dT N , on non-relativistic DM particles χ in the Galactic halo. For a DM mass m χ and density profile ρ χ (r), this induces a relativistic CRDM flux incident on Earth of [22,46] Here r denotes the Galactic position, and dσ χN /dT χ is the differential elastic scattering cross section for accelerating a DM particle to a kinetic recoil energy T χ . For DM particles initially at rest, this requires a minimal CR energy T min Furthermore, in the second line of Eq. (2.2), we have introduced an effective distance D eff that allows us to express the CRDM flux in the solar system in terms of the relatively well measured local interstellar CR flux, dΦ N LIS /dT N , and the local DM density, for which we adopt ρ local χ = 0.3 GeV/cm 3 [65] (noting that our final limits are independent of this choice). The advantage of this parameterisation is that uncertainties deriving from the integration over the volume relevant for CRDM production, dΩ d , are captured in a single phenomenological parameter D eff . Indeed, despite the complicated underlying physics, this parameter is surprisingly well constrained, with uncertainties dominated by the vertical extent of the confinement zone of Galactic CRs. In what follows, we will use a fiducial value of D eff = 10 kpc. 1 We note that our final limits only depend logarithmically on this quantity, for large interaction rates, or scale as D −1/2 eff when attenuation in the soil or atmosphere is inefficient, respectively.
When computing the CRDM flux in Eq. (2.2), we take into account the four most abundant CR species, N = {p, He, C, O}, for which high-quality determinations of the local interstellar fluxes exist [68]. The fluxes of heavier nuclei are subject to significant uncertainties for the energies of interest to us, see e.g. the discussion in Ref. [69], not least due to apparent discrepancies between AMS-02 data [70][71][72] and earlier measurements. We also note that the CRDM flux contribution from these heavier elements is strongly formfactor suppressed at large T χ , see section 3, and hence anyway not relevant for constraining DM with masses m χ 0.1 GeV.
Attenuation: On its way to the detector, the CRDM flux given by Eq. (2.2) is attenuated due to scattering of the CRDM particles with nuclei in the atmosphere and soil (overburden) above the experimental location. This effect can be well modelled by the energy loss equation which can be used to relate the average kinetic energy at depth z, T z χ , to an initial energy T χ at the top of the atmosphere. Here, the sum runs over the nuclei N in the overburden, i.e. no longer over the CR species, and ω χ is the energy loss of a DM particle in a single collision. For elastic scattering, ω χ is equal to the nuclear recoil energy T N . In that case, the maximal energy loss of a DM particle with initial kinetic energy T z χ is given by is the (squared) CMS energy of the process. For inelastic scattering on the other hand, which we will discuss in more detail in section 4, the energy loss can in principle be as high as ω max χ = T z χ . For the purpose of this work we will mostly be interested in the Xenon-1T detector, located at a depth of z = 1.4 km in the Gran Sasso laboratory. In this case the limestone overburden has a density of 2.71 g/cm 3 [73], mostly consisting of an admixture of CaCO 3 and MgCO 3 , and attenuation in the atmosphere can be neglected; in terms of weight percentages the dominant elements are O (47.91%), Ca (30.29%), C (11.88%), Mg (5.58%), Si (1.27%), Al (1.03%) and K (1.03%) [74]. We note that Eq. (2.4) only provides an approximate description of the stopping effect of the overburden, which is nonetheless sufficiently accurate for our purposes. For a detailed comparison of this approach with Monte Carlo simulations of individual particle trajectories, see Refs. [16,63,[75][76][77] Detection: The elastic scattering rate of relativistic CRDM particles arriving at underground detectors like the Xenon-1T experiment is given by Note that the above integral is over the energy of the DM particles before entering the atmosphere. On the other hand, the elastic scattering cross section dσ χN /dT N must still be evaluated at the actual DM energy, T z χ , at the detector location, which requires numerically solving Eq. (2.4) for T z χ (T χ ). The lower bound on the integral then represents the minimal initial CRDM energy that is needed to induce a nuclear recoil of energy T N at depth z, i.e. T min ). This can be obtained by inverting the solution of Eq. (2.4), where T z,min χ is given by the right-hand side of Eq. (2.3) under the replacement (T χ , m χ , m N ) → (T N , m N , m χ ). In general, the elastic nuclear scattering cross section dσ χN /dT N is a function of both s and the (spatial) momentum transfer, If the dependence on s can be neglected or the (dominant) dependence on Q 2 factorizesas in the case of standard form factors -then the rate in the detector given in Eq. (2.7) will have an identical Q 2 -dependence as compared to the corresponding rate expected from the standard population of non-relativistic halo DM. As pointed out in Ref. [22], this salient feature makes it possible to directly re-interpret published limits on the latter (conventionally expressed as limits on the scattering cross section with protons) into limits on the former. Otherwise, for an accurate determination of the expected count rate in a given analysis window, one would in principle have to also model the detector response in the evaluation of Eq. (2.7) and then infer limits based on the full detector likelihood (e.g. with a tool like DDCalc [78,79]).

Nuclear form factors
The target nuclei used in direct detection experiments are typically larger than the de Broglie wavelength of DM with standard Galactic velocities, at least for heavy nuclei, implying that the incoming DM particles only 'see' part of the nucleus. Since the elastic scattering process is fundamentally induced by a coupling between DM and the constituents of these nuclei, this means that it should be suppressed by a nuclear form factor, G 2 (Q 2 ), compared to the naive expectation that the nuclear cross section is merely a coherent sum of the cross sections of all the constituents (for recent pedagogic accounts of conventional direct DM searches, see e.g. Refs. [80,81]). 2 For CRDM, this effect is amplified, given the smaller de Broglie wavelengths associated to the faster moving upscattered DM particles. These nuclear form factors are essentially Fourier transforms of the number density of nucleons inside the nucleus, usually approximated by the experimentally easier accessible charge density. A common parameterization is the one suggested by Helm [82], which is based on modelling the nucleus as a hard sphere with a Gaussian smearing (in configuration space). For heavy nuclei we follow instead a slightly more accurate approach and implement model-independent form factors [83], based on elastic electron scattering data. Concretely, we implement their Fourier-Bessel (FB) expansion approach, with parameters taken from Ref. [84]. For nuclei where the FB parameters are not available, notably Mg and K, we use model-independent Sum of Gaussians (SOG) form factors instead. 2 We focus here on spin-independent elastic scattering. For spin-dependent scattering, the sum would not be coherent and hence generally result in much smaller cross sections. This prevents standard DM from being stopped in the overburden before reaching the experimental location -unless the scattering cross section per nucleon is so large that it becomes incompatible with other astrophysical constraints. A detailed treatment of attenuation in the Earth's crust is, hence, less relevant in this case. For Q 2 (0.1 GeV) 2 one starts to resolve the inner structure of the nucleons themselves, which we discuss in more detail in section 4. Let us however briefly mention that in the case of He, this effect is already largely captured by the above description in that we take the SOG form factors from Ref. [84] (thus improving on the simple dipole prescription used, e.g., in Ref. [22]). For the proton, we adopt the usual dipole nucleon form factor, noting that the nuclear form factor would formally equal unity, with Λ p = 0.843 GeV. This provides a very good fit to experimental data up to momentum transfers of at least Q 2 ∼ 1 GeV 2 , with an agreement of better than 10% for Q 2 ≤ 10 GeV 2 [85,86]. We note that our final results are highly insensitive to such large momenta.
In the rest of the section, we will briefly describe the impact of nuclear form factors on the CRDM flux and the attenuation of this flux on its way to the detector. In both cases the effect is sizeable, motivating the need for a precise modelling of G 2 (Q 2 ).

Impact on production
The solid lines in Fig. 1 show the expected CRDM flux before attenuation, cf. Eq. (2.2), for a range of DM masses. For the purpose of this figure, we have assumed a constant elastic scattering cross section σ p SI = σ n SI on nucleons, i.e. a nuclear cross section given by Here, describes the usual coherent enhancement, in this case proportional to the square of the atomic number A of nucleus N . In the rest of the expression, µ χN (µ χp ) is the reduced mass of the DM/nucleus (DM/nucleon) system and the maximal DM energy T max χ that can result from a CR nucleus with energy T N is given by the right-hand side of Eq. (2.5) after replacing T z χ → T N and m χ ↔ m N . In the left panel of the figure, we show that neglecting nuclear form factors (dashed lines) would lead to a significant overestimate of the CRDM flux at high energies. For m χ 0.1 GeV, the form factor suppression even becomes the dominant effect to determine the overall normalization of the flux, while for lower DM masses, the peak of the distribution is entirely determined by the fact that the CR flux itself peaks at GeV energies. This suppression in the flux leads to a rapid deterioration of CRDM limits. Modelling form factors correctly is thus particularly important for the highest DM masses that can be probed by cosmic-ray upscattering, i.e. for m χ ∼ 1 − 10 GeV.
In the right panel of Fig. 1, the contributions from the individual CR nuclei to the CRDM flux are shown. At low energies the dominant contribution is always from Helium, closely followed by the one from protons. The high-energy part of the CRDM flux, on the other hand, is almost exclusively due to CR protons because the contribution from heavier CR nuclei is heavily form-factor suppressed. In addition, for m χ 1 GeV, the peak amplitude of the CRDM flux -which typically has the most constraining power in direct detection experiments -is almost exclusively determined by CR p and He nuclei (see also Fig. 2 below to better gauge the relevant range of energies after attenuation in the overburden). For lower DM masses, on the other hand, including further high-Z CR species than those taken into account here could in principle increase the relevant part of the CRDM flux by up to ∼ 50 % [63]. In what follows, we conservatively neglect these contributions, in view of both the larger uncertainties in the underlying CR fluxes and the fact that we are mainly interested in DM masses around the GeV scale.

Impact on attenuation
We now turn our attention to assessing the effect that the form factor suppression has on the attenuation of DM particles on their way to the detector in a direct detection experiment. For concreteness we will again focus on the case of Xenon-1T, where Xe nuclei recoiling with an energy of at least T Xe = 4.9 keV trigger a detectable signal [5]. In Fig. 2, we show the minimal initial DM energy that is required to kinematically allow for this, after penetrating through the Gran Sasso rock. In practice this is done by numerically solving Eq. (2.4) with DarkSUSY. Dash-dotted lines indicate the result when conservatively assuming that the Figure 2. Minimal kinetic energy T χ that a DM particle must have at the surface of the Earth (z = 0) in order to trigger a signal in the Xenon-1T experiment, as a function of a (constant) spin-independent scattering cross section σ p,n SI on nucleons. Different colors correspond to different DM masses, as in Fig. 1. Dash-dotted lines show the kinetic energies that would be necessary when computing the attenuation in the zero momentum transfer limit. Dashed lines illustrate the effect of adding the expected form factor suppression, cf. section 3, while solid lines show the result of our full treatment, including also inelastic scattering events (discussed in section 4). stopping power in the overburden is as efficient as in the zero-momentum transfer limit (as in Ref. [22]), while dashed lines show the effect of adding the additional form factor suppression for high Q 2 (as in Refs. [38,63]). Solid lines, finally, demonstrate the effect of also adding the attenuation power of inelastic scattering events, as described in detail below in Section 4.
For small cross sections, attenuation is inefficient and, as expected, the three approaches give the same answer. In this limit, the difference in the required DM energy is entirely due to the well-known kinematic effect, cf. Eq. (2.3), that lighter particles require a higher energy to induce a given recoil of much heavier particles (up to a minimum energy of T χ ≥ m Xe T Xe /2 = 17.3 MeV in the limiting case where m χ → 0). Correspondingly, this also means that the CRDM fluxes cannot actually be probed by Xenon-1T for the entire range of T χ shown in Fig. 1; unless m χ 10 MeV, however, the lowest detectable energy is always smaller than the energy at which the CRDM flux peaks.
For large cross sections, on the other hand, Fig. 2 shows a pronounced difference between the three approaches: while in the case of a constant cross section (dash-dotted lines) the energy loss equation results in an exponential attenuation, adding form factors (dashed lines) implies that the required initial DM energy only rises as the square root of the scattering cross section in the Q 2 = 0 limit. In fact, we note that this is exactly the behaviour one would expect from Eq. (2.4) for a cross section that falls off very rapidly at large momentum transfers. Comparing again to Fig. 1, this correspondingly enlarged range of kinetic energies that becomes kinematically accessible to Xenon-1T will inevitably lead to significantly larger rates in the detector -which, indeed, is exactly the conclusion reached in Refs. [38,63]. However, such a strong suppression of the physical stopping power of the Gran Sasso rock for a relativistic particle is highly unphysical. As we discuss in the next section, this is simply because the DM particles will start to scatter off the constituent nucleons themselves, albeit not coherently across the whole nucleus. Adding this effect (solid lines), results again in exponential attenuation in the overburden -though only at significantly larger cross sections than what would be expected when adopting a constant cross section for simplicity.

Inelastic Scattering
Our discussion so far has largely neglected the impact of inelastic scattering events of relativistic DM particles incident on nuclei at rest, or vice versa. Physically, the inclusion of inelastic scattering processes is non-negotiable and should be considered in a full treatment. This is because, whilst the form factor suppression described above is the relevant feature in the transition from coherently scattering off the whole nucleus to only parts of it, once the DM or nucleus transfers a sufficiently large amount of energy ω, the scattering will probe individual nucleon-, or even quark-level processes. The result is an additional contribution to the total scattering cross section that can easily dominate in the large energy transfer regime. As far as CRDM limits are concerned, the most important effect that the inclusion of inelastic scattering modifies is the attenuation of the flux through the Earth or atmosphere. Not including it, therefore, will lead to an overly optimistic estimate as to the amount of parameter space that is ruled out via this mechanism. 3 Let us note that inelastic scattering of non-relativistic DM, resulting in the excitation of low-lying states in the target nuclei, was previously both studied theoretically [19,[87][88][89] and searched for experimentally [90][91][92][93]. Here we concentrate on different types of inelastic processes that are only accessible to nuclei scattering off high-energy DM particles.
The rest of this section is organised as follows: firstly we give a qualitative description of the most important inelastic scattering processes, such as the excitation of hadronic resonances or quasi-elastic scattering off individual nucleons. Secondly, we explain how we obtain a quantitative estimate of these complicated nuclear interactions by making a direct analogy to the case of neutrino-nucleus scattering. In this regard, we make use of the public code GiBUU [94,95]. Finally, we will explain how to build this into the formalism described in section 2 in terms of the DM energy loss, see Eq. (2.4).

Scattering processes and associated energy scales
There are a number of relevant contributions to scattering cross sections on nuclei that are associated to certain characteristic energies or nuclear length scales. In the highly nonrelativistic limit, as described above, coherently enhanced elastic scattering dominates. At somewhat higher energies, more specifically momentum transfers corresponding to (inverse) length scales smaller than the size of the nucleus, the elastic scattering becomes form factor suppressed -a description which physically assumes a smooth distribution of scattering centres throughout the nucleus. The main characteristic of elastic scattering in both of these regimes is that the energy loss of the incident DM particle is uniquely related to the momentum transfer by ω = Q 2 /(2m N ).
This relation no longer holds for inelastic scattering processes, which are expected to become relevant at even higher energies. For our purposes, these inelastic processes can be broadly split up into three scattering regimes, depending on the energy that is transferred (see also Fig. 3 below, as well as a review [96] for the discussion of the analogous situation in the case of neutrino-nucleus scattering): • Quasi-Elastic Scattering (ω 10 −2 GeV): At suitably large energy transfers, the form factor suppression cannot be totally physical. This is because the incident DM particles will probe directly the constituent nucleons, which are inherently not smoothly distributed. Quasi-elastic scattering (QE) dominates for 10 −2 GeV ω 1 GeV, and describes this situation, i.e. where the dominant scattering is directly off individual protons (and neutrons) inside the nucleus, χ p(n) → χ p(n).
• Excitation of Hadronic Resonances (ω 0.2 GeV): At higher energies still, DM-nucleon scattering can excite nuclear resonances such as χ p → χ (∆ → pπ 0 ) etc., leading to a wide variety of hadronic final states. Often, the contribution due to the lowest lying ∆ resonances (DR) is distinguished from contributions from higher resonances (HR) since the former can be well resolved and starts playing role at considerably smaller transferred energies. In a complicated nucleus such as 16 O, both the QE and resonance contributions to the scattering cross section must be resolved numerically, taking into account effects such as the nuclear potential and spin statistics.
• Deep Inelastic Scattering (ω 1 GeV): Most DM couplings to nuclei and nucleons result from more fundamental couplings to quarks or gluons. As such, once the energy transfer is large enough to probe the inner structure of the nucleons (ω 1 GeV), then deep inelastic scattering (DIS) of DM with partons inside the nucleons can occur. Again, this should be resolved numerically to give an accurate estimate of the impact at the level of the scattering cross section.

Computation of the inelastic cross section for neutrinos
Due to the complicated nuclear structure of the relevant atomic targets in the Earth, or in the composition of cosmic rays, it is typically not possible to analytically compute all the contributions to DM-nucleus scattering described above. Instead, to estimate their impact on our conclusions and limits, we will make a direct connection with the physics of neutrino-nucleus scattering for which numerical codes -such as GiBUU [94] -are capable of generating the relevant differential cross sections.
In more detail, we draw the analogy between neutral current neutrino-nucleon scattering via processes such as ν p → ν p and DM-nucleon scattering. Numerically modelling the neutral current quasi-elastic scattering, resonances and deep inelastic scattering as a function of the energy transferred to the nucleus, ω, allows us to understand the relative importance of these processes as a function of the incoming neutrino energy (or DM kinetic energy T χ ). Of course, since these codes are tuned for neutrino physics, simply outputting the differential cross sections such as dσ νN /dω is not sufficient. To map the results onto DM, see section 4.3 below for further details, we should re-scale the results so as to respect both the relative interaction strengths and model dependences such as e.g. the mediator mass. In general, we expect this approach to provide a good estimate of the DM-nucleus cross section (at least) for contact interactions and scattering processes dominated by mediators in the t-channel.
At the level of implementation, we choose the settings in the GiBUU code described in Tab. 1 (see end of text). Since we are interested in quantifying the effect of inelastic scattering on the attenuation of the CRDM flux as it passes through the Earth, we mostly focus on the total inelastic scattering cross section, i.e. the sum over all the processes described in the previous section. We numerically calculate this for the most abundant nuclei in the Gran Sasso rock, N = {O, Ca, C, Mg, Si, Al, K}. Fundamentally, inelastic cross sections are expressed in terms of double-differential cross sections like d 2 σ νN /dQ 2 dω, since for inelastic scattering Q 2 and ω are independent variables. For integrating the energy loss equation, Eq. (2.4), however, it suffices to compute On the other hand, the full information about the Q 2 -dependence of d 2 σ νN /dQ 2 dω provided by GiBUU still remains a highly useful input to our analysis. This is because the double-differential cross sections of the individual inelastic processes turn out to sharply peak at values of Q 2 that have simple relations to ω. For example, the peak position for the QE contribution corresponds to the 'elastic' relation (2.8) for nucleons. As described below, this information will be used for setting realistic reference values of Q 2 to capture the model-dependence of the DM cross sections.

Mapping to the dark matter case
Having described the technical details of how we obtain the neutrino-nucleus inelastic cross sections using GiBUU, we now turn our attention to the mapping of these quantities onto DM models. This is a necessary step for two broad reasons: (a) the interaction strength governing the DM-nucleus interactions is typically very different from the neutrinonucleus SM value, and (b) the way the interaction proceeds via e.g. a contact interaction or mediator exchange can lead to substantially different kinematics and non-trivial Q 2 -or s-dependences. The total scattering cross section dσ χN /dω consists of the coherent elastic scattering contribution that we compute analytically for each of the models considered in this work, and the inelastic scattering cross section that we want to estimate based on the GiBUU output: Here dσ SI /dω| el is the differential DM-nucleon elastic cross section, excluding nucleon form factors such as the one given in Eq. (3.1). The sum runs over the various individual processes, i ∈(QE, DR, HR, DIS), which all have characteristic reference values of where the respective inelastic cross section peaks. In the second step above, we thus choose to rescale the inelastic scattering events to the elastic scattering off a point-like nucleon. This rescaling is motivated by the fact that for inelastic contributions like QE, the underlying process is much better described by scattering on individual nucleons than on the entire nucleus. The factor thus quantifies the ratio of the inelastic scattering process on a nucleus to the elastic scattering on an individual nucleon.
We now make the simplifying assumption that this ratio is to a certain degree modelindependent, based on the expectation that DM should probe the inner structure of nucleons in a similar way as neutrinos do when only neutral current interactions are involved. Physically, indeed, this closely resembles the situation both for contact interactions and t-channel mediators. The model dependence thus dominantly comes from the structure of the term dσ SI /dω| el , and we approximate Here, the inelastic neutrino-nucleus cross section dσ i νN /dω inel (E ν , ω) can be obtained using the GiBUU code, as described in section 4.2, and we evaluate it at the incoming DM kinetic energy, E ν = T χ . On the other hand, a possible estimate for the denominatorthe elastic neutral current neutrino-nucleon cross section without the form factor -is the average of the proton and neutron cross sections in the ω → 0 limit [96]: Here τ p 3 = 1 and τ n 3 = −1, θ W is the weak mixing angle and G F is the Fermi constant. The axial vector and strange quark contributions are encoded in the parameters ∆ S ≈ −0.15 (see, e.g., Ref. [97] for a discussion) and g A = 1.267 [98], respectively. Numerically the square bracket evaluates to a factor of ∼ 2.24 (2.01) for neutrons (protons). Let us stress, however, that this formula is valid only for energies relevant for inelastic scattering, 0.1 GeV E ν 10 GeV. At much smaller energies, only the valence quarks contribute to the scattering, and we would instead have for neutrons, while the scattering on protons is strongly suppressed by a factor of It is worth noting that in principle, we could improve the assumption made in Eq. (4.4) for the quasi-elastic process, because there is a well-controlled understanding of the analytic QE cross section via the Llewellyn-Smith formalism (see section V of Ref. [96]). For clarity, we choose to take a consistent prescription across all inelastic processes, and we have checked that including the full QE cross section would only introduce an additional O(1) factor in the DM QE cross section. For the numerical implementation in DarkSUSY, we pre-tabulate I ν,i from T χ = 0.01 GeV up to energies of T χ = 10 GeV, with 200 (101) equally log-spaced bins in T χ (ω) and a normalization as given by Eq. (4.5), and then interpolate between these values. 4 We also must choose the reference values for the transferred momentum Q 2 i,ref , which allows us to account for e.g. mediators that may be much lighter than the electroweak scale. Importantly, each process (quasi-elastic, ∆-resonance,...) is expected to have a different characteristic Q 2 -ω dependence that takes into account the relevant binding energies and kinematic scaling. For example, in the case of elastic scattering, the relation Q 2 = 2m N ω holds, whilst for quasi-elastic processes, the relevant scattering component is a nucleon such that the cross section is peaked around Q 2 ∼ 2 m ω, where m ≡ (m n + m p )/2. The resonance of a particle with mass m res can be accounted for by noting that part of the transferred kinetic energy is used to excite the resonance, such that the cross section peaks around Q 2 ∼ 2 m (ω − (m res − m)). We have confirmed these expectations numerically by comparing directly to the doubly-differential cross section extracted from GiBUU. From this 4 For significantly higher energies, GiBUU is no longer numerically stable. Furthermore, the underlying equations that describe the interaction processes begin to fall outside their ranges of validity as the Z boson mass starts to get resolved. At higher energies, where anyway only the DIS contribution is nonnegligible, a reasonable estimate can still be obtained by a simple extrapolation GeV. By running GiBUU up to Eν ∼ 30 GeV, we checked that this prescription traces the peak location (in ω) of the DIS contribution very well, independently of the exact choice of T ref χ . We also confirmed that the peak value of I becomes roughly constant for such large energies. On the other hand, higher-order inelastic processes are expected to become increasingly important at very large energies, not covered in GiBUU. We therefore only add the above extrapolation as an option in DarkSUSY, and instead completely cut the incoming CRDM flux at 10 GeV in the default implementation. As a result, our bounds on the interaction strength may be overly conservative for small DM masses mχ 0.1 GeV. numerical comparison we further extract that Q 2 ∼ 0.6 m (ω−ω DIS ), with ω DIS = 1.0 GeV, constitutes a very good fit to the peak location of the DIS cross section. In summary, we take the following reference values across the four inelastic processes: Here, ∆m ∆ = 0.29 GeV is the mass difference between the ∆ baryon and an average nucleon, and ∆m res = 0.40 GeV is an estimate for the corresponding average mass difference of the higher resonances (we checked that our final limits are insensitive to the exact value taken here).
To illustrate this procedure concretely, we consider the simple case of a contact interaction where, cf. Eq. (3.2), dσ SI /dω| el. = σ SI /ω max and ω max = 2 m(T 2 χ + 2χT χ )/((m + m χ ) 2 +2mT χ ). The results for the rescaled inelastic cross section (blue) are shown in Fig. 3 for a DM mass m χ = 1 GeV incident on a 16 O nucleus. In this figure, we also compare to the coherent elastic contribution (green) and highlight the balance between the relative contributions to the total (integrated) cross section σ tot χN . In particular, we see that above kinetic energies T χ 0.2 GeV, the inelastic contribution dominates, clearly motivating the necessity of its inclusion. This is consistent with the picture previously encountered in Fig. 2, where we could see the impact of inelastic scattering on the energy loss. More concretely, the result lies in some intermediate regime between the G(Q 2 ) = 1 and G(Q 2 ) = 1 cases, the former/latter leading to conservative/overly optimistic limits respectively. In the next section we will derive the relevant CRDM limits in the σ SI − m χ plane for a number of models to make this point quantitatively.
Let us conclude this section by briefly returning to the implicit assumption of isospinconserving DM interactions that we made above, with σ SI = σ p SI = σ n SI . Interestingly, neutral-current induced inelastic scatterings between neutrinos and nucleons hardly distinguish between protons and neutrons [96], such that the factor I χ,i ≈ I ν,i indeed becomes, by construction, largely independent of the nucleon nature. Naively, one would thus conclude that isospin-violating DM couplings can easily be incorporated in our treatment of inelastic scattering by replacing σ SI → (1/A) × (Zσ p SI + (A − Z)σ n SI ) in Eq. (4.2). When doing so, however, it is important to keep in mind that the nucleon cross sections should be evaluated at energies that are relevant for inelastic scattering, not in the highly non-relativistic limit. At these high energies, isospin symmetry is typically largely restored because the nucleon couplings are no longer exclusively determined by the valence quarks, and instead receive corrections from a large number of sea quarks (and, in principle, gluons). As pointed out above, the example of neutrino scattering illustrates this effect very clearly: even though isospin is almost maximally violated at low energies, the effective neutrino couplings to neutrons and protons agree within ∼ 5 % at energies around 0.1 GeV, cf. Eqs. (4.5) and (4.6). In practice, however, a possible complication often arises in that the nucleon couplings g n and g p are only provided in the highly non-relativistic limit. In that case, an educated guess for σ SI in the second term of Eq. (4.2) is to anyway take the leading order (Born) expression -but to adopt (effective) values for both nucleon couplings that correspond to Figure 3. Comparison between the elastic (green, lower energies) and inelastic (blue, higher energies) contributions to the DM-nucleus differential cross section dσ χN /dω, where ω is the DM energy loss. This figure shows these contributions for a constant isospin-conserving DM-nucleus cross section, with m χ = 1 GeV and N = 16 O. The small colorbar on the inset of the plots, along with the stated numerical ratio, indicates the balance between elastic and inelastic scattering in terms of the contribution to the integrated cross section σ tot χN .
the maximum of |g p | and |g n | in the non-relativistic limit. This induces a model-dependent uncertainty in the normalization of the inelastic contribution that can in principle only be avoided by fully implementing the concrete interaction model in a code like GiBUU. On the other hand, the neutrino example illustrates that this error should generally not be expected to be larger than a factor of ∼ 2, implying that for most applications such a more sophisticated treatment is not warranted.

Contact interactions and beyond
In sections 3 and 4 we have discussed in detail the Q 2 -dependence that arises due to both form factor suppression and inelastic scattering, as well as the impact this has on the production and attenuation of the CRDM flux. This does not yet take into account, however, the possible angular and energy dependence of the elastic scattering cross section itself. In fact, for (sub-)GeV DM, a significant dependence of this type is actually expected in view of null searches for new light particles at colliders. For example, it has been demonstrated in a recent global analysis [99] that it is impossible to satisfy all relevant constraints simultaneously (even well above GeV DM masses) and at the same time maintain the validity of an effective field theory description at LHC energies. Of course, this necessarily introduces a model-dependent element to the discussion, and in this section, the aim will be to analyse the most generic situations that can appear when considering models beyond simple contact interactions. Concretely, in section 5.2 we will study the case of a light scalar mediator, a light vector mediator in section 5.3, and the scenario where DM particles have a finite extent in section 5.4. In all these cases, we will re-interpret the published Xenon-1T limits and assess whether there is a remaining unconstrained window of large scattering cross sections for GeV-scale DM. Just before this, however, in section 5.1 we will briefly revisit the (physically less motivated) case of a constant cross section, which can be viewed as the highly non-relativistic limit of a contact interaction. This will allow us to illustrate how the resulting CRDM constraints compare with established bounds from both surface and astrophysical experiments, as well as provide a more direct comparison with the existing literature.

Constant cross section
For the discussion of a constant cross section, we will again consider the case of spinindependent scattering with isospin conserving nucleon couplings, cf. Eq. (3.2). In the left panel of Fig. 4, we show our improved constraints from a re-interpretation of the Xenon-1T limits in this case. Broadly, these updated and refined CRDM limits cover the mass range up to m χ 10 GeV for cross sections 10 −31 cm 2 σ SI 2 × 10 −28 cm 2 .
For comparison, we also indicate (with dash-dotted lines) the limits that result when neglecting both form-factor dependence of the cross section and inelastic scatterings in the attenuation part. As expected, this leads to a shape of the excluded region very similar to that originally derived in Ref. [22], where the same simplifying assumptions were made. As a result of our improved treatment of CR fluxes and form factors, however, the limits indicated with dash-dotted lines are overall slightly more stringent than what is reported in that analysis. We find that for very light DM, with m χ 10 MeV, this simplistic treatment actually leads to rather realistic limits, the reason being that for highly relativistic particles the typical momentum transfer is always so large that efficient inelastic scattering becomes relevant. For heavier DM masses, on the other hand, this treatment clearly overestimates the stopping power because it neglects the form factor suppression relevant for semi-relativistic DM scattering on nuclei. Limits on a constant spin-independent DM-nucleon scattering cross section as a function of the DM mass, based on a re-interpretation of Xenon-1T limits on non-relativistic DM [5] for the CRDM component studied in this work (solid lines). Dash-dotted lines show the excluded region that results when assuming a constant cross section in the attenuation part (as in Ref. [22]). Dashed lines show the effects of adding form factors in the attenuation part, but no inelastic scattering, resulting in limits similar to those derived in Ref. [63]. For the latter case, for comparison, we also show the effect of artificially cutting the incoming CRDM flux at the indicated energies. Right panel. Updated CRDM limits (coinciding with the solid lines from the left panel) in comparison to limits from the Lyman-α forest [100], the Milky Way satellite population [101], gas clouds in the Galactic Centre region [102], the XQC experiment [76,103], and a recently analysed storage dewar experiment [104,105]. We also show upper limits on the cross section as published by the CRESST collaboration [6] (solid green lines), based on a surface run of their experiment, along with the maximal cross section where attenuation does not prevent DM from leaving a signal in the detector [16]. Alternative limits are indicated by green dashed [76] and dash-dotted lines [106], based on the assumption of a thermalization efficiency of th = 2 % and th = 1 %, respectively, which is significantly worse than the one adopted in the CRESST analysis.
Dashed lines furthermore show the effect of adding the form factor suppression during the attenuation in the soil, as done in Ref. [63], but still not including inelastic scattering. Clearly, this vastly underestimates the actual attenuation taking place and therefore appears to exclude very large cross sections. 5 In order to gain a better intuitive understanding for the shape and strength of our final limits, finally, we also indicate the effect of neglecting inelastic scattering and instead artificially cutting the CRDM flux (prior to entering the soil) above some given energy. The resulting upper limit on the cross section that can be probed in this fiducial setup strongly suggests that inelastic scattering events very efficiently stop the incident CRDM flux in the overburden as soon as they become relevant compared to elastic scattering events. From Fig. 4, and well in accordance with the expectations from section 4, this happens at CRDM energies T χ 0.2 GeV.
In the right panel of Fig. 4 we show our improved constraints from a re-interpretation of the Xenon-1T limits in comparison with complementary limits from direct probes of the DM-nucleon scattering cross section. At small DM masses the dominant constraint results from analysing the distribution of large-scale structures as traced by the Lyman-α forest. This is based on the fact that protons scattering too strongly off DM would accelerate the latter and thereby suppress the matter power spectrum at sub-Mpc scales. Such limits have recently been significantly tightened [100], utilizing state-of-the-art cosmological hydrodynamical simulations of the intergalactic medium at redshifts 2 z 6. Similar bounds from the CMB (not shown here) are generally weaker by up to three orders of magnitude [100,107,108], while the Milky Way satellite population [101] -as inferred from the Dark Energy Survey and PanSTARRS-1 [109] -places bounds that are roughly one order of magnitude weaker. Beyond cosmological bounds, cold gas clouds near the Galactic Center provide an interesting complementary testbed, in particular at high DM masses, where halo DM particles scattering too efficiently on the much colder baryon population would heat up the latter [110]. Here we show updated constraints [102] based on the cloud G357.8-4.7-55, noting that these constraints might be improved by more than one order of magnitude if G1.4-1.8+87 is indeed as cold as T ≤ 22 K (as reported in Refs. [111,112] but disputed in Ref. [113]). We also display the limits [76] that result from the ten minutes' flight of the X-ray Calorimetry Rocket (XQC) [103], based on the observation that ambient DM particles scattering off the silicon nuclei in the quantum calorimeter would deposit (part of) their energy in the process [14,114,115]. In deriving these XQC limits, one must take into account that the recoil energy of a silicon nucleus potentially thermalizes much less efficiently in the calorimeter than the e ± pairs produced from an incoming X-ray photon, such that a nuclear recoil energy T N will leave a signal equivalent to a photon with a reduced 'thermal' recoil energy T T = th T N . Concretely, the limits shown in the plot are based on the very conservative assumption of a thermalization efficiency factor of th = 0.02. 6 Furthermore, in order to directly probe sub-GeV DM with very large cross sections, the CRESST collaboration has performed a dedicated surface run of their experiment [6], deliberately avoiding the shielding of the Gran Sasso rock used in the standard run [116]. The result of this search is the exclusion region indicated by the solid green line in Fig. 4. Here, upper bounds on the cross section correspond to the published limits, obtained under the assumption that any attenuation in the overburden can be neglected. Modelling the effect of attenuation with detailed numerical simulations also results in the exclusion 6 When the scattering is mediated by a Yukawa-like interaction, a perturbative description of the scattering process may no longer be adequate. In that case the constraints shown here, in particular for XQC, receive corrections due to non-perturbative effects leading to resonances or anti-resonances in the scattering cross section [106]. Here, we will not consider this possibility further, noting that a variation of the relatively uncertain value of th anyway has a larger impact on the XQC limits [76]. region limited from above [16], coming from the fact that one must have a sufficiently large flux of DM particles at the detector location. In a series of papers, Farrar et al. have claimed that the CRESST thermalization efficiency adopted in the official analysis is too optimistic [76,105,106,117], challenging the general ability of the experiment to probe sub-GeV DM. We indicate the resulting alternatives to the published CRESST limits in the same figure, albeit noting that the underlying assumption of an efficiency as low as th ∼ 1 % is not supported by data or simulations. For example, no indication for such a dramatic loss of efficiency at low energies is observed for neutrons from an AmBe neutron calibration source [118].
To summarise, Fig. 4 illustrates the fact that the existence of the CRDM component provides an important probe of strongly interacting light DM. In particular, below m χ 100 MeV, it restricts parameter space that is otherwise either unconstrained or only testable with cosmological probes (which -at least to some degree -are subject to modelling caveats regarding the Lyman-α forest and the non-linear evolution of density perturbations at small scales; see, e.g., Refs. [119,120]). The CRDM component also leads to highly relevant complementary constraints up to DM masses of a few GeV, especially when noting that these constraints are independent of the thermalization efficiency discussion above.

Scalar mediators
As our first example beyond a constant scattering cross section we consider the case where a new light scalar particle φ mediates the interaction between DM and nucleons. We thus consider the interaction Lagrangian L int = −g χ φχχ − g p φpp − g n φnn , (5.1) and assume, for simplicity, isospin conservation (g p = g n ). At the level of the effective nuclear interaction Lagrangian, the dominant interaction terms with scalar (N 0 ) and fermionic (N 1/2 ) nuclei are thus given by 7 Here, the dimensionful coupling to scalar nuclei has been normalized such that both terms in the above expression result in the same scattering cross section in the highly nonrelativistic limit. In addition, the coupling to individual nucleons is coherently enhanced across the nucleus, resulting in an effective coupling to both scalar and fermionic nuclei given by  where G N is the same form-factor as in the case of a 'constant' cross section. For the resulting elastic scattering cross section for DM incident on nuclei at rest we find (5.4) where µ χp is the reduced mass of the DM/nucleon system and is the spin-independent scattering cross section per nucleon in the ultra non-relativistic limit. For reference, the kinematic quantities T max N , s and Q 2 are given by Eqs. (2.5), (2.6) and (2.8), respectively. For the production part of the process, where CR nuclei collide with DM at rest, one simply has to exchange T N ↔ T χ and m χ ↔ m N in these expressions for kinematic variables -but not in the rest of Eq. (5.4) -in order to obtain dσ χN /dT χ .
In the left panel of Fig. 5 we show the resulting CRDM fluxes for this model. For small kinetic energies these fluxes are, as expected, identical to those shown in Fig. 1 for the case of a constant cross section. This is the regime where Q 2 = 2m χ T χ is smaller than the masses of both the mediator and CR nuclei, such that Eq. (5.4) reduces to Eq. (3.2). For Q 2 m 2 φ , on the other hand, the presence of a light mediator clearly suppresses the fluxes.
Note that the matrix element also contains a factor of (Q 2 + 4m 2 χ ), which additionally leads to a flux enhancement for fully relativistic DM particles, T χ 2m χ . In the figure, this latter effect is clearly visible for the case of m χ = 10 MeV and a heavy mediator. In general, the appearance of such model-dependent features demonstrates the need to use the full matrix element for the relativistic cross section. This is in contrast to the non-relativistic case, where a model-independent rescaling of the cross section by a factor of (1 + Q 2 /m 2 φ ) −2 is usually sufficient to model the effect of a light mediator (see, e.g., Refs. [121][122][123]).
In the right panel of Fig. 5, we explore the minimal CRDM energy T χ that is needed to induce a detectable nuclear recoil. Compared to the situation of a constant scattering cross section (depicted by the solid lines for easy comparison), the attenuation is as expected rather strongly suppressed when light scalar mediators are present (with the exception of the case with m χ = 10 MeV and m φ = 100 MeV, where the cross section is enhanced due to the (Q 2 +4m 2 χ ) factor in the squared matrix element). In order to understand the qualitative behaviour of T min χ (z = 0) better, we recall from the discussion of Fig. 2 that there are two generic scaling regimes for solutions of the energy loss equation. Firstly, for cross sections with no -or only a mild -dependence on the momentum transfer, T min χ (z = 0) grows exponentially with increasing σ NR SI . Secondly, in the presence of an effective cutoff in the cross section (like when form factors or light mediators are introduced), T min χ (z = 0) ∝ σ NR SI for large energies T χ . These different regimes are clearly visible in the figure. For the green dot-dashed curve (m χ = 1 GeV, m φ = 100 MeV), for example, one observes as expected an initial steep rise at the smallest DM energies -until the form factor and mediator suppression of the cross section cause a scaling with σ NR SI for kinetic energies above a few MeV. At roughly T χ 0.1 GeV, inelastic scattering kicks in, leading again to an exponential suppression of the flux. For even higher energies, finally, the scattering cross section falls off so rapidly that the required initial DM energy once again only grows as σ NR SI . Turning our attention to the resulting CRDM limits, it is worth stressing here that σ NR SI , as introduced in Eq. (5.5), is a somewhat artificial object that only describes the cross section for physical processes restricted to Q 2 m 2 φ . In a direct detection experiment like Xenon-1T this is necessarily violated for m φ 2m N T thr N ∼ 35 MeV, given that T thr N = 4.9 keV is the minimal recoil energy needed to generate a signal. A natural consequence of this is that making a straight-forward comparison to the σ SI appearing in the 'constant cross section' case discussed in section 5.1 is challenging. Instead, the best we can achieve in terms of a meaningful comparison is to define a reference cross sectioñ  2), and the fact that s ≈ (m χ + m N ) 2 for the energies of interest here, thatσ p Xe,SI should be interpreted as the effective CRDM cross section per nucleon that is dominantly seen in the Xenon-1T analysis window. It is thus this quantity, not the σ NR SI from Eq. (5.5), that should be compared to the published Xenon-1T limits on the DM-nucleon cross section.
This also allows us to address the question of how the limits on the DM-nucleon coupling coming from the CRDM component compare to the complementary constraints introduced in section 5.1 (cf. the right panel of Fig. 4). In order to do so, one first needs to realize that all of those limits are derived under the assumption of non-relativistic DM and a constant cross section. In reality, however, they probe very different physical environments and typical momentum transfers. In order to allow for a direct comparison, therefore, they also need to be re-scaled to a common reference cross section. Assuming that the DM energies in Eq. (5.4) are non-relativistic, a reported limit on the DM-nucleon cross section σ p SI from an experiment probing typical momentum transfers of the order Q 2 ref would correspond to a cross section of in the Xenon-1T detector. As an example, consider the CRESST surface run [6], where a threshold energy of ∼ 20 eV for the sapphire detector would imply Q 2 ref ∼ (0.98 MeV) 2 / th . Similarly, a thermal recoil energy of of 29 eV in XQC corresponds to Q 2 ref ∼ (8.7 MeV) 2 for the nuclear recoil on Si nuclei (assuming th = 0.02 as for the unscaled limits). Turning to cosmological limits, a baryon velocity of v rms b ∼ 33km/s at the times relevant for the emission of Lyman-α photons [124] implies typical momentum transfers from the Helium nuclei to DM of Q 2 ref ∼ 4µ 2 χHe × 10 −8 . This means that, for the range of DM and mediator masses considered here, the cross section at these times becomes roughly constant and we can approximate Q 2 ref ≈ 0 in Eq. (5.7). The same goes for the constraints stemming from the MW satellite abundance, which are sensitive to even lower redshifts and thus smaller momentum transfers [101,125].
In Fig. 6 we show a subset of these correspondingly rescaled constraints 8 -for mediator masses m φ = 1 MeV, 10 MeV, 100 MeV and 1 GeV -along with the full CRDM constraints derived here. We also indicate, for comparison, with dotted black lines where non-perturbative couplings would be needed in this model to realize the stated cross section. This line is only visible for the case of m φ = 1 GeV, which demonstrates that it is generically challenging to realize large cross sections without invoking light mediators. The presence of an abundant species with a mass below a few MeV, furthermore, would affect how light elements are produced during big bang nucleosynthesis (BBN). For a 1 MeV particle with one degree of freedom, like φ, this can be formulated as a constraint of τ > 0.43 s [130] on the lifetime of such a particle. Physically, this constraint derives 8 Upper bounds on the excluded cross section, due to attenuation effects, cannot simply be rescaled as in Eq. (5.7). For the sake of Fig. 6, we instead adopt a rather simplistic approach [16,[126][127][128] to estimate these limits by requiring that the most energetic halo DM particles, with an assumed velocity vmax, can trigger nuclear recoils above the CRESST threshold of 19.7 eV/ th after attenuation in the Earth's atmosphere. For the average density and distribution of elements in the atmosphere, we follow Ref. [129]. By treating vmax and the effective height of the atmosphere, ha, as free parameters, we can then rather accurately fit the results of more detailed calculations [16,76] for the case of a constant cross section -with numerical values in reasonable agreement with the physical expectation in such a heuristic approach. Finally, we adopt those values of vmax and ha to derive the corresponding limits for the case of a scalar mediator, as displayed in Fig. 6. from freeze-in production of φ via the inverse decay process. Since φ → γγ (apart from φ →νν) is the only kinematically possible SM decay channel, the translation of this bound to a constraint on the SM coupling g p is somewhat model-dependent. For concreteness we consider the Higgs portal model, where τ > 1 s at m φ = 1 MeV corresponds to a squared mixing angle sin 2 θ = (8.62 × 10 2 g p ) 2 > 3.8 × 10 −4 [131]. The area above the dashed line in the top left panel of Fig. 6 requires either a larger value of g p than what is given by this bound, or a non-perturbative coupling g 2 χ > 4π. This confirms the generic expectation that for very light particles BBN constraints are more stringent than those stemming from the CRDM component [46,132].
Our results demonstrate that in the presence of light mediators the largest DM mass that can be constrained due to CR upscattering is reduced from about 10 GeV, cf. Fig. 4, to just above 1 GeV (for m φ ∼ 1 MeV). This is a direct consequence of the suppressed CRDM production rate discussed above. On the other hand, the reduction of the cross section also implies a smaller attenuation effect, thus closing parameter space at larger cross sections. More importantly, complementary constraints from cosmology and dedicated surface experiments become more stringent in the presence of light mediators, once they are translated to a common reference cross section. To put this in context, let us first recall that in the constant cross section case, Fig. 4 tells us that cross sections σ SI 2 · 10 −31 cm 2 are safely excluded across the entire DM mass range (or σ SI 6 · 10 −31 cm 2 when assuming that the thermalization efficiency of CRESST is indeed as low as 2 %). From Fig. 6 we infer that these limits can be somewhat weakened for sub-GeV DM, when considering light meditators in the mass range 10 MeV m φ 100 MeV (as we will see further down, the situation of a vector mediator is not appreciably different from that of the scalar mediator shown here). Concretely, the upper bound on the cross section now becomesσ SI 3 · 10 −31 cm 2 , independently of the DM and mediator mass. For a 2 % thermalization efficiency of CRESST [76] and a narrow range of mediator masses, 10 MeV m φ 100 MeV, a small window opens up above the maximal cross section that can be probed with CRESST. The reason is the gap between Lyman-α bounds and the weakened CRESST limits from Ref. [76] that is visible in the figure, for m φ 10 MeV, and which is closed by the CRDM limits only for mediator masses of m φ 30 MeV. Nominally, for m χ ∼ 2 GeV and m φ ∼ 30 MeV, this would allow for cross sections as large asσ SI ∼ 4 · 10 −29 cm 2 . In either case, the conclusion remains that CRDM leads to highly complementary limits, and that this relativistic component of the DM flux is in fact crucial for excluding the possibility of very large DM-nucleon interactions.

Vector mediators
We next consider the general case of a massive vector mediator V , with interactions given by L = V µ (g χ χγ µ χ + g p pγ µ p + g n nγ µ n) . (5.8) We will again assume g n = g p for simplicity, noting that smaller values of the ratio g n /g p can lead to significantly smaller cross sections (see, e.g., Refs. [123,133]); in our context this would mostly imply that the attenuation in the overburden becomes less relevant, Solid purple lines show the updated CRDM limits studied in this work. We further show limits from the Lyman-α forest [100], the XQC experiment [76,103], the CRESST surface run [6,16] and an alternative analysis of the CRESST limits [76]. All these limits are rescaled to match the situation of a light mediator, as explained in the text. The parameter region above the dotted black line in the bottom right panel requires non-perturbative couplings, while the area above the dotted line in the top left panel is excluded by BBN.
leading to more stringent constraints. In analogy to Eq. (5.2), this implies the following dominant interaction terms with scalar and fermionic nuclei, respectively: where the effective mediator coupling to nuclei, g N , is again given by the coherent enhancement stated in Eq. (5.3). For the elastic scattering cross section on nuclei we find Here, the cross section in the ultra-nonrelativistic limit, i.e. for Q 2 → 0 and s → (m N + m χ ) 2 , agrees exactly with the result obtained for the scalar case, as expected. For large energies and momentum transfers, on the other hand, the behaviour is different. The resulting CRDM fluxes are nonetheless so similar to the scalar case shown in the left panel of Fig. 5 that we refrain from plotting them separately. Differences do exist, however, for the stopping power in the overburden. In the left panel of Fig. 7 we therefore show the minimal initial kinetic energy needed by a CRDM particle to induce detectable nuclear recoils in Xenon-1T. Compared to the scalar case, cf. the right panel of Fig. 5, the attenuation is more efficient for highly relativistic DM particles due to the s-dependence of the terms in the second line of Eq. (5.10). As before, the effect of these model-dependent terms from the scattering amplitude is most visible for highly relativistic particles, with small m χ , and large mediator masses, where the suppression due to the factor (1 + Q 2 /m 2 A ) −2 is less significant. In the right panel of Fig. 7 we compare the final exclusion regions for the situations considered so far, i.e. for a contact interaction, scalar mediators and vector mediators, respectively. For the sake of comparison in one single figure, we plot here the cross section in the ultra-nonrelativistic limit. For an interpretation of these limits in comparison to complementary constraints on DM-nucleon interactions we thus refer to the discussion of Fig. 6, noting that the rescaling prescriptions for vector and scalar mediators are qualitatively the same. The first thing to take away from Fig. 7 is that, as expected, the exclusion regions for heavy mediators resemble those obtained for the constant cross section case. The figure further demonstrates that the only significant difference between scalar and vector mediators appears at smaller mediator masses, where the former are somewhat less efficiently stopped in the overburden. It is worth noting, however, that this region of parameter space where the vector and scalar cases differ substantially is nonetheless excluded by Lyman-α bounds. The general discussion and conclusions from the scalar mediator case explored in the previous subsection thus also applies to interactions mediated by vector particles.

Finite-size dark matter
As a final generic example of a Q 2 -suppressed cross section let us consider the situation where the DM particle itself has a finite size that is larger than its Compton wavelength. Various models of such composite DM have been extensively studied in the literature [134][135][136][137][138][139][140][141]. In fact, Ref. [142] even suggests that DM with masses above 1 GeV cannot be pointlike for DM-nucleon cross section 10 −25 cm 2 . The corresponding scattering cross section then takes the same form as in the point-like case, multiplied by another factor G χ (Q 2 ) 2 that reflects the spatial extent of χ [143][144][145]. Specifically, just as for nuclear form factors, we have where ρ χ (x) is the distribution of the effective charge density that the interaction couples to. For simplicity we will choose a dipole form factor of the form 9 G χ (Q 2 ) = 1 + r 2 with r χ being the r.m.s. radius of the DM particle, r 2 χ = d 3 x x 2 ρ χ (x). We then multiply G 2 χ (Q 2 ) with Eq. (3.2) in order to obtain dσ χN /dT N , thus describing an effective scalar 9 The exact choice of the form factor does not significantly affect our results, as long as Gχ(Q 2 ) < Gχ(0) = 1. An interesting, qualitatively different situation occurs when Gχ(0) = 0, i.e. for a form factor that grows with Q 2 . This is, e.g., realized if the scattering is mediated by a dark U (1) under which χ is neutral [143,145]. We will not consider this class of models in this work. interaction with the usual coherent enhancement inside the nucleus -but where each of the nucleons only 'sees' some fraction of the entire DM particle.
In a very similar fashion to what happens in the presence of a light mediator φ, such a cross section features a sharp suppression for momentum transfers exceeding a 'mass' scale m φ ∼ √ 12/r χ . Sharper than in that case, in fact, as the suppression scales with a power of Q −8 rather than just Q −4 . This is clearly visible in the left panel of Fig. 8, where we plot the expected CRDM flux for DM with a finite size, for various values of m χ and r χ . For example, for r χ = 10 fm, we have √ 12/r χ ∼ 68 MeV and the cutoff indeed appears at only slightly smaller values of T χ than in the case of the 100 MeV mediator displayed in Fig. 5 (for m χ = 1 GeV). The slope above the cutoff, however, is twice as steep -as expected from the Q −8 suppression.
In the right panel of Fig. 8 we show how the constraints on a constant DM-nucleon cross section weaken when considering the situation where the DM particles themselves have a finite extent. Concretely, for a DM radius of r χ = 1 fm (r χ = 10 fm) the maximal DM mass that can be probed decreases from ∼ 10 GeV to about 4.5 GeV (1.1 GeV). The reduced CRDM flux for extended DM, cf. the left panel of the figure, also visibly weakens the lower bound on the exclusion region. At the same time, attenuation is also less efficient for a given cross section in the non-relativistic limit (inelastic scattering still effectively cuts off the incoming CRDM flux above ∼0.2 GeV, explaining e.g. the upper, almost horizontal boundary of the exclusion region in the r χ = 10 fm case). For r χ 1 fm, this starts to significantly enlarge the excluded region to higher cross sections. On the other hand, it should be noted that for composite DM particles the interaction cross section may not actually continue to drop as Q −8 for very large momentum transfers, as would be implied by Eq. (5.13). At some point, instead, inelastic scattering events on the DM constituents will take over, in analogy to what we discussed for nuclei in section 4. This is particularly relevant if the DM constituents are themselves finite in size, in which case the upper boundaries of the exclusion regions shown in Fig. 8 would be overly optimistic for very large r χ .
Similar to the discussion in section 5.2, a comparison of the limits shown in Fig. 8 with complementary limits requires a re-scaling of σ SI to a common reference cross section. Due to the strong form factor suppression, this rescaling has an even larger effect than in the light mediator case; concretely, instead of Eq. (5.7), the rescaling of reported limits, σ p SI , to those relevant for the Xenon-1T detector now takes the form (5.14) Qualitatively, however, this does not change the lesson learned in the light mediator case: while limits from the CRDM component can be weakened by increasing r χ , this will inevitably strengthen complementary bounds from cosmology. As a result, we find once again an absolute upper bound on the cross section of aboutσ SI ∼ 3 · 10 −31 cm 2 , independently of the DM mass and size. Also in this case there is a small loophole to this statement if one is willing to assume that the thermalization efficiency of CRESST is as small as 2 %: when tuning the size of the DM particles to r χ 10 fm, we find that cross sections two orders of magnitude larger may in that case be viable for DM masses in a narrow range between around 1 GeV and 2 GeV.
In section 5 we discussed various generic situations where the amplitude for elastic scattering shows a significant dependence on the momentum transfer, and how this impacts the conclusions about whether a window of large scattering cross sections remains open or not. In this section we complement those more model-independent considerations by taking a closer look at a specific DM candidate in the GeV range, with relatively large nuclear interactions. Concretely, it has been conjectured that a neutral (color-flavor-spin-singlet) bound state of six light quarks uuddss may exist, and provide a plausible DM candidate that would evade all current constraints despite its baryonic nature [14,[146][147][148][149]. In particular, this sexaquark S (to be distinguished from a generic 6-quark state, often referred to as hexaquark) would form early enough to behave like standard cold DM during both big bang nucleosynthesis and recombination. It would thus not be in conflict with the independent, and rather precise, measurements [150,151] of the cosmological baryon density during these epochs.
Compared to the H-dibaryon that was suggested earlier [152] and thoroughly studied both theoretically and experimentally (see Refs. [153,154] for reviews), furthermore, the S should be much more tightly bound, leading to weaker interactions with ordinary baryons and thus evading direct searches. Such a particle would be absolutely stable for m S < m D + m e 1.88 GeV, and decay with a lifetime exceeding the age of the Universe for m S 2 GeV [148]. Determining its expected mass exactly, however, is challenging; lattice simulations, for example, remain somewhat inconclusive (see, e.g., Refs. [155][156][157][158] where the results for binding energies of the H-dibaryon state range from ∼17 MeV to ∼75 MeV relying, however, on unrealistically large quark masses). Even if the sexaquark is stable on cosmological timescales, its relic abundance would generally be much smaller than the observed DM abundance if one assumes that its interactions in the early universe are of the order of the strong force [159,160]. If instead, one postulates much weaker interactions due to the assumed compactness of the sexaquark, thermal equilibrium with the SM heat bath would not be possible to maintain after the QCD phase transition and the correct DM abundance might be achieved -in a region of parameter space claimed to evade all existing constraints [148].
Motivated by this intriguing possibility, for simplicity we will adopt the description of sexaquark interactions from Ref. [148], i.e. we model the interaction with nucleons by the exchange of a vector meson. In particular, the relevant interaction terms with the flavour-neutral mixture of φ and ω, denoted by V , are given by and we adopt the value m V = 1 GeV used in Ref. [148] for our calculations. The value of g n = g p ∼ 2.6 √ 4π can be inferred from the literature on the one-boson-exchange model [161] although O(1) uncertainties can be expected here. 10 The coupling g S is largely unknown, though simple scaling arguments suggest that is very roughly of the order of ∼ 0.1 [148]. Following that reference, we will treat α SN as a free parameter that we will generously vary in the interval (10 −3 , 10). Importantly however -at least in this parameter range -the DM relic abundance is independent of α SN . Instead, the final abundance of S is set by an independent coupling constantg [148] that describes the (much weaker) sexaquark-breaking interactions within the effective description. This coupling does not directly enter the analysis presented here.
We treat the interaction of V with nuclei similarly to that in section 5.3, i.e. we describe it by the effective Lagrangian (5.9) with the coherently enhanced, effective coupling g N 10 In particular, we note that modern analyses of low-energy baryon-baryon scattering consider processes beyond single meson exchange [162], and that baryon-baryon interactions can also be treated within the more systematic approach of chiral perturbation theory [163]. However, given the significant uncertainties on the sexaquark couplings we consider the one-boson-exchange approximation to be sufficient for our purposes.
given by Eq. (5.3). For the elastic scattering cross section on nuclei we thus find Here, is the scattering cross section on nucleons in the non-relativistic limit and µ Sp (µ SN ) is the reduced mass of the sexaquark-nucleon (nucleus) system. Compared to the treatment in section 5.3, we introduce an additional form factor G V related to the cutoff in the one-boson-exchange models. In this context, exponential cutoffs are mostly used and the cutoff mass Λ V is fitted to data (and can in principle differ for different meson exchange channels). For example, within the fit to data taking into account hyperon-nucleon interactions [161], these cutoff masses were found to range between 820 MeV and 1270 MeV. Since yet lower cutoff masses appear in related literature (e.g., down to 590 MeV in [164]), we generously vary Λ V between 500 and 1500 MeV. We note that for Λ V 1500 MeV, CRDM limits become in fact independent of the cutoff scale. In Fig. 9 we show the parameter space in the α SN vs. m S plane where sexaquark DM is excluded because of the irreducible CRDM component. For a better direct comparison, we also indicate the preferred mass range according to Ref. [148], along with the complementary limits presented in that analysis. From this figure, it is clear that our new limits close a significant part of the viable parameter region where sexaquarks could be the dominant DM component -even without taking into account the CRESST results. In particular, we note that the Lyman-α limits [100] shown in figures 4 and 6 were presented subsequent to the analysis of Ref. [148] and are significantly stronger than the CMB limits indicated in Fig. 9. The apparently open window at α SN ∼ 0.3 is thus also robustly excluded. On the other hand, a small open window remains for α SN 4 · 10 −3 . While not being in conflict with the DM abundance, as explained above, we recall that such values of α SN are somewhat smaller than intrinsically expected.
Let us, finally, briefly comment on the fact that the DM-nucleon scattering cross section can, strictly speaking, only be calculated perturbatively in the Born limit, α SN µ χN m V . Outside this regime, non-relativistic scattering in a Yukawa potential exhibits parametric resonances where the scattering amplitude is significantly enhanced or suppressed. This non-perturbative effect is well-known from the self-scattering of cold DM in the presence of light mediators [167], and it is the origin of the resonant structure in the complementary limits from Ref. [148] that is visible in Fig. 9. For our CRDM limits, on the other hand, this additional complication does not arise because such non-perturbative corrections are Figure 9. Effective sexaquark coupling α SN vs. sexaquark mass m S . The purple region shows the parameter range that is excluded by the analysis in this work, assuming that sexaquarks make up all of the cosmologically observed DM; different line styles correspond, as indicated, to cutoff masses Λ V /GeV ∈ {0.5, 1, 1.5} in the one-boson exchange approximation. All other constraints are, for easier comparison, directly reproduced from Fig. 10 of Ref. [148], conservatively assuming an attractive Yukawa force between S and nuclei. The thin vertical stripe corresponds to the mass range where, according to that analysis, the sexaquark would be a viable DM candidate without being in conflict with other particle physics observation, in particular the stability of deuterons based on SNO data [165]. The upper end of that mass range may increase from 1890 MeV to up to 2054 MeV if sexaquark DM does not accumulate in the Earth at the level claimed in Ref. [166]. largely irrelevant for relativistic scattering; in fact, already for the typical velocities during the freeze-out process of thermally produced DM, v χ ∼ 0.3, the impact is strongly suppressed [167]. The CRDM limits are thus also robust w.r.t. underlying model assumptions such as whether the force mediated by the Yukawa potential is attractive or repulsive.

Summary and Conclusions
For sizeable elastic scattering rates between DM and nuclei there is an irreducible relativistic component of the flux of DM particles arriving at Earth. This extends the sensitivity of conventional direct detection experiments both to sub-GeV masses and to scattering cross sections above the limit set by a too efficient attenuation of the DM flux on the way to the detector. While such large scattering cross sections are also constrained by complementary probes from astrophysics and cosmology, it has repeatedly been pointed out that there might be an open window of relatively strongly interacting DM with a mass in the ballpark of ∼ 1 GeV.
We find that the CRDM component in the DM flux generically closes this window, under rather minimal assumptions. In order to arrive at this conclusion, we included in our analysis a detailed treatment of the inelastic scattering of DM off nuclei (section 4). We demonstrate that this provides an important additional stopping channel for CRDM particles on their way to direct detection facilities -unlike for non-relativistic DM, where only elastic scattering is relevant. We also investigated the extent to which a possible energy or momentum-transfer dependence of the cross section could weaken our general conclusions. For this purpose, we considered i) a class of simplified models where the scattering with nuclei is mediated by a light scalar (section 5.2) or vector (section 5.3) particle, as well as ii) situations where DM particles cannot be described as being point-like (section 5.4). In all these cases, the additional momentum-transfer dependence indeed weakens the limits from direct detection -which however is compensated for by a corresponding strengthening of complementary limits, in particular from cosmology. In combination, these limits stringently constrain the possibility of cross sections larger than a few times 10 −31 cm 2 , over a wide range of DM masses. Interestingly, this is largely independent of underlying modelling assumptions such as the mass of new mediator particles or the DM particles' radius.
Finally, an exotic QCD bound state that is produced well before BBN, has repeatedly been put forward as a potential DM candidate. While it is theoretically unclear whether such states could actually exist, adding to significant experimental constraints, it is certainly an intriguing idea to have a 'baryonic' DM candidate that would in fact evade the strong evidence from BBN and CMB against this possibility. However, cosmic-ray upscattering of such particles leads to stringent new constraints that have not previously been pointed out in this context. For the concrete case of stable sexaquark DM, as discussed in section 6, we find that the parameter space giving the correct cosmological abundance is strongly pressured.
For the analysis performed in this work we used the numerical tool DarkSUSY [64] to compute CRDM fluxes and limits. In doing so we significantly expanded the general numerical routines provided therein, adding in particular inelastic scattering, the contribution from CRs beyond p and He, and an updated treatment of nuclear form factors in the context of CRDM attenuation. These updates will be included in the next public release of the code.