Hadrophilic Light Dark Matter from the Atmosphere

Light sub-GeV dark matter (DM) constitutes an underexplored target, beyond the optimized sensitivity of typical direct DM detection experiments. We comprehensively investigate hadrophilic light DM produced from cosmic-ray collisions with the atmosphere. The resulting relativistic DM, originating from meson decays, can be efficiently observed in variety of experiments, such as XENON1T. We include for the first time decays of $\eta$, $\eta^{\prime}$ and $K^+$ mesons, leading to improved limits for DM masses above few hundred MeV. We incorporate an exact treatment of the DM attenuation in Earth and demonstrate that nuclear form factor effects can significantly impact the resulting testable DM parameter space. Further, we establish projections for upcoming experiments, such as DARWIN, over a wide range of DM masses below the GeV scale.


I. INTRODUCTION
Extensive astrophysical observations point towards the existence of dark matter (DM), currently constituting ∼ 85% of matter in the Universe (see [1,2] for a review). Despite decades of searching for its possible nongravitational interactions, the nature of the DM remains elusive. One extensively explored paradigm is that of Weakly Interacting Massive Particles (WIMPs), where the DM has a typical mass in the range of 1 GeV-100 TeV, and often appears in theories aiming to address the hierarchy problem. A less studied scenario, which could also appear in variety of well-motivated models (see e.g. [3][4][5][6][7]), is sub-GeV (light) DM. This later case is the focus of this article.
A cornerstone of DM searches are direct DM detection experiments, which looks for energy depositions on a target material associated with interactions within the experiment of traversing DM from the Milky Way halo. Conventional ton-scale direct DM detection experiments with keV-level energy thresholds (see e.g. [9,10]) face significant challenges to search for light DM, since it is not expected to produce sufficient nuclear recoil when DM interacts with the target material.
Observational signatures of DM could be drastically distinct from nominal expectations if the DM speed distribution contains a significant relativistic component. Such component can be naturally realized when DM is produced in cosmic-rays collisions with the atmosphere. The "atmospheric beamdump" has been historically exploited as an essential instrument for ex-ploring the properties of neutrinos, contributing to the discovery of neutrino oscillations [11]. Recent studies have highlighted that the atmospheric beamdump is a unique tool to explore new physics beyond the Standard Model (SM), including long-lived particles [12], heavy neutral leptons [13], millicharge particles [14][15][16], particle DM [8,17], supersymmetric particles [18], and magnetic monopoles [19]. This natural beam dump produces a steady flux that can be exploited by terrestrial experiments. Other possible sources of a relativistic DM component include boosted DM (see e.g. [20]), DM accelerated in astrophysical sources [21][22][23], the result of cosmicray DM scattering [24][25][26][27][28][29] or evaporation of primordial black holes [30]. These DM components are distinct from conventional halo DM analyses (e.g. [31,32]).
In this article, we comprehensively analyze light hadrophilic DM produced in the atmosphere by cosmicray collisions. In particular, we focus on hadrophilic dark sectors, which have been discussed within variety of scenarios (see e.g. [33][34][35][36]). We improve upon earlier work [8] of hadrophilic DM producing in an atmospheric beamdump in many significant aspects which include: considering a more realistic model of the atmospheric flux, including of additional mesons whose deay produces DM, and perform a proper treatment of the attenuation of DM through the Earth. We also present projections of this scenario for upcoming and proposed experiments.

II. HADROPHILIC LIGHT DARK MATTER
Collisions of the almost isotropic cosmic-ray flux, primarily composed of protons, with the atmosphere results in copious production of Standard Model particles from GeV to well over TeV energies [37]. This is directly analogous to a fixed target beam dump experiment. The atmospheric beam dump can produce particles associated arXiv:2203.12630v2 [hep-ph] 16 Aug 2022 with physics beyond the SM.
To explore sub-GeV DM, we consider a minimal realization within a framework of flavor-specific, hadrophilic, scalar mediator [33,34]. The effective Lagrangian includes where S is a singlet scalar mediator with mass m S , χ is a singlet Dirac fermion DM with mass m χ stabilized by a Z 2 symmetry, and g u and g χ are effective couplings.
Here, the scalar mediates the interactions between the dark sector and the visible sector, in our work, the up quarks.
In the medium where the DM propagated and in the experiments instrumented volume, DM will scatter against a nucleus with differential cross-section given by where A, T r , and m N are the nucleon number, recoil energy, and nucleus mass, respectively. The term F H is the nuclear Helm form factor [38] (see also e.g. [39]) where q = √ 2m N T r is the momentum transfer, j 1 (x) is a spherical Bessel function of the first kind, R 1 = R 2 A − 5s 2 with R A 1.2A 1/3 fm and s 1 fm. The g χ is the scalar-DM coupling, which we fix at 1.0, while the effective scalar-nucleon couplings, y spp , y npp are given by [34] with m p , m n , and m u the masses of the proton, neutron, and up-quark. The nuclear interaction cross-section can be related to the conventional spin-independent cross-section, σ SI , by

