Closing the window on GeV Dark Matter with moderate (b) interaction with nucleons
Abstract
We improve limits on the spinindependent scattering cross section of Dark Matter on nucleons, for DM in the 300 MeV – 100 GeV mass range, based on the DAMIC and XQC experiments. Our results close the window which previously existed in 1 – 8 GeV mass range, for a DMnucleon cross section of order b, assuming the standard velocity distribution.
Prepared for submission to JCAP
Closing the window on GeV Dark Matter with moderate (b) interaction with nucleons

Center for Cosmology and Particle Physics, Department of Physics, New York University,
4 Washington Place, New York, NY 10003, USA
1 Introduction
As constraints on Weakly Interacting Massive Particle (WIMP) DM have become more and more stringent, and most of the parameter space for WIMP DM compatible with LHC bounds has been ruled out, interest in DM in the mass range below 10 GeV has increased; see [1] for a recent review.
The existence of a window for moderatelyinteracting DM (b where is DMnucleon cross section) in the 0.3 – 100 GeV mass range was emphasized long ago [2]. Since that time both the lower and upper limits of the allowed cross section regime have moved. First, the reanalysis of XQC by Erickcek et al. [3] concluded that the upper edge of the allowed region was actually higher than determined by [2], such that a window remained in spite of the constraint from below based on DAMIC [4], which superseded the limit from CRESST [5] used in [2]. Here, we reexamine the limits from XQC and DAMIC pertinent to the existence of a window for b. We find in fact the window is closed by these experiments.
XQC — the Xray Quantum Calorimeter [6] — was an Xray detector aboard a soundingrocket. Due to its small overburden it was sensitive to DM with relatively strong interactions. The XQC data was first used to place a rough upper bound on the DMproton cross section in [7], then the bound was refined in [2]. Later, [3] claimed the bound of [2] to be overly stringent. Checking the analysis in [3], caused us to suspect that those authors assumed that the XQC detector is shielded by the body of the soundingrocket; they have confirmed this was indeed their assumption (A. Erickcek private communication). However for the lower boundary of the region that can be excluded by the XQC, the DM has negligible interaction probability in the rocket body^{1}^{1}1We only examine here the “window region” at b and have not considered other boundaries of the excluded region reported by Erickcek et al., which were also not examined in [2].. Thus our result for the lowcrosssection reach of XQC is almost a factor 10 stronger than that of [3]. While our result updates and tightens in some mass regions the bounds of [2], it is largely consistent with the original analysis in [2].
The DAMIC — Dark Matter In CCDs [8] — experiment was operated at a depth of 106.7 meters underground. As a result, DAMIC is shielded both by the Earth and the lead shield around it. If the DMnucleon cross section is large enough, the DM particles lose too much energy in collisions prior to reaching the detector, to be able to make an energy deposit above threshold for the detector^{2}^{2}2We denote as capable, those particles reaching the detector with enough energy to potentially trigger it.. Thus DAMIC only places a constraint on DM, for cross sections below some limit. Prior to the work reported here, the only analysis of the DAMIC data to address the moderatelyinteracting window was that of Kouvaris and Shoemaker [4], KS below, who used a mean scatteringlength and energyloss approximation, and assumed purely vertical propagation, following Starkman et al. [9], SGED below.
We report here the results of our detailed MonteCarlo simulation of DM particles propagating in the Earth’s crust to reach DAMIC’s depth. We developed an importance sampling technique to make the otherwise very numerically intensive simulations tractable; details of this method are given in [10]. We find that as much as four ordersofmagnitude more events actually would be seen, for the limiting cross section of the analysis found in [4] using the SGED/KS approximation, than observed by DAMIC (see figure 14 and section 4.5). This discrepancy translates to 1.8 – 5.6 times stronger limit on the DMproton cross section from DAMIC depending on the DM mass. The DMATIS code, which we developed for this study, will be publicly available online [11].
Our improved limits from reanalyzing the XQC and DAMIC data, close the window for moderatelyinteracting DM in the mass range 1 – 8 GeV^{3}^{3}3Rich et al. [12] reports exclusion of b DM candidate heavier than 2.1 GeV. Unfortunately, a detailed description of their silicon detector, exact time of the flight, exposure time, location and orientation of the detector is not provided in [12], so we cannot evaluate that analysis. To be more conservative, we restrict to limits using XQC and DAMIC. , assuming the standard DM velocity distribution. We note the strong dependence of the upper bound on the DMnucleon cross section implied by the XQC data on the DM velocity distribution. We show how the limit depends on the local DM velocity distribution, for some representative DM masses.
2 Preliminaries
The number of DM particles with velocities in the interval , which pass through a surface element of the target in a detector during a time interval , is
(2.1) 
where is the DM local mass density, is the mass of each DM particle, is the unitnormalized velocity distribution of DM particles in the detector rest frame, and is the angle between the DM velocity vector and the normal vector to the surface element.
The probability of a DM particle scattering off a nucleus of mass number in the path length element of the target along the trajectory of this DM particle, denoted by , is
(2.2) 
where is the projection of the path length element along the normal vector of the target’s surface element. is the mean interaction length of the DM particle passing through a target of nuclear mass number and number density , for DMnucleus cross section .
Using eqs. (2.1) and (2.2) and allowing for a mixed composition target, the number of nuclei in a target of volume scattered by DM particles with velocities in during the time interval in the limit of single scattering () is
(2.3) 
where is the number of nuclei of mass number in the target. Eq.(2.3) shows that only the total mass of nuclei each of mass number in the target determines the number of scattered nuclei in each bin of the DM’s velocity.
A DM particle with speed and energy in the detector rest frame transfers part of its initial energy to the target nucleus in the detector. This nuclear recoil energy is
(2.4) 
where is the scattering angle in the CM frame and is the DMnucleus reduced mass.
The momentum transfer is small for a light ( GeV) DM particle that scatters off a nucleus. Therefore, the corresponding wavelength is larger than the nuclear radius. This has two consequences:

