Accelerating Earth-Bound Dark Matter

A fraction of the dark matter may consist of a particle species that interacts much more strongly with the Standard Model than a typical weakly interacting massive particle (WIMP) of similar mass. Such a strongly interacting dark matter component could have avoided detection in searches for WIMP-like dark matter through its interactions with the material in the atmosphere and the Earth that slow it down significantly before reaching detectors underground. These same interactions can also enhance the density of a strongly interacting dark matter species near the Earth's surface to well above the local galactic dark matter density. In this work we propose two new methods of detecting strongly interacting dark matter based on accelerating the enhanced population expected in the Earth through scattering. The first approach is to use underground nuclear accelerator beams to upscatter the ambient dark matter population into a WIMP-style detector located downstream. In the second technique, dark matter is upscattered with an intense thermal source and detected with a low-threshold dark matter detector. We also discuss potential candidates for strongly interacting dark matter and we show that the scenario can be naturally realized with a hidden fermion coupled to a sub-GeV dark photon.


I. INTRODUCTION
Dark matter (DM) makes up 25% of the energy budget of the Universe and manifests itself on many distance scales, from the cosmological horizon down to individual satellite galaxies [1,2]. Despite a plethora of evidence for DM based on its gravitational influence on visible matter, the identity of DM remains a mystery. The simplicity and success of the lambda-cold dark matter (ΛCDM) model of the early Universe [1], where DM evidence is seen prior to the CMB decoupling, motivates "elementary" DM candidates. Many particle physics candidates have been suggested as prospective particle DM [3][4][5], and a diverse set of experimental programs have been devised to detect them in a laboratory setting [6].
Approaches to DM include "bottom-up" phenomenological models of DM as well "top-down" scenarios where DM is part of a more complete framework that addresses other shortcomings of the Standard Model (SM) [3][4][5]. A common element in nearly all these theories is nongravitational interactions between DM and Standard Model particles. Such interactions can mediate the production of DM in the early Universe from collisions within the SM plasma, providing an explanation for the large DM abundance observed today. These interactions also form the basis of most attempts to detect DM in the laboratory. * mckeen@triumf.ca † mamoore@mit.edu ‡ dmorri@triumf.ca § pospelov@umn.edu ¶ hramani@stanford.edu An extremely robust experimental effort has been developed to detect DM in the form of weakly interacting massive particles (WIMPs), where DM connects to the SM through the weak force. The standard detection technique is to search for WIMPs scattering off atomic nuclei or electrons in the target and imparting some of their kinetic energy to the detector. Large-scale detectors with energy thresholds down to E thr keV have put very strong bounds on weak-scale DM and significantly constrain or exclude many WIMP models of DM [7][8][9]. New detector developments have achieved even lower detection energy thresholds down to E thr ∼ eV, albeit with smaller detector masses and volumes [10][11][12].
While searches for WIMPs have strong motivation, it is also important to consider how to find other DM candidates that would not show up in standard direct detection experiments. In this paper, we investigate scenarios where some or all of the DM interacts "strongly" with visible matter, corresponding to the scattering length of DM in the Earth being small compared to the overburden of a typical direct detection experiment [13][14][15][16]. 1 By exchanging momentum with particles in the Earth, strongly interacting DM particles quickly thermalize with the surrounding matter with kinetic energies E 0.05 eV, much lower than the detection threshold of nearly all underground WIMP detectors [21]. Consequently, strongly interacting DM does not deposit sufficient energy to be detected by scattering in standard direct detection experiments [14,[22][23][24][25][26][27][28].
Various strategies have been employed to test such strongly interacting DM. Indirect bounds have been derived by considering the impact of strongly interacting DM on primordial nucleosynthesis [29], the cosmic microwave background [30,31], Milky Way satellite galaxies [30,31], and large scale structure [32]. More direct limits have been obtained from DM searches at the surface of the Earth such as the CRESST [33] and EDEL-WEISS [34] surface runs as well as reanalyses [35][36][37][38] of the rocket-based X-Ray Quantum Calorimeter (XQC) detector [39]. In this case the column depth of material that the DM passes through on its way to the detector is far smaller so that it can impart considerable energy and be detected. However, the mitigation of the cosmicray-created backgrounds proves to be a considerable (but not insurmountable [40]) challenge. Lighter strongly interacting DM with mass m χ GeV can also be upscattered by cosmic rays to produce a detectable signal in near-surface neutrino detectors [41][42][43] and bounds have been determined by PROSPECT [44] and PandaX-II [45]. Limits on the terrestrial population of such DM have been obtained from the heating it can cause in cryogenic Dewars [46][47][48] as well as the deexcitation of certain long-lived isomers of tantalum (Ta) it can induce by scattering [49,50]. In Fig. 1 we summarize existing bounds on a strongly interacting DM candidate χ that makes up a fraction f χ = ρ χ /ρ DM of the total cosmological DM density and that scatters with nuclei through a spin-independent (SI) interaction. More details on these exclusions and their experimental basis can be found in Appendix A.
In this paper, we propose an alternative strategy for detecting strongly interacting DM based on upscattering that leverages precisely the same tendency of this type of DM to thermalize that makes it so challenging to detect deep underground. Specifically, the robust interactions between DM and the Earth lead to efficient capture, inhibited evaporation, and slow propagation through the crust. Together, this generates an enhanced abundance of strongly interacting DM within the Earth, both in terms of an equilibrated thermal population [46] as well as a transient sinking population [49,51]. These large DM densities open the possibility of an active "two-step" detection process where SM projectiles are used to upscatter DM within the Earth to energies that are observable in nearby detectors.
We consider two realizations of this approach in the present work. First, we investigate DM searches using existing or planned nuclear accelerators installed in deep underground labs for precision studies of nuclear reactions in low radiation environments. These include the LUNA accelerator currently deployed at the Laboratori Nazionale di Gran Sasso (LNGS) featuring proton and helium ion beams with current I b = 1 mA and terminal accelerating beam voltage V b = 0.4 MV [52,53]. Upgrade plans include LUNA-MV at LNGS with V b = 0.5-3.5 MV, I b = 0.1-1.0 mA and beams of protons, helium ions, and carbon ions [54,55], JUNA at the China Jin-Ping underground Laboratory (CJPL) with V b = 0.4 MV, I b = 2.5-12 mA and beams of protons and helium ions [56], CASPAR at the Sanford Underground Research Facility (SURF) with V b = 1.1 MV, I b = 0.25 mA and beams of protons and helium ions [57], and possibly even higher power, accelerators with gradients above V b 10 MV dedicated to neutrino physics [58]. These ion beams can upscatter DM to energies in the keV range that standard WIMP DM detectors are most sensitive to; therefore a detector placed downstream could potentially detect this DM. Given the fact that the same labs usually host a variety of DM experiments, an implementation of such a detection scheme appears to be feasible.
The second approach that we study is the upscattering of DM by gas atoms in a hot thermal source. DM that is thermalized near the Earth's surface has kinetic energies near 0.025 eV, well below the threshold of any nearfuture direct detection experiment. However, in the presence of a thermal source containing gas with temperature T ∼ 3000 K, 2 DM near the Earth's surface can be upscattered to kinetic energies around a eV, potentially within reach of planned detectors. While this presents a challenging signal, it is not outside the realm of possibility, with existing experiments already probing electron recoils with ∆E ∼ eV [10][11][12] and future proposals extending down to energy depositions below ∆E ∼ 10 meV [59][60][61][62][63]. This technique could also be further improved by optimizing the thermal sources.
Following this introduction, we discuss the capture, slowing, and resulting density distribution of strongly interacting DM in the Earth in Sec. II. Next, in Sec. III we compute the upscattering and detector event rates of strongly interacting DM in deep underground accelerator facilities for a hypothetical detector placed downstream of the beam. In Sec. IV, we present estimates for event rates in low-threshold detectors from DM upscattered by a thermal source. In Sec. V, we discuss specific realizations of strongly interacting DM. Finally, Sec. VI is reserved for our conclusions. A review of bounds on strongly interacting DM is given in Appendix A, while technical details about thermal upscattering are presented in Appendix B.

II. OVERVIEW OF DARK MATTER ACCUMULATION
As the Earth passes through the local DM halo it can interact with and trap DM particles [64][65][66][67]. This occurs when DM scatters with material in the Earth such that its velocity falls below the local escape velocity. Once Lyman-α XQC Tantalum isomer Dewar heating CDMS surface EDELWEISS surface CRESST surface Deep underground FIG. 1. Current bounds on a strongly interacting dark matter component χ that scatters with nuclei through a spin-independent interaction and that makes up a fraction fχ = 1, 10 −3 , 10 −6 , 10 −9 of the total dark matter density as a function of its mass and effective per-nucleon cross section. Details of how these bounds are obtained and a list of the original experimental sources are given in Appendix A.
captured, a strongly interacting DM particle will undergo further scatterings to thermalize and sink until it reaches an equilibrium distribution within the Earth. Over the history of the Earth, this leads to a thermal (or Jeans) population of DM within it. In addition to the thermal population, an itinerant traffic jam population can arise when the sinking rate of DM is slow relative to its capture rate [49][50][51]. We present our estimates for both populations in this section.

A. Thermal density
The thermal density is the population of DM captured by the Earth over its history that has thermalized [66,67]. When DM interacts strongly with the Earth, several new qualitative features arise relative to more weakly interacting scenarios, as emphasized in Ref. [46]. In particular, DM is expected to reach local thermal equilibrium with the baryonic matter around it, and thermal evapo-ration of DM from the Earth is impeded by scattering. Together, these features can lead to a significant thermal population of captured DM near the surface of the Earth [46].
For a DM particle χ that does not annihilate significantly, the total number of χ particles captured in the Earth over its history is [65] where t ⊕ 4.54 × 10 9 yr is the age of the Earth, C is the total capture rate, and E is the loss rate due to evaporation. We apply the calculations of Ref. [46] specific to strongly interacting DM to estimate the rates of capture and evaporation.
The capture rate for the large cross sections that we study here is well approximated by the Earth geometric cross section times the local DM flux up to small correc-tions [46], Here, R ⊕ 6370 km is the Earth radius, v eff 250 km/s is the effective local DM flux velocity, and f cap is the fraction of DM that is captured when it strikes the Earth. The capture fraction is less than unity because some of the DM that encounters the Earth will be reflected and is estimated to be [46] f cap = 4 π where v esc (R ⊕ ) 11.2 km/s is the surface escape velocity and β χ = 2m χ m N /(m χ + m N ) 2 is the mean fractional energy loss per collision with nucleus N . Up to the factor of f cap , this capture rate matches well with other approaches [68,69]. Evaporation of DM arises from thermal collisions that upscatter DM particles beyond the local Earth escape velocity [65,66]. When the interaction length of DM is much smaller than the radius of the Earth, this process is strongly impeded by the rescattering of DM before it exits the Earth or the atmosphere. Assuming a DM number density profile of n χ (r) and local velocity v th,χ (r) in the Earth, the evaporation rate from Ref. [46] is where all quantities are to be evaluated at the lastscattering radius r LSS . This is obtained from the condition with n N (r) the local number density of species N . Determining the DM evaporation rate requires the radial distributions of its mean velocity v th,χ (r) and its number density n χ (r) within the Earth and the atmosphere. The mean velocity is just the local thermal one, v th,χ (r) = 8 T (r)/π m χ where T (r) is the local temperature. For the number density, we follow Ref. [46] and obtain it by balancing the local DM thermal pressure p χ = n χ T with gravity, where T = T (r), and M r is the total mass enclosed at radius r from the center of the Earth. Note that this expression is based on the assumption of local thermal equilibration, which is expected to occur for the range of large cross sections of interest. We integrate this expression using the PREM Earth density and temperature profiles from Refs. [70,71] and the NRLMSISE-00 model of the atmosphere [72] and set n χ (∞) = n halo χ .
Our estimates for the thermal density of captured DM are subject to a number of uncertainties whose total impact we expect to be small in terms of the regions in the m χ -σ χn plane where the density is strongly enhanced. The most significant one is our simplified treatment of evaporation following Ref. [46], which neglects the possibility of DM escaping from inside the last scattering surface, r < r LSS . This would tend to push the masses for which evaporation strongly depletes the thermal DM density to slightly higher values. However, we find that the depletion due to evaporation varies very rapidly with mass, and thus, the shift in mass will be mild. A further approximation is our assumption of thermal equilibration. We find that this is a good approximation when the DM density is significantly enhanced near the Earth's surface over the halo value since this requires large cross sections and, therefore, frequent scattering with material in the Earth, to avoid major losses to evaporation.

B. Traffic jam density
In addition to the thermal population of accumulated dark matter, large DM cross sections also give rise to an itinerant traffic jam population consisting of infalling DM that has not yet reached its equilibrium distribution [49,51]. The local density of this component can be obtained from the steady-state conservation of DM flux in the form where v z v eff /2 is the mean inward flux velocity of local DM beyond the Earth and v z (r) is the average downward DM velocity at radius r. Obtaining n χ (r) is therefore just a matter of computing v z (r).
DM entering the Earth is pulled inward by gravity while being slowed by collisions with the surrounding material. To describe these effects, it is convenient to define the inward coordinate z = R ⊕ − r and split the evolution of v z into two regimes: v > v th,χ and v < v th,χ , where v ≥ v z is the DM speed and v th,χ 8 T /π m χ is the local thermal DM speed. In both regimes, we find that the downward velocity is well described by the relation with the sum running over all relevant targets N with mass m N , µ N = m N m χ /(m N +m χ ) is the reduced mass, n N is the target number density, σ T,N is the transfer cross section, andṽ N = max{v, v th,N } where v th,N 8 T /π m N is the relevant thermal velocity of the target. Evaluating Eq. (8) also requires that we specify the DM speed v. For v > v th,χ , the speed is reduced by collisions and can be tracked through the rate of energy loss [73] such that This expression neglects gravity, which is a good approximation for the parameters of interest in this work. Note as well that v and v z have different kinematic dependencies on the masses for m χ m N . In this limit, multiple collisions are needed to reduce v by a significant amount even though a single collision can reflect the DM particle and change v z by order unity. We apply Eq. (9) to estimate the slowing of v down to v th,χ . Once v → v th,χ , further collisions are not expected to change its speed on average, and thus, we fix v = v th,χ from this point on.
As for the calculation of the enhanced thermal density, our results for the traffic jam density are subject to several uncertainties. The main sources here are our neglect of reflection of incoming DM and evaporation once it reaches a steady state. These only become significant at lower masses m χ 3 GeV, and for this reason we do no extend our traffic jam calculation below this value.

C. Enhanced densities
In Fig. 2 we show the enhanced density of DM at a depth of z = 1.4 km under the Earth, corresponding to the overburden at LNGS, as a function of mass m χ and per-nucleon cross section σ χn . Both the thermal and traffic jam contributions to the enhanced density are included. To connect the per-nucleon cross section to cross sections on nuclei, we take σ T,N = min{A 2 (µ χN /µ χn ) 2 σ χn , 4 π r 2 N }, where A is the atomic mass of nucleus N , and r N (1.2 fm) A 1/3 is the nuclear radius. This form corresponds to a SI point interaction with a nuclear form factor of unity together with a saturation of the cross section at the geometric area of the nucleus [49,74]. Since most of the scatterings leading to capture and accumulations have a low momentum transfer relative to the inverse nuclear radii 1/r N , we expect that setting the nuclear form factor to unity should be a good approximation.
The DM densities shown in Fig. 2 are much larger than the local halo density, particularly for larger cross sections. This enhancement has two primary features corresponding to the thermal and traffic jam components, respectively. The greatest enhancement between m χ ∼ 1-10 GeV comes from thermal accumulation and coincides with that found in Ref. [46]. The density enhancement from this component at the locations of other underground accelerators is very similar. Evaporation depletes this population for m χ 1 GeV, while for m χ 10 GeV, the thermal population is mainly located deeper within the Earth. 3 Instead, the dominant enhancement at larger 3 Since we do not consider evaporation effects in our traffic jam masses m χ 10 GeV comes from the traffic jam population. If χ makes up only a fraction f χ of the total DM energy density, the densities shown in Fig. 2 are reduced by the same factor.

III. UPSCATTERING OF DARK MATTER BY ACCELERATOR BEAMS
In this section we investigate the upscattering of strongly interacting dark matter by the beams of deep underground accelerators such as LUNA [52,53], LUNA-MV [54,55], JUNA [56], and CASPAR [57]. We compute the upscattering rates as well as the detection rates through elastic nuclear scattering in a xenon detector of modest size.

A. Dark Matter Upscattering by Accelerator Beams
Consider a beam of nuclei of mass m b and kinetic energy E b m b incident on a cloud of DM particles χ effectively at rest. If a beam nucleus collides with a DM particle in the cloud, the DM will be upscattered to a velocity where θ is the angle of the outgoing DM relative to the beam direction. Should the upscattered DM particle collide with a target nucleus N = (A, Z) in a detector, the calculations, we only include this component of the enhanced density for mχ ≥ 3 GeV.
nucleus will recoil with kinetic energy where α is the angle between the recoiling nucleus and the incident DM direction. We note that all the factors multiplying E b in this expression are less than unity and represent the combined kinematic suppression from the two scattering reactions involved. In Fig. 3 we show the maximum nuclear recoil energies E max R setting cos θ = 1 for DM upscattered by beams of protons (left) or carbon (right) with kinetic energies E b = 0.4 MeV (solid) and E b = 1.0 MeV (dashed) on targets of hydrogen (H), helium (He), germanium (Ge), and xenon (Xe). These recoil energies fall within the sensitivity windows of typical underground nuclear recoil DM detectors.
Given an accelerator beam of particles with energy E b , total current I b , and charge per particle Q b , the differen-tial rate of DM upscattering per unit beam travel length is where c θ = cos θ corresponds to the outgoing DM angle relative to the beam, z ∈ [−L/2, L/2] ranges over the beam travel region after full acceleration, and n χ is the local χ DM number density. From this, we see that total rate of upscattered DM is proportional to the quantity where L is the total length over which the fully accelerated beam travels.

B. Detection of upscattered dark matter
For a detector placed downstream of the beam, the measured rate of DM scattering in the detector is where = (θ, z) is the path length in the detector for a DM particle upscattered at point z through angle θ, n N is the number density of the target nucleus, σ χN is the total DM-nucleus cross section, P thr (θ, E thr ) is the probability that the scattering will yield a recoil energy above the detector threshold E thr , and P sh (θ, z) is the probability for DM to scatter in material before reaching the detector. The exponential factor is the probability for a DM particle from (z, θ) to scatter at least once in the detector; it reduces to σ χN n N when this combination is much less than unity. In the second line of Eq. (14), we have factored the expression into a total upscattering rate times a dimensionless acceptance factor for scattering above threshold in the detector.
The result of Eq. (14) is very general, and it is instructive to evaluate its components for a specific detector geometry. We consider a spherical detector of radius r located along the beam axis a distance d from the end of the beam pipe, as illustrated in Fig. 4. For this configuration, the DM path length is where D = L/2 − z + d + r and θ s = sin −1 (r/D).

C. Application to spin-independent pointlike DM
If we specialize further to DM that scatters primarily through a spin-independent point interaction, the differential DM-nucleus cross section is where is a nuclear form factor for SI scattering [75,76], andσ N is given in the Born approximation bȳ for an effective per-nucleon cross section σ χn . Since we study very large cross sections in this work, we also consider the possibility that the Born approximation on which Eq. (17) is based might break down. While the way in which this occurs depends on the detailed interactions between DM and nucleons, there exists a simple prescription based on geometric saturation that provides a reasonable approximation to calculations in a wide range of models [49,74]. Specifically, we bound from above the total nuclear cross section derived from Eq. (16) by the geometric cross section σ χN ≤ 4πr 2 N with r N 1.2 fm A 1/3 . This is equivalent to the replacement ofσ N in Eq. (16) byσ N,eff defined bȳ With this form, we can express the nuclear cross section portion of Eq. (14) by where x thr = E thr /E max R . We can also specify the upscattering rate more precisely if we specialize to a SI interaction. For a low-energy beam of protons, If the interaction connects DM to protons and neutrons with equal strength, we can identify σ χp = σ χn defined in Eq. (17). This result can also be generalized to lowenergy beams of nuclei. Using the saturation prescription described above, we find with and σ b,eff defined as forσ N,eff in Eq. (18) Fig. 5 we show the estimated detector rates of beam upscattered DM as a function of mass m χ and per nucleon cross section σ χn assuming a pointlike SI interaction for a potential beam and detector apparatus. We take beam parameters motivated by the LUNA accelera- In both cases, the detector is assumed to be a sphere of radius r = 10 cm located a distance d = 50 cm from the end of the beam pipe with a lower energy threshold of E thr = 5 keV.
tor [52,53] with an accelerated beam section of L = 5 m and a kinetic energy per particle of E b = 0.4 MeV for proton beams with current I b = 1.0 mA (left) and carbon 12 C + beams with current I b = 0.1 mA (right). For both beam types, we assume a detector consisting of a sphere containing liquid xenon of radius r = 10 cm located along the beam axis at a distance d = 50 cm from the end of the beam pipe with a lower detection energy threshold of E thr = 5 keV. See Fig. 4 for details of the setup.
The detector scattering rates shown in Fig. 5 are significant and suggest that this method could be used to test strongly interacting DM even for fractional densities f χ 1. These rates trace the DM density enhancements shown in Fig. 2 to a large degree. They are largest for m χ ∼ 1-10 GeV, corresponding to the enhanced thermal DM population discussed in Sec. II, although there is also a shoulder at larger masses from the traffic jam population. For masses below m χ ∼ 1 GeV, the detection rates are reduced by the lower DM population due to evaporation as well as the energy threshold we assume for the detector. This is most clearly visible in the right panel of Fig. 5 where the kinematic mismatch between the masses of the beam carbon nuclei and the m χ ∼ 1 GeV DM leads to reduced upscattering energies and a sharp cutoff for DM masses below a GeV. We also note that the scattering rates fall off more quickly at smaller cross sections than the DM density. This arises from the dependence of the detector acceptance on the cross section since smaller cross sections produce a lower probability for the DM to scatter in the detector as it passes through.
In Fig. 6, we compare the potential reach of the accelerator upscattering method to existing bounds on DM with SI scattering on nuclei for DM fractions f χ = 1, 10 −3 , 10 −6 , 10 −9 as a function of the DM mass m χ and per nucleon cross section σ χn . The green (proton beam) and yellow (C + beam) contours encircle the regions that would yield at least 10 events per year for beam and (liquid xenon) detector parameters as described above for  Fig. 1 and discussed in Appendix A. These plots show that accelerator upscattering can potentially test fractional DM components beyond what has been achieved so far.
While our analysis neglects several potentially important aspects such as backgrounds and attenuation of upscattered DM due to shielding, we argue that this technique is a promising and realistic method for testing strongly interacting DM. The detector configuration studied here is relatively modest and in line with smaller liquid xenon detectors such as the XENON10 [77] and ZEPLIN-III [78] experiments. Moreover, these experiments achieved near-zero backgrounds, albeit within less active environments.
A potential source of backgrounds relative to these more traditional dark matter search experiments are those from the beam itself as well as the nuclear targets the beam is directed to. Some of these can be mitigated or avoided by choosing a suitable detector location. For LUNA-MV, the primary central beam is split and directed to a pair of target stations located at large angles [54,55]. The current setup would allow a detector to be seated directly behind the beam splitter and away from the beam path and target stations. Similarly, JUNA [56] and CASPAR [57] feature a bend in the beam after the initial accelerating segment and a section of evacuated beam pipe, and could similarly accommodate a detector directly behind the bending magnets. Residual beam backgrounds produced with the low beam energies of these accelerators have very small penetrating power in the case of protons and charged nuclei, while secondary neutrons and photons may require a layer of high Z material to mitigate. Much thicker shielding placed around the detector away from a beam aperture would reduce backgrounds from the nuclear reactions studied in the target stations, although it may be necessary to avoid searching for DM while neutron-producing reactions are being studied. These background sources would be relatively straightforward to measure in situ before staging a DM search.
For very large DM cross sections, the upscattered DM beam can be attenuated by accelerator components and detector shielding between the beam pipe and the detector material. We estimate that the attenuation of DM from a reasonable thickness < 5 cm of steel or lead is minimal but can be significant for thicknesses greater than this. To avoid losing potential signal to shielding, a mostly unshielded beam aperture can be left for upscattered DM to pass through combined with thicker shielding elsewhere.
Despite these important practical challenges, the large DM rates we find suggest that they could be overcome. Of course, a more detailed study beyond the scope of the present work would be needed before installing an actual detector.

IV. UPSCATTERING OF DARK MATTER BY THERMAL SOURCES
We seek next to determine the upscattering rate of strongly interacting dark matter using very hot thermal sources. This second, complimentary method is presented in the context of hot gas surrounding the filament of a conventional incandescent lightbulb. Given that such filaments can heat gas to temperatures on the order of O(10 3 K), DM could be scattered to kinetic energies near O(eV) by the high energy Boltzmann tail. Such energies are large enough to be potentially detectable at nearfuture experimental setups.
A. Thermal upscattering of DM Consider a gas of SM particles heated to a temperature T g . This gas can scatter on a population of DM thermalized with the Earth at a temperature T χ T g , accelerating the DM and potentially making it more easily detectable. For example, the filaments in typical commer-cial lightbulbs are heated to around 3000 K, along with a substantial number of gas particles surrounding the filament. To estimate this number, we assume that the volume occupied by these hot gas particles is π 2 coll fil with coll the scattering length of the gas and fil the length of the filament. Taking coll = 0.5 mm, fil = 20 cm, and a gas pressure of 0.7 atm leads to the crude estimate of N g ∼ 10 17 hot gas particles. In this case the rate of hot gas particles a unit length of DM particles encounters is where we have normalized the mass of the gas particles, m g , to that of 4 He. Comparing this to Eq. (13), we see that it corresponds to a substantially larger rate of DM upscattering, albeit at much lower energies.
Assuming that the gas nuclei scatter on DM through a spin-and momentum-independent interaction with cross section σ χg , the flux of DM a distance d from the region of hot gas is In this expression we have assumed that the cross section for DM to scatter on the gas nuclei scales as A 4 g and has not reached geometric saturation. The kinetic energy distribution of the upscattered DM, in the limit that T χ T g and assuming that the interaction matrix element is momentum independent, is 4 µ χg is the reduced mass of the DM and gas. Thus, the typical kinetic energy of DM upscattered by the hot gas is of order T g unless there is a large hierarchy between m χ and m g . With T g = 3000 K = 0.26 eV, a detector with an energy threshold of order 0.1 to 1 eV placed near the hot gas can potentially detect the upscattered DM.

B. Detection of thermally upscattered DM
Gas-based detectors may offer the best prospects for sensitivity to eV-scale energy depositions from nuclear 4 We show the expression for finite Tχ in Appendix B. scattering in the near future. NEWS-G [79] has demonstrated energy thresholds on the order of hundreds of eV. Pushing even lower, exciting rotational and vibrational molecular modes in a gas could offer a pathway to detecting O(0.1 eV) energy depositions [80]. In the following, we roughly estimate the sensitivity of an idealized future detector, assuming that the detector is sensitive to all DM with kinetic energy above E min χ with a probability of interacting given by L det / scatt where L det is the length of the detector and scatt is the DM scattering length in the detector material. We leave estimates folding in realistic detection efficiencies for future work.
We will crudely assume that the scattering length in the detector is dominated by elastic scattering of DM on detector nuclei with a cross section σ χN so that scatt ∝ 1/σ χN . The rate of scattering near a thermal source producing the flux of DM given in Eq. (24) is then In this expression, N det is the number of nuclei in the detector with atomic number A det , and f det is the fraction of DM with kinetic energy above E min χ that can be detected, (26) This approximation is valid when the thermalized DM temperature can be neglected, T χ T g . The numerical estimate in Eq. (25) uses the scaling σ χN ∝ A 4 det σ χn to relate the cross section to scatter on a nucleus to that on a nucleon, assuming a point interaction with neutrons and protons of equal strength. In the results that follow, however, we make use of the more sophisticated relation between σ χN and σ χn described in Sec. III C.
In Fig. 7, we show the per-gram detector rates as a function of mass m χ and per-nucleon cross section, assuming a pointlike SI interaction using the DM densities computed in Sec II. We fix N g = 10 17 with T g = 3000 K, a distance d = 10 cm away from the detector which should allow for sufficient thermal insulation. We assume that the hot gas is helium or argon and that the accelerated DM scatters on carbon nuclei in a detector filled with carbon monoxide gas. Optimistically, we set E min χ = 0.8 eV; detailed studies to more accurately model the detector response are beyond the scope of this work.
Another possibility for low-threshold detectors that has been extensively discussed in the DM direct detection literature involves DM scattering on electrons [81]. Unfortunately, this does not obviously provide a clear pathway for the detection of DM upscattered by nuclei in a hot gas. This is largely due to the fact that DM that is most efficiently captured in the Earth and subsequently upscattered by nuclei has m χ ∼ 1-100 GeV (cf. Fig. 2) which is not well kinematically matched to maximize the momentum transfer when scattering on electrons.

V. REALIZATIONS OF STRONGLY INTERACTING DARK MATTER
Having presented two promising search techniques for a strongly interacting dark matter component χ with mass in the range m χ ∼ 1-100 GeV, we turn next to an examination of potential candidates for such a species. We focus on two possibilities: DM made up partially or entirely from particles charged under the strong force; and DM connecting to the SM through a new dark mediator particle. As a specific realization of the second class of candidates, we perform a detailed study of asymmetric DM coupled to a dark photon.

A. Dark matter with QCD interactions
Proposals for stable relics with constituents charged under quantum chromodynamics (QCD) have a long history as well as more recent interest. Early examples include stable di-baryons consisting of six quarks [17] and tightly bound quark nuggets of macroscopic size [82]. Recent studies of DM connected to QCD include candidates composed of the standard quarks and gluons, both microscopic [15,17,83,84] and macroscopic [82,85,86], as well as DM candidates containing exotic particles charged under the strong force [18,87]. Among these proposals, the methods presented in this work are most applicable to microscopic relics consisting only of quarks and gluons. They could also be applied to bound states involving new QCD-charged states, but these are very strongly constrained by direct [88][89][90] and indirect [91,92] searches at the Large Hadron Collider (LHC) and must be heavier than at least a few hundred GeV to have avoided detection, putting them beyond our primary region of sensitivity.
The best-studied proposals for a stable QCD relic are the sexaquark (S) and the H-dibaryon with quark quantum numbers uuddss [15,17,83,84]. The state S must lie in the mass range m S ∼ 1860-1890 MeV to be stable while also not destabilizing nuclei [84,93,94], while the H-dibaryon would be slightly heavier at around 2150 MeV [17]. Direct searches for a stable di-baryon have been inconclusive [95,96] and its study on the lattice is very challenging [97]. There has also been some debate over the astronomical properties of a stable dibaryon, with Ref. [94] finding that it would obtain a very small density fraction n S /n DM ∼ 10 −11 through thermal freeze-out, and Ref. [98] arguing that the existence of this state would drastically change the evolution and properties of neutron stars. Both of these works rely on reasonable estimates for di-baryon interaction cross sections, while Refs. [99][100][101] counter that these can be significantly suppressed or avoided. We do not attempt to resolve this debate here, but we do note that a stable QCD relic, possibly making up only a small fraction of the total DM abundance, would be expected to have a mass and nucleon interaction cross sections within the range our proposals are most sensitive to.

B. Asymmetric dark matter coupled to a dark photon
Dark matter can also have a large nucleon scattering cross section if it connects to the SM through a relatively light mediator particle. In general, both lighter mediators and larger couplings increase the resulting cross section in the limit of zero momentum transfer. At the same time, these properties also tend to make the mediator easier to observe, and there is a generic tension between large DM-nucleon cross sections and direct bounds on the mediators [74,102,103].
As a concrete example of a DM candidate with a large nucleon cross section that also illustrates the tension with mediator searches, we investigate a Dirac fermion χ coupled to a massive dark photon that connects to the SM through gauge kinetic mixing. A similar analysis could be applied to other dark vectors that couple to baryons such as B [104,105] or B − L [106][107][108]. We focus on an asymmetric density of the dark relic χ to maintain a large enhanced density in the Earth.
The simple theory we consider consists of Dirac fermion χ coupled to a sub-GeV massive dark photon A with strength g [20,109,110]. The corresponding Lagrangian is where the last term is the kinetic mixing with SM hypercharge. The dark photon mass can arise from a Higgs [20,109,110] or Stueckelberg [111,112] mechanism. For m χ ∼ 1-100 GeV, m A ∼ 10-100 MeV, and within the currently allowed range, we find that small relic fractions f χ 1 and large per-nucleon cross sections σ χn ∼ 10 −32 -10 −27 cm 2 are obtained. In this range, we show that accelerator upscattering can test the theory beyond what has been possible through other methods, particularly when the relic density of χ is set in part by a charge asymmetry relative toχ.
Annihilation of the χ fermion in the early Universe is dominated by χχ → A A for m χ > m A and α = g 2 /4π 2 α [20]. In this limit, the Born-level annihila- tion cross section at zero velocity is For larger values of α , the total cross section during annihilation in the early Universe can be enhanced by the Sommerfeld effect [110,[113][114][115] and bound state formation [109,[116][117][118]. We follow the methods of Ref. [119] to compute the relic density of χ while taking these effects into account and allowing for a possible charge asymmetry. In this approach, the mass of the A mediator is neglected, which we have checked to be a good approximation during freeze-out for the parameter range we study by comparing the massless mediator result to an estimate for the massive result based on analytic results for the Hulthèn potential [120,121]. The resulting total relic fraction of χ andχ can be written in the form [122] where r = nχ/n χ is the ratio ofχ to χ relic densities, andf χ is the DM fraction in the completely asymmetric limit (r → 0) as determined by the charge asymmetry whose origin we do not specify. The ratio r is obtained from a freeze-out calculation and decreases exponentially with an increasing annihilation cross section [122].
The kinetic mixing interaction of Eq. (27) leads to SI scattering between relic χ particles and nucleons mediated by the dark photon. For m A m Z , the scattering is almost entirely with protons with a zero-momentum cross section of [20] Extending this to scattering on nuclei N = (A, Z) at where σ χN = (µ N /µ p ) 2 Z 2 σ p , F N is the same nuclear form factor as before, and we have introduced a DM form factor It is straightforward to generalize our prescription for saturation of the total DM-nuclear cross section at σ χN ≤ 4π r 2 N by replacing |F N | 2 → |F N | 2 |F χ | 2 in Eq. (18). In Fig. 8, we show detector event rates for accelerator upscattered DM in this theory for benchmark dark photon masses of m A = 15 MeV (top) and m A = 50 MeV (bottom). We assume the same accelerator and xenon detector parameters as in Sec. III, and we compute rates for a beam of energy of E b = 0.4 m and currents I b = 1.0 mA for protons (left) and I b = 0.1 mA for 12 C + (right). For two reasons, the rates in this scenario are reduced relative to the SI point interaction considered in Sec. III and shown in Fig. 5. First, the nuclear cross section scales as Z 2 instead of A 2 prior to saturation. And second, for m A = 15 MeV, there is a significant additional suppression from the DM form factor of Eq. (32) corresponding to m 2 A < |q 2 | in the beam or detector scattering. For m A = 50 MeV, we find that the DM form factor is typically close to unity.
In Fig. 9, we show the event rates per gram of gaseous CO detector for dark photon DM upscattered by a thermal source. As in Sec. IV, we assume a source at temperature T = 3000 K consisting of 10 17 atoms, either argon (left) or helium (right). The detector is taken to consist of 1 g of CO gas that is sensitive to DM kinetic energies above E min χ = 0.8 eV. As for accelerator upscattered DM, these rates are reduced relative to SI DM with equal couplings to protons and neutron due to the scaling with Z 2 rather than A 2 prior to (nuclear scattering saturation). In contrast, however, this method In both cases, the detector is assumed to be a sphere of xenon of radius r = 10 cm located a distance d = 50 cm from the end of the beam pipe with a lower energy threshold of E thr = 5 keV. The solid (dashed) red lines show the boundary above which α > 1 is needed to achieve the corresponding proton cross section for a 15 (50) MeV dark photon. is not impacted by the benchmark dark photon masses (m A = 15, 50 MeV) since the momentum transfers in this thermal source method are always much smaller, By working within a specific theory, we are able to de-rive consistency conditions on the range of nucleon scattering cross sections σ χp . Fixing m A and demanding α ≤ 1 to maintain perturbativity, the largest possible σ χp corresponds to the greatest value of allowed by searches for dark photons [123]. Assuming the dark pho-ton decays visibly, the limit is approximately 8×10 −4 from NA48/2 [124] and BaBar [125] for both m A values. This upper bound on σ χp is shown by the solid red lines in Fig. 8 and agrees with Refs. [74,102,103]. A lower bound on σ χp can also be derived if we make the further assumption of a standard cosmological history and demand that the χ relic density be primarily asymmetric, which implies a lower bound on α for a given DM fraction f χ . Combined with experimental lower bounds on [126][127][128][129], this provides a lower limit on the proton cross section given m A and f χ . We find that the condition is not relevant to the parameter space tested by our methods.
In Fig. 10, we compare the sensitivity of our proposed accelerator upscattering method to this DM candidate to other bounds on it from direct detection and cosmology. The accelerator and detector parameters are the same as above. Each of the four panels in the figure correspond to different values of the χ DM fraction f χ = 1, 10 −3 , 10 −6 , 10 −9 . The contour lines in Fig. 10 correspond to event rates of at least 10 per year for proton (p, green) or carbon (C + , yellow) beams, or a lightbulb filled with argon gas (blue), and a dark photon mass of 15 MeV (solid) or 50 MeV (dashed). The gray regions show the collected exclusions on χ discussed in Appendix A.

VI. CONCLUSIONS
In this paper, we have proposed two new methods to detect strongly interacting dark matter (DM) that has accumulated in the Earth. Both methods leverage the tendency of strongly interacting DM in the local halo to scatter with nuclei in the Earth to slow down and thermalize. On one hand, this leads to a DM population with density greatly enhanced compared to the local halo density. On the other hand, this DM population has a reduced kinetic energy, making it difficult to detect in standard direct detection searches. The two methods we propose use energetic sources, in the form of underground nuclear accelerators or thermal sources, to upscatter the DM population accumulated in the Earth to make it easier to detect. We analyzed this scenario in a model-independent way, assuming only a point interaction between DM and nucleons, as well as in the context of a well-motivated model where DM interacts with the Standard Model through the exchange of a dark photon.
The first method we investigated for detecting strongly interacting DM is based on nuclear beams produced by underground accelerators, such as the LUNA [55] facility at LNGS. These accelerators feature beams of protons or light nuclei with kinetic energies in the E b ∼ MeV range and currents near I b ∼ 1 mA. When a beam nucleus strikes an ambient DM particle with mass m χ ∼ 1-100 GeV, the resulting DM energy typically falls into the sensitivity range of standard WIMP DM detectors. Based on realistic accelerator beam parameters and a modest detector proposal, we find potentially observable DM event rates for masses in the range m χ ∼ 1-100 GeV for both a pointlike spin-independent interaction and for scattering mediated by a sub-GeV mass dark photon. This technique is able to probe parts of parameter space unconstrained by other searches, as seen in Figs. 6 and 10.
In the second method, ambient DM is upscattered by hot gas in a thermal source. This has the virtue of producing locally a very large flux of upscattered DM, but only with kinetic energies at the eV scale, much lower than the sensitivity range of standard WIMP detectors. However, proposed ultralow threshold detectors could soon be sensitive to DM in this energy range [10-12, 59-63, 79, 80]. Combined with the enormous effective fluxes of thermal sources and enhanced local DM densities in the Earth, this second approach could be a promising way to test strongly interacting DM in the near future.
The two methods we have proposed can be characterized by: (i) the kinetic energy imparted to the upscattered DM; (ii) the effective rate of upscattering, which we define to be the flux of beam particles times the volume of DM impacted, corresponding to Eqs. (13) and (24). The ranges of these parameters is sketched in Fig. 11 for the two methods we have studied. These two methods of DM acceleration overlap well with two major thrusts in DM direct detection: searches for WIMP-like DM through nuclear scattering with energy deposition at the keV scale and searches for the scattering or absorption of light DM with energy transfers near the eV scale. Figure 11 offers a useful comparison of these techniques and also motivates studying other potential ways of accelerating and detecting strongly interacting DM for other combinations of energies and upscattering rates. For comparison, we also show in Fig. 11 a similar quantity for deep underground detectors surrounded by thermalized DM, where now the beam is identified with the detector volume moving through the surrounding DM at thermal velocity.
It is important to point out that the existing beam sources that we have considered are not optimized for maximizing the effective DM upscattering rate shown on the y axis of Fig. 11. It would be interesting to explore the feasibility of increasing this parameter at the expense of beam energy, collimation, and monochromaticity using other sources. A natural extension of beamlines and storage rings to higher densities are devices with plasma confinement with the main application being a nuclear fusion reactor. Temperatures up to 160 MK (≈ 10 keV) have been reached and sustained for O (1000) seconds at the EAST reactor [130] with ITER not too far behind [131]. The typical number densities required for fusion are in the range of 10 20 m −3 with plasma volume near 840 m 3 [132]. The viability of a low-threshold detector close to such a reactor and the impact of reactor shielding need to be analyzed in order to understand whether the large number of energetic particles confined in the plasma would help to offset much larger cosmogenic backgrounds existing on the Earth's surface to make it an efficient way of probing strongly interacting dark matter. One could also explore whether much more compact plasma retaining devices can be installed in underground laboratories with the main purpose of serving as a source of DM acceleration.
Going to higher energy beams, ambient DM could be upscattered to higher kinetic energies, making it potentially easier to detect or distinguish from background. This could lead to signals in accelerator neutrino experiments with larger energy thresholds, albeit with larger volumes compared to ones considered in this work, in analogy to cosmic ray upscattered DM [41]. Examples of such setups include LSND [133,134], ISODAR [135], and DUNE [136]. These setups could lead to sensitivity to larger DM masses owing to the superior ratio of beam energy to detection energy compared to accelerators considered in this work. The ISODAR proposal [135] looks especially promising in this respect, as it is planned to be installed at an underground laboratory with intrinsically low cosmogenic background. In the analysis above, we have presented exclusions on strongly interacting dark matter from direct detection experiments on the surface of the Earth and underground, measurements by satellite detectors, and cosmological observations. These bounds are typically presented in terms of SI DM with equal couplings to protons and neutrons that makes up the full DM abundance. In this appendix, we explain briefly how we generalize these bounds to a species χ that makes up only a fraction f χ of the total DM density or that couples primarily to protons through a dark photon mediator. We also generalize previous analyses to include saturation of the nuclear cross section.
Direct DM searches typically present their (spin-independent) results as a lower exclusion limit on the per-nucleon cross section σ χn . As the DM interaction with nuclei becomes stronger, DM will interact significantly with nuclei in the Earth, detector shielding, and atmosphere, causing it to slow down. Since realistic experiments have finite energy sensitivity, this leads to an upper exclusion limit on the interaction strength that they are sensitive to that depends on their energy threshold and detector location. A full treatment of this effect typically requires a detailed simulation of DM scattering in the experimental overburden [26,37,143] that goes beyond the scope of this work. However, a simple estimate of overburden scattering based on a set of simplifying assumptions [14] (that are not always justified) has been found to give a reasonable approximation of the simulation results [26], and we make use of it in our analysis. Specifically, in this estimate, DM particles are treated as traveling in a straight line from the halo to the detector, they are assumed to lose an average amount of energy in each scattering with nuclei along this path, and nuclear form factors are neglected. The condition for DM to be undetectable in a given experiment is then where E th χ is the lowest DM kinetic energy to which the detector is sensitive and depends on the detector energy threshold of the signal looked for (e.g. elastic nuclear scattering or the Migdal effect [144,145]), E i,max χ = m χ (v esc + v eff ) 2 /2 is the maximum DM kinetic energy in the halo, d is the detector depth, and where the integral runs over all nuclear species N = (A, Z) in the overburden with local number density n N (z). Applying the condition of Eq. (A1) and using Eqs. (16,17) to relate nuclear cross sections to a pernucleon value then leads to an upper exclusion of σ χn . When possible, we take upper exclusion limits from published experimental results, but when they are not given, we estimate them using the condition of Eq. (A1).
We apply a further rescaling to generalize limits on DM that makes up the full cosmic abundance to a species that is only a fraction f χ ≤ 1 of the total density. Given a lower exclusion of the per-nucleon cross section σ lo nχ , reducing the DM fraction just decreases the expected number of DM events and the rescaling is simply For upper exclusions of the per-nucleon cross σ hi χn , the dependence of the event rate on the DM-nucleon cross section is much steeper than linear [26]. Indeed, within the threshold approximation of Eq. (A1), the local DM density does not enter at all. As a result, we do not apply a rescaling by f χ to upper exclusions σ hi nχ . The collected exclusions derived using these methods for SI DM are summarized in Fig. 1 for f χ = 1, 10 −3 , 10 −6 , 10 −9 .
The prescription we use here can also be generalized to DM for which the connections between the nuclear and nucleon cross sections of DM do not follow exactly Eqs. (16,17). Given a bound on the per-nucleon cross section for minimal SI DM, it can be converted into a limit on the relevant nuclear cross sections and then recast as a limit on an effective per nucleon cross section for other scenarios. We use this approach to add the effect saturation of the nuclear cross section as described in Eq. (18), which turns out to be relevant only for some upper bounds from surface experiments and XQC. We also apply this method to derive limits on dark photonmediated DM that couples almost exclusively to protons. The limit on σ χp for such DP DM given a bound on σ χn for minimal elastic nuclear scattering of SI DM is then where η = η(v min ) is the usual halo function [4], F N is the nuclear SI form factor [75,76], and F χ is the DM mediator form factor given in Eq. (32). The exclusions on dark photon mediated DM derived with these methods are shown in Fig. 12 for DM fractions of This appendix contains our calculations of the DM energy distribution produced by scattering with hot SM particles in a thermal source that we use in Sec. IV. While we show the approximate distribution in Eq. (24), we use the more exact treatment below to generate the event rates that go into Fig. 7.
As a starting point, we consider strongly interacting DM that has thermalized with the local environment on Earth. We take it to be nonrelativistic with a Maxwell-Boltzmann phase space distribution, and a temperature T χ ∼ 300 K determined by the local thermodynamic temperature. Some of this DM encounters SM gas particles near a thermal source at a high temperature T g T χ . For our purposes we assume that these gas particles achieve thermal equilibrium with the thermal source. While they are comparatively hot, these particles' phase space distribution, f g , can also be approximated by a Maxwell-Boltzmann distribution as in (B1) with m χ → m g , T χ → T g . Given a matrix element for DM to elastically scatter on the gas of M, the rate of DM scattering is where n g is the number density of gas particles. k and p label the DM and gas momenta and primes denote final state quantities.
Using the delta functions to integrate over the unobserved final state gas momentum and reexpressing the initial momenta in terms of the center-of-mass velocity V and relative velocity v, allows us to write the rate as Γ = n g µ χg m g 8(2π) 5 V dV v 2 dv d cos θ k dk d cos θ × f χ (k)f gas (p)δ cos θ − k 2 + m 2 where θ (θ ) is the angle between V and v (k ). If we make the simplifying assumption that the matrix element is independent of the momenta then we can perform these integrations to find the total rate: The normalized flux distribution of the upscattered DM flux with respect to its kinetic energy, E χ = k 2 /(2m χ ), is then In the limit T χ → 0 this reduces to the expression in Eq. (24).