Saturated Overburden Scattering and the Multiscatter Frontier: Discovering Dark Matter at the Planck Mass and Beyond

We show that underground experiments like LUX/LZ, PandaX-II, XENON, and PICO could discover dark matter up to the Planck mass and beyond, with new searches for dark matter that scatters multiple times in these detectors. This opens up significant discovery potential via re-analysis of existing and future data. We also identify a new effect which substantially enhances experimental sensitivity to large dark matter scattering cross-sections: while passing through atmospheric or solid overburden, there is a maximum number of scatters that dark matter undergoes, determined by the total number of scattering sites it passes, such as nuclei and electrons. This extends the reach of some published limits and future analyses to arbitrarily large dark matter scattering cross-sections.


I. INTRODUCTION
While the presence of dark matter has been inferred from astrophysical observations and cosmological data, its nature remains enigmatic. Dark matter searches in the last few decades have sought out weakly-interacting massive particles, wimps, which scatter at most once as they pass through underground detectors. We show how ongoing underground experiments could find dark matter which scatters multiple times as it travels through these detectors. We refer to these multiply interacting particles as mimps, rather than wimps. We also show that a hitherto-neglected effect in dark matter studies, Saturated Overburden Scattering (sos), enhances published and prospective sensitivities to large dark matter scattering cross-sections.
Dark matter more massive than the so-called unitarity limit of ∼ 100 TeV arises naturally in grand unified theories which predict stable colored and electroweak states [1,2]. It has long been appreciated that supermassive, strongly-interacting, stable particles arise in supersymmetric models [3], and can be produced out of equilibrium in the early universe [4][5][6][7]. Depending on the reheating temperature, the equation of state in the primordial epoch, and any substantial increases of entropy in the early universe [8][9][10][11][12], super-heavy dark matter could make up the bulk of "missing mass" observed in galaxies. For low reheating temperatures, heavy and strongly-interacting dark matter could be a subdominant fraction of mass in the dark sector. Altogether, the discovery of super-heavy relic particles may provide the first evidence for supersymmetry, grand unified theories, and new particle dynamics prior to Big Bang Nucleosynthesis. This motivates the search for such particles in the data of existing experiments and with future detectors.
Here we investigate detection of such heavy, stronglyinteracting dark matter. In particular, in Section II, we detail the saturated overburden scattering effect, wherein dark matter scatters with every target along its path through the overburden. In Section III, we show the reach obtainable by existing dark matter direct detection experiments looking for multiply interacting dark matter, and evaluate the signatures of such mimps. Section IV shows how, in much of the mimp parameter space, the angle of entry into the detector can be used to discriminate backgrounds and validate signals, and to determine the mass, cross-section, and local density of mimps with a single experiment.

II. SATURATED OVERBURDEN SCATTERING
On its path to a detector, dark matter may be slowed by scattering with the atmosphere, Earth overburden, and detector shielding. Usually, as the dark matternucleon or dark matter-electron scattering cross-section is increased, the kinetic energy of dark matter arriving at a detector is decreased, due to more frequent scatters en route. Because experiments require the dark matter kinetic energy to exceed some minimum energy threshold to observe scattering, dark matter particles arriving too slowly become undetectable. This results in an upper limit on the cross-section for the sensitivity of such detectors. However, above a certain cross-section σ sos , the dark matter scatters with every nucleus or electron along its path. Increasing the cross-section beyond σ sos does not further dampen the particle's kinetic energy.
We thus find that the optical depth τ od , i.e. the approximate total number of recoils from a single dark matter particle passing through an overburden of length D, is correctly given as where n is the number density and n 1/3 the inverse distance between nuclei or electrons, and σ is the dark matter-nucleon or dark matter-electron scattering crosssection. The "min" function accounts for the fact that dark matter scatters at most once with each of the  Figure 1. The bound from the dama experiment on the spinindependent dark matter-nucleon scattering cross-section, before [13] (dashed blue) and after (solid yellow) accounting for saturated overburden scattering in Earth with density ρcrust = 2.7 g/cm 3 . In the region above the dotted (red) line, dark matter scatters with every nucleus in the atmosphere and Earth's crust on its way to the detector. Therefore, increasing the dark matter-nucleon cross-section beyond this threshold does not increase the stopping power of the detector overburden. This excludes previously-allowed parameter space, up to a limit from the cooling of gas clouds in the Milky Way [14]. n 1/3 D nuclei or electrons it passes. This treatment is valid so long as the typical momentum transfer is larger than the inverse distance between scattering sites. For super-heavy dark matter in the Milky Way with velocity ∼ 0.001c scattering elastically with nuclei or electrons, this assumption holds. Setting the arguments of the "min" function equal to each other, the cross-section for which overburden scattering is saturated is manifestly σ sos = n −2/3 . Note that this saturation cross-section is independent of the length of the overburden, and depends only on its density.
Typical dark matter direct detection experiments are conducted ∼ 1 − 2 km underground, where the Earth's crust constitutes the main overburden. We find that for spin-independent scattering σ sos = 1.3 ×10 −20 cm 2 , assuming a crust density of 2.7 g/cm 3 and a crust composition as given in Section IV. Spin-dependent scatters occur mainly off 27 Al, giving σ sos ∼ 10 −18 cm 2 . As illustrated in Figure 1, the effect of saturating the overburden by scattering with all nuclei along the dark matter's path modifies prior bounds [13], extending the range of crosssections probed to arbitrarily large values for very heavy dark matter. We leave the re-visitation of other bounds vis-a-vis sos to future work, which include those from the rrs balloon [15], the xqc calorimeter [16,17], cresst-i Al 2 O 3 experiments [18,19], as well as indirect bounds, e.g. from Earth heating [20]. We note that for a large enough nuclear scattering cross-section, dark matter prohibits the cooling of interstellar molecular gas clouds in the Milky Way [14], which however only becomes relevant for σ nχ 10 −10 cm 2 (m χ /(10 13 GeV)) in the spin-independent case, where m χ is the dark matter mass.

III. DARK MATTER SCATTERING AT HIGH MULTIPLICITY
For the high dark matter masses considered here, the probed cross-sections are so large that particles may scatter with higher multiplicity than the customarilyassumed single scatter. Numerically, this can be understood by examining the number of dark matter particles found by a detector, where Φ = (ρ χ /m χ )A det v χ t exp is the integrated flux of dark matter particles passing through the detector, A det the area of the detector, v χ the average dark matter velocity, and t exp the exposure (or observation) time. As dark matter passes through the detector, it will scatter a number of times roughly equal to the optical depth τ = n det σL det where L det and n det are the detector length and number density of scattering sites in the detector.
Limits on σ set by conventional searches are obtained in the single-scatter limit τ 1. Comparing Φ to τ , it is clear that bounds on the cross-section weaken with increasing dark matter mass, σ bound ∝ m −1 χ . For high enough dark matter mass and hence cross-section, this scaling breaks down as τ → 1, since transiting dark matter particles begin to scatter multiple times in the detector. In a typical single-scatter search, such events will have been cut away in the data analysis. As we will see, this breakdown occurs at a special point in (m χ , σ Tχ ) space where Φ ∼ 1 and τ ∼ 1. In the multiscatter limit, where all transiting dark matter scatters at least once (τ 1), it is the area of the detectors which determines mimp sensitivity. This can be compared to the singlescatter case, where the detector sensitivity is given by its volume or equivalently its total number of scattering sites.
If τ 1, so that all dark matter particles scatter at least once in a detector, by setting the flux of throughgoing dark matter particles to one, Φ ∼ 1, we obtain the maximum mass to which an experiment is sensitive as m max For higher masses, the experiment lacks enough exposure to see any transiting dark matter. Typically, every generation of dark matter experiments increases the area A det of the detector as well as the exposure time t exp , thereby increasing m max χ . Thus for each new detector, there will be high-mass dark matter candidates that can only be uncovered with a multiscatter search. This multiscatter frontier encompasses a large portion of parameter space that needs be considered by such experiments.

A. Prospects for Multiscatter Detection
There is a substantial but presently unexplored region of dark matter parameter space for which dark matter scatters multiple times while traversing terrestrial detectors. A number of prior studies have considered dark matter scattering multiple times while transiting detectors [13,[15][16][17][18][19][21][22][23][24][25][26][27] and astronomical bodies [20,21,28,29]. Since our findings can be easily applied to a given experiment with adequate resolution and sensitivity to identify multiple scatters, we simply consider two example technologies for dark matter detection: liquid xenon time projection chambers and bubble chambers. Multiply interacting dark matter could be found by analyzing data already collected at these experiments, such as the Xe detectors xenont [30], pandax-ii [31], and lux [32], and the bubble chamber experiment pico- [33] containing C 3 F 8 . We also show the improvement in sensitivity at the future 50-ton detector darwin [34] and the future pico upgrade, pico- [35]. These experiments are poised to uncover mimps in regions that are orders of magnitude beyond past efforts such as dama (using NaI crystals) [13], and reanalyses of edelweiss (Ge) and cdms (Ge, Si) [24], which were able to search for dark matter up to m χ 10 16 GeV.
For convenience, we define a mimp scattering threshold σ mimp = (L det n det ) −1 as the cross-section at which the Poisson expectation value for the number of recoils from a single dark matter particle in a given detector τ = 1. Given Poisson statistics, Poiss(multiplicity, τ ) = (e −τ τ multiplicity )/(multiplicity!), this happens when dark matter scatters at least twice during more than 25% of its detector transits. We adopt canonical definitions for the dark matter-nucleon cross-section in the spinindependent (si) and spin-dependent (sd) cases [36]. For a target with atomic number A, dark matter-nucleus reduced mass µ Tχ , and dark matter-nucleon reduced mass µ nχ , the si and sd per-nucleon cross-sections in terms of the per-nucleus cross-section are given by where the proportionality is in the m χ m T limit. For simplicity we will assume isospin-independent scattering and thus set a n = a p = 1. Figure 2 shows the order-of-magnitude sensitivities achievable by mimp searches for events with high multiplicity, assuming a local dark matter density of 0.3 GeV/cm 3 . In practice, regions covered by some searches where τ 1 overlap due to Poisson fluctuations. Regions to the left of the solid green (dashed green) lines can be probed with 1-year (10-year) run-time. This mass reach is determined by requiring that the integrated flux Φ over the exposure time t exp is at least 2.3 events (corresponding to the 90% confidence upper limit for zero observed events).
Single-scatter sensitivities σ single are shown for comparison, and are separated from mimp regions by horizontal green lines at the threshold cross-section σ MIMP . These single-scatter sensitivities are taken from Refs. [30][31][32][33]37] and are rescaled using the following parameters: for {xenont, darwin, pico-, pico-}, we simplify the detector as having an effective length L det = {96 cm, 250 cm, 50 cm, 100 cm}. We assume densities of 2.94 g/cm 3 for Xe and 1.36 g/cm 3 for C 3 F 8 , and isotopic abundances of 27% 129 Xe, 4% 130 Xe, 21% 131 Xe, 27% 132 Xe, 10% 134 Xe, and 9% 136 Xe, and 100% of C 3 F 8 . Following [38] we take (J A , S p , S n ) = ( The relative behavior of these curves may be understood from their scalings with detector length, where for simplicity we here take a spherical geometry for the detector. The multiscatter threshold cross-section then scales inversely with the effective diameter of the detector, σ MIMP ∝ L −1 det . The maximum reachable mass scales with the effective target area and the exposure time, m max χ ∝ L 2 det t exp . As usual, the single-scatter crosssection sensitivity σ single ∝ (L 3 det t exp ) −1 . We point out that depending on their actual exposure times, detectors such as xenonnt [39], lz [40], darwin and pico- are capable of probing dark matter masses of the order 10 19 GeV, i.e. up to and beyond the Planck scale.

B. Signatures at High Multiplicity
Perhaps the most salient feature of a super-heavy, multiply interacting, stable particle is that it will leave a mostly collinear track of nuclear recoils as it passes through a meter-scale detector. Comparing the maximum detector-frame scattering angle for mimps, sin α max = m T /m χ , with the maximum number of recoils in the detector n 1/3 det L det , we find that in the m χ m N limit, the maximum total deflection angle of a dark matter particle is Ω max n 1/3 L det sin α max , or Thus, the nuclear recoils of a transiting mimp are typically collinear, although for m χ 10 13 GeV, degree-level deflection along the length of a detector becomes plausible. We remark that, since solar and atmospheric neutrinos will not leave tracks in detectors, the "neutrino floor" [41,42] is less of a concern for mimp searches. Detectors such as lux, pandax-ii and xenont are liquid xenon time projection chambers (tpcs) where, for  Figure 2. Prospects for probing dark matter scattering at high mass and high multiplicity, given the examples of liquid xenon or bubble chamber experiments. The (diagonal, blue) lines indicate the usual sensitivity from a zero-background single-scatter search. However, above the (horizontal, green) lines labeled σMIMP, more than 25% of the dark matter traversing the detector scatters more than once, requiring a dedicated multiscatter analysis to probe this parameter range. The ultimate high-mass reach is given by the (vertical, green) lines from the requirement that at least 2.3 dark matter particles traverse the detector for the stated effective detector area and exposure time. For detectors such as darwin, this limit can lie beyond the Planck mass ∼ 10 19 GeV. each interaction, two signals are observed with photomultiplier tubes: an O(10 ns) pulse ("S 1 ") from scintillation in the scattering target, followed by an O(µs) pulse ("S 2 ") from electrolumiscence of electrons that have drifted into the gas above the target liquid. The drift time is O(1 ms), allowing for clear separation of S 1 and S 2 . In comparison to these timescales, dark matter transits a 1 m detector length in ∼ 5 µs. The relative strength of S 1 vs S 2 helps distinguish dark matterinduced nuclear recoils from electronic recoils from β and γ radiation, which comprises the main background.
A mimp transiting such a detector would produce multiple S 1 's and S 2 's, each characteristic of a nuclear recoil of relatively high energy in the range of 10's of keV. Whether the pulses would appear individually or merged is determined by the timing between successive scatters, which is typically the transit time divided by the number of recoils τ . To understand the basic signatures let us neglect the fact that tpcs have a much better resolution along their symmetry axis compared to the horizontal. Then, for τ 5 (τ 500) the S 2 (S 1 ) pulses merge into elongated pulses S 2 (S 1 ). There are thus three qualitatively distinct mimp signatures: (1) a series of S 1 s followed by a series of S 2 s, for 1 τ 5; (2) a series of S 1 s followed by a merged S 2 , for 5 τ 500; and (3) an elongated S 1 followed by an S 2 , for τ 500, where the S 1 and S 2 will overlap at least partially. Backgrounds to mimp scattering will be exceedingly small. Signature (1) at small multiplicity can be mimicked by the pile-up of individual single-scatter background events, which however will happen predominantly at the surface of the detector and thus can be fiducialized. Another potential background in this regime is from fast decays such as the 214 BiPo coincidence which occurs in the 222 Rn decay chain [43,44], but will be of little concern since the alpha-decay usually deposits much more energy than expected from mimps. At intermediate multiplicity, radiogenic neutrons might mimic the expected mimp signature, but they do not usually travel at non-relativistic speeds nor scatter along a straight line. A background to elongated S 2 events comes from instrumental sources of drifting electrons [45,46], which will not usually conspire with multiple S 1 events to mimic the required mimp signature. At large multiplicity, tracks from through-going muons will display much shorter S 1 pulses than expected from mimps, deposit much more energy than expected from dark matter, and can typically be vetoed by means of active shielding.
At pico-, the energy deposited by dark matter scattering off superheated C 3 F 8 nucleates a bubble, identified visually with stereo cameras and acoustically with piezoelectric transducers. A mimp nucleates multiple bubbles during its ∼2.5 µs-long transit through the detector volume. The acoustics, sampled at a rate of 0.4 µs −1 , can identify at most 2.5/0.4 ∼ 6 individual bubble nucleations, whereas the cameras could image up to 100 bubbles. With the background from radioactive neutrons limited [33], the visual position reconstruction makes pico well-suited for a mimp search.

IV. DETERMINING THE MASS, CROSS-SECTION, AND LOCAL DENSITY OF MIMPS
In some regions of parameter space it is possible to determine the mass, scattering cross-section, and local density of mimps using the directions of their arrival. To explain this, we must first detail our overburden calculations. The depletion of dark matter's kinetic energy after passing through a path D, along which lie nuclei with number densities n i and masses m i is where E i and E f are the respectively initial and final dark matter kinetic energies and zβ i = z 4m i m χ /(m χ + m i ) 2 is the fraction of kinetic energy lost in each scatter, where z is a kinematic factor set to 1/2 for evenly distributed center-of-mass scattering angles [29]. To calculate this approximate upper bound on a detector's sensitivity to dark matter-nucleon scattering, we require that the final kinetic energy E f of a mimp, which travels at typical Galactic velocities of 10 −3 c, still has enough kinetic energy to deposit 1 keV when recoiling against the target nucleus in an elastic collision.
To determine the number density of nuclei in the path of the dark matter as it transits either from above through the crust, or from below through the crustmantle-core system, we model the interior of Earth using elemental abundances given in [47,48]. Specifically, we model the Earth as a sphere, and take Earth's crust to  Figure 3 shows the spin-independent dark matter-nucleon crosssections above which the dark matter's kinetic energy is diminished so much that it becomes undetectable, for a path through 2 km of the crust or for the opposing radial path through the Earth, passing through 62 km of crust, The solid lines indicate the boundaries above which the mimp has lost too much of its kinetic energy propagating through the overburden, either from above through 2 km of Earth crust or from below through 12756 km of crust, mantle and core. In the parameter space between these lines, it is possible to use the acceptance cone of the detector to pinpoint the mass, scattering cross-section, and local density of mimps. Dashed lines indicate the typical number of recoils from a single dark matter particle passing through xenont or pico-. The saturated overburden scattering for Earth with density 2.7 g/cm 3 is indicated with a dotted line. 5734 km of mantle and 6960 km of core. Spin-dependent scattering occurs off Al, Na, Ca, and 29 Si, with nuclear spin expectation values ( S p , S n ) given in [38].
With the overburden specified, we now show how mimp directionality can be used to determine the mass, crosssection, and local density of mimps with a detection in just a single detector. As a preliminary caution, we note that while typical directional dark matter detection effects rely on a spatial anisotropy of incoming dark matter velocities (the dark matter velocity towards Earth is peaked in a direction pointing toward Cygnus), the effect we introduce here does not depend on dark matter velocity anisotropy. This directional effect would occur even if the dark matter velocity were perfectly isotropic at the surface of Earth, and arises simply because the overburden above and below a direct detection experiment is different.
The directional analysis we propose proceeds as follows: first, the number of recoils from a single mimp in the detector gives a measurement of the optical depth τ = n det σ T χ L det , where the detector length L det and nuclear density n det are known. Thus the mimp scattering cross-section can be determined directly, and independently of the local dark matter density, by measuring the number of recoils per mimp event along the length of a detector. Next, for parameter space between the two diagonal overburden lines in Figure 3, some fraction of the mimps will not arrive at the detector. In fact, the arriving flux will have an angular dependence, since upward-going mimps will pass through more Earth overburden. More precisely, the mimp mass and cross-section uniquely determine the "angle of acceptance" at a direct detection experiment, defined as the angle between a vector pointing in the normal direction from the detector to the Earth's surface, and a vector pointing from the detector towards where the dark matter flux through the Earth is halved. This halving of flux will occur near an angle where all dark matter particles coming up through the Earth are blocked, which depends on dark matter's mass and velocity distribution; see Equation (5). After appropriately modeling the dark matter velocity distribution and the Earth's interior, it would be possible, using Equation (5), to pinpoint the dark matter mass. (This is in contrast to usual single-scatter wimp analyses where at high masses the spectral shape carries no information about the wimp mass.) With the mimp flux measured by the experiment (Φ = (ρ χ /m χ )A det v χ t exp ) this is enough to determine the local mimp mass density.
A more detailed simulation would also incorporate the orientation of the dark matter speed, which is roughly maximal towards the Cygnus constellation, relative to the dark matter detector. Intriguingly, since the acceptance angle depends slightly on the dark matter velocity distribution, it is larger when the detector is so oriented that the higher-velocity mimps from Cygnus are traveling through more of Earth. A full analysis thus requires keeping track of the detector's varying orientation with respect to the dark matter wind due to Earth's rotation and orbit around the Sun, as well as the arrival times of observed mimp events.

V. OUTLOOK
Super-massive relic particles with per-nucleon crosssections in the 10 −40 − 10 −20 cm 2 range are motivated by grand unified, supersymmetric, and other theories which contain stable heavy particles charged under Standard Model gauge symmetries. We have demonstrated that multiply scattering massive particles, with masses in the 10 10 − 10 20 GeV range, require dedicated analyses, and may be discovered in data already collected by dark matter direct detection experiments, along with future experimental searches. As explained in Section III, it is no coincidence that the multiscatter frontier opens up just as the flux limit is approached for a given detector: as dark matter particles become rarefied with increasing mass, a detector has to be able to interact with every transiting particle for a meaningful signal. For cross-sections larger than this, dark matter will scatter multiple times in the detector.
The new phenomenology presented in this article reveals many avenues for future research. Previously published limits on dark matter at large cross-sections must be re-derived after taking saturated overburden scattering into account. The sensitivity of dark matter detectors to leptophillic mimps scattering multiple times with target electrons must be explored. A directional analysis for validating mimp signals at direct detection experiments incorporating a more detailed model of the Earth's interior can be undertaken. Dropping the assumption of simple spherical detectors and taking into account the direction of drift in liquid xenon tpcs will result in diurnal modulation signatures. mimp searches at detectors such as Borexino [49], sno+ [50], juno [51], and mathusla [52] should be investigated. Finally, the links between grand unified theories with stable colored and electroweak states, the cosmological production of the same, and terrestrial searches for mimps should be sought out. Altogether, the search for dark matter in multiple-scatter events has the potential to revolutionize our understanding of fundamental physics and the early universe.