The scatterings in the CM frame are isotropic.

The spinindependent scattering amplitudes of individual nucleons add up coherently:
(2.5) where and are the neutron and proton spinindependent amplitudes; and are the mass number and the atomic number respectively.
In a general discussion, Kurylov and Kamionkowski [13] showed that in the small momentum transfer limit only the scalar and the axial vector terms of the interaction Lagrangian survive. Using this result, the DMproton spinindependent cross section is
(2.6) 
Specializing to the case of equal neutron and proton spinindependent couplings in eq. (2.5), the spinindependent DMnucleus cross section in the small momentum transfer limit is
(2.7) 
In real collisions with nonzero momentum transfer, not all DMnucleon scattering amplitudes add up coherently. When the momentum transfer, defined by , increases such that the corresponding wavelength is no longer large compared to the nuclear radius, the cross section begins to fall with the increase of . This effect can be taken into account by a form factor, , where is the effective nuclear radius. The effective nuclear radius can be approximately found by fitting the muon scattering data to a Fermi distribution [14]
(2.8) 
with parameters: fm, fm, and fm. We use the following analytical expression which is proposed by Helm [15] to calculate the nuclear form factor
(2.9) 
The DMnucleus cross section in the nonzero momentum transfer limit and the DMproton cross section in the zero momentum transfer limit are related by:
(2.10) 
Due to the isotropy of the scatterings in the CM frame, DM particles with a given speed produce a uniform distribution of recoil energies, i.e. . So, the differential energydeposit rate is [16]
(2.11) 
where the lower limit of integration, , is the minimum speed that a DM particle needs in order to deposit the recoil energy into the detector:
(2.12) 
3 XQC bounds
To find the upper bound on the DMproton cross section for b DM, we use the observed energydeposit spectrum of the Xray Quantum Calorimetry (XQC) experiment, a rocket experiment launched on March 28, 1999. The Earth shielded the XQC detector from part of the DM flux. The body of the rocket carrying the XQC detector might also shield the detector, depending on the DM interaction length. In the following subsections, we will consider both of these shielding scenarios to calculate the expected energydeposit spectrum of XQC and the corresponding DMproton cross section upper bound.
3.1 Velocity distribution and geometrical shielding factors
Assuming an isothermal spherical density profile, the DM velocity distribution in the Galactic rest frame, , is characterized by a Maxwellian distribution with a velocity dispersion of truncated at an escape velocity . Correspondingly, the velocity of DM particles in the detector frame can be calculated by subtracting the velocity of the detector in the Galactic rest frame, , from the DM velocity in the Galactic rest frame . The black histogram in figure 0(a) shows the DM speed distribution seen by a detector that moves with the Earth’s velocity, , in the Galactic rest frame and Galactic coordinate system, for velocity dispersion and escape velocity [3].
XQC flight 27.041UG was launched 1999 March 28 at 09:00UT from White Sands, New Mexico [6]. The normal vector to the Earth’s surface at the position and time of flight launch was in the direction in the Galactic coordinate system. The data was collected a few minutes after flight launch at a timeaveraged height of 201 km above the Earth’s surface, when the normal vector of the detector was pointing toward . Figure 0(b) shows the zenith angle distribution of DM particles in various coordinate systems centered on the XQC detector. 1 corresponds to a DM particle that reaches the detector from above along the axis in each coordinate system. The cyan histogram shows the zenith angle distribution of DM particles in the coordinate system whose axis is normal to the Earth’s surface at the position and time of flight launch. The orange histogram is the zenith angle distribution in the coordinate system whose axis is aligned with the normal vector to the surface of the XQC detector during the detection period. Figure 2 shows a schematic view of these two coordinate system and their corresponding axis.
i  