III. ATMOSPHERIC DARK MATTER PRODUCTION
The scalar mediator S can be abundantly produced in meson decays in cosmic-ray air showers. The production is modeled by the cascade equation (e.g. [40]), which, omitting the interaction terms, contains the following source term for DM production where the sum M runs over all possible mesons that can decay to S, ρ is the density of the atmosphere, X is the column depth (i.e. the depth in g/cm 2 measured from the top of the atmosphere to production point), and λ M ≡ γ M β M cτ M is the laboratory decay length of the meson that includes the boost factor γ M β M and its proper lifetime τ M . Here, the differential production rate of mesons in the shower per unit of solid angle is given by In Eq. (6), the number of scalar mediators with energies between E S and E S + dE S produced in the decay of the meson M is given by (dn/dE S ). The explicit expression for the last quantity depends on whether the scalar is produced in two-body or three-body decays. We focus on the simplest case involving only two-body decays of mesons. In this case, the energy distribution in the laboratory frame is given by where Br(M → S +P ) is the meson branching fraction to the scalar and an additional particle P , K is the Källen function [41], and p M is the meson momentum. For scalar mediators sufficiently strongly coupled to the DM, they will promptly decay to a pair of χχ in the particle shower. The resulting DM production rate can be obtained using mediator production of Eq. (6) via Eq. (8) implies that the DM is produced at the same point (height) as the scalar mediator, while accounting for the energy distribution of the DM in the decay. If the coupling between the scalar and the dark matter is not strong enough, the mediator could be made a long-lived particle and the production of DM would be suppressed. Integrating Eq. (8) over the column depth yields the expected DM flux at the surface of the Earth per unit of solid angle and per energy bin.
In this work, we include for the first time the production of scalar mediator from the decays of η, η as well as charged kaons in the atmosphere, via η → π 0 S, η → π 0 S as well as K + → π + S channels 1 . We review the branching ratios for these processes in Appendix A. We numerically solve for the production rate of the mesons in the atmospheric cascade using the Matrix Cascade Equation (MCEq) software package [44,45], which includes several models for the cosmic-ray spectrum, hadronic interactions, and atmospheric density profiles. For our base results we employed the SYBILL-2.3 hadronic interaction model [46], the Hillas-Gaisser cosmic-ray model H3a [47], and the NRLMSISE-00 atmospheric model [48] to obtain the production rate of mesons, which have been extracted following the procedure of Ref. [12,13].
The resulting DM flux at the surface of the Earth is shown in Fig. 1. We show each of the progenitor mesons η, η and K + , considering g u = 10 −3 , g χ = 1, m S = 0.2 GeV, and m χ = m S /3. With a mass of m η = 957.8 MeV, η allows for additional kinematic regimes compared to η and K + , since these have smaller masses; namely m η = 547.9 MeV and m K + = 493.7 MeV, respectively.

IV. DARK MATTER FLUX ATTENUATION
In analogy to SM particles such as muons, the flux of DM with sizable interaction strength is attenuated as it traverses the medium from the top of the atmosphere towards the experiment. At large cross-sections the medium becomes "optically thick," causing direct DM detection experiments to dramatically lose their sensitivity in this part of the DM parameter space (e.g. [49,50]).
The DM energy loss is described by where T r is the energy lost by a DM particle in a collision with nucleus N , and To account for DM flux attenuation, we numerically solve Eq. (9).
Our analysis takes into account the relevant effects of the nuclear form factor of Eq. (3), which in contrast have been neglected in previous work employing analytic approximation of Eq. (9) to treat attenuation [8]. The effects of the nuclear form factor in attenuation can significantly impact which DM parameter space regimes are direct DM detection experiments are sensitive to, as has been recently highlighted in Ref. [51].
In Fig. 1 we display the DM flux attenuation using the exact numerical treatment of Eq. (9) as well as analytic approximation given in Ref. [8]. We highlight that the approximate treatment did not include the effects arising from the nuclear form factor, which suppress large momentum transfer regimes resulting in an overestimation of the DM energy losses for underground experiments lying at a depth below ∼ 1 km of rock overburden. Having found T χ (z) with the procedure outlined above, we can determine the resulting atmospheric DM flux as a function of depth dΦ χ /dT (z).

