Earth-Scattering of super-heavy Dark Matter: updated constraints from detectors old and new

Direct searches for Dark Matter (DM) are continuously improving, probing down to lower and lower DM-nucleon interaction cross sections. For strongly-interacting massive particle (SIMP) Dark Matter, however, the accessible cross section is bounded from above due to the stopping effect of the atmosphere, Earth and detector shielding. We present a careful calculation of the SIMP signal rate, focusing on super-heavy DM ($m_\chi \gtrsim 10^5 \,\,\mathrm{GeV}$) for which the standard nuclear-stopping formalism is applicable, and provide code for implementing this calculation numerically. With recent results from the low-threshold CRESST 2017 surface run, we improve the maximum cross section reach of direct detection searches by a factor of around 5000, for DM masses up to $10^8 \,\,\mathrm{GeV}$. A reanalysis of the longer-exposure, sub-surface CDMS-I results (published in 2002) improves the previous cross section reach by two orders of magnitude, for masses up to $10^{15} \,\,\mathrm{GeV}$. Along with complementary constraints from SIMP capture and annihilation in the Earth and Sun, these improved limits from direct nuclear scattering searches close a number of windows in the SIMP parameter space in the mass range $10^6$ GeV to $10^{13}$ GeV, of particular interest for heavy DM produced gravitationally at the end of inflation.


I. INTRODUCTION
Direct detection experiments aim to detect Dark Matter (DM) by measuring the energy of nuclear recoils, induced by the interaction of DM in the Galactic halo with nuclei in the detector [1,2]. As standard direct detection experiments probe smaller and smaller cross sections for weak-scale DM masses [3,4], low-threshold experiments [5,6] and DM-electron scattering searches [7] are beginning to search for DM with GeV-scale masses and below. With no confirmed detection thus far, we are driven to search as much of the DM parameter space as possible. Indeed, there are well-motivated particle DM candidates ranging from the very light [8,9] to the very heavy [10][11][12][13], interacting with the Standard Model very weakly [14][15][16] or very strongly [17][18][19]. In this paper, we revisit direct detection constraints on such strongly-interacting DM, with a focus on super-heavy candidates.
There has been much recent interest in signals and constraints for DM which interacts strongly with nuclei and whose distribution may therefore be affected by scattering in the Earth [20][21][22]. Much of this has been focused on light (sub-GeV) DM where lower limits from conventional direct detection experiments weaken, meaning that moderate DM-nucleon scattering cross sections have not yet been excluded [23][24][25][26][27][28][29][30]. However, care must be taken in the treatment of DM scattering in the Earth. Light DM may be deflected substantially when it scatters with nuclei, drastically changing the incoming direction and path length [31][32][33]. However, the recently developed DaMaSCUS Monte Carlo code [34] accounts for these effects and can be used to reliably predict signals of strongly-interacting light DM over a range of cross sections [35]. * b.j.kavanagh@uva.nl Heavier DM, on the other hand, is already strongly constrained by direct detection experiments. But constraints typically cut off at large values of the DMnucleon cross section, where scattering in the Earth slows the DM particles before they reach the detector [36,37], rendering them undetectable. At even larger cross sections, airborne experiments provide constraints, as the less dense atmosphere is not as effective at stopping DM particles [30]. A number of small windows remain open between these constraints for DM masses larger than ∼ 10 5 GeV, in the mass range favoured by gravitational production at the end of inflation [38][39][40]. Though these regions are disfavoured by complementary constraints from considering DM capture and annihilation in the Earth and Sun [41][42][43], it is important to cross-check using direct searches, which rely on an entirely different set of assumptions. In this work, we demonstrate that the maximum cross section excluded by existing direct detection experiments can be increased by more than three orders of magnitude, excluding these remaining windows using nuclear scattering experiments alone. We also provide a numerical code, verne [44], which can be used to calculate the Earth-stopping effect for heavy DM and derive lower limits on the DM-nucleon cross section.
We focus here on two experiments: the CRESST 2017 surface run [6], operated at the Max-Planck-Institute for Physics (MPI) Munich in 2017, and CDMS-I [45,46], operated at the Stanford Underground Facility (SUF) in 1998-1999. The CRESST 2017 surface run provides excellent sensitivity to strongly interacting massive particles (SIMPs) because of its low threshold of 20 eV, as well as the fact that it was operated in a surface building, meaning that there is limited shielding by the Earth. We extend the SIMP limits presented in Ref. [30] to large masses. However, the short (∼ few hours) exposure of the CRESST 2017 surface run means that it is not sensitive to DM masses larger than ∼ 10 8 GeV, when the arXiv:1712.04901v4 [hep-ph] 30 Jun 2018 number of DM particles crossing the detector over the course of the exposure becomes too small to generate an appreciable signal. Instead, the CDMS-I experiment had a longer (∼ 1 year) exposure and therefore probes masses up to ∼ 10 15 GeV. Limits on SIMPs from CDMS-I were derived in Ref. [37] and we present a refinement of those limits, taking into account form factor suppression of the SIMP scattering underground, as well as modelling the full (anisotropic) velocity distribution of incoming SIMPs.
In Sec. II, we present the general framework of direct detection, then in Sec. III, we describe the formalism for calculating the final DM velocity at the detector, after propagation through the atmosphere, Earth and detector shielding. We focus on heavy DM (m χ 10 5 GeV) with large scattering cross sections, which allows us to make two simplifications. First, DM with a cross section ∼ 10 −28 cm 2 would typically scatter N scat ∼ O(5000) times before reaching the CDMS-I detector at SUF. Thus, we can approximate the energy losses as continuous over the path of the particles. In addition, for each scatter with a nucleus of mass m N , the DM is deflected by an angle δθ ∼ m N /m χ 10 −3 [31]. The total deflection of the DM particle is then expected to be ∆θ ∼ δθ/ √ N scat O(10 −5 ). We can thus approximate the trajectories of DM particles as straight lines, allowing us to use the standard nuclear stopping formalism and avoiding the need for complicated Monte Carlo calculations.
With this formalism in hand, we derive the final velocity distribution at the detectors in Sec. IV, followed by the expected signal rate and resulting constraints in Sec. V. The final constraints are shown in Fig. 5, for DM masses in the range 1 GeV to 10 15 GeV, with the caveat that at low masses these results need to be supplemented with dedicated Monte Carlo simulations. All code for performing the calculations in this paper, along with all results and plots, is made freely available online, along with the numerical code verne, at https://github.com/bradkav/verne [44].

II. DIRECT DETECTION OF SIMPS
In order to determine the maximum cross section probed by a particular experiment, we must compare the number of observed events with the number of expected SIMP signal events. For DM scattering with a nucleus N , the rate of nuclear recoils of energy E R per unit detector mass is given by [47], Here, ρ χ is the local DM density, which we fix to the benchmark value of 0.3 GeV cm −3 [48], noting that there is some uncertainty on the true value [49]. For the differential DM-nucleus cross section dσ χN /dE R we assume standard spin-independent (SI) scattering. We integrate over the DM velocity distribution f (v), including the contribution of all particles with sufficient DM speed to excite a recoil of energy E R . That is, we include all particles with v > v min = m N E R /(2µ 2 χN ), where µ χN is the DM-nucleus reduced mass.
For weakly interacting massive particles (WIMPs), the velocity distribution is typically assumed to have the form: 1 This velocity distribution defines the so-called Standard Halo Model (SHM), for which we assume an escape speed in the Galactic frame v esc = 533 km/s [50] and dispersion [48]. The mean DM velocity is given by v χ = −v lab (t), the velocity of the lab with respect to the Galactic rest frame. The lab velocity is time-varying, arising from the motion of the Sun around the Galactic centre, the orbit of the Earth around the Sun as well as the daily revolution of the Earth [51,52].
Equation (2) corresponds to the initial velocity distribution of DM particles arriving at Earth and will be modified if DM undergoes substantial interactions in the atmosphere or in the Earth itself. This modification depends on the path travelled by the DM particles or, equivalently, the position of the detector on Earth with respect to the DM flux. It is useful to define an angle γ, defined as the angle between the position vector of the detector r det (as measured from the centre of the Earth) and the direction of the mean DM flux: With this definition, the mean DM flux appears to come from directly overhead for γ = 180 • while the flux appears to come from below (through most of the Earth) for γ = 0 • . Explicit expressions for the velocity distribution in terms of the angle γ are given in Appendix A. In Fig. 1, we plot the value of γ for the Stanford Underground Facility (SUF, latitude 37.4 • N, longitude 122.2 • W) for a 1-year period covering the CDMS-I exposure presented in Ref. [45,46]. The slow variation due to the Earth's orbit and the fast oscillations from the Earth's rotation are both visible. On average, the flux of DM appears to come from above, at an angle ∼ 53 • off the vertical. This suggests that the typical path length for DM particles arriving at SUF will be relatively short, as the particles do not need to cross the entire Earth to reach the detector. The variation of γ during the exposure means that DM particles will travel different path lengths through the atmosphere and Earth at different times. This in turn leads to a variation in the stopping effect and a time variation in the velocity distribution f (v, γ(t)) of particles reaching the detector. The (time-varying) rate of recoils in the detector is then: where M is the fiducial mass of the experiment, (E R ) is the nuclear recoil detection efficiency 2 and the analysis window of the experiment is E R ∈ [E min , E max ]. In the case of the CRESST 2017 surface run, we calculate the expected number of signal events by integrating over the full 5.31 hour exposure on the night of 16th February 2017 (see the ancillary material to Ref. [6] for more details). For CDMS-I, we do not know the exact times of the exposure, so we integrate the event rate over the course of a typical day, then multiply by the total number of exposure days: We select this typical day as 21st April 1999 (somewhere mid-way through the full exposure).
While γ = 180 • leads to the shortest path length for incident DM particles and therefore the largest rate, the detectors will generally not spend a long time oriented such that γ = 180 • . It is therefore important to calculate the speed distributions for a range of values of γ and take all of these into account. In the next two sections, we explain how to calculate the DM velocity at the detector for a given trajectory, followed by the procedure for evaluating the speed distribution at the detector, accounting for all possible trajectories.

III. NUCLEAR STOPPING OF DARK MATTER
In order to determine the final velocity of DM particles at the detector, we use the 'nuclear stopping' approach, first presented in Ref. [36]. Our treatment follows closely a number of references, including Ref. [30]. However, there are a number of refinements to the standard picture which we will make particular note of. For now, we keep the DM mass m χ general, before focusing on the case of super-heavy DM (m χ m A ). As DM particles traverse a medium, consisting of different nuclear species which we label with the index i, they may interact and lose energy. As we argued in the introduction, we can assume that this scattering is continuous, meaning that the rate of change of the average DM energy is given by: where n i (r) is the number density of nuclei of species i at position r, σ i (v) is the total DM-nucleus scattering cross section and v is the DM speed. The average recoil energy (or equivalently the average change in DM energy) for scattering with nucleus i is given by: Here, E max i = 2µ 2 χi v 2 /m i is the maximum recoil energy which can be transmitted to nucleus i through elastic scattering.
For standard spin-independent (SI) scattering, the differential cross section is given by [53,54]: The strength of the interaction is parametrised by the DM-proton cross section at zero momentum transfer σ SI p and the total cross section is coherently enhanced by the number of nucleons in the target nucleus, A 2 i . The form factor F 2 i (E R ) (which we take to have the Helm form [55,56]) arises due to the finite size of the nucleus and suppesses the cross section at large recoil energies. For the SI interaction, the mean energy loss is then given by: The coherence factor C i (m χ , v) reflects the suppression of the mean recoil energy (relative to that expected for a point-like nucleus) owing to the nuclear form factor. For very light DM (and for point-like nuclei) C i (m χ , v) → 1 meaning that this factor is typically neglected in studies of the Earth-stopping of sub-GeV mass DM [30]. For heavy DM and heavy target nuclei, however, the effects of the form factor can be large. For a DM particle of mass 10 5 GeV scattering off lead, for example, C Pb (m χ , v) drops below 1% for DM speeds larger than around 200 km/s. This reduces the shielding efficiency of lead (and the stopping efficiency of Earth elements) for fast, heavy DM. It is therefore imperative that we include it here.
We now examine in more detail the argument given in the introduction that the scattering can be described as continuous and that ultra-heavy DM particles should travel in a straight line. The typical number of scattering events which a DM particle experiences is: where D is the distance travelled. This corresponds to a mean free path of λ ≈ 2 mm. We have taken here typical number densities in the Earth's crust and assumed that the DM is much heavier than the most abundant elements in the crust. For a detector 10 m underground, DM particles with a cross section of σ SI p = 10 −28 cm 2 will scatter around 5000 times with a roughly gaussian error of √ N scat ∼ 70. Treating the stopping as continuous (rather than as a series of discrete scatters) should therefore introduce an error of O(1%). As we will see, we are able to exclude DM particles with even larger cross sections that this, in which case the approximation only improves.
For very light DM particles, the angle α (measured in the lab frame) by which their path is deflected after a scatter is uniform. As we increase the DM mass, however, α becomes increasingly peaked in the forward direction (see Sec. 3.3.3 of Ref. [31] for a more detailed discussion). In the heavy DM limit, we find that α ∈ [0, m i /m χ ]. Over many scatters, the typical total deflection of the particle is For the case of a 10 5 GeV particle, scattering 5000 times, this corresponds to a typical angular deflection of O(10 −5 ) radians (or a 0.1 mm deflection over 10 m of travel). With increasing DM mass, the typical angular deflection becomes smaller and smaller.
In light of this, we assume from now on that the DM particles travel in straight line trajectories.
Writing the DM energy as E χ = m χ v 2 /2, we can rewrite Eq. (6) in terms of the DM speed v and the straight-line distance travelled by the DM particle D: where on the second line we have taken the limit m χ m i , for all target nuclei of interest. For the numerical results of the paper, we use the full expression on the top line of Eq. (11), but it is useful to note that for very heavy DM, the stopping effect scales as σ SI p /m χ . This means that once we have calculated the final velocity distribution at the detector for a given (large) DM mass and cross section, the results for another DM mass can be obtained trivially by rescaling the cross section appropriately.
Because of the non-trivial dependence of the coherence factor C i (m χ , v) on the DM speed, it is not possible to solve Eq. (11) analytically. Instead, it must be solved numerically, starting from some initial speed v i and integrating forward along the trajectory of the DM particle to obtain the final speed at the detector v f . Each DM particle is propagated through three regions, in order: • Atmosphere -particles are propagated along a straight line, towards the detector, from the top of the atmosphere to the surface of the Earth. We include stopping by Oxygen and Nitrogen nuclei and take the atmospheric density profiles from the International Standard Atmosphere [57], which extends up to a height of 80 km.
• Earth -particles then propagate in a straight line from the surface of the Earth to the detector. We include 8 different Earth elements -O, Si, Mg, Fe, Ca, Na, S, Al -using density profiles tabulated in Ref. [58] (based on data from Refs. [59,60]).
• Shielding -finally, the particles propagate through any shielding which surrounds the detector. In the case of the atmosphere and the Earth, the path length depends on the angle of incidence, while for the shielding we assume a constant path length, regardless of incoming direction.
Explicit expressions for the path length through the atmosphere and Earth are given in Appendix A.
For the CDMS-I detector as SUF, we assume that the detector is located 10.6m underground and is surrounded by 16cm of lead shielding 3 . For the CRESST 2017 surface run at MPI, the only substantial shielding comes from the walls of the surface building itself, consisting of 30cm of concrete. We model this by assuming that the detector is positioned at a depth of 30 cm underground. In addition, we include shielding from 1mm of copper surrounding the detector.
In Fig. 2, we show an example of the final DM speed v final as a function of the initial DM speed v initial for a particular trajectory, namely particles travelling vertically downwards from the top of the atmosphere to the detector. In the left panel, we show results for a detector at SUF and in the right panel for a detector in a surface building at MPI. In both cases, we assume a DM mass of 10 5 GeV, but note that in the left panel the DM-nucleon cross section is roughly 20 times larger than in the right.
For a detector at MPI (left panel), we note that the atmospheric crossing (blue line) has the largest slowing effect, as described in Ref. [30]. Crossing the 30 cm of concrete which makes up the walls of the surface building (green line) also gives a reduction in speed. Neglecting the coherence factor C i (m χ , v) in Eq. (11), we would expect v final ∝ v initial , as is usually assumed when considering nuclear stopping. Above around 400 km/s, however, the final speed after crossing the Earth is higher than expected from this simple proportionality. This is due to the coherence factor C i (m χ , v). More energetic particles typically impart much less than the maximum possible recoil energy, especially when scattering with heavy nuclei. Finally, the particles must cross 1mm of copper (orange line), which gives only a small reduction in the speed.
For the case of a detector at SUF (right panel), we consider a smaller DM-nucleon cross section, meaning that the atmosphere (blue line) has only a small effect on the slowing of DM. Instead, the DM particles must cross 10.6m of earth to reach the detector, leading to substantial slowing at low speeds (green line). Again, high speed particles are slowed less and this is particularly noticeable after the particles have also crossed the lead shielding (orange line), where only the fastest-moving DM particles survive. Neglecting the coherence factor, we would expect that these fast-moving DM particles would be effectively stopped and would not be detectable in an underground experiment.
In this section, we have assumed that the energy losses of the particles are continuous and are well-described by the average number of scatters. Indeed, some particles will undergo more or fewer scatters than this average, introducing some variation into the final velocity of the particles. For a particle scattering around 5000 times (see Eq. (10)), the typical error on the number of scatters is O(1%). We can accommodate this into our formalism as an uncertainty on the distance travelled by the particle D ∼ N scat λ. In the right panel of Fig. 2, we include a shaded green band, showing the effect of varying the detector depth by ±4% (or roughly a 3σ variance in the number of underground scatters). The resulting error on the DM velocity is around 5%, further supporting the assumption that the energy losses can be treated as continuous in this many-scatter regime.

IV. VELOCITY DISTRIBUTION AT THE DETECTOR
We have so far discussed the propagation of a given DM particle to the detector. Now we examine the impact on the total population of DM particles to obtain the final speed distribution at the detector 4 . We assume that the initial velocity distribution f (v i ) (at the upper atmosphere, before scattering) has the Maxwell-Boltzmann form given in Eq. (2). The final distribution of velocities at the detectorf (v f ) is obtained by a change of variables: where in passing to the last line we have used thatv f = v i ≡v, because the DM particles are assumed to travel in straight lines.
The speed distribution at the detector is then: Here, we understand v i to be the initial speed required to obtain a final speed of v f for particles travelling along the trajectory specified by the directionv. We note also that v i and the derivative dv i /dv f depend on the incoming DM direction (and so must be included inside the integral over incoming angles). Evaluating the final speed distributionf (v f ) at the detector then requires us to integrate Eq. (11) backwards from the detector to the top of the atmosphere. An alternative approach would be to evaluatef (v f ) by a full Monte Carlo simulation, following the trajectories of a large number of particles between individual scattering events. However, given the large number of scatters expected, O(1000) or more, such simulations could be prohibitively slow. A simplified Monte Carlo simulation, as performed in Ref. [30], is also possible, in which the final distribution of velocities is found from generating a sample of initial velocities and solving Eq. (11) for each. From Eq. (13), we see thatf (v f ) is typically dominated by trajectories where dv i /dv f is large. The 'backpropagation' method we use automatically takes care of 4 We could also explore the full 3-dimensional velocity distribution at the detector, which would be relevant for directional detectors. However, we leave this for future work. the proper weighting by the derivative dv i /dv f (which is calculated numerically), thus minimising the number of trajectories which must be evaluated.
For given values of m χ , σ SI p and γ, we evaluatef (v f ) over a range of values for v f , from 1 km/s up to v max , the maximum speed which a DM particle can have when it reaches the detector (maximised over all incoming directions). As a cross-check, we can verify that this procedure maintains the normalisation of the distribution function.
We integratef (v f ) in the range v f ∈ [1 km/s, v max ] and add to this the fraction of the initial particles which have a speed lower than 1 km/s when they reach the detector 5 . We find that the distribution is always normalised to unity, typically to better than 1 part in 10 3 . Figure 3 shows the final DM speed distribution at a detector at SUF for a range of values of γ, the direction of the average DM flux. The attenuation effect becomes more pronounced as we go from γ = 180 • to γ = 0 • . In the former case, DM particles have a shorter path on average to reach the detector, while in the latter, the typical DM particle must cross the entire Earth to reach the detector. We note, however, that γ denotes only the average incoming DM direction; the DM particles have a distribution of incoming directions. This is why the γ = 0 • case does not lead to complete attenuation: some particles still arrive at the detector from roughly 5 The initial velocity which gives a final speed v f less than 1 km/s is obtained by back-propagating particles of v f = 1 km/s for a range of incoming directions.
overhead (though these typically have smaller speeds in the first place).
In Fig. 4, we show the final speed distribution as a function of the interaction cross section for detectors at both MPI (left) and SUF (right), fixing the mean incoming DM direction to γ = 126 • (a typical value for both detector sites). As expected, increasing the DM-nucleon cross section reduces the maximum speed of particles which arrive at the detector and increasingly populates the low-speed tail of the distribution. For a detector near the surface (left), the final speed distribution typically resembles a 'compressed' version of the initial distribution. This is because the final speed scales roughly linearly with the initial speed as demonstrated in the left panel of Fig. 2. The result is a roughly exponential cut off in the final speed distribution, with the maximum DM speed decreasing with increasing cross section.
In the case of SUF (right panel of Fig. 4), the final speed distribution is not simply a rescaled Maxwell-Boltzmann. Instead, slower particles tend to be stopped almost completely (populating the low-speed tail) while faster particles experience comparatively less stopping (due to the coherence factor C(m χ , v)), leading to a bump at higher speeds. Because the speed distribution does not drop as rapidly as for a detector at MPI, we expect that a large rate should be observable by a detector at SUF, so long as the maximum final speed lies above the threshold speed of the experiment (shown as a grey band).

V. CONSTRAINTS
In order to determine the maximum cross section which can be probed by the CRESST 2017 surface run and the CDMS-I run at SUF, we must compare the observed number of recoil events with the predicted signal, given by Eq. (5), evaluated by performing a grid scan over m χ and σ SI p . Assuming (conservatively) that all observed events could be signal, we set a single-bin Poisson limit at the 90% confidence level on the minimum value of the cross section. Given the density of the grid scan, we expect these lower limits to be accurate at the level of 2-3%.
For CDMS-I, we consider an analysis window for nuclear recoil energies of E R ∈ [10, 100] keV and use an estimated Gaussian energy resolution of 2.4 keV [46]. The Germanium detector modules at CDMS-I observed 27 veto anti-coincident events, consistent with single-scatter and multiple-scatter nuclear recoils. The total exposure was M T = 15.8 kg days, obtained over roughly 1 year in 1998-1999. In CDMS-I, the signal rate varies rapidly with cross section; the 'bump' in the final speed distribution (shown in the right panel of Fig. 4) means that the signal rate is very large as long as the maximum DM speed lies above the experimental threshold. For a DM mass of 10 5 GeV, we limit the cross section to be larger than σ SI p = 1.55 × 10 −28 cm 2 . We have also checked that in this case, the total energy deposited by a single particle in a 1cm-thick Ge detector module (due to multiplescattering) is O(20 keV), meaning that such a particle should fall within the analysis window of the experiment and be in principle detectable.
For the CRESST 2017 surface run, using an Al 2 O 3 target, we assume a total exposure of 4.6 × 10 −5 kg days, an energy window of E R ∈ [19.7, 600] eV and a Gaussian energy resolution of 3.74 eV at threshold. In this exposure, 511 events were observed in the region of interest. Because of the exponential cut off of the speed distribution in this case (see the left panel of Fig. 4), the number of signal events varies more slowly with cross section. For a DM mass of 10 5 GeV, we limit the cross section to be larger than σ SI p = 4.5 × 10 −27 cm 2 . As for CDMS-I, the limit for larger masses can be obtained by simply rescaling this cross section by a factor m χ /(10 5 GeV). Figure 5 shows the resulting constraints from CDMS-I in blue and from the CRESST 2017 surface run in red. For CDMS-I, we extend the limit up to DM masses of O(10 15 ) GeV. At such high masses, the number density of DM particles is very low and the rate of DM particles crossing a 10 cm-scale detector is O(1) per year. At larger masses than this, the limit from CDMS-I disappears, as we would not expect any substantial number of DM particles to impinge on the detector over the course of the 1 year exposure. In the case of the CRESST 2017 surface run, we truncate the limit at a mass of O(10 8 ) GeV. For a 5 mm-scale crystal exposed over 2.27 hours, the number of 10 8 GeV DM particles crossing the detector is O(200) and above this mass, the flux of DM particles is too small to explain the O(500) events in the CRESST 2017 surface run dataset. 6 We also plot the limit obtained using our procedure in the mass range 1 GeV to 10 5 GeV. For DM particles heavier than around 50 GeV (the mass of an iron nucleus), the assumptions we have made here should not be too strongly violated. In particular, as long as m χ > m i the deflection of DM particles will be predominantly in the forward direction and so, with a large number of scatters, the typical deflection should not be too large.
Instead, for light DM the limits should be treated with caution because the deflection of DM particles in a single scatter can become large and therefore the assumption of straight-line trajectories is violated. We note, however, that while a number of particles are deflected away from the detector, some of these should be replaced by those particles which are deflected towards the detector (c.f. Ref. [31]). For isotropic deflection (in the limit of low DM masses), these two effects should balance. The 6 We note also that the right hand limit of the grey shaded region in Fig. 5 should not be taken as exact. While Xenon1T has a large enough area and exposure to probe up to DM masses of around 10 17 GeV, more detailed calculations (including stopping effects) are required to estimate where the limit should be truncated. dashed bounds plotted in Fig. 5 therefore correspond to such a scenario, where the net DM flux at the detector is unchanged and the typical path length of particles is not substantially increased by the deflections. As discussed in Refs. [32,33], those particles reaching the detector (and contributing predominantly to the limit) may not correspond to the average particles but to those which scattered less often than expected, with a resulting shorter path length than average. In the end, dedicated Monte Carlo simulations (using codes such as DaMaS-CUS [34,35]) will be required to obtain a precise limit in this region.
The results of Fig. 5 indicate that constraints from the CDMS-I run at the Stanford Underground Facility extend up to larger cross sections than previously thought, by as much as 2 orders of magnitude over a wide range of DM masses. In addition, the constraints from CRESST 2017 surface run extend the limits by an additional factor of ∼ 50 for DM masses below ∼ 10 8 GeV. In particular, these revised constraints now close a number of gaps in the parameter space which persisted between underground experiments and high-altitude experiments in the mass range between 10 7 and 10 13 GeV. Though these parameter windows are also disfavoured by limits on anomalous heating in the Earth [41,67] and neutrino telescope constraints [42], we have demonstrated that they are also ruled out independently by dedicated direct searches for DM-nuclear recoils.
Where does the improvement in the CDMS-I constraints come from, compared to the previous analysis of Ref. [37]? The first contribution comes from accounting for the form factor suppression of the scattering rate, encoded in the coherence factor C(m χ , v) in Eq. (11). The impact of form factor suppression was demonstrated in Fig. 2, where we observed that fast-moving particles are slowed less than would be expected from scattering off point-like nuclei in the Earth and shielding. While initially slow-moving particles are effectively stopped, a fast-moving population is preserved and can reach the detector.
The second contribution comes from considering the full (lab-frame anisotropic) 3-dimensional distribution of DM velocities. In Fig. 1, we showed that at SUF, the flux of DM particles appears to come mostly from above. This means that the typical DM particle must travel a much shorter path length to reach the detector than would be expected from an isotropic DM flux (as assumed in Ref. [37]). The impact of this is clear in the left panel of Fig. 4: particles impinging from above experience a much smaller stopping effect and arrive at the detector with a larger speed. The overall effect is that DM particles with much larger cross sections than previously thought can arrive at the detector with appreciable speed, increasing the maximum cross section reach of the experiments.
How robust are these constraints? Here we have assumed that the free, asymptotic velocity distribution matches the Standard Halo Model, which is typically assumed for WIMP Dark Matter. Of course, strong interactions could alter this distribution (see e.g. Ref. [68]) before the particles ever reach the Earth. In addition, there are also uncertainties on the WIMP velocity distribution from N-body simulations [69]. We could of course study specific velocity distributions using the methods presented here, but unless these lead to a substantial change in the maximum velocity of SIMPs arriving at Earth, we do not expect the limits to change by a large amount. Indeed, the largest effect would perhaps come from uncertainties in the Galactic escape speed itself, The orange shaded region is excluded by high-altitude experiments such as RSS [61], XQC [62,63], IMP 7/8 [64], IMAX [65] and SKYLAB [66] (collected in Ref. [41]). The grey shaded regions are excluded by previous analyses of direct detection experiments. This region is bounded from below by constraints from Xenon1T [3] and from the left by CRESST-III [5] and the CRESST 2017 surface run [6,30]. The solid black line denotes the previous cross section reach from Ref. [37] (reported in Ref. [41]). The shaded blue area shows the additional region of parameter space excluded by the present reanalysis of the CDMS-I (SUF) limits [45,46], while the shaded red region shows the new limits from the CRESST 2017 surface run [6]. We also exclude the area below the dashed blue and red lines, though we note that all constraints on SIMPs in this mass range come with caveats, as discussed in the text.
v esc ∈ [492, 586] km/s (90% confidence) [50]. A larger escape speed would increase the maximum cross section which can be probed by underground direct detection experiments. As can be seen in Fig. 4, increasing the cross section by a few tens of percent can reduce the maximum velocity of DM particles arriving at the detector by over 100 km/s. We would therefore expect uncertainties in v esc to contribute at most an O(30%) correction to the bounds we have derived here. As we argued in the introduction, the angular deflection of DM particles (per scatter) is δθ ∼ m N /m χ . As we move to larger DM masses and cross sections, the deflection per scatter decreases while the number of scatterings increases. This means that more and more DM particles will follow trajectories close to straight lines and the average energy loss (Eq. (6)) will become an increasingly accurate description of the entire population of particles. Of course, DM particles which undergo fewer scatters than the average will typically arrive at the detector with a larger speed, meaning that there could still be a (highly suppressed) high speed tail of particles not captured by our description. Refs. [32,33] suggest an importance sampling method to capture these effects, though the rel-  [41,67] and DM-cosmic ray interactions (green, solid) [71] are also shown. Note that these complementary constraints are model-dependent and typically require further assumptions about Dark Matter, beyond its nucleon scattering cross section (e.g. self-annihilation into Standard Model states). Large areas of the SIMP Dark Matter parameter space are ruled out by at least two constraints, both direct searches for DMnuclear scattering and indirect constraints. evant code has not yet been released. In any case, given the large number of scattering events we consider here, along with the rapid change in the CDMS-I event rate with cross section, we do not expect this high speed tail to significantly alter our results.

VI. CONCLUSIONS
In this work, we have presented limits on stronglyinteracting massive particle (SIMP) Dark Matter from direct detection experiments. We take into account the stopping effect of the atmosphere, Earth and detector shielding using a nuclear-stopping formalism. The DM particles are assumed to scatter continuously and travel along straight line trajectories, an approach we argue is valid for very heavy DM (m χ 10 5 GeV). We account for the full 3-dimensional velocity distribution of the incoming DM, as well as form factor suppression of the DMnucleus scattering, which leads to smaller energy losses for high speed particles than previously assumed.
Recent results from the CRESST 2017 surface run [6] allow us probe up to roughly 5000 times larger cross sections than was previously possible. For masses larger than about 10 8 GeV, the rate of DM particles cross-ing the detector during the CRESST 2017 surface run is small and a detector with a larger exposure time is required. Our reanalysis of the CDMS-I experiment [45,46], operated at a depth of around 10m, extends previous limits [37] by two orders of magnitude in cross section, up to DM masses O(10 15 ) GeV. This highlights the importance of careful modelling when studying the Earth-scattering of Dark Matter. With this in mind, we note that while we also present limits on DM masses down to 1 GeV, in that case our method is not strictly valid, and such limits must be confirmed with dedicated Monte Carlo simulations.
These updated constraints close a number of windows in the parameter space in the mass range 10 6 GeV to 10 13 GeV, between air-borne and high-altitude experiments on the one hand and ground-based and underground experiments on the other. For super-heavy DM produced gravitationally at the end of inflation [38,39], the DM mass is typically around the inflaton mass scale, which lies precisely in this range [70]. These constraints, which come solely from direct searches for DM-nucleus scattering, are complemented by independent indirect constraints from DM-cosmic ray interactions [71], from Ice-Cube searches for the annihilation products of SIMPs captured in the Sun [42] and from limits on the heat flux from the Earth (which would be affected by SIMP capture and annihilation) [41,67]. We emphasise however that several of these indirect constraints rely on further assumptions about the Dark Matter particle beyond its DM-nucleon scattering properties (such as its self-annihilation into Standard Model particles).
The current constraints are summarised in Fig. 6. With the updated constraints presented here, SIMP Dark Matter is now ruled out in the approximate range σ SI p /m χ ∈ [10 −37 , 10 −25 ] cm 2 /GeV, in many cases using multiple independent probes, relying on different sets of assumptions.

ACKNOWLEDGMENTS
The author thanks Laura Baudis, Jonathan Davis and Timon Emken for helpful discussions and comments on an early version of this paper. The author also thanks Ciaran O'Hare for providing useful code for calculating the lab velocity as a function of time. The author gratefully acknowledges support from the Netherlands Organisation for Scientific Research (NWO) and from the European Research Council (Erc) under the EU Seventh Framework Programme (FP7/2007-2013)/Erc Starting Grant (agreement n. 278234 -'NewDark' project). This research was also supported by the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe". Part of this work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.
Detector Earth Atmosphere FIG. 7. Coordinate system for Earth scattering. Illustration of the coordinate system used to calculate the impact of DM scattering in the Earth. The incoming direction of a single DM particle towards the detector is denoted by v, while the direction of the mean DM flux is denoted vχ = −v lab (t).
In the associated text, we give expressions for the DM path length to the detector and the velocity distribution in this coordinate system.

Appendix A: Coordinate system
Here, we give explicit expressions for the DM path length and velocity distribution in terms of the angular coordinates specifying the DM trajectory and the detector position (with respect to the mean DM flux direction). This coordinate system is illustrated in Fig. 7. We denote the detector position as r det (as measured from the centre of the Earth) and we assume that the detector is at a vertical depth d from the surface of the Earth. We denote the Earth radius as R E ≈ 6371 km and the vertical height of the atmosphere as h A ≈ 80 km. With these definitions, we have |r det | = r det = R E − d.
We assume that a specific DM particle arrives at the detector in directionv, which makes an angle θ with the vertical. 7 The straight-line path length from the top of the atmosphere to the detector for such a particle is given by: In order to correctly compute the scattering rate for DM particles, we must account for the radial dependence of the nuclear density in the Earth and atmosphere. For a DM particle which has travelled a distance D along its trajectory, as measured from the top of the atmosphere, the radial distance from the Earth's centre is: r = (R E + h A ) 2 + D 2 + 2D(r det cos θ − ) . (A2) 7 Note that θ = 0 corresponds to an upward-going DM particle.
With this, we can calculate the final speed v f from the initial speed v i (at the top of the atmosphere), along a given trajectory, using We emphasise that dv/dD (given explicitly in Eq. (11)) depends on both r and v and so Eq. (A3) must be solved numerically using an ODE solver. The DM velocity distribution is given in vectorial form in Eq. (2). Here, we rewrite this in terms of angular coordinates. We describe the incoming DM direction using θ, the polar angle measured from the downward vertical at the detector, and φ, the azimuthal angle (which is suppressed in Fig. 7). With these definitions, we can write: where cos δ = sin γ sin θ cos φ + cos γ cos θ .
In calculating the velocity distribution, we fix the magnitude of v lab = 220 km/s, with the time dependence arising through the angle γ = γ(t): We finally note that path length of the DM does not depend on the azimuthal angle φ, meaning that the final velocity at the detector depends only on the initial velocity and the polar angle, v f = v f (v i , cos θ). This means that in evaluating the final speed distribution at the detector (see Sec. IV), the two angular integrations can be performed separately.
In terms of the angular coordinates, the final speed distribution is then, where we define The values of f (v i , cos θ, γ) are easily tabulated as a function of v i and cos θ cos γ, meaning that only a single integration over θ is required (at a given value of v f and γ) in order to evaluate the final velocity distribution (Eq. A8).