1  2936  0.38  0 
2  36128  0.5  11 
3  128300  1  129 
4  300540  1  80 
5  540700  1  90 
6  700800  1  32 
7  800945  1  48 
8  9451100  1  31 
9  11001310  1  30 
10  13101500  1  29 
11  15001810  1  32 
12  18102505  1  15 
13  4000  1  60 
In the case that the body of the rocket is essentially transparent to DM particles, only the Earth shields the XQC detector. The horizontal displacement of the rocket is negligible compared to its vertical displacement and the Earth’s radius. To find the geometrical shielding factor in this case, we use the velocity distribution in the frame with the axis pointing toward the flight zenith, . Using the timeaveraged height of the XQC measurement, km, DM particles must have at least the following zenith angle to be able to reach the detector
(3.1) 
where is the Earth’s radius. The purple dashed line in figure 0(b) corresponds to this angle. The geometrical shielding factor in this case is
(3.2) 
where is the DM velocity distribution in the frame with its normal axis pointing toward vertical direction at the launch of the rocket carrying XQC and is the normalization factor.
If the body of the rocket carrying XQC completely shields the detector, only DM particles coming through the 1 steradian opening angle centered on the normal to the detector, , would reach the detector; this was the implicit assumption of Erickcek et al. [3] (private communication). The red dashdotted line at cos in figure 0(b) corresponds to the minimum zenith angle of a DM particle that reaches the XQC detector. The Earth was not blocking any part of the 1 steradian opening angle as expected since the purpose of the mission was to detect astrophysical diffuse Xray sources. The geometrical shielding factor in the bodyshielding scenario is
(3.3) 
3.2 Cross section limits
Each of the 34 XQC calorimeters was composed of a 0.96 m film of HgTe mounted on a 14 m substrate of Si [6]. The total mass of each component of the XQC detector is
(3.4) 
The exposure time of the experiment depended on the recoil energy threshold, , so we write seconds, where is the fraction of the exposure time that XQC was sensitive to recoil energy . Table 1 contains the sensitivity factors and the observed number of events in each bin taken from Erickcek et al. [3].
Substituting the geometrical shielding factors introduced in the previous subsection into eq. (2.13), the XQC differential detection rate for each shielding scenario, labeled by for the Earthshielding and bodyshielding scenarios respectively, is
(3.5) 
where and are the unitnormalized speed distributions of DM particles which are not shielded by the Earth and rocket body respectively (refer to figure 0(a)).
The differential detection rate can be calculated for a given DM mass, , a trial DMproton cross section, , and taking the local DM density 0.3 . The crudest approach is to pick an interval from the observed energydeposit spectrum whose signal to expected background ratio is the most sensitive and find the 90% CL DMproton cross section. This method was used by Zaharijas and Farrar [2] to analyze the XQC data. Erickcek et al. [3] used a more sophisticated approach which exploits the shape of the observed energydeposit spectrum by binning this spectrum and calculating a figure of merit to find the 90% CL DMproton cross section upper bound. Yellin [17, 18] argues that still a better method to find the cross section upper bound for experiments with an uncertain background is to use the maximum gap method or the optimal interval method depending on the number of events in a depositenergy interval. In the following, we adopt the Erickcek et al. method, because the Yellin method requires individual event energies to which we do not have access. While the Erickcek et al. method may not be optimal according to [17, 18], using it enables us to directly compare our result with theirs. Moreover, given that we close the window by increasing the DAMIC exclusion regime, having the most stringent possible limit from XQC is not essential.
The basic strategy of Erickcek et al. [3] is to step through trial choices of cross section and for each of them define a quantity
(3.6) 
where is the observed number of events and is the expected number of events for the given trial cross section in the ith energy bin, and is the Heaviside function. Only bins in which the predicted number is larger than the observed number are included in the sum. Thus as the trial cross section is decreased from some large value, both the values decrease and also bins which contribute at larger become excluded from the sum. Therefore, is a monotonically decreasing function of , eventually becoming zero when in every bin. The 90% CL upper limit on the cross section is defined to be the value of for which 90% of Poissondistributed datasets for the expected values, have a value smaller than in the data. Specifically:

from Poisson’s distribution with predicted mean value , the observed number of events in each bin is sampled to find the mock ’s.

The value is calculated for each mock , using eq. (3.6).

The number of cases in which the value of is smaller than the original value of is counted, and the confidence level is calculated using CL counts/(total number of trials).
The inferred DMproton cross section limit is taken to be the trial cross section for which this fraction is 90% CL.
Figure 3 shows the DMproton cross section upper bound as a function of the DM mass for the two different shielding scenarios. The bodyshielding bounds are bigger than the bounds in the Earthshielding scenario by a factor of 8.9 on average. The difference is primarily due to the ratio of the transmission factors, i.e. , with a small difference due to the slightly different speed distributions of DM particles in the Earthshielding and the bodyshielding scenarios (refer to figure 0(a)).
The aluminum body of the rocket carrying XQC produces an overburden of ([6] and Erickcek private communication [19]), corresponding to 3.7 cm aluminum. The green dotted line in figure 3 shows the DMproton cross section corresponding to the DM interaction length for an aluminum target of 3.7 cm thickness. The red dotted line shows the crosssection corresponding a mean energyloss length of 3.7 cm in aluminum. The mean energyloss length is
(3.7) 
where is the average nuclear recoil energy of a DM particle scattering off a nucleus of mass number .
For a 100 GeV DM particle, the probability of interacting with an aluminum nucleus in the rocket body is (in figure 3, for 100 GeV, compare the value of solid orange line, 3.3, and the value of green dotted line, 1.5). Therefore even for the heaviest DM candidate considered here, the body of the XQC rocket needs to be at least 50 times denser or thicker than 3.7 cm aluminum to begin to shield the detector significantly. We conclude that the rocket body did not shield the XQC detector for the smallest values of that can be ruled out by the XQC experiment. Thus the appropriate scenario is the onlyEarthshielding one.
Finally, we check the validity of our simplification in ignoring the geometry of the XQC detector itself and using the total mass of each type of nuclei in the XQC detector to calculate the depositenergy rate in eq. (3.5). The mean interaction length of DM particles in a target with a mix of nuclei can be calculated by
(3.8) 
where for the XQC detector. Figure 4 shows this mean interaction length using the DMproton cross section upper bound in the two shielding scenarios. The mean interaction length in the bodyshielding and the Earthshielding scenarios are respectively at least 10 and 100 times bigger than the typical size of the XQC detector ( mm). Thus DM particles do not have multiple scatterings in the XQC detector for the smallest values of that are ruled out by the XQC experiment and eq. (3.5) is accurate.
3.3 Sensitivity to the DM velocity distribution
In the above discussions, we have only considered the usual DM velocity distribution. However, as discussed in [20], for a large enough cross section, DM interactions with gas in the Galactic disk can in principle significantly alter the local DM velocity distribution, causing the DM to corotate with the solar system to some degree. Given that the energy available to deposit in XQC depends quadratically on the velocity, reducing the DM relative velocity distribution reduces the flux of particles capable of producing an energy deposit above threshold. We illustrate this effect by determining the cross section bounds obtained from XQC, under the assumption that the DM velocity relative to XQC is given by a Maxwellian with peak velocity . Figure 5 shows how, as the relative speed of DM particles decreases, the upper bound on DMproton cross section increases and XQC’s bounds become weaker. Note that for crosssections above , the bodyshielding limit applies. For example, if 1.7 GeV DM particles corotate with the baryonic matter such that their speed distribution is a truncated Maxwellian in the XQC rest frame a window opens for peak velocity , so that both the XQC and DAMIC bounds can be evaded.
4 DAMIC bounds
With the weakly interacting massive particle (WIMP) DM paradigm in mind, most DM direct detection experiments are built underground. The layer of the Earth’s crust acts as a shield that filters out unwanted cosmicray and radioactive backgrounds while not affecting the DM flux. For a b DM candidate, this layer of the Earth’s crust also shields DM particles.^{4}^{4}4We also considered the effect of the overburden of the Earth’s atmosphere. For the DAMIC experiment at 106.7 m, the Earth’s atmosphere has a 27 smaller column depth ( ) in comparison to the Earth’s crust. Due to the exponential dependence of the DM particle’s final energy to the overburden column depth, the overburden from the Earth’s atmosphere has a negligible effect in slowing down DM particles before reaching DAMIC’s detector, for cross sections in the DAMIC exclusion region.
In this moderately strongly interacting regime, where the DM particle mean interaction length in the Earth’s crust is much smaller than the detector depth , DM particles scatter off nuclei in the Earth’s crust at least once during their passage. In each scattering, the DM particle loses a fraction of its initial energy. After a few scatterings in the Earth’s crust, three scenarios are imaginable for each DM particle:

It scatters back to the Earth’s atmosphere.

Its energy falls below , the minimum energy to potentially trigger the detector, before reaching the detector.

It reaches the detector at the depth of with .
The minimum energy of a DM particle to scatter a silicon target nucleus in the DAMIC detector and produce the threshold recoil energy is
(4.1) 
from Equation (4.3) for backward scattering in the CM frame, i.e. , and substituting for the reduced mass.
The DAMIC apparatus, installed at the rather shallow underground depth of 106.7 meters at Fermilab [8], is an ideal experiment for constraining b DM candidates. The DAMIC detector, made of silicon, accumulated a total exposure of 107 gdays from June 2010 to May 2011, and observed a total of 106 events with an ionization signal between 40 (eV electron equivalent energy) and 2 k.
4.1 MonteCarlo simulation
We start this subsection by discussing the steps of a bruteforce MonteCarlo simulation which implements the actual path length and scattering angle distributions to find the velocity distribution of DM particles at the DAMIC detector’s depth. Here are the steps of the bruteforce MonteCarlo simulation:
Step 0: For the long exposure time (from June 2010 to May 2011) of DAMIC, we considered the timeaveraged initial zenith angle of DM particles to be isotropic, i.e. has a uniform distribution in . (We changed the definition of hereafter, so that 1 corresponds to a DM particle that travels along the direction from above the detector). Also, we considered the shielding layer of the Earth’s crust above the DAMIC detector as a slab of material due to the smallness of the detector depth in comparison to the Earth’s radius. Therefore, only the vertical displacement of DM particles in the Earth’s crust which we denote by is tracked. Knowing the zenith angle , the depth of a DM particle at the first scattering is selected from the scattering length distribution (details below, eq. (4.8)). The DM particle’s initial speed is sampled from the speed distribution on the Earth’s surface (figure 0(a)). If its initial energy is smaller than the minimum required energy to trigger the DAMIC detector, we count this DM particle as one of the particles which does not give a signal. If its energy is above the minimum energy and the particle is already at the DAMIC depth, , it is counted as a capable DM particle. Otherwise, the DM particle enters the first scattering iteration.
Step 1: Choose the target nucleus. To calculate the energy loss, we first choose the target nucleus for the given iteration of scattering. The probability that a given DM particle scatters off a nucleus of mass number is
(4.2) 
The nuclei O (46.5), Si (28.9), Al (8.9), and Fe (4.8) make up of the Earth’s crust mass constituents. We scale the mass abundances of these elements to sum up to to calculate the number density of each element in the Earth’s crust. The relative probability is a function of the DM particle’s mass and the target element’s mass, but it is independent of the DMproton cross section.
Step 2: Choose the CM scattering angle, , fixing the final energy. The scattering is isotropic in the CM frame, thus has a uniform distribution. Given the target element mass number , the ratio of the DM particle’s energy after scattering to its energy before scattering is:
(4.3) 
where is the scattering angle in the CM frame.
The scattering angle in the lab frame, , and the scattering angle in the CM frame, , are related by
(4.4) 
Step 3: Choose the azimuthal scattering angle, , fixing the new direction. Figure 6 shows a schematic view of the scattering parameters in the lab frame. To find the zenith angle after scattering , it is sufficient to rotate in the primed coordinate system, , around the axis with zenith angle . This gives a recursive relation for the zenith angle which is also a function of angles and :
(4.5) 
Generate for a random number in , to calculate the zenith angle after scattering using eqs. (4.4) and (4.5).
Step 4: Find the position of the next interaction. The probability of a DM particle scattering exactly once, off of a nucleus of specified mass number in , is the product of the probability that the DM particle has not scattered in the distance and the probability that it scatters in the distance (see e.g. [21])
(4.6) 
The path length distribution of a DM particle passing through the Earth’s crust with different constituents is the sum of all elements’ probability distributions
(4.7) 
Using eqs. (4.5) and (4.7), the vertical displacement between scatterings and 1 can be calculated
(4.8) 
where is a random number representing the path length’s cumulative probability distribution.
Steps 1 – 4 are repeated until each DM particle is categorized as follows

If at any stage the depth of a DM particle is less than zero (positive direction of taken to be downward with at the Earth surface), that particle is not tracked any more and is counted as one of the particles which scatters back to the Earth’s atmosphere.

If the energy of a DM particle becomes smaller than the minimum energy, , before reaching the DAMIC detector, that DM particle is no longer tracked and is counted as one of the events with energy below the detector’s threshold.

If a DM particle reaches the DAMIC detector’s depth, , with potentially enough energy to trigger the detector, , that particle is counted as a capable DM particle which can potentially trigger the detector. Such particles are called capable DM particles in this paper.
There is a 6inch lead shielding layer around the DAMIC detector. Even though the thickness of this shielding layer is much smaller than the Earth’s crust shielding layer of 106.7 meters, the lead shield can significantly reduce the energy of heavier DM particles for which the Earth’s crust nuclei are not so well massmatched. So, DM particles which end up in the third category go through more iterations of scattering in the lead shield. There is no need for step 1 as there is only one type of nucleus in the lead shield and replaces in eqs. (4.7) and (4.8).
4.2 Importance sampling
Only an extremely small fraction of DM particles (as small as for DM mass 4 GeV at the limiting cross section) reaches the DAMIC detector with enough energy to potentially produce an energy deposit above threshold (see figure 6(a)). This means, for a DM mass of 4 GeV, the propagation of particles must be simulated using the actual path length distribution to achieve one DM particle with enough energy at the DAMIC surface to potentially trigger the detector. To make the determination of the DM velocity distribution at the DAMIC detector’s depth computationally feasible, we developed an importance sampling technique which reduces the computational cost of the bruteforce simulation by increasing the mean path length that each DM particle travels. We account for this modification of the path length by assigning a weight factor to each simulated DM particle [10].
To increase the sample size of capable DM particles, the path length distribution which is used in eq. (4.7) is modified to
(4.9) 
where shifts the average path length from to while keeping the general shape of the path length distribution the same.
To correct for not sampling from the true path length distribution for the th DM particle, we assign a weight for the th scattering iteration, . The total weighting factor for the th particle (having scatterings in the Earth’s crust and lead shield) is
(4.10) 
This method, called importance sampling, enables our MonteCarlo simulation to reconstruct the DM velocity distribution on the surface of the DAMIC detector using 100 to 1000 times less computational resources for 0.4 – 0.8, by wasting less computational power in simulating the events which end up losing too much energy in the Earth’s crust or the lead shield. In [10], we discuss how to calculate the distribution of any physical quantity using the importance sampling method and also validate our implementation of the strategy by comparing numerous physical distributions obtained with the bruteforce and importance sampling methods.
4.3 Cross section limits
The DAMIC detector accumulated a total exposure of 107 gdays from June 2010 to May 2011, and observed a total of 106 events with an ionization signal between 40 and 2 k [8]. Interpreting a spectrum of ionization signal in terms of the underlying nuclear recoil energy, requires understanding the relationship between the two. For high energies, the Lindhard model [22] can be used to find the nuclear recoil energy spectrum which produces the measured ionization signal spectrum. The Lindhard model gives that the nuclear recoil energy corresponding to the maximum ionization signal 2 k is 7 keV. However, for low energies the Lindhard model fails. The DAMIC team measured the ionization signal of silicon nuclei down to 60 , and found a significant deviation from the Lindhard model below 1 keV nuclear recoil energy [23]. To find the threshold nuclear recoil energy corresponding to the DAMIC electronic threshold, we extrapolate the DAMIC team’s measurement of ionization factor from 60 to 40 and find the equivalent threshold recoil energy of 550 eV. We could more conservatively limit our analysis to the ionization signal above 60 , where the corresponding nuclear recoil energy 700 eV has been determined experimentally. That changes the DMproton crosssection bounds by less than one percent, but increases the minimum DM mass to which the bound applies from 1 GeV to 1.16 GeV. The functional dependence seems clear and the extrapolation small, so we adopt 40 and 550 eV as the threshold.
We use the total number events observed by DAMIC, in the 40 – 2 range, to find the DMproton cross section lower bounds. Considering the observed spectral shape does not improve the cross section bounds, due to the number of DM particles with sufficient energy being exponentially dependent on the number of scatterings and hence on the cross section.
For any assumed cross section, eq. (2.13) can be integrated to calculate the expected total number of events for DAMIC’s silicon target ():
(4.11) 
where is the normalized speed distribution of capable DM particles and is the fraction of capable DM particles at the Earth’s surface. is the attenuation parameter (the ratio of the number of capable DM particles at the detector to the number of capable DM particles at the Earth’s surface) which is normalized to one at the Earth’s surface. The extra factor of two in the denominator of eq. (4.11) is due the Earthshielding of DM particles entering the Earth from below the horizon.
The DMATIS code [11], a library written in Python and MATHEMATICA, implemented the importance sampling MonteCarlo simulation. For DM mass and trial DMproton cross section , the DMATIS code calculates the attenuation parameter and .
Figure 6(a) shows , the attenuation parameters of capable DM particles reaching the lead shield and reaching the DAMIC detector, as a function of DM mass, corresponding to the 90% CL DMproton crosssection lower bounds that we calculated from DAMIC’s data.
Figure 6(b) shows the attenuation parameter for the number of events in the DAMIC energy range (the ratio of the number of events in DAMIC at its actual depth, to the number of events if the detector were on the Earth’s surface), as a function of mass, corresponding to the 90% CL cross section that we have calculated from the DAMIC data. The difference between the Earthonly and the total attenuation parameters of the number of events increases by several orders of magnitude over the mass range considered. This demonstrates the importance of considering the thin (6 inches) lead shield around the DAMIC detector in obtaining the correct limits.
As the cross section increases, a smaller and smaller fraction of DM particles can reach the detector. As the trial cross section is decreased from some large value, the attenuation parameter increases faster than . Therefore, the expected total number of events is a monotonically decreasing function of . The 90% limit on the allowed (by DAMIC) cross section is defined to be the value of for which the expected total number of events is 123. (The equivalent 90% CL upper limit on 106 observed total number of events is 123.) The solid line in figure 8 shows the 90 CL DMproton cross section sensitivity at detector depth of 106.7 meters underground. The result from the SGED/KS method (dashed line), which we discuss briefly in the next subsection and in detail in [10], is added for comparison.
Using our 90% CL cross section limits as a function of mass that we have calculated from the DAMIC data, the mean interaction lengths in the Earth’s crust, in a target of copper, and in a target of lead are given in figure 9. The mean interaction lengths in a target of copper for these cross sections are at least ten times bigger than the copper shield thickness ( cm) around the DAMIC detector. This justifies our shielding model which ignores the copper shield around the DAMIC detector.
4.4 Properties of DM particles that successfully triggered DAMIC
It is enlightening to contemplate some properties of the DM particles that produce events in the DAMIC detector, for representative DM masses 1.7, 10, and 100 GeV. In each case, we use the 90% CL sensitivitylimit cross section for these DM masses, i.e. , 7.2, and 9.3.
Figure 10 shows the distributions of the final energies and the recoil energies of successful DM particles, i.e. DM particles that reach the DAMIC detector and scatter the detector giving a depositenergy above the DAMIC threshold, 550 eV. In this mass range, as the DM mass increases, DM particles can lose a larger fraction of their initial energy while passing through the Earth’s crust and the lead shield, and deposit a smaller fraction of their initial energy to the silicon nuclei of the DAMIC detector, and still produce an above threshold energydeposit.
Figure 11 shows the zenith angle distribution of successful DM particles. Although successful DM particles mainly enter and exit the Earth’s crust vertically for any mass, the heavier DM particles exit the lead shield with a much broader zenith angle distribution. This broadening is due to the similarity of the DM and target nuclei masses in the Earth’s crust and the lead shield.
Figures 11(a) and 11(b) show distributions of the zenith angles in all scatterings of successful DM particles in the Earth’s crust and the lead shield. Depending on the mass of DM particles and the material that these DM particles pass through, the distribution of the averaged zenith angle changes. Clearly, DM particles that trigger the detector successfully do not move vertically with zero deflection angle.
4.5 Deficiencies of the SGED/KS method
Starkman, Gould, Esmailzadeh and Dimopolous [9] proposed a continuous energyloss and vertical propagation approximation, to estimate the maximum cross section for which a given experiment is sensitive to strongly interacting DM particles. They took the inverse energylosslength, to be where is the mean fractional energy loss for isotropic scattering and is the scattering length in the material. This estimation method was applied by Kouvaris and Shoemaker [4], to analyze DAMIC’s data.
Figure 8 shows that the SGED/KS method underestimates the true sensitivity range of DAMIC. Figures 12(a) and 12(b) show that DM particles that trigger the detector are exploiting the tails of the path length and scattering angle distributions to lose less energy. Using mean values of these quantities is very inaccurate. The problem is most acute when the DM energyloss per interaction can be significant due to the DM and target nucleus masses being comparable. Figures 11 and 12 show that the assumption of purely vertical propagation is inaccurate. For a more detailed discussion of the problems with the SGED/KS method see [10].
Although our lower bounds on the DMproton cross section are larger by only a factor of 1.8 – 5.6 than KS, this difference corresponds to a four orderofmagnitude larger expected total number of events for the KS limiting cross section (compare orangesolid and cyansolid lines in figure 14). This reveals how badly the SGED/KS method fails to calculate the energy loss of DM particles propagating in the Earth and the lead shield. Nonetheless, as explained in [10], the SGED approximation work to get an orderof magnitude estimate of the cross section sensitivity which was the original SGED purpose; the problem is to apply it for extracting a limit as done in KS [4].
5 Conclusion
In this work, we have derived the constraints on moderately interacting dark matter which follow from the XQC and DAMIC experiments. Previous analyses made shortcuts which significantly reduced the apparent constrainingpower of these experiments and indicated a window for 0.1 – 1 b, for dark matter mass in the 1 – 100 GeV range. Our new, more stringent limits, close the window for 1 – 8 GeV and exclude the possibility that Dark Matter having a mass 1 – 100 GeV, has moderate interactions with nucleons, when combined with limits from other experiments.
Our upper limit on the DMnucleon scattering cross section from the XQC soundingrocket experiment [6] is a factor 8.9 lower (more restrictive) than reported by Erickcek et al. [3]; where analysis made the assumption that the rocket body shields the detector which is not the case in the the relevant cross section range. Our results refine but generally confirm, the earlier analysis of XQC of [2].
The DAMIC experiment operating 106.7m underground, is shielded by the Earth and lead shielding. Due to energy loss in the overburden, only a tiny fraction of the flux of DM particles at the surface of the Earth, has sufficient energy to trigger the detector. This means that there is a limiting cross section, above which DAMIC is insensitive to DM. In order to determine what crosssections DAMIC can exclude, one must accurately determine the attenuation of DM particles capable of producing an above threshold energydeposit, as a function of cross section. The analysis of Kouvaris and Shoemaker [4], used a poor approximation to the energy loss and overestimated the attenuation enormously (see figure 14). With the limiting cross section we find, the KS approximation implies 100% attenuation, and for their reported limit, they overestimate the attenuation by a factor . With a correct treatment of energy loss, we find that DAMIC is sensitive to cross sections a factor 2 or more higher than KS reports, extending the exclusion range of DAMIC to higher cross sections and overlapping the new XQC upper bound.
Besides the exclusion of this type of Dark Matter, our results shed light on the process of energy loss in the Earth or another overburden. We showed that the commonlyused SGED [9] approximation, based on the mean number of interactions corresponding to an assumed cross section and vertical propagation to the detector, and the average energyloss per interaction, is grossly misleading in some circumstances. The problem is particularly severe when the mass of nuclei in the overburden is within a factor fewto10 of that of the DM particle, so that a substantial fraction of the DM’s energy can be lost in each collision. In this case, there is a population of trajectories with lowerthanaverage scattering angles and longerthanaverage path lengths between interactions, which can arrive at the detector with sufficient energy to be capable of triggering the detector. Such a population is completely missed in the SGED approximation used by [4] to analyze the DAMIC experiment, resulting in a four orderofmagnitude larger number of events expected in DAMIC over a large mass range (see figure 14). We present more details of this phenomenon in a companion paper [10], where we also present an importance sampling technique to reduce computational time by a factor 100 – 1000. This makes the MonteCarlo simulation of trajectories and energy loss feasible, even when the attenuation parameter is smaller than . The DMATIS code, which we developed for this study, will be publicly available online [11].
Figure 15 summarizes our 90% CL DMproton spinindependent cross section bounds for DM masses of 0.3 – 100 GeV, using the XQC and DAMIC energydeposit spectra. This closes the previouslyexisting window for b) cross sections for DM masses in the 1 – 8 GeV range. All of the limits presented here assume the standard DM velocity distribution. For completeness, we performed simulations assuming other Maxwellian velocity distributions; the results are presented in Figure 5.
Acknowledgements: We thank A. Erickcek and D. McCammon for information about XQC and the analysis of [3]. MSM acknowledges support from James Arthur Graduate Assistantship; the research of GRF has been supported by NSFPHY1212538 and AST1517319.
References
 [1] G. B. Gelmini, Light WIMPs, 1612.09137.
 [2] G. Zaharijas and G. R. Farrar, Window in the dark matter exclusion limits, Phys. Rev. D 72 (Oct, 2005) 083502.
 [3] A. L. Erickcek, P. J. Steinhardt, D. McCammon and P. C. McGuire, Constraints on the interactions between dark matter and baryons from the xray quantum calorimetry experiment, Phys. Rev. D 76 (Aug, 2007) 042007.
 [4] C. Kouvaris and I. M. Shoemaker, Daily modulation as a smoking gun of dark matter with significant stopping rate, Phys. Rev. D90 (2014) 095011, [1405.1729].
 [5] G. Angloher, M. Bruckmayer, C. Bucci, M. BÃ¼hler, S. Cooper, C. Cozzini et al., Limits on {WIMP} dark matter using sapphire cryogenic detectors, Astroparticle Physics 18 (2002) 43 – 55.
 [6] D. McCammon, R. Almy, E. Apodaca, W. B. Tiest, W. Cui, S. Deiker et al., A high spectral resolution observation of the soft xray diffuse background with thermal detectors, The Astrophysical Journal 576 (2002) 188.
 [7] B. D. Wandelt, R. Dave, G. R. Farrar, P. C. McGuire, D. N. Spergel and P. J. Steinhardt, Selfinteracting dark matter, in Sources and detection of dark matter and dark energy in the universe. Proceedings, 4th International Symposium, DM 2000, Marina del Rey, USA, February 2325, 2000, pp. 263–274, 2000. astroph/0006344.
 [8] J. Barreto, H. Cease, H. Diehl, J. Estrada, B. Flaugher, N. Harrison et al., Direct search for low mass dark matter particles with {CCDs}, Physics Letters B 711 (2012) 264 – 269.
 [9] G. D. Starkman, A. Gould, R. Esmailzadeh and S. Dimopoulos, Opening the window on strongly interacting dark matter, Phys. Rev. D 41 (Jun, 1990) 3594–3603.
 [10] M. S. Mahdawi and G. R. Farrar, Energy loss during Dark Matter propagation in an overburden, 1712.01170.
 [11] M. S. Mahdawi and G. R. Farrar, DMATIS, Dark Matter ATtenuation using Importance Sampling, will be made available at https://github.com/msm550/DMATIS (2017) .
 [12] J. Rich, R. Rocchia and M. Spiro, A search for strongly interacting dark matter, Physics Letters B 194 (1987) 173 – 176.
 [13] A. Kurylov and M. Kamionkowski, Generalized analysis of the direct weakly interacting massive particle searches, Phys. Rev. D 69 (Mar, 2004) 063503.
 [14] G. Fricke, C. Bernhardt, K. Heilig, L. Schaller, L. Schellenberg, E. Shera et al., Nuclear ground state charge radii from electromagnetic interactions, Atomic Data and Nuclear Data Tables 60 (1995) 177 – 285.
 [15] R. H. Helm, Inelastic and elastic scattering of 187mev electrons from selected eveneven nuclei, Phys. Rev. 104 (Dec, 1956) 1466–1475.
 [16] J. Lewin and P. Smith, Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil, Astroparticle Physics 6 (1996) 87 – 112.
 [17] S. Yellin, Finding an upper limit in the presence of an unknown background, Phys. Rev. D 66 (Aug, 2002) 032005.
 [18] S. Yellin, Extending the optimum interval method, ArXiv eprints (Sept., 2007) , [0709.2701].
 [19] A. L. Erickcek, A new constraint on stronglyinteracting dark matter from xray quantum calorimetry, Senior Thesis (2003) .
 [20] G. R. Farrar, 6quark Dark Matter, 2017. 1711.10971.
 [21] S. Tavernier, Experimental Techniques in Nuclear and Particle Physics. 2010, 10.1007/9783642008290.
 [22] J. Lindhard, V. Nielsen, M. Scharff and P. Thomsen, Integral equations governing radiation effects. (notes on atomic collisions, iii), Kgl. Danske Videnskab., Selskab. Mat. Fys. Medd. Vol: 33: No. 10 (Jan, 1963) .
 [23] A. E. Chavarria et al., Measurement of the ionization produced by subkeV silicon nuclear recoils in a CCD dark matter detector, Phys. Rev. D94 (2016) 082007, [1608.00957].