V. EXPERIMENTAL CONSTRAINTS
The differential event rate for a detector with N t nucleons is given by where (T r ) is the detection efficiency taken from reference [52]. In Fig. 2 we show the event spectra as a function of the nuclear recoil energy T r for different target materials, including xenon (Xe), oxygen (O) and hydrogen (H). The bumps in the spectrum follow the shape of the nuclear Helm factor in the interaction cross-section. We also display in Fig. 2 the total integrated number of nuclear recoil events above experimental threshold T r , for different materials and considering characteristic model parameter values of g u = 10 −3 and m χ = 0.067 GeV. These results highlight the importance of low O(1)keV experimental thresholds available for large DM experiments, in accord with the literature results for other physics analyses and complementarity studies (e.g. [53][54][55][56]). We note that hadrophilic DM could also induce interactions in large neutrino experiments. However, such interactions require significantly higher thresholds and the resulting event rates are expected to be suppressed even taking into account the larger available fiducial mass.
To set limits on the model up-quark coupling g u from XENON1T, we integrate Eq. (11) in the recoil energy window with a threshold of 5 keV and a maximum recoil value of 40.9 keV. Following Ref. [52], we set a limit by requiring to observe more than 3.56 events for an exposure of 278.8 live-days. The limit in the plane g u −m S is shown in Fig. 3. For comparison, we also display constraints on the model from E787/E949 and MiniBooNE experiments following Ref. [34]. We also display the sensitivity reach projections for the future upgrade of XENON1T, the xenon-based DARWIN experiment [57], which is expected to have a size that is about 50 times larger with a fiducial mass of ∼ 50 tons. The effects of implementing an exact numerical treatment of attenuation that also includes the nuclear form factor are shown in Fig. 3 with the line "absorption," above which DM flux would be attenuated if the approximate treatment has been used instead. This further highlights the importance of attenuation for identifying and exploring the optimal DM parameter space with direct DM experiments. Inclusion of additional mesons allows us to set new limits above few hundred MeV, as shown on Fig. 3. This is in contrast to earlier studies of Ref. [8] that focused solely on η and hence were not sensitive to this regime. Depending on the details of the full model as well as cosmological history, additional constraints could become relevant in the parameter space of interest, such as those from Big Bang nucleosynthesis and cosmic microwave background (see e.g. [34,58]). We stress, however, that our analysis results are independent of cosmology and the mechanism of DM production.
The limits found for the parameters on the model and up-quark coupling g u can be translated into a stringent constraint for the spin-independent DM cross-section using Eq. (5). The results are shown in Fig. 4 for DM masses below the GeV scale. This highlights the significance of atmospheric DM production for exploration of light sub-GeV DM.

VI. CONCLUSIONS
Light sub-GeV DM is a promising target for exploration beyond the traditional WIMP paradigm. However, the associated DM parameter space is challenging to probe with conventional direct DM detection experiments due to their low nuclear recoils. We have comprehensively explored production of hadrophilic light DM from cosmic ray collisions with the atmosphere, whose relativistic flux leads to observable experimental signatures. The decays of mesons resulting from collisions copiously source an irreducible flux of light DM for all terrestrial experiments. In this sense, our work applies to any hadrophilic scalars which couple to a new light neutral state, irrespective of whether or not that state accounts for the DM density.
Additionally, for the first time we analyze DM production associated with η, η as well as K + mesons. We implement an exact treatment of DM flux attenuation as it traverses the Earth towards an underground detector, finding that nuclear form factor effects can dramatically impact the viable DM parameter space. We set new limits on hadrophilic light DM from XENON1T experiment and establish sensitivity projections for future experiments, such as DARWIN which will cover new low mass DM and low mass mediator parameter space. Future extensions of this work may include extending the DM model to examine inelastic DM transitions, which will impact both the recoil spectra and the DM attenuation. Another future extension would be to examine DM flux contributions beyond meson decay, such as scalar bremsstrahlung [34] which may provide access to higher mass DM.