Constraining light dark matter upscattered by ultrahigh-energy cosmic rays

Light halo dark matter (DM) particles upscattered by high-energy cosmic rays (CRs) can be energetic, and become detectable by conventional direct detection experiments. The current constraints derived from space-based direct CR measurements can reach $\mathcal{O}(10^{-31})\text{ cm}^{2}$ for a constant DM-nucleon scattering cross section. We show that if the CR energy spectrum follows a power law of type $\sim E^{-3}$, the derived constraints on the scattering cross section will be highly insensitive to DM particle mass. This suggests that ultrahigh-energy CRs (UHECRs) indirectly measured by ground-based detectors can be used to place constraints on ultralight DM particles, as $E^{-3}$ is a very good approximation of the UHECR energy spectrum up to energy $\sim10^{20}\text{ eV}$. Using the recent UHECR flux data, we show that the current constraints derived from space-based CR measurements can in principle be extended to ultralight DM particles far below eV scale.


I. INTRODUCTION
Although compelling astrophysical evidence supports the existence of dark matter (DM) in the Universe, whether or not DM participates non-gravitational interactions is still an open question. The majority of the current DM direct detection (DD) experiments search for nuclear recoil signals from the scatterings between the halo DM particle and target nucleus. As the typical detection threshold of the current experiments is of O(keV), searching for light halo DM below GeV is in general challenging. The reason is that for lighter halo DM particles the kinetic energy is lower, and the energy transferred to the target nuclei is suppressed. For instance, for a DM particle with mass m χ ∼ 1 GeV and a typical DM escape velocity ∼ 540 km s −1 , the elastic scattering off a target nucleus with mass ∼ 100 GeV leads to a maximal recoil energy ∼ 0.06 keV which is significantly lower than the typical detection threshold. Several physical processes have been considered to lower the detection threshold such as using additional photon emission in the inelastic scattering process [1] and the Migdal effect [2, 3], etc.. The same DM-nucleus scattering process may leave imprints in some cosmological and astrophysical observables, which can be used to place constraints on the scattering cross section. The resulting constraints are in general much weaker but can be applied to lower DM particle masses unreachable to the DD experiments. For instance, from the spectral distortion of the cosmic microwave background (CMB), a constraint of σ χp 5 × 10 −27 cm 2 for m χ at 1 keV-TeV can be obtained [4]; the constraints from the population of Milky Way (MW) satellite galaxies can reach σ χp 6 × 10 −30 cm 2 for m χ 10 keV [5]; and the measurement of the gas cooling rate of the Leo T dwarf galaxy can also lead to a constraint of σ χp 3 × 10 −25 cm 2 for m χ 1 MeV [6,7].
Recently, it was shown that important constraints can be derived from the scattering between cosmic ray (CR) particles and DM particles in the local Universe. High-energy CR particles in the Galaxy can scatter off halo DM particles, which results in the energyloss of CRs [8], the production of γ-rays [9,10] and energy-boost of DM particles [11][12][13], etc.. In the last process, a small but irreducible component of DM (referred to as CRDM) can obtain very high kinetic energies. These energetic CRDM particles can scatter again with the target nuclei in underground detectors, and deposit sufficient energy to cross the detection threshold, which greatly extend the sensitivity of the current DD experiments to sub-GeV DM particles [11][12][13][14][15][16][17][18]. It has been shown that in this approach the constraints on constant DM-nucleon (DM-electron) spin-independent scattering cross section can reach ∼ 10 −31 (10 −34 ) cm 2 for DM particle mass down to at least ∼ 0.1 MeV (∼ 1 eV) [11,12] (for constraints on energy-dependent cross sections, see e.g. [14,15]).
It is of interest to explore the potential of this approach in constraining even lighter DM particles far below MeV scale, as some well-motivated DM candidates such as QCD axions and axion-like particles can be extremely light. Constraining lighter DM requires better information on the CR spectra at higher energies. Note, however, that all the current analysis [11][12][13][14][15][16][17][18] on the detection of CRDM adopted the CR fluxes from either the parametrizations in [19,20] or the GALPROP code [21,22], which are inferred from the space-based direct CR measurements (e.g. PAMELA [23,24], AMS-02 [25,26] and CREAM-I [27] etc.). For current space-based experiments the CR fluxes which can be measured with reasonable precision are typically with energy 200 TeV (see also the data of CREAM-III [28], CALET [29] and DAMPE [30]). Towards higher energies, the statistic uncertainties increase rapidly due to the limited acceptance of space-based experiments [31][32][33]. Naively extrapolating these analyses to higher energies will leads to incorrect conclusions, as the spectral feature of the CR flux start to change above 1 PeV. Alternatively, high-energy CR can be measured indirectly by ground-based air-shower detectors. Despite larger uncertainties in energy scale and mass resolution, this approach can measure the CR flux to much higher energy due to the huge acceptances. For detecting lighter DM particles, the local DM number density is higher. However, the energy transfer from the scattering process become less efficient, and the CR flux is known to decrease rapidly towards higher energies. Whether or not ultrahigh-energy CRs (UHECRs, defined as CR with total energy E > PeV) can be used to place useful constraints on ultralight DM will depend strongly on the spectral feature of the UHECR flux.
In this work, we show that as long as the energy spectrum of CR flux follows a power law ∼ E −α with α 3, the derived constraints on the DM-nucleon scattering cross section will not decrease towards lower DM mass. In the limit of α = 3, the constraints will DM mass independent. This justifies using UHECR to place stringent constraints on ultralight DM particles, as the UHECR all-particle spectrum above the "knee" structure (at ∼ 3 PeV) can be well-approximated by a power law with α ≈ 3. From the recent UHECR nucleus flux data, we obtain the following results: the constraints on the spin-independent DM-nucleon scattering cross section can be 10 −(32−31) cm 2 for DM particle mass down to extremely small value ∼ 10 −12 eV, which expands the currently known constraints derived from spacebased direct CR measurements by around ten orders of magnitude in DM mass, and close a large previously unconstrained parameter space; the most stringent constraints are found to be at DM mass ∼ 10 −5 eV and ∼ 10 −11 eV, due to the "knee" and "toe" structure in the UHECR flux, respectively; this CRDM approach will completely loss sensitivity for DM mass below 10 −14 eV as the UHECR flux is highly suppressed above ∼ 10 20 eV, a phenomena possibly related to the scatterings between UHECRs and cosmic microwave background (CMB) photons [34,35]. The constraints obtained in this work are highly model-independent and conservative, as only the elastic scattering process is required and very conservative choices of parameters are adpoted. The constraints obtained in this work are derived based on the observables of the present-day local Universe, which are complementary to other constraints derived from the data of earlier epochs of the Universe. This paper is organized as follows: In section II, we discuss the spectral feature of the DM flux upscattered by UHECRs. In section III, we discuss the effect of earth attenuation of the CRDM kinetic energy. The nuclear recoil sepectrum and the constraints from direct detection experiment Xenon-1T is discussed in section IV. We summarize the work and give some remarks in section V.

II. CR-UPSCATTERED DARK MATTER FLUX
A. Single CR component case In the generic process of elastic scattering between an incident particle A with kinetic energy T A and a target particle B at rest, the recoil energy of particle B in the laboratory frame is given by T B = T max B (1 − cos θ)/2, where θ is the scattering angle of particle B in the center-of-mass (CM) frame. The maximal recoil energy of particle B is given by where m A(B) is the mass of particle A(B). We assume that the scattering is isotropic in the CM frame, such that the differential cross section dσ AB /dT B in the laboratory frame is simply related to the total cross section σ AB as dσ AB /dT B = σ AB /T max B . In the case of CR-DM scattering, if the CR particle i (i = H, He, . . . ) is highly relativistic, i.e., the Lorentz factor γ i ≈ T i /m i 1, but m χ is small enough such that γ i m i /2m χ , the maximal recoil energy of the CRDM can be approximated as T max χ ≈ 2m χ γ 2 i . The CRDM particle with kinetic energy T χ can scatter again with the nucleus N (with mass m N ) in either the outer crust of Earth or the detector of the underground DM direct detection experiments. The maximal recoil energy T max N of the nucleus which is also the maximal energy-loss of CRDM particle can be well approximated as T max N ≈ 2T 2 χ /m N . Note that T max N is independent of DM particle mass.
After being upscattered, the CRDM particles travel through the Galaxy in straight lines as they are not deflected by the interstellar magnetic fields. The observed flux (number of particles per unit area, time and solid angle, dN/dAdtdΩ) of CRDM at the surface of Earth can be approximately written as where dΦ LIS i /dγ i is the local interstellar CR flux measured at Earth. The integration lower limit γ min i ≈ (T χ /2m χ ) 1/2 is the minimal Lorentz factor required to produce T χ . The form factor F (Q 2 χ ) is evaluated at the momentum transfer Q 2 χ = 2m χ T χ . For very light DM, F (Q 2 χ ) ≈ 1 is an excellent approximation. In the above expression we have assumed that the CR energy spectrum in the Galactic halo is not significantly different from that in the local interstellar (LIS) region, i.e., dΦ i (r)/dγ i ≈ dΦ LIS i /dγ i . In this case, the information of halo DM density distribution can be parameterized into a single parameter D eff where ρ loc χ ≈ 0.3 GeV/cm 3 is the local DM density in the Solar system, and the integration is performed along the line-of-sight (l.o.s). For typical DM profiles, the value of D eff is around a few kpc. In this work, we make a conservative choice of D eff = 1 kpc as a benchmark value.
Let us start with a simple case where the flux of a CR species i can be parametrized by a single power law with index α i and a cutoff at a Lorentz factor γ i,cut as follows where Φ 0 i is a normalization factor. The power-law behavior is expected if CRs are accelerated by the diffusive shock waves of the Galactic supernova-remnants (SNRs) and the pulsar wind, etc., and the cutoff represents the maximal energy that can be achieved by the acceleration process. If m χ is sufficiently small such that γ i,cut m i /2m χ , which is easily justified for sub-eV CRDM, the approximation of T max χ ≈ 2m χ γ 2 i can be used in the whole integration range of Eq. (2), and the corresponding CRDM flux can be obtained analytically as follows where Γ is the incomplete Γ-function, t = (T χ /T max χ,cut ) 1/2 with T max χ,cut = 2m χ γ 2 i,cut the maximal energy of CRDM upscattered by UHECR at the cutoff γ i,cut .
In the region far below the cutoff, i.e., T χ T max χ,cut , which corresponds to the case where the CR flux is essentially a single power law ∼ T −α i i . Using the asymptotic behavior of the incomplete Γ-function Γ(a, z) → −z a /a for z 1, the CRDM flux can be approximated as As F (Q 2 χ ) ≈ 1 is a very good approximation for ultralight DM, the above expression shows that the CRDM flux follows a power law ∼ T −(1+α i )/2 χ . The mass dependence of the CRDM flux is proportional to m (α i −3)/2 χ , which shows that as long as the CR flux is hard enough, namely, α i 3, the resulting CRDM flux will not decrease with decreasing DM mass.
Note that in a wide energy range the CR flux is close to the case of α i ≈ 3. Direct and indirect measurements show that from a few GeV up to the "knee" (at ∼ 3 PeV), the primary CR all-particle spectrum approximately follows a single power law with index α i ≈ 2.7. Above the "knee" the spectrum softens to α i ≈ 3.1. Before reaching the highest observed energy ∼ 10 20 eV, there are several minor spectral structures such as the "second knee" at ∼ 10 17 eV, the "ankle" at ∼ 8 × 10 18 eV and the "toe" at ∼ 3 × 10 19 eV. The corresponding power-law indices vary around the α i ≈ 3 case. Consequently, the DM upscattered by UHECR in this ultrahigh-energy region should fluctuate around the power law T −2 χ , and is highly insensitive to DM mass; and the recoil event rate and the derived bounds on σ χi should be independent of m χ as well, as the recoil energy or the energy loss in the χN scattering is almost independent of DM mass. Eq. (6) also suggests that the CR electrons are less efficient in constraining ultralight DM particles, as the CR electron flux follows a power law with power index α e ≈ 4 after ∼ 0.9 TeV [36][37][38].
In a different region where T χ is close to the cutoff T max χ,cut , the effect of cutoff in CR flux will be significant. Using the asymptotic behavior of Γ(a, z) → z a−1 e −z for large z, the CRDM flux is given by Since T max χ,cut is proportional to m χ , lighter CRDM particles will have an earlier cutoff. A final cutoff in the CR flux around ∼ 10 20 eV is expected from the inelastic scattering between UHECR particles and CMB photons as predicted by Greisen, Zatsepin and Kuzmin [34,35], which is supported by recent observations [39][40][41]. The cutoff in the UHECR flux essentially sets the scale of the minimal m χ that can be constrained by this approach.

B. Multiple CR component case
The primary CR flux in the ultrahigh-energy region may receive contributions from different sources such as SNRs, pulsar winds and active galactic nuclei (AGN), etc. (for recent reviews see e.g. [42][43][44]). The multi-source description is also essential to reproduce the observed various spectral structures of UHECRs. Thus a realistic description of the UHECR flux necessarily contains multiple components, which can be written as Φ i = j Φ ij with j = 1, . . . , n. For each component j, the flux Φ ij takes the form of Eq. (4) with the power index α i and cutoff γ i,cut replaced by α ij and γ ij,cut , respectively. Thus we adopt the following form of the primary CR flux [45] dΦ LIS where Φ 0 ij and α ij are the normalization factors and power indices, respectively, for a CR species i in the component j. Following the reasoning of Peters [46], the CR species in each component j should share a common cutoff in rigidity R j , which leads to γ ij, where Z i is the electric-charge of the CR species i. In Ref. [45] four different parametrization are found to be in good agreement with the current UHECR data [39,[47][48][49][50][51][52][53][54][55][56][57]. We choose one of the "Global-Fit" parametrization with n = 3. The best-fit values of the rigidity cutoffs are R 1,2,3 = 120 TV, 4 PV and 1.3 EV, respectively [45]. Compared with other parametrizations, this one is the most economic and conservative as the final cutoff of R 3 is the lowest, which leads to the lowest CRDM flux at high energy region. The details of the parametrizations are summarized in Appendix-A.  (1) and (2) without using approximations. In the calculation we take the dipole form factors for light CR species H and He [58], and the Helm form factor for heavier species [59,60]. In the energy region where T χ is far below the lowest cut off, since α ≈ 2.7 the CRDM fluxes follow an approximate power law ∼ T −1.85 χ and scale with DM mass as m −0.15 χ which is a rather weak m χ -dependence. Thus lighter CRDM particle have slightly larger flux. In the cutoff dominated region, due to the superposition of various cutoffs γ ij,cut in the three components, the CRDM fluxes fluctuate around the power law case of T −2 χ over many orders of magnitude in kinetic energy before reaching the last cutoff, and are insensitive to m χ , which is expected from Eq. (6) and can be clearly seen in Fig. 1. Above the final cutoff R 3 , the CRDM flux drops rapidly; the flux of lighter CRDM particle drops faster, as expected from Eq. (7). Taking the case of m χ = 10 −10 eV as an example, the CR protons in the three components lead to the induced CRDM flux with three cut off at T max χ,cut = 3.3 × 10 −9 , 3.6 × 10 −6 and 0.38 GeV, respectively, as can be seen from the inset of Fig. 1. If the cutoff is too low, the CRDM cannot be energetic enough to produce enough recoil energy to be detected by the DD experiments. Thus a lower limit on m χ for a given DD experiment exists.

III. EARTH ATTENUATION
Before arriving at the underground detectors, CRDM may loss a non-negligible fraction of energy due to the same elastic scatterings with nucleus within the outer crust of Earth. We adopt an analytic approach for Earth attenuation based on average energy loss [61,62]. For simplicity, we only consider elastic scatterings which is an irreducible process and neglect the effect of form factor. The decrease of T χ with depth z due to the elastic scattering with the nucleus N in Earth's crust is given by where ρ N is the mass density of nucleus N in the crust, and T N stands for the nucleus recoil energy which equals the energy loss of the incident CRDM particle. Using the expression of T max N the energy loss can be approximated as where κ = N ρ N σ χN /m 2 N . We consider the case where the scattering is isospin conserving, namely, σ χn ≈ σ χp such that the cross sections at nucleus and nucleon level are simply related by σ χN ≈ A 2 N σ χp with A N the nucleus mass number of N . We further adopt the relation m N ≈ A N m p which is a very good approximation. Under these two simplifications, the factor A N cancels out in the expression of κ, which gives κ ≈ σ χp ρ ⊕ /m 2 p with ρ ⊕ ≈ 2.7 g·cm −3 the average mass density of Earth's outer crust [63]. After integrating Eq. (9), the CRDM kinetic energy T z χ at depth z is related to that at surface as Thus the effect of energy loss is independent of both m χ and the chemical composition of the crust. For a small enough cross section σ χp m 2 p /(zρ ⊕ T z χ ), one obtains T z χ ≈ T χ . In the opposite case where σ χp is large enough, T z χ will stop tracking T χ , and reach a maximal value T z,max The appearance of T z,max χ is due to the rapid energy loss proportional to T 2 χ for relativistic incident particles, which leads to a sharp cutoff in the CRDM flux at depth z. The CRDM flux dΦ χ /dT z χ at depth z can be evaluated from that at surface dΦ χ /dT χ through the relation dΦ χ /dT z χ = (dΦ χ /dT χ )(dT χ /dT z χ ) .

IV. DIRECT DETECTION
A. Recoil event spectrum The differential nuclear recoil event rate per target nucleus mass Γ = dN/dM N dtdT N of the χN scattering at depth z is given by where Q 2 N = 2m N T N . The scale of T z χ relevant to direct detection is set by the lower limit of the integration T z,min χ (T N ) ≈ (T N m N /2) 1/2 . For a detector located at depth z ∼ 1 km and a target nucleus mass m N ∼ 100 GeV, in order to produce a recoil energy T N which is close to the threshold T thr N ∼ 1 keV, the required minimal T z χ is ∼ 10 MeV. Thus the condition of T z χ ≈ T χ leads to a typical requirement of σ χp O(10 −27 ) cm 2 . Since the lower bounds of the excluded cross section is expected to be much smaller σ χp 10 −31 cm 2 [11], in deriving the lower bound of the excluded σ χp , the effect of Earth attenuation can be safely neglected. Using the CRDM flux from Eq. (5), the recoil event rate can be written as where t = [m N T N /(8m 2 χ γ 4 i,cut )] 1/4 . In the case where the CR flux follows a power law with index α i , the recoil event rate can be obtained analytically from Eqs. (6) and (12) as follows which explicitly shows that if α i 3, the recoil event rate is proportional to T −3/2 N and does not decrease with a decreasing m χ . For the case α i ≈ 3, the derived upper limit on σ χp will be insensitive to m χ . In Fig. 2, we show the full numerical results of recoil event rate of the scattering between the CRDM particles and xenon nuclei. The approximate power-law behavior of the recoil spectrum can be clearly seen for T N keV. Above ∼ 10 keV, the suppression due to the Helm form factor is visible. Due to the power-law like spectrum of the recoil event rate, the experiment with lower threshold T thr N will be more sensitive to lighter DM particles.
For large enough cross sections, the appearance of T z,max χ leads to a cutoff in the recoil event spectrum. If the corresponding recoil energy is below the threshold T thr N , it will form a blind spot for direct detection. This possibility is illustrated in the inset of Fig. 2, which will lead to m χ -independent upper bounds on the excluded value of σ χp .

B. Xenon-1T detector response and data analysis
We numerically calculate excluded regions in the (m χ , σ χp ) plane at 90% C.L. for CRDM from the data of Xenon-1T experiment located at depth z ≈ 1.4 km [64,65]. The Xenon-1T experiment utilizes the liquid xenon time projection chambers to detect the recoil energy of the target nuclei from their scattering with DM particles. The deposited energy can produce a prompt scintillation signal (S1) and ionization electrons which are extracted into gaseous xenon and produce proportional scintillation light (S2). The S2/S1 signal size ratio allows for discrimination between nuclear recoils and electron recoils. Since the nuclear recoil event rate from the collisions with CRDM is quite different from that with the nonrelativistic halo DM, to be more accurate on the effect of threshold, we derive the limits directly from the distribution of the signals of S1 and S2, rather than naively rescaling the reported experimental limits from halo WIMP searches [11,[13][14][15].
For the calculations from the deposited recoil energy T N to the position-corrected signals cS1 and cS2 b , we closely follow Ref. [66]. For a deposited recoil energy T N , the averaged number of photons N γ and charges N e are given by where W is the average energy required to create either an excitation or ion-electron pair in the liquid xenon, L is the Lindhard factor, N ex /N i is the excitation-to-ion ratio and r is the recombination probability. The prompt and scintillation light hitting the PMT photocathode will produce photoelectrons (PEs). The average number of PEs observed per emitted photon is described by the gain factor g 1 (x, y, z), and the amplification factor for charge signal is described by the parameter g 2 (x, y). The S1 and S2 signal are constructed from N pe and N prop . The bias and fluctuations are modeled by Gaussian distribution with Gauss(δ s1 , ∆δ s1 ) and Gauss(δ s2 , ∆δ s2 ) respectively. Finally, the spatial dependences of S1 and S2 will be corrected which leads to the signal cS1 and cS2 b . The Xenon-1T collaboration has performed the S1+S2 analysis based with an effective exposure of 1 ton-yr. In the analysis the DM search was performed between 3 < cS1 < 70 PE, corresponding to an average keV nr of 4.9 − 40.9 keV with an effective exposure of one ton-year [64] . The Xenon-1T collaboration also formed the S2-only analysis with cS2 b > 150 PE, corresponding to a threshold of 0.7 keV nr with an effective exposure of 22 ton-days [65].
We calculate the signal distributions of the scattering between CRDM particles and target Xenon nuclei, and derive the excluded regions in (m χ , σ χp ) plane for the Xenon-1T data (S1+S2) using the binned Poisson statistic approach [67,68]. The distribution of the background events are taken from the Xenon-1T analysis. The calculation procedure, main parameters and extended results with different parametrizations of CR flux are summarized in appendix-B.

C. Results
In Fig. 3, we show the final constraints from analysing the Xenon-1T data. It can be seen from the figure that the lower bound of the excluded region reaches σ χp 10 −(32−31) cm 2 can be extended down to DM particle mass ∼ 10 −12 eV. In the sub-eV region, the shape of the excluded region is directly related to the structures in the UHECR flux. The most stringent limits of σ χp 3×10 −32 cm 2 at m χ ∼ 10 −11 eV and 10 −5 eV correspond to the "knee" and the "toe" structures of the UHECR flux. The exclusion region closes at m χ ∼ 10 −14 eV, which corresponds to the observed suppression of the UHECR flux at ∼ 10 20 eV possibly due to the GZK cutoff. For comparison purpose, in Fig. 3 we also show the results using the CR proton and Helium fluxes in [19] which are based on the space-based direct CR measurements and are only available for energy below 200 TeV [25,27,73]. It can be seen that using the groundbased indirect UHECR measurements can extend the previous constraints by about ten orders of magnitude towards lower DM particle mass. The constraints are conservative, as we adopted the "Global-Fit" parametrization of the UHECR flux. For other parametrizations such as "H3a", "H4a" and "Global-Fit4", the resulting constraints are even stronger towards lower DM mass, as it is shown in Fig. 6 of Appendex-C.
CRDM mass is very sensitive to the detector threshold as the recoil event spectrum from CRDM is quite different from halo DM. Due to the lower threshold ∼ 0.7 keV nr of the S2only data [65], the constraints from the S2-only data which has an exposure about an order of magnitude smaller than that of the S1-S2 data (with a threshold of ∼ 4.9 keV nr [64]) turns out to be more sensitive to lighter CRDM below 10 −12 eV.
In the simple analytic approach adopted in this work, the upper bound of the excluded region due to the Earth attenuation is found to be σ χp 8 × 10 −28 cm 2 , and is almost insensitive to m χ , as expected from the fast energy loss proportional to T 2 χ in Earth attenuation discussed in Sec. III.

V. DISCUSSIONS/CONCLUSIONS
In summary, we have derived novel constraints on ultralight DM boosted by UHECRs. The constraints obtained in this work are highly model independent, as only the DM-nucleon scattering cross setion is relevant. The approach only requires the information of the presentday local Universe, even the standard cosmology is not assumed. Thus the constraints are complementary to that derived from different epochs of the Universe, such as the observa-tions of CMB [4,72,74,75], Lyman-α forest [76] and 21 cm radiations [74], etc.. If the ultralight DM particles reached chemical and kinetic equilibrium in the early Universe, very stringent constraints will arise from BBN and CMB. However, the conditions for reaching thermal equilibrium require more information on the cross section of DM annihilation and production process, which are in general model dependent. For ultralight DM particles, annihilating into most standard model particles are kinematically forbidden (for the analysis try to connect DM scattering and annihilation cross sections using crossing symmetry, see, e.g. [77]).
In this work, we considered the simplest case where the scattering cross section is a constant, i.e., energy and momentum-transfer independent. For relativistic scatterings, it more likely that the differential scattering cross section depends on the energy of both the incident and outgoing particles. For some simple models, such as the fermionic DM with a scalar mediator, it has been shown that the differential scattering cross section can be greatly enhanced at high momentum transfer. Consequently, the constraints on the total cross section at the zero momentum transfer can be many orders of magnitude stronger for lighter DM particles [15]. The approach proposed in this work can be extended to the case with energy-dependent cross sections in a straight forward manner. CR particles with energy above a few hundred TeV are mainly measured by the groundbased air-shower arrays which detect the cascades of secondary particles from the interactions of primary CR particles in the Earth atmosphere. In such indirect experiments, the information about the chemical composition is limited to determining the relative abundance of the main species or groups. Thus, it is likely that there exists different parametrizations which can describe the data equally well. We take the n-component power-law parametrization in which the CR total energy spectrum of the CR species i has the following form [45] dΦ LIS where j = (1, . . . , n) is the component index, E i (in unit of GeV) is the total energy of CR species i. The normalization constants c ij are related to Φ 0 ij in Eq. (5) of the main text by is the mass of CR species i. A global analysis to the recent UHECR data has been performed in [45]. We adopt one of the three-component "Global-Fit" model as the benchmark model with the parameters listed in Table I. In this parametrization, the first rigidity cutoff R i is around 100 TV which is the typical maximal energy from the acceleration of SNR with magnetic field around a few µ Gauss. It also well reproduce the observed hardening in the CR all-particle spectrum above 200 GeV [23,78]. In the figure, we also list a slightly extended four-component parametrization (referred to as "Global-Fit4").
TABLE I. Normalization constants c ij , power indexes α ij , and rigidity cutoffs R j in the parametrization of "Global-Fit" and "Global-Fit4" in [45]. The parameters of "Global-Fit4" which are different from those of the "Global-Fit" are shown in parentheses.
Two alternative parametrizations [79] based on the Hillas model [80] are also found in good agreement with data, which are labled as "H3a" and "H4a" in [45]. The major difference from the "Global-Fit" parametrization is that the first rigidity cutoff is quite high about 4 PV, which is responsible for the "knee" structure. In this type of parametrization the "ankle" represent the transition between the galactic and extra-galactic contributions. The corresponding parameters are listed in Table II. The CR all-particle fluxes of the four parametrizations are shown in Figure 4 together with the recent experiments data.

Appendix B: Xenon-1T data analysis
For the data analysis of the Xenon-1T experiment, we adopt the signal response model described by the Xenon-1T collaboration in [66]. The Xenon-1T experiment utilizes the liquid xenon time projection chambers to detect the recoil energy of the target nuclei from the scattering with DM particles. The deposited energy can produce a prompt scintillation light signal (S1), and the ionization electrons extracted from liquid xenon into gaseous xenon can produce proportional scintillation light (S2). The S2/S1 signal size ratio allows for discrim-  Table I, but for parametrizations of "H3a" and "H4a" in [45]. The parameters of "H4a" which are different from those of "H3a" are shown in parentheses.  4. (Left) CR all-particle spectra for four type of parametrizations in [45] together with the current experimental data [39,[47][48][49][50][51][52][53][54][55][56][57]. Right) Contributions from the three individual components in the "Global-Fit" parametrization in [45] ination between nuclear recoil (NR) and electron recoil (ER) events. For a deposited recoil energy T N , the produced total number of quantum N q is the sum of the number of excitons N ex and ion-electron pairs N i , which follows a binomial distribution N q ∼ Binom(T N /W, L) where W = 13.8 eV is the average energy required to create either an exciton or ion-electron pair in the liquid xenon, and L is the Lindhard factor. In the case of NR, it is given by [81] where k is a normalization factor, the function g( ) is proportional to the ratio of electric and nuclear stopping power, which can be parametrized as g( ) = 3 0.15 + 0.7 0.6 + , where = 11.5(T N /keV)Z −7/3 with Z = 54 the atomic number of the xenon nucleus. The distribution of N i is described by a binomial distribution N i ∼ Binom(N q , 1/(1 + N ex /N i )), where N ex /N i is the averaged exciton-to-ion ratio. The number of excitons is given by N ex = N q − N i . The excitons contribute to scintillation photon signals through de-excitation process. The ionized electrons have a probability of r to be recombined into xenon atoms and produce scintillation photons, and a probability of (1 − r) to escape the ion-electron pair. Thus the number of escaped electrons is given by N e ∼ Binom(N i , 1 − r) and the total number of photons is N γ = N ex + N i − N e . The recombination probability r is modeled by a Gaussian distribution r ∼ Gauss( r , ∆r), where r is the mean value and ∆r is the variance. The mean value is calculated using the Thomas-Imel box model [82,83] where ς = 0.057F −0.12 with F = 81 V/cm the electric field. The variance is described by ∆r = q 2 (1 − e −T N /q 3 ) with q 2 = 0.034 and q 3 = 1.7. In summary, the averaged number of photons N γ and charges N e are given by The photon and charge signals are converted into photoelectron (PE) emission of the photomultiplier tube (PMT) photocathode. The corresponding gain factors are: g 1 (x, y, z) the average number of PEs observed per emitted photon, and g 2 (x, y) the amplification factor for charge signals. Both are spatial dependent.
The spatial-dependence of the signals S1 and S2 are corrected, which results in the corrected signals cS1 and cS2 b (corresponding to the S2 signals from the bottom PMTs). These two quantities can be understood as spatial-averaged signals. For simplicity we use the spatial-averaged gain factors of g 1 = 0.142 and g 2 = 11.4 for cS1 and cS2 b , respectively. Thus in this case the number of PE is given by N pe ∼ Binom(N γ , g 1 ) and that of proportional signal is given by N prop ∼ Gauss(N e g 2 , √ N e ∆g 2 ), with ∆g 2 /g 2 = 0.25. The cS1 and cS2 b signals are constructed from N pe and N prop . The biases and fluctuations in the construction process is modeled as where we adopt mean values of δ s1(s2) = −0.065 (0.034), and variances ∆δ s1(s2) = 0.110 (0.030).
In the left panel of Fig. 5, we show the Monte Carlo (MC) simulation of the signal distributions from a typical scattering between non-relativistic halo DM and xenon nucleus for m χ = 200 GeV and σ χp = 4.7 × 10 −47 cm 2 . We adopt the Maxwellian distribution of the standard halo model (SHM) for DM with the most probable velocity v 0 = 220 km/s, the escape velocity v esc = 544 km/s, the local DM density ρ loc χ = 0.3 GeV/cm 3 and the Earth velocity v E = 232 km/s [84]. We assume the scattering process is elastic, spin-independent and isospin-conserving, and adopt the Helm form factor [60]. We find that the figure is in reasonable agreement with the result of the Xenon-1T collaboration [64]. For deriving the constraints on σ χp from halo DM, we consider two different two statistic methods. The first one is the Binned Poisson (BP) method [67,68]. First, Let us consider the single-bin case. Given an expectation value of λ = b + s events with s the theoretical prediction and b the expected background, the probability of observing N obs events is given by the Poisson distribution P (N obs |λ) = Poiss(N obs |λ).
The value of λ p at which P (λ > λ p |N obs ) is smaller than α is excluded at 1 − α confidence level (C.L.). For example, the 90% C.L. exclusion limit means α = 0.1. The required λ p can be obtained from P (λ > λ p |N obs ) = P (N < N obs |λ p ) < α. In the case of multiple bins, if (1 − α bin ) is the probability of seeing λ < λ p in that bin, the possibility (1 − α) of seeing λ < λ p in any of the bins is given by the binomial distribution where N bin is the number of bins. For a desired exclusion level of 1 − α, we then use this relation to determine α bin , and find the value of λ p . The second one is the Maximal Likelihood (ML) method. In this method, the joint likelihood is obtained by the product of individual likelihoods in each bin i, i.e., L = i L i where L i = Poiss(N obs |λ i ) is the Poisson distribution. The theoretical prediction of the event number depends on DM parameters, e.g. λ i = λ i (m χ , σ χp ) The test statistics is defined as wherem χ andσ χp are the best-fit DM parameters which maximize the likelihood. For a given value of m χ , the TS should approximately follow a χ 2 -distribution with one degree-offreedom [85]. The value of σ χ for which TS > 2.7 are excluded at 90% C.L.. For the S1-S2 combined data analysis, the Xenon-1T collaboration adopted the energy regions of interest (ROI) for cS1 as 3 < cS1 < 70 PE, corresponding to an average energy of 4.9-40.9 keV nr . The ROI for cS2 b is 50.1 < cS2 b < 7940 PE. The selection of ROIs affects the total acceptence. We take the total acceptance due to the data selection, reconstruction, noise rejection, S1-S2 correlation and single scattering, etc.. from Fig. 14 of [86]. In deriving the constraints on CRDM, we use the data of cS2 b distribution shown in Fig. 4 of [64]. In the figure the distribution is shown with respect to the rescaled quantity (cS2 b − µ ER )/σ ER , where µ ER and σ ER are the ER mean and 1σ quantile, respectively. We take µ ER = 1958 PE and σ ER = 408 PE from the cS2 b distribution shown in Fig. 3 of [64]. The number of signal counts s i in a given bin of cS2 b signal is given by where E X = 1 tone · year is the total exposure of the XENON-1T data, (cS2 b ) low(up) i is the lower (upper) endpoints of the i-th bin of the corresponding cS2 b signal. 2 = A 1 (cS1)A 2 (cS2 b ) is the total efficiency of cS2 b , where A 1,2 are the acceptance within the ROI of cS1 and cS2 b , respectively, and are vanishing outside the ROIs. The value of cS1 and cS2 b are related to each other through Eq. B6, so the total efficiency of S2 can be written as a function of cS2 b only. For the background event number b i , we directly adopt the overall expected background given by XENON-1T, which include electric recoils, neutron, surface, accidental coincidence (AC), and coherent elastic neutrino-nucleus scatters (CEνNS). In the right panel of Fig. 5 we show the upper limits on σ χp at 90% C.L. for SHM DM using the BP and ML statistic methods. The DM local density and velocity distribution are the same as that for Fig. 5. In the figure, the limits obtained by the Xenon-1T collaboration using a full profile likelihood analysis is also shown for comparison. It can be seen that the limits obtained in BP and ML approaches are quite conservative. As we are interested in conservative constraints on CRDM properties, we adopt the BP statistic approach in the main text.
Using the relation between the S2 signal and the averaged recoil energy for NR and ER process shown in Fig. 6 of the supplementary material of [65], we convert the Xenon-1T S2-only data in Fig. 4 Ref. [65] into the NR recoil energy distribution dN/dT N . The total number of events in the i-th energy bin is given by For s i , we have where (T N ) low(up) i is the lower (upper) endpoints of the i-th energy bin, and ex (T N ) is the effective exposure which is a function of recoil energy T N shown in Fig. 1 of Ref. [65]. For the background event number, we use the expected background given by XENON-1T, which include cathode, CEvNS, flat ER background for S2-only data.
Appendix C: Excluded regions for CRDM Making use of the Xenon-1T S1-S2 data, we derive the excluded regions in (m χ , σ χp ) plane using the BP statistic approach for the four different parametrizations of the primary CR flux. The results are shown in Figure 6. It can be seen that the constraints based on the "Global-Fit" parametrization is quite conservative compared with other parametrizations. Global-Fit Global-Fit4 H3a H4a FIG. 6. Exclusion regions in the (m χ , σ χp ) plane from the Xenon-1T data on the S1-S2 signals and four CR parameterizations.