DaMaSCUS: The Impact of Underground Scatterings on Direct Detection of Light Dark Matter

TL;DR

This paper introduces DaMaSCUS, a Monte Carlo simulation tool that models how underground scatterings affect the detection of light, strongly interacting dark matter, impacting recoil spectra and diurnal modulation signals.

Contribution

The paper presents a novel publicly available simulation code that accurately models underground scatterings of dark matter, accounting for their effects on detection signals and diurnal modulation.

Findings

Underground scatterings significantly alter dark matter recoil spectra.

The simulation predicts diurnal modulation patterns in detection signals.

Earth shadowing effects depend on detector depth and location.

Abstract

Conventional dark matter direct detection experiments set stringent constraints on dark matter by looking for elastic scattering events between dark matter particles and nuclei in underground detectors. However these constraints weaken significantly in the sub-GeV mass region, simply because light dark matter does not have enough energy to trigger detectors regardless of the dark matter-nucleon scattering cross section. Even if future experiments lower their energy thresholds, they will still be blind to parameter space where dark matter particles interact with nuclei strongly enough that they lose enough energy and become unable to cause a signal above the experimental threshold by the time they reach the underground detector. Therefore in case dark matter is in the sub-GeV region and strongly interacting, possible underground scatterings of dark matter with terrestrial nuclei must be taken into account because they affect significantly the recoil spectra and event rates, regardless of whether the experiment probes DM via DM-nucleus or DM-electron interaction. To quantify this effect we present the publicly available Dark Matter Simulation Code for Underground Scatterings (DaMaSCUS), a Monte Carlo simulator of DM trajectories through the Earth taking underground scatterings into account. Our simulation allows the precise calculation of the density and velocity distribution of dark matter at any detector of given depth and location on Earth. The simulation can also provide the accurate recoil spectrum in underground detectors as well as the phase and amplitude of the diurnal modulation caused by this shadowing effect of the Earth, ultimately relating the modulations expected in different detectors, which is important to decisively conclude if a diurnal modulation is due to dark matter or an irrelevant background.

Tables

Table 1. Table 1 : Benchmark points for the MC simulations with m χ = 500 subscript 𝑚 𝜒 500 m_{\chi}=500 MeV. Here ⟨ N sc ⟩ delimited-⟨⟩ subscript 𝑁 sc \langle N_{\rm sc}\rangle is the average number of underground scatterings of the simulated trajectories and the sample size N sample subscript 𝑁 sample N_{\rm sample} is the number of recorded velocity data points per isodetection ring.

Simulation ID	${\sigma }_{\chi n}^{\mathrm{SI}}\left(0\right)$[pb]	$⟨N_{\mathrm{sc}}⟩$	$N_{\mathrm{sample}}$
‘SS’	0.521	0.12	${10}^7$
‘MS1’	4.26	1.0	${10}^7$
‘MS10’	41.2	10.0	${10}^6$
‘MS50’	300.0	$\gtrsim$50.0	$5\cdot {10}^5$

Table 2. Table 2 : Relative Element Abundances in the Earth Core and Mantle  [ 54 ] . For each element we use the most abundant isotope.

Element	Core[%]	Mantle[%]		Element	Core[%]	Mantle[%]
56Fe	85.5	6.26		32S	1.9	0.03
16O	0	44		52Cr	0.9	0.26
28Si	6	21		23Na	0	0.27
24Mg	0	22.8		31P	0.2	0.009
58Ni	5.2	0.2		55Mn	0.3	0.1
40Ca	0	2.53		12C	0.2	0.01
27Al	0	2.35		1H	0.06	0.01
				Total	100.26	99.83

Equations

$$\frac{{\sigma }_{\chi p}}{{\sigma }_{\chi e}}\simeq {\left(\frac{{\mu }_{\chi p}}{{\mu }_{\chi e}}\right)}^2,$$

$$P(L)$$

$${\lambda }_{\text{MFP}}^{-1}(\vec{x},\vec{v})$$

$$P(\text{scattering on j})$$

$${\sigma }_{\chi A}(q^2)={\sigma }_{\chi A}(0)F_A^2(q^2),\text{where }{F}_A(0)=1.$$

$${\sigma }_{\chi A}^{\mathrm{total}}$$

$${\sigma }_{\chi A}^{\text{SI}}(0)={\sigma }_{\chi n}^{\text{SI}}(0)\frac{{\mu }_{\chi A}^2}{{\mu }_{\chi n}^2}{A}^2,$$

$$P(L)=\xi \in (0,1),$$

$$\vec{\Delta }(\vec{x},\vec{v})=L{\vec{e}}_v,$$

$${\vec{v}}^{\prime }=\frac{m_A\left|\vec{v}\right|\vec{n}+m_{\chi }\vec{v}}{m_A+m_{\chi }}.$$

$$\frac{\underset{0}{\overset{q^2}{\int }}\mathrm{d}{q}^2{F}_A^2(q^2)}{\underset{0}{\overset{q_{\max }^2}{\int }}\mathrm{d}{q}^2{F}_A^2(q^2)}=\xi \in (0,1)$$

$$\cos \alpha =1-2\frac{q^2}{q_{\max }^2}.$$

$${\vec{v}}_{\text{ini}}={\vec{v}}_{\text{halo}}-{\vec{v}}_{\oplus }(t).$$

$$f_{\text{halo}}(\vec{v})$$

$${\vec{r}}_{\text{ini}}$$

$$\Theta (t)=\arccos \left[\frac{{\vec{v}}_{\oplus }(t)\cdot {\vec{x}}_{\mathrm{lab}}(t)}{v_{\oplus }(t)(r_{\oplus }-d_{\mathrm{lab}})}\right],$$

$$A_k=2\pi (r_{\oplus }-d_{\mathrm{lab}}{)}^2\left[\cos ({\Theta }_k)-\cos ({\Theta }_k+\Delta \Theta )\right].$$

$$\delta ({\Phi }_{\mathrm{lab}})=100\frac{R_{\max }-R_{\min }}{R_{\max }}\%,$$

$$n_{\text{J2000.0}}$$

$$+\frac{h}{24}+\frac{m}{24\times 60}+\frac{s}{24\times 60^2}-730563.5,$$

$$\tilde{Y}$$

$$\tilde{M}$$

$$T_{\text{J2000.0}}\equiv \frac{n_{\text{J2000.0}}}{36525}.$$

$$\displaystyle=86\,400s\,\big{[}0.7790\,5727\,32640+n_{\text{J2000.0}}\text{mod }1+0.0027\,3781\,1911\,35448\;n_{\text{J2000.0}}\big{]}$$

$$+$$

$$=\text{GMST}+E_e(T_{\text{J2000.0}}),$$

$$E_e(T_{\text{J2000.0}})$$

$$\Delta \psi$$

$$\Omega$$

$$L$$

$${ϵ}_A$$

$$\text{LAST}(\lambda )=\text{GAST}+\frac{\lambda }{360^{\circ }}86400s.$$

$${\vec{x}}^{\text{(equat)}}=\mathcal{N}{\vec{x}}^{\text{(lab)}},\text{ with }\mathcal{N}=\left(\begin{array}{ccc}-\sin \varphi & -\cos \theta \cos \varphi & \sin \theta \cos \varphi \\ \cos \varphi & -\cos \theta \sin \varphi & \sin \theta \sin \varphi \\ 0& \sin \theta & \cos \theta \end{array}\right),$$

$${\vec{x}}^{\text{(geo-ecl)}}=-{\vec{x}}^{\text{(equat)}}.$$

Full text

Preprint: CP3-Origins-2017-20 DNRF90

1 Introduction

The direct detection of dark matter (DM) particles from the galactic halo has been an ongoing endeavour for the last three decades [1, 2]. The various direct detection experiments have been the most straight forward strategy to shed light on one of the most intriguing questions in physics, the nature of dark matter. Yet no conclusive experimental evidence for DM has been found so far and we are still confronted with the discrepancy between the virtually conclusive gravitational evidence on all large scales, from galactic to cosmological [3], and the absence of any Earth-based experimental evidence.

The conventional direct detection approach is based on elastic DM-nucleus collisions and the subsequent observation of the nuclear recoil. Experiments such as LUX [4, 5, 6] have been successful in constraining the standard ‘WIMP paradigm’, putting severe bounds on interaction cross-sections for DM masses above several GeV. The continuing null results have therefore led to a shift towards the relaxation of underlying assumptions and a new focus on DM models beyond the classic WIMP.

One aspect of this shift is the redirection of experimental efforts towards lighter DM. Even though we do not make assumptions about UV-completions in this paper, there exists a series of models accommodating light DM, e.g. asymmetric DM [7, 8, 9, 10, 11, 12, 13]. Sub-GeV weakly interacting particles evade the common direct searches due to their soft nuclear recoil energies falling below experimental recoil energy thresholds, typically of the order of keV. Experiments like DAMIC [14], CRESST-II [15], EDELWEISS-III [16] or CDMSlite [17] pushed the limits of conventional detectors by using low-mass target nuclei and realizing recoil thresholds as low as $\mathcal{O}$(100eV), such that the sensitivity reaches to DM masses down to $m_{\chi}\approx 500$ MeV for CRESST-II. For even lower DM masses however new detection strategies are required, since the discrimination between soft nuclear recoils and background becomes a serious problem.

A very promising idea is the search for DM-electron scatterings, which lead to new detection signatures sensitive to masses below $\sim$GeV. Different approaches have been suggested, such as electron ionizations and excitations in atoms or semiconductors [18, 19, 20, 21]. To trigger this kind of signals the kinetic energy of the DM particles needs to exceed binding energies of the order of only $\mathcal{O}$(10eV) for atoms and $\mathcal{O}$(eV) for semiconductors, rendering possible the discovery of light DM, provided that the detector is sensitive to such low energy deposits. Using this channel first limits on sub-GeV DM using DM-electron scatterings have been presented for XENON10 [22] and XENON100 [23]. Other proposed targets for DM-electron scattering experiments were scintillators [24], superconductors [25, 26], two-dimensional targets [26] and superfluid helium [27, 28], some of which potentially sensitive to light DM with masses as low as $\mathcal{O}$(keV).

Apart from utilizing DM-electron scatterings, further search strategies include the use of the Bremsstrahlung of the nuclear recoil in a conventional detector [29], allowing to constrain MeV-scale DM e.g. with xenon detectors [29, 30]. Additional methods exploit the disintegration of chemical bonds [31], or other effects [32, 33].

In light of these new ideas concerning the direct searches for sub-GeV DM, there is a crucial aspect that should not be neglected. While it is true that elastic DM-nucleon scatterings of light DM are not directly observable below a certain DM mass, it also means that the corresponding cross-section is still allowed to be large and scatterings on nuclei inside the Earth may occur nevertheless [20]. A prime example in the context of DM-electron scattering experiments are models which include a dark photon mixing kinetically with the SM photon [34, 35]. In the heavy mediator limit the DM-electron and the DM-proton scattering cross-section are related via

[TABLE]

where $\mu_{XY}$ refers to the reduced mass of two particle species $X$ and $Y$. For DM masses of the order $\mathcal{O}$(100MeV) this leads to a interaction strength hierarchy of about $(\mu_{\chi p}/\mu_{\chi e})^{2}\sim(m_{\chi}/m_{e})^{2}\sim\mathcal{O}(10^{4}-10^{5})$. Hence DM-electron scattering cross-sections being tested by new experiments are accompanied in this model by much stronger, yet unobservable, DM-proton interactions.

In cases like this, sub-GeV DM particles scatter elastically on terrestrial nuclei while they travel through the Earth towards the detector, leading to deformed DM density and velocity distributions due to deflection and deceleration, which in turn has an impact on any direct detection experiment regardless of its specific detection channel or search strategy.

The typical signature of this deformation are diurnal modulations of the detection signal rate. The average distance a DM particle moves through the Earth’s bulk mass to reach a detector varies as the Earth rotates. Therefore the pre-detection scattering probability changes periodically over a sidereal day, and with it the modification of the DM density and velocity distribution. This is especially true for detectors in the southern hemisphere, which are much more sensitive to this ‘Earth shadowing’ effect. This effect has been quantified in early MC simulations in the context of the classic WIMP [36, 37, 38], and further studied in the context of hidden sector DM [39, 40, 41] and DM-electron scattering experiments [20]. In the most extreme case underground detectors might lose detection sensitivity altogether, since the rock of the Earth crust above the laboratory, meant to serve as a background shield, effectively blinds the experiment via the nuclear stopping and screening of the incoming DM particles. This possibility has been studied in [42, 43, 44] and more recently in the context of DM-electron scatterings [45]. Experimental efforts to observe a possible diurnal modulation have been carried out in the early 90s by the COSME-II detector [36] and later by the DAMA collaboration [46, 47]. Both experiments were located in the northern hemisphere and failed to find any evidence for diurnal modulations. A promising future experiment in the southern hemisphere is the SABRE experiment at the Stawell Underground Physics laboratory [48, 49], which is designed to test the DAMA/LIBRA discovery claim [50] and whose location would be significantly more sensitive to diurnal modulation due to Earth scatterings.

The detection signature of Earth scatterings has recently been investigated using analytic methods [51]. Therein the authors published the EarthShadow code, which allows to quantify the lab-frame DM distribution distortions in the case, where the DM particles scatter at most once before passing through the detector. Therefore this approach is restricted to the single-scattering regime. For the case of multiple scatterings numerical methods become necessary. In this paper we present the Dark Matter Simulation Code for Underground Scatterings (DaMaSCUS), a Monte Carlo simulation code for individual DM trajectories, which allows to calculate local distortions of density and velocity distributions for any number of scatterings. In this sense DaMaSCUS complements and generalizes the EarthShadow code, which we use as a crucial consistency check for the new MC simulations. The simulations take the Earth’s orientation in the galactic frame into account, as well as its composition and layer structure and time-dependent velocity through the halo, while it orbits the Sun. By simulating billions of particles we investigate how underground DM-nucleus scatterings affect the local DM velocity distributions by statistical means. This information may then be used to compute time dependent distortions of recoil spectra and diurnal modulations of event rates for any experiment while precisely accounting for its location, underground depth and search strategy. It is even possible, although computationally more expensive, to investigate the Earth’s crust screening effect for very strong DM-nucleus interactions.

We should stress that this is an important effect that can play an important role if DM is light and sufficiently strongly interacting. Underground scatterings could render light DM completely unable to produce detectable nuclear recoils at the usual $\sim$ 1 km depths of most current detectors, leaving a hole in the DM parameter space which will not be covered neither with larger exposure nor with lower energy thresholds. Shallow-site or surface detectors looking for a diurnal signal might be the only working strategy for discovering DM in this part of the parameter space [44].

This paper is organized as follows. After a brief review of basic DM-nucleus scatterings in section 2 we describe the MC simulation algorithm of DaMaSCUS in section 3. Our main results are presented in section 4, before we conclude and give an outlook on future steps in section 5. In addition we provide a set of extensive appendices for the interested reader, which contain short reviews of necessary astronomical relations and computational details of our simulations.

Together with this paper we also make the DaMaSCUS v1.0 code publicly available, together with documentation and some illustrative videos [52].

2 DM Scatterings on Terrestrial Nuclei

The central physical process of interest are elastic collisions of DM particles on nuclei of the Earth’s bulk mass. As the particles pass through the Earth’s mantle and core they may interact with matter depending on their interaction cross-section. In this section we review the underlying dynamics and probabilities. First of all, in order to describe the DM particle’s underground motion we need to model the Earth and its layer structure. The Earth’s density increases with the depth, such that we have more nuclei near the core for the DM particle to scatter on. We implement the mass density profile of the Preliminary Reference Earth Model (PREM) [53], which separates the Earth into 10 distinct layers. Furthermore the chemical composition changes depending on the underground depth as well. We distinguish two different compositional layers, the core and the mantle, for each of which we implement the 9 and 14 most abundant nucleus species respectively [54]. The details are summarized in appendix B.

We consider a DM particle moving through matter. The probability for a particle of velocity $\vec{v}$ to scatter on some nucleus, after freely travelling a distance $L$, is

[TABLE]

In general the mean free path $\lambda_{\rm MFP}$ is a local and speed dependent property and given by

[TABLE]

where we sum over all nuclear isotopes present at $\vec{x}$. $n_{A_{i}}(\vec{x})$ is the number density of the nucleus with atomic mass number $A_{i}$. Substituting the Earth’s local chemical abundances, or the mass fractions of the different atomic species $f_{A_{i}}(\vec{x})$, as well as the density profile of the Earth $\rho_{\oplus}(\vec{x})$, we can rewrite Eq. (2.3)

[TABLE]

In figure 1 we show the average underground mean-free-path together with experimental constraints, which show that for sub-GeV DM the possibility of multiple Earth scatterings is still very viable.

Suppose the DM particle scatters at $\vec{x}$, the probability to scatter on a certain nucleus species $j$ is given by

[TABLE]

The scattering cross-section $\sigma_{\chi A}$ between a DM particle and a nucleus of mass number $A$ generally depends on the momentum transfer $q$ at least through a form factor that accounts for loss of coherence. For high momentum transfers with de Broglie wavelengths below the nucleus’s size, DM cannot scatter coherently with all the nucleons composing the nucleus [56]. This loss of coherence is taken care of by a form factor $F_{A}(q^{2})$ via

[TABLE]

The total scattering cross-section is obtained, by averaging over all possible momentum transfers/recoil energies,

[TABLE]

For sub-GeV DM it will not be necessary to take the loss of coherence into account111However an approximative Helm form factor is implemented in DaMaSCUS and may be used for the simulation of heavier DM. since $F_{A}(q^{2})\approx 1$, and hence $\sigma^{\rm total}_{\chi A}\approx\sigma_{\chi A}(0)$. But it is interesting to note that the form factor makes even the spin-independent total scattering cross-section velocity dependent via $q_{\rm max}$.

In modelling the DM-matter interactions we take the bottom-up framework of non-relativistic effective theory [57]. However in this work we only present results for the first operator, better known as spin-independent DM-nucleon interactions. A MC exploration of the other operators in analogy with [51] will follow in a later publication. Hence for now we exclusively consider isospin non-violating spin-independent interactions, for which the zero momentum transfer cross-section is given as

[TABLE]

with the DM-nucleon cross-section $\sigma^{\text{SI}}_{\chi n}$ and the corresponding reduced mass $\mu_{\chi n}$.

3 Monte-Carlo Simulations with DaMaSCUS

Having covered the basics we introduce the Dark Matter Simulation Code for Underground Scatterings (DaMaSCUS). DaMaSCUS performs simulations of individual particles traversing through the Earth’s mantle and core undergoing scatterings on terrestrial nuclei, which deflect and decelerate the particle. It accounts for the changing composition and density throughout the Earth as well as its motion in the DM halo and its orientation in the galactic frame. A statistical sample of trajectories can be analysed to give precise estimates of the modified local DM number density and velocity distribution, which for any given direct detection experiment will be time-dependent. This allows to compute the local signal rate and its diurnal modulation for any specific experiment, similarly to the recent EarthShadow code [51]. There diurnal modulations have been computed under the assumption that the DM particle scatters at most once on terrestrial nuclei before reaching a detector. Using MC simulations we are not restricted to the single scattering regime and can generalize these findings by simulating particles with any number of underground scatterings. We compare our MC simulation results with the ones of the EarthShadow in the single-scattering regime where the latter is valid as an extra consistency check of DaMaSCUS.

The basic idea of the MC simulation is to follow individual DM particle on their journey through Earth’s interior, as they scatter on terrestrial nuclei resulting in trajectories not unlike a random walk. By recording how the scatterings diffuse DM particles underground, we can derive precise estimates of the DM density and velocity distributions at the location of detectors of interest.

As we will demonstrate below, the particles are sent underground with appropriate initial conditions $(t_{\rm ini},\vec{x}_{\rm ini},\vec{v}_{\rm ini})$ and assumed to move on a straight line until hitting a nucleus. In order to find the distance $L$ a particular particle travels freely we employ (2.1) and solve the equation

[TABLE]

where $\xi$ is a uniformly distributed random number. Then we define the displacement vector $\vec{\Delta}(\vec{x},\vec{v})$ as

[TABLE]

where $\vec{e}_{v}$ is the unit vector in the direction of $\vec{v}$. This vector points from the particle’s original position $\vec{x}$ to the nucleus on which it scatters. The solution of (3.1) is found by an algorithm which combines analytic and numerical methods, for details we refer to appendix C.

Now that the location of the first scattering event is known, the particular nucleus $A$ involved in the scattering can be inferred from (2.4). The particle will deflect and decelerate and its resulting velocity after the scattering is given by a simple relation for elastic collisions,

[TABLE]

Here the only unknown part is $\vec{n}$, the unit vector pointing into the direction of the DM particle’s velocity after the scattering in the center-of-mass-frame. We define the scattering angle $\alpha=\sphericalangle(\vec{v}_{\chi},\vec{n})\in[0,\pi]$ as the angle between the incoming and outgoing direction of the DM particle in the CMS-frame. For spin-independent cross-sections there is no preferred value for $\alpha$. However we also note that by including the form factor, the scattering angle $\alpha$ will no longer be uniformly distributed even for the SI case. Instead we solve

[TABLE]

for $q$, where $\xi$ again is a uniformly distributed random number. The scattering angle is then given by

[TABLE]

In the case of sub-GeV DM, $q^{2}\approx\xi q_{\rm max}^{2}$, and we obtain a uniform distribution in $\cos\alpha$.

The procedure repeats itself as the DM particle continues with its new velocity from the position of the scattering. Again we have to solve (3.1) and find the next scattering location and nucleus, and so on. This is repeated until the particle reaches the Earth surface again or until the velocity of the particle drops below a given speed threshold. The algorithm is summarized by a flow chart in figure 2. It returns a list of events, which make up the DM particle’s trajectory.

The choice of initial conditions for the simulated DM particles is critical. To find the initial time is trivial, we can set $t_{\text{ini}}$ to a random value or just start at $t_{\text{ini}}=0$. The initial velocity is straight forward as well and has two components,

[TABLE]

The first term is the velocity component in the galactic rest frame $\vec{v}_{\text{halo}}$, for which we choose the Standard Halo Model (SHM),

[TABLE]

where $N_{\text{esc}}=\pi v_{0}^{2}\left(\sqrt{\pi}v_{0}\text{Erf}\left(\frac{v_{\text{esc}}}{v_{0}}\right)-2v_{\text{esc}}\exp\left(-\frac{v_{\text{esc}}^{2}}{v_{0}^{2}}\right)\right)$ is the normalization constant and $\Theta(x)$ the Heaviside step function. We use the standard parameter $v_{0}=220$km/s and $v_{\rm esc}=544$km/sec. The second component is the Earth’s velocity $\vec{v}_{\oplus}(t)$ relative to the galactic rest frame, causing the “DM wind”. This velocity is time dependent and changes over the course of one year, giving rise to annual modulation of detection signals [2, 58]. For details on its determination we refer to appendix A.

Once the initial time and velocity are determined, the question of the initial position is more delicate. For one of course, the particle should start its trajectory outside the Earth and move towards the planet. But more importantly, the initial positions have to be distributed uniformly in space. They can not be chosen arbitrarily, e.g. simply on top the Earth’s surface as done in [36, 59, 37, 60, 38]. This point is subtle but crucial. Simulating particles with initial velocities given by (3.6), which exclusively start their trajectory on the surface of the Earth, creates a finite volume bias at shallow depths, i.e. exactly where detectors are located. Over proportionally many particles are sent into the Earth with narrow angles. We have checked that in the case of a transparent Earth, the aforementioned choice of initial conditions results in a DM over density close to the surface, in contrast to the expected uniform DM density. Apart from the numerical verification, it can also be shown analytically that the aforementioned choice of initial conditions does not describe the real situation.

Instead an effectively uniform distribution of the initial positions is realized by choosing a random point on a circular disk of radius $r_{\oplus}$ at a distance $R$ from the Earth center and perpendicular to $\vec{v}_{\text{ini}}$, see figure 3.

[TABLE]

where $\xi\in(0,1)$ and $\phi\in[0,2\pi)$ are uniformly distributed random numbers, and $\vec{e}_{x}$ and $\vec{e}_{y}$ span the disc. Together with a random starting time $t_{\text{ini}}$, this is equivalent to the choice of a random point inside a cylinder with radius $r_{\oplus}$, orientated parallel to $\vec{v}_{\text{ini}}$. This ensures the effectively uniform distribution of initial positions in space.

The halo DM velocity distribution (3.7) is isotropic. Boosting our system into the frame of the Earth via (3.6) breaks the isotropy, but nonetheless our system still has a rotational symmetry around the axis parallel to the Earth’s velocity $\vec{v}_{\oplus}(t)$. We call the polar angle of this axis isodetection angle $\Theta$. Note that in [51] and in the EarthShadow code the authors use an equivalent angle defined as $\gamma=180^{\circ}-\Theta$. As the name suggests, along a constant isodetection angle the DM particle’s velocity distribution as well as direct detection event rates will also be constant. We exploit this symmetry for our MC simulations to define small but finite-sized isodetection rings as done in [36, 37, 38], see figure 4. However we place these rings not at the Earth’s surface, but at a finite underground depth, which can be adjusted for different experiments e.g. to 1400 meter for the LNGS. This way the MC simulation can also be used to investigate the Earth’s crust shielding effect for very large interaction cross-sections. The effect of nuclear stopping on the sensitivity of underground direct detection experiments has been studied in [43, 44, 61] and lately in the context of DM-electron scattering experiments using MC simulations in [45].

We embed a spherical surface inside the Earth with radius $r_{\oplus}-d_{\rm lab}$, where $d_{\rm lab}$ is the depth at which the direct detection experiment of interest is placed underground. A detector with his fixed coordinates will travel through the isodetection rings in a non-trivial way, see appendix A.4. The position of a detector in terms of the isodetection angle $\Theta$ is given by

[TABLE]

where we have to use (A.22) and (A.27). Note that $\vec{v}_{\oplus}(t)$ only changes marginally over a matter of days. We show the evolution of (3.9) for different laboratories around the globe in figure 5. It illustrates how different experiments move through the isodetection rings and already hints at, which experiments will be more sensitive to the Earth’s shadowing effect and the corresponding diurnal modulations.

More concretely, we divide up the detector sphere into 180 isodetection rings with width $\Delta\Theta=1^{\circ}$ symmetric about the axis defined by $\vec{v}_{\oplus}(t)$. We label the rings with $\Theta_{k}$ where $k\in[0,179]$. The area of ring $\Theta_{k}$ is then

[TABLE]

We choose $\Delta\Theta$ sufficiently small, such that the velocity distribution will be approximately constant over a single isodetection ring’s surface. The final goal of our MC simulations is to derive a precise estimate for the local speed distribution function for each of these rings based on the individual simulated particles including statistical uncertainties.

For each simulated trajectory we record if, how often, where (in terms of isodetection ring $\Theta_{k}$), and with what velocity the particle passes through the detection sphere. The different isodetection rings have varying surface areas, see (3.10). Furthermore the particle flux will also vary with $\Theta$ because of the DM wind. As a consequence more particles pass through certain rings than others. Yet we want to have the same statistics for each ring. We therefore repeat simulating DM trajectories until we have collected the same velocity data sample size $N_{\rm sample}$ for each ring. This value must of course be chosen sufficiently large. This way we accumulate $180\times N_{\rm sample}$ DM velocity vectors. For each isodetection ring we can determine the local distribution functions, DM densities and event rates independently.

In order to get an accurate estimate for the speed distribution function we need a non-parametric density estimation method. We employ histograms, weighting the data properly, in order to estimate the distribution function. For details on the distribution estimates we refer to appendix D. There we also show in detail how direct detection rates can be derived from the MC data. In this work we take a CRESST-II like detector as an illustrative example. Our framework can however easily be extended to other detector types or detection strategies.

Early MC simulations of a similar kind have been performed in the 90s [36, 37, 38]. The focus of these works laid on the standard WIMP model with DM masses of 50 GeV and higher and direct detection using nuclear recoils, whereas we focus on sub-GeV DM and have new detection techniques in mind. Apart from this there is a number of essential differences in the implementation of DaMaSCUS, two of which we emphasize here. For one we employ the corrected method of finding appropriate initial conditions for the DM particles, as has been described above. This ensures a spatially uniform distribution inside the Earth and avoids density overestimation close to the surface. Secondly we embedded the isodetection rings underground instead of directly on the Earth’s surface, so we can simulate cases, where the rock above an experiment acts as a DM shield. Naturally a very strong interaction between DM and nuclei is required for this to occur. Another minor difference is a more refined modelling of the Earth.

4 Results

This study investigates the terrestrial effect of DM-nucleus underground scatterings in the regime of sub-GeV mass and strong enough (yet unconstrained) DM-nucleus cross sections that allow single or multiple scatterings before the DM particle arrives at the detector. We start by presenting simulation results in the case where DM can scatter at most once, because in this case we can directly compare the MC results to the analytic ones obtained in [51] with the EarthShadow code. As we will demonstrate our MC results agree perfectly with the ones of  [51], reinforcing our confidence in our simulation code as we use it afterwards in cases of multiple DM scattering where the EarthShadow code is no longer valid.

We choose four benchmark points for a DM mass of 500MeV, which lies directly at the boundary of the detection sensitivity of CRESST-II. The lowest cross-section used corresponds to an underground scattering probability of $10\%$. The MC simulation for this point will be compared directly to the analytic results of [51], given that a $10\%$ probability results at most to a single scattering (double scattering has a probability of $\sim 1\%$). The other three cross-sections are higher and tuned such that the average number of underground scatterings is 1, 10 and 50 respectively. For a summary see table 1.

The resulting DM density and velocity distribution functions can be used to compute direct detection event rates for any kind of experiment and search strategy.

The computations were performed on the Abacus 2.0, a 14.016 core supercomputer of the DeIC National HPC Center at the University of Southern Denmark and typically involve the simulation of up to $10^{11}$ DM particle trajectories.

4.1 Single Scattering: Comparison to EarthShadow

In order to perform a vital consistency check, we compare the results of the analytic methods of the EarthShadow code with our MC simulation. Since the analytic approach can only cover the single-scattering regime, we choose a DM-nucleon cross-section corresponding to an underground scattering probability of $10\%$. We can employ the respective EarthShadow routine and find $\sigma_{\chi n}^{\rm SI}(0)\approx 0.5$pb for a DM mass of 500 MeV. For these parameters we find good agreement with the MC simulations, where $\sim 90\%$ of the simulated particles cross the Earth freely, while the remaining scatter at least once.

As mentioned earlier, when comparing the two approaches we should keep in mind that deviations of the order of a percent should be expected, since about a percent of the particles will scatter not just once but twice, which is only accounted for by DaMaSCUS. These small deviations will affect the large $\Theta$ region, where DM particles have travelled the longest distances underground.

The distortions of the DM distribution result from deceleration and deflections of DM particles due to underground scatterings, which redistribute the DM inside the Earth. In the left panel of figure 6 we see the speed distribution functions for $\Theta=0^{\circ}$ and $180^{\circ}$, which are normalized to a DM density of $0.3\text{ GeV cm}^{-3}$. Already in the single-scattering regime it is clear that for high values of $\Theta$ the DM density gets reduced, because DM particles travel longer distances underground and scatter away from their original path more likely. Particles entering the Earth at low values of $\Theta$ on the other hand have to cross only short underground distances of the order of the detector depth $d_{\rm det}$ and will most likely not scatter before reaching this depth. However other particles which originally were not on a path towards this region may be deflected and still end up here, leading to an overall increase in the local DM density. Overall the distortion of the DM distribution is due to deflections rather than the deceleration of DM.

The main outcome of this comparison however is the excellent agreement and consistency of the local DM distribution between DaMaSCUS and EarthShadow. Therefore it doesn’t surprise that this agreement also translates to the direct detection event rate shown in the right panel of figure 6. It clearly shows the Earth’s ‘shadow’, the decrease of the signal rate for large values of $\Theta$. The decrease is mostly caused by the locally depleted DM population due to deflections. The effect of DM deceleration is less severe in this case, but might still be crucial for detection strategies relying on DM particles from the tail of the velocity distribution.

Detectors, whose revolution around the Earth involve high isodetection angles would be able to measure this decrease as a diurnal modulation. As seen in figure 5 this is a signature pronounced more in laboratories of the southern hemisphere. The experimental facilities beneath Gran Sasso for example only cover roughly $\Theta\in(0^{\circ},90^{\circ})$, where the signal rate is slightely increased but effectively constant.

One should mention here that there are two more sources of diurnal modulation in the DM signal that can compete or exceed the shadowing effect in low enough DM-nucleus cross sections. One is related to the rotational velocity of the Earth around its own axis ($\sim 1\text{km}/\text{s}$), which is superimposed with the total velocity of the Earth in the rest frame of the Galaxy, causing a small daily fluctuation of the DM flux of the order $\sim 10^{-4}$ at the detector. The second source is due to gravitational focusing, i.e., the fact that the Earth works as a gravitational lens causing small fluctuations on the DM density depending on the relative position of the detector with the center of the Earth and the direction of the DM wind [62]. All three sources of modulation have different amplitudes and phases, producing a final modulation which is the superposition of all.

4.2 Multiple Scatterings

Now that consistency between the analytic and the MC results is established in the single-scattering regime we can confidently turn to higher cross sections of the multiple-scatterings or diffusion regime. We investigate three benchmark points, i.e., cross sections tuned to result in 1, 10 and 50 underground scatterings on average, see table 1.

We start off with the local DM speed distribution distortions depicted in figure 7. Compared to the single-scattering regime results, the two main observations are the much more severe deceleration and depletion of the DM population for higher values of $\Theta$, i.e. deep in the Earth’s shadow. Especially for the benchmark point ‘MS50’ we see clearly how underground scatterings deplete regions with high speeds and increase the slow DM population significantly. The effect is less pronounced in the other benchmark points, but it can be inferred by observing a horizontal shift in the peak of the distributions as $\Theta$ changes. On the other hand the overall height of the DM speed distributions functions show how the deflections of DM particles enhance the DM density for locations facing the DM wind, while simultaneously deplete the DM population for large values of $\Theta$. We again see the same effects as in the single-scattering regime but in a more severe way. A new feature for cross sections of the order of $\mathcal{O}(100\text{pb})$ is the second peak in the lower right panel on figure 7, populated by very slow particles. Although current detector thresholds cannot probe this second peak, more sensitive future detectors could observe this peak as a bump in the low energy recoil spectrum. Something similar has been claimed to be produced by DM gravitationally bound to the Earth [63, 64].

The local DM velocity distribution functions in figure 7 are the central output of DaMaSCUS, since they allow us to compute direct detection rates of any kind. They also encode the local average speed and DM density, which we plot separately in the first two panels of figure 8 for the four benchmark points. Depending on the interaction strength the average speed is decreased through nuclear stopping with a clear minimum for large $\Theta$. We observe a slight but noticeable increase for the isodetection rings between 60 and $90^{\circ}$. It is most notable for the ‘MS50’ benchmark point, but it is also visible in ‘MS1’, where we see that the average speed is increased slightly even beyond the expectation for free particles. It could be explained by particles from the fast velocity tail of our initial distribution which typically enter the Earth parallel to the DM wind at low values of $\Theta$ and deflect towards high-$\Theta$ isodetection rings.

In addition another impact of the elastic DM-nucleus collisions is the redistribution of DM particles. The local DM density, while being increased by up to 100% at $\Theta=0^{\circ}$ for the benchmark point with the strongest interaction, drops to values below the halo density as we move further into the Earth’s DM shadow. The decrease of both speed and density results to a decrease of direct detection events, as we see in the left lower panel of figure 8, where we again took a CRESST-II like detector as a n example. The event rate drops significantly more than in the single-scattering regime. Therefore the overall local DM density is higher in the northern hemisphere compared to the southern, with an increasing difference as we move to stronger DM-nucleon interactions.

The thereby induced diurnal modulation of the signal rate for a detector of fixed location on Earth traversing the isodetection rings in accordance with Eq. (3.9) is shown in the lower right panel of figure 8, where we selected the SUPL ($37.07^{\circ}$S,$142.81^{\circ}$E) detector in Australia as its location in the southern hemisphere is very sensitive to diurnal modulations of this kind. However such modulations are not a signature reserved only for experiments in the South. To see the dependence on an experiment’s latitude, we define the percentile signal modulation $\delta$ observable at a given laboratory via

[TABLE]

where $\Phi_{\rm lab}$ is the laboratory’s latitude. We show the modulation as a function of the latitude in figure 9 for the four benchmark points.

It is generally true that the diurnal modulation caused by underground scatterings is maximal for experiments in the southern hemisphere. In the multiple-scattering regime however we see that such modulations can be significant almost anywhere on the globe with the exception of the poles’ neighbourhood. Even for experiments at Gran Sasso we find a diurnal signal modulation of 10%, 45% and almost 60% for $\langle N_{\rm sc}\rangle=1$, 10 and 50 respectively. Nevertheless a location such as the SUPL is strongly favoured with modulations of $\sim 18\%$, $\sim 65\%$ and more than $95\%$ for the same set of cross sections. For the two benchmark points with the largest cross-section we find that a potential sub-GeV direct detection experiment in the southern hemisphere should expect a vastly reduced signal rate during a significant part of the day while the planet’s bulk mass shields off the sought particles. In this case we have significant sensitivity to DM during certain hours of the day only.

5 Conclusions and Outlook

In this paper we study the effect of underground DM-nucleus scatterings on the direct detection prospects of sub-GeV DM with sufficiently large cross sections that could accommodate such pre-detector underground scatterings. This effect can be important in various contexts. One is the possibility of a subdominant strongly interacting component of DM not abiding by direct detection constraints. Additionally, DM in the sub-GeV region is practically unconstrained by current direct detection underground experiments. Since in this region, DM can have sufficiently strong DM-nucleus interactions, underground scatterings before DM particles reach the detector must be taken into account. At best, these underground scatterings can distort the nuclear recoil spectrum making hard to establish beyond any doubt DM discovery. In the worst case scenario, detectors at current deep sites might be completely blind to this part of DM parameter space regardless of exposure and even if they lower significantly their energy thresholds. Similarly for DM constraints based on DM-electron interactions in the same mass region, one should consider carefully the terrestrial effect from nuclear stopping before reaching the detector.

In order to study this effect in full generality, we developed the DaMaSCUS code where we can simulate the trajectory of halo DM particles of given cross section and mass that cross the Earth, potentially scatter underground and eventually scatter in the detector. In order to do this, we use a state-of-the-art density and composition profile of the Earth and improve the generation of initial conditions from simulations done in the past. The output of our simulation is the local DM density and velocity distribution of DM at the location of a given underground detector, time and day of a year.

For DM-nucleon cross section where DM particles can scatter at most once before reaching the detector, the results of our MC simulation are in excellent agreement with the analytical calculation of  [51]. However our simulations have the big advantage that they can be used for higher yet experimentally allowed DM-nucleon cross sections, where DM particles can possibly scatter more than once. In fact we presented several benchmark points where the average number of underground scatterings is 1, 10 and 50. As expected, we found that with increasing cross section the local DM density and velocity distribution deviate more and more from the transparent Earth, where the effect of underground scatterings is ignored completely. We found that this affects strongly the number of potential events in an underground detector. This is due to two effects: DM deceleration and DM deflection due to underground scatterings. We found that the effect can be so strong, that a potential DM signal in that parameter space will have a large diurnal modulation simply because as the Earth rotates around its own axis, DM particles travel different distances in order to reach the detector, thus increasing the probability/number of underground scatterings. We predicted the amount of diurnal modulation for different labs and we verified that the southern hemisphere has larger daily fluctuations in the DM signal. Note that for sufficiently strong cross section, the diurnal modulation can be next to $\sim 100\%$. Our simulation can also provide real time correlations between DM signals of detectors in different latitudes, facilitating the task of discriminating DM from potentially other backgrounds with daily modulation. Our study revealed also another interesting feature. For sufficiently large cross section, the velocity distribution acquires a second sharp peak at low energies. Although this is currently out of experimental reach, with the advent of detectors with lower thresholds, this peak could create a very distinct feature in the low energy recoil spectrum, identifying DM beyond any doubt.

We leave a lot of things for future work. We plan to extend our work by studying the shadowing effect for different types of DM-nuclei interactions beside the spin-independent one as well as for DM-electron interactions. The latter is extremely interesting since the core of the Earth is assumed to be metallic. This means that one should treat some of the atomic electrons as free particles. It should be stressed here that free electrons do not require a minimum energy to excite. They can absorb even tiny amounts of energy, thus decelerating even low energetic DM particles. In addition we intend to investigate the shadowing effect in the case of long range forces between DM and nuclei, e.g. in models where the interaction is mediated by a light dark photon. In this case DM particles might interact collectively with multiple atoms as they pass through, losing energy perhaps in a similar fashion as massive objects lose energy through gravitational dynamical friction. We also plan to include in our simulation the other two sources of diurnal modulation, i.e., due to gravitational focusing and the rotational velocity of the Earth. Moreover we plan to study the shadowing effect on directional detectors. A first attempt was presented in [61]. Finally we plan to use the simulation in order to make a precise estimate of the effect of DM gravitationally bound to the Earth on the nuclear recoil spectrum of underground detectors. Although there is a promising analytical estimate [63, 64] that can be explored by future detectors, a precise MC simulation will give a more accurate estimate about the density and velocity distribution of bound DM in the Earth.

Acknowledgments

This work is partially funded by the Danish National Research Foundation, grant number DNRF90 and by the Danish Council for Independent Research, grant number DFF – 4181-00055. Computation/simulation for the work described in this paper was supported by the DeIC National HPC Centre, SDU.

Appendix A Astronomical Conventions and Coordinates

For the sake of completeness and in order to serve as a reference we review the astronomical basics required for our simulations.

A.1 Sidereal Time

In order to keep track of the Earth’s rotation, we use a time unit based on that very rotation called sidereal time. A sidereal day is the time interval of one rotation relative to vernal equinox $\Upsilon$. It is slightly shorter than a mean solar day: $23.9344699$ hours [65]. Note that the sidereal time is often given as an angle instead of a time unit.

The local apparent sidereal time (LAST) is the time since the local meridian passed $\Upsilon$. In this section we will show how to compute the LAST for any time and location. But first, we have to introduce a reference time. The time in our problem are always given relative to 01.01.2000 12:00 TT (or GMT), a commonly used reference time denoted as J2000.0. In order to calculate the fractional number of days $n_{\text{J2000.0}}$ relative to J2000.0 for a given date $D.M.Y$ and time $h:m:s$ (UT) we use the following relation [66],

[TABLE]

and $\lfloor\cdot\rfloor$ is the floor function. For example the 15.02.2016, 8:00am corresponds to $n_{\text{J2000.0}}=3332.83$. We also define the epoch,

[TABLE]

Next we show in detail how to compute the Local Apparent Sidereal Time (LAST) anywhere on Earth starting from the Universal Time. We start with the formula for the Greenwich mean sidereal time (GMST) in seconds [65],

[TABLE]

For the Greenwich apparent sidereal time (GAST) we’ll have to add the equation of equinoxes,

[TABLE]

To obtain the local apparent sidereal time (LAST) at a location with latitude and longitude $(\Phi,\lambda)$, we just add the longitude,

[TABLE]

Note that for western longitudes, $\lambda$ is negative. It should always be made sure that $\text{LAST}\in(0,86400)$. Comparing to the tables of [65], the errors of these approximations are of the order of tens of milliseconds.

A.2 Coordinate Systems

We introduce the relevant coordinate systems and how to transform in between them [66, 65]. All coordinate systems are rectangular and right-handed.

the galactic frame ‘(gal)’: A heliocentric coordinate system, the $x$-axis points towards the galactic center, the $z$-axis points towards the galactic north pole. The $x$- and $y$-axis span the galactic plane. 2. 2.

the heliocentric, ecliptic frame ‘(hel-ecl.)’: A heliocentric coordinate system, the $x$-axis points towards vernal equinox $\Upsilon$, the $z$-axis points towards the ecliptic north pole. The $x$- and $y$-axis span the ecliptic plane. See figure 11. 3. 3.

the geocentric, ecliptic frame ‘(geo-ecl.)’: A geocentric coordinate system, the $x$-axis points towards vernal equinox $\Upsilon$, the $z$-axis points towards the ecliptic north pole. The $x$- and $y$-axis span the ecliptic plane. 4. 4.

the geocentric, equatorial frame ‘(equat)’: A geocentric coordinate system, the $x$-axis points towards vernal equinox $\Upsilon$, the $z$-axis points towards the Earth north pole. The $x$- and $y$-axis span the equatorial plane. See figure 11. 5. 5.

the laboratory frame ‘(lab)’: A coordinate system with the detector in the origin. The $x$-axis points towards east, the $y$-axis towards north and the $z$-axis to the sky.

Since all our calculations will in the end be done in the galactic frame, we need the transformation matrices, which are time-dependent.

1.) (lab) $\longleftrightarrow$ (equat):

The rotation from the laboratory frame to the equatorial frame is done by

[TABLE]

where $\theta=\frac{\pi}{2}-\Phi$ and $\phi=2\pi\frac{\text{LAST}(\Phi,\lambda)}{86400\text{s}}$.

2.) (hel-ecl.)$\longleftrightarrow$(geo-ecl.):

The simplest transformation is the one between the two ecliptic frames,

[TABLE]

3.) (geo-ecl.) $\longleftrightarrow$ (equat):

To transform a vector $\vec{x}^{\text{(geo-ecl)}}$ to equatorial coordinates, the necessary rotation is

[TABLE]

where $\epsilon=23.4393^{\circ}-0.0130^{\circ}T_{\text{J2000.0}}$ is the obliquity or axial tilt of the ecliptic.

4.) (equat)$\leftrightarrow$(gal):

The equatorial frame at J2000.0 can be related to the galactic frame,

[TABLE]

The three angles, namely the J2000.0 right ascension of the north galactic pole $\alpha_{\text{GP}}$, the J2000.0 declination of the north galactic pole $\delta_{\text{GP}}$ and the longitude of the north celestial pole in J2000.0 galactic coordinates $l_{\text{CP}}$, are

[TABLE]

To rotate a vector in equatorial coordinates at J2000.0 to any time epoch $T_{\text{J2000.0}}$, we use

[TABLE]

The equatorial precession angles are

[TABLE]

Summary:

In order to transform from any frame to any other frame, we just multiply the corresponding rotation matrices. We can follow this flow chart, where we take the inverse matrix if opposing the arrow’s direction.

(hel-ecl.)(geo-ecl.)

(equat)

at $T_{\text{J2000.0}}$

(lab)

(equat)

at J2000.0

(gal)-\mathds{1}$$\mathcal{R}$$\mathcal{N}$$\mathcal{P}$$\mathcal{M}

For example, to go from equatorial to galactic coordinates for any given epoch $T_{\text{J2000.0}}$ we use

[TABLE]

Now it is easy to write the axis vectors $\vec{e}_{x}=(1,0,0)^{T}$, $\vec{e}_{y}=(0,1,0)^{T}$ and $\vec{e}_{z}=(0,0,1)^{T}$ of the equatorial and heliocentric-ecliptic frame transformed into the galactic one,

[TABLE]

The axis vectors of the heliocentric-ecliptic frame are

[TABLE]

This completes our review on coordinate systems.

A.3 Earth’s Velocity in the Galactic Frame

The standard reference for the Earth velocity in the context of direct detection has long been the review by Smith and Lewin [56]. However as pointed out first in [67] and confirmed in [66] there is an error in the first order correction due to the orbit’s eccentricity. We take the full expression for the Earth velocity vectors from [66]. The Earth’s velocity with respect to the galactic rest frame is the sum of three components,

[TABLE]

These components are

the galactic rotation,

[TABLE] 2. 2.

the sun’s motion relative to nearby stars,

[TABLE] 3. 3.

the Earth’s orbital velocity $\vec{v}_{e}(t)$ relative to the sun, visualized in figure 11,

[TABLE]

where the unit vectors $\vec{e}_{i,\text{hel-ecl}}^{\text{(gal)}}$ are given in Eqs. (A.19) and (A.20).

Finally we list the mean velocity $\langle v_{\oplus}\rangle$, the Earth’s eccentricity $e$, the mean longitude $L$ as well as the longitude of the perihelion $\varpi$,

[TABLE]

The inclusion of the correction due to the orbit’s eccentricity may not be relevant for the results reported in this paper. We include them for the sake of completeness and potential future applications of the simulation code.

A.4 Laboratory Position and Velocity

As a first step we find the spherical coordinate angles $(\theta_{\text{lab}},\phi_{\text{lab}})$ of the detector’s position in the geocentric equatorial coordinate system. We specify the location of a detector in the Earth through the latitude and longitude $(\Phi_{\text{lab}},\lambda_{\text{lab}})$ and the underground depth $d_{\text{lab}}$ of the laboratory.

The x-axis of the equatorial coordinate system points towards the vernal equinox or March equinox. Therefore the spherical coordinate angles are given by

[TABLE]

Since we use sidereal seconds, the rotation frequency is simply $\omega_{\text{rot}}=\frac{2\pi}{86400s}$. Now that we have the lab’s spherical coordinates we can transform the position vector into galactic coordinates,

[TABLE]

The velocity component of the laboratory due to Earth’s rotation is given by

[TABLE]

The unit vectors are given by (A.16) and (A.17). The rotation velocity at the equator is $v_{\rm eq}\approx 0.465\text{ km s}^{-1}$.

Appendix B Earth Model

We model planet Earth by dividing it into two sets of layers. On the one hand we have the compositional layers, which differ by their chemical element abundances. Based on [54] we implement two compositional layers, the core and mantle. On the other hand we have ten mechanical layers, which, for our purpose, only differ by their density profile. We adopt the density profile of the Preliminary Reference Earth Model (PREM) [53]. The chemical abundances of the two compositional layers are listed in table 2. The core has a radius of about 3480km.

For the mechanical layers the Earth’s density for each layer $l$ is parametrized as

[TABLE]

The coefficients are given in [53] and listed in table 3, the density profile is plotted in figure 12.

Appendix C Hybrid-Algorithm for the Displacement Vector

In this section we describe in detail how we solve Eq. (3.1) for $L$ using a mixture of analytic and numerical methods. First we can rewrite the equation as

[TABLE]

The straight forward approach to find $L$ and the displacement vector $\vec{\Delta}$ of a DM particle, e.g. in between two scatterings, is to solve (C.1) by numerical integration along the particle’s path. Since it is most likely that a particle, before it scatters or leaves the Earth, passes several layers, some maybe twice, without scattering, it will be reasonable to use the fact that we can solve the integral analytically inside a fixed layer $l$, as we will see now.

The first step towards the analytic solution for a given Earth layer $l$ and DM velocity $\vec{v}$ is to separate the space and velocity dependence of the mean free path in (2.3),

[TABLE]

where

[TABLE]

The factors $g_{l}$ are constant inside each of the two compositional layers, since they only depend on the nuclear composition of the respective layer. As we have seen, in the case of SI interactions and light DM, we can neglect the nuclear form factor. In this case the factors $g_{l}$ lose their dependence on the velocity and only have to be computed once in the beginning of the MC simulations. Otherwise they have to be updated each time the DM particle changes its velocity.

With the parametrisation of the density (B.1) we can solve the integral analytically for a constant layer $l$.

[TABLE]

We used $\cos\alpha=\frac{\vec{x}_{0}\cdot\vec{v}}{x_{0}v}$ and $\tilde{L}\equiv\sqrt{L^{2}+x_{0}^{2}+2Lx_{0}\cos\alpha}$.

In conclusion, it is possible to calculate the probability of scattering after a travelled distance $L$ through a fixed layer $l$,

[TABLE]

where $\Lambda_{l}$ is an analytic function. This fact should be exploited since it will occur fairly often, that a DM particle passes an Earth layer without scattering. In that case we add up the individual terms from each layer, the particle passes through without interacting with the term of the layer where the scattering takes place,

[TABLE]

Here, $t^{l}_{\text{exit}}$ it the time a particle spends inside layer $l$ before leaving. Before a scattering event the particle moves from layer-boundary to layer-boundary and at each layer change we add up the layer’s contribution to $\Lambda_{\text{total}}$. Only if the particle scatters will a numerical method be used in that very layer to calculate the right term. More specifically, in a given layer $l$ there are two possibilities:

The particle passes through the layer without interaction and we jump to the point of layer exit and add the layer’s contribution to $\Lambda_{\text{total}}$,

[TABLE] 2. 2.

The particle does scatter inside layer $l_{s}$, we have to solve

[TABLE]

for $L_{s}$. Looking at (C.4) it’s obvious we have to do that numerically. Since it is a monotonous function, the Newton Method of finding roots will converge quickly.

Using this hybrid algorithm of analytic and numerical methods, which is also depicted in figure 13, will save computation time, especially if the interaction cross-section is chosen small, such that most layers are passed without collisions.

Appendix D From MC to Direct Detection Event Rates

In this appendix we present details of the data analysis and the statistical uncertainties. The DaMaSCUS simulations generate velocity data, which need to be processed to finally obtain local properties like the DM density, velocity distribution or direct detection event rates further down the line.

All of the following steps are performed independently for each isodetection ring.

D.1 Local DM Speed Distribution

We start by describing the procedure to obtain a histogram estimate of the local velocity distribution function for each of the isodetection rings. Even though our simulations in principle provide us with the full velocity data, including the directional information, we focus on the speed distribution for now. It will be sufficient for the non-directional direct detection rate computations of this study, and it is straight-forward to extend the procedure to the full velocity distribution, which we leave for a future treatment of directional detection. Hence the goal here is to estimate the speed distribution,

[TABLE]

As the simulation code tracks a DM particle on its path through the Earth’s bulk mass, it records and saves its velocity $\vec{v}_{i}$ the moment the particle passes through one of the isodetection rings at position $\vec{r}_{i}$. Simulations continue until we accumulate the same statistical sample size of velocity data points $N_{\rm sample}$ for every isodetection ring. This resulting number of simulated particles is large enough that they accurately track the true underlying distribution in the underground phase space inside the Earth and we can obtain good histogram density estimates close to the isodetection ring surfaces.

The next question is how to relate the speed data of the particles passing through the isodetection ring with the local speed distribution in the immediate neighborhood of this surface. By immediate neighborhood we refer to the surrounding volume in sufficient proximity, such that the local distribution function and density therein can be considered constant and particles will not scatter before crossing the boundary. In other word the volume size is assumed to be significantly smaller than the one defined by the local mean free path of the DM particles. We do not actually measure the particle flux in the MC simulation, i.e. the particles are not sent into the Earth continuously in time and their crossing of the rings is not time-tagged, which would be necessary to estimate the particle crossing rate. Instead we fire a number of particles and wait until each particle has run its course, waiting long enough to collect all particles passing the isodetection ring. Obviously faster particles will arrive before slow ones, but since we wait long enough we collect them all. Therefore the speed distribution constructed from these collected particles does not track the flux $\Phi$ but rather $\Phi/v$. Given that the flux is related to the distribution function (see e.g. [68, p.272]) via

[TABLE]

where $n_{\chi}$ is the DM number density, it is clear that the speed distribution of collected particles tracks $f(\vec{v})\cos\gamma$, where $\gamma$ is the angle between the velocity and the normal of the surface at the point of crossing. Hence in order to estimate $f(\vec{v})$ or $g(v)$ each data point needs to be weighted by the reciprocal cosine of the crossing angle $\gamma$, i.e.

[TABLE]

It is clear why we have to weigh our distribution like this in order to get the actual speed distribution. For a patch of the isodetection ring with area $dA$, a given particle sees an effective area $dA\cos\gamma$.

We estimate the distribution function underlying the data in a non-parametric way using histograms. At the beginning of the data processing we divide up the distribution’s domain $(0,v_{\rm esc}+v_{\oplus})$ into $N_{\rm bins}=\lceil\frac{v_{\rm esc}+v_{\oplus}}{\Delta v}\rceil$ histogram bins $B_{1}=[0,\Delta v)$, $B_{2}=[\Delta v,2\Delta v)$,…, $B_{N_{\rm bins}}=[(N_{\rm bins}-1)\Delta v,N_{\rm bins}\Delta v]$. To find a suitable bin width $\Delta v$ we use Scott’s normal reference rule [69],

[TABLE]

The height of the weighted histogram bin $i$ is given by

[TABLE]

where $N=\Delta v\sum\limits_{j=1}^{N_{\rm sample}}w_{j}$ normalizes the histogram, so that $\hat{g}(v)$ truly estimates the probability density function $g(v)$. Furthermore we determine the variance of the bin height $W_{i}$ based on Poisson statistics,

[TABLE]

The average speed of all particles passing a certain isodetection ring is nothing but the weighted mean,

[TABLE]

We also use the standard error (SE) approximation by Cochran (1977) [70],

[TABLE]

We performed a consistency check by running the simulations for a transparent Earth, thus without any DM-nuclei interactions. We retrieved the correct average speed for each isodetection ring, and the histogram estimates of $g(v)$ were statistically stable and successfully reproduced the Maxwell-Boltzmann distribution in perfect agreement with the Standard Halo Model, which was used to generate the initial conditions. We have also checked that our results are robust to changes in the number of simulated particles $N_{\rm sample}$.

D.2 Local DM Density

Next we want to extract an estimate $\hat{\rho}_{\chi}$ of the local DM number density. For this purpose we make use of two observations. The first is the fact that for free DM particles, i.e. $\sigma_{\chi n}=0$, the number density is constant throughout space and simply given by $n^{(0)}_{\chi}=\rho^{(0)}_{\chi}/m_{\chi}$ with $\rho^{(0)}_{\chi}=0.3\,\text{GeV}/\text{cm}^{3}$. Secondly we utilize that without the normalization the area of the histogram (D.6) is directly proportional to the local number density in the MC simulation.

We can therefore perform an initial run of free trajectory simulations without any scatterings and relate this to the main run including scatterings giving us the local density. For both the initial and the main simulation we sum up the weights of all particles passing through a given isodetection ring. The ratio of these two sums is proportional to the ratio of the two densities.

[TABLE]

The superscript $(0)$ denotes the initial simulation run with no scatterings. As mentioned before, the total number of simulated particles $N_{\rm total}$ of the main run is determined by the demand of a common data sample size for all isodetection rings. Hence it will differ from the number of particles simulated in the initial run $N_{\rm total}^{(0)}$. This is taken into consideration via

[TABLE]

The standard deviation is obtained by propagating the errors of $W_{\rm total}$,

[TABLE]

where $\sigma^{2}_{W_{\rm total}}\equiv\sum\limits_{j=1}^{N_{\rm total}}w_{j}^{2}$. Again we point out, that these steps are done independently for each isodetection ring.

D.3 Direct Detection Rates

Our study is relevant in case DM-nuclei cross sections are large enough that pre-detection underground scatterings can occur. These scatterings will change the DM distribution and the predictions for any direct detection experiment and detection channel. To quantify this effect we will have to give precise estimates of the local event rates based on our MC estimates described in the last section. In this work we present the case of nuclear recoil detectors as a first application.

The recoil spectrum for conventional detectors is given by

[TABLE]

Here $X_{A}$ is the target mass fraction of atoms with mass number $A$. The function $\eta(v_{\rm min})$ can be calculated analytically for the Standard Halo Model, see e.g. the preprint of [71]. We already found the local DM density. The second necessary ingredient entering the event rate calculation is a histogram estimate of the $\eta$-function in (D.13), for which we just add up the bin areas222The bin width $\Delta v$ is the same as in (D.4).

[TABLE]

where $W_{j}$ is given in (D.5). This way we obtain a histogram estimate $\hat{\eta}$ for the true $\eta$-function,

[TABLE]

Together with the local DM density (D.11) we are ready to compute the MC recoil spectrum by substituting $\hat{\rho}_{\chi}$ and $\hat{\eta}$ into (D.13). The residual steps towards the total event rate for different experiments does not deviate from the analytic case. This study focusses on sub-GeV DM and we choose a CRESST-II type detector as an example. However we also implement a simplified computation event rate for a LUX-type detector in DaMaSCUS in the case one wants to study heavier DM. Both methods are described in appendix B of [51]. The various integrations are done numerically using the trapezoidal rule. The uncertainty of the resulting integral is obtained by integrating the uncertainty of the integrand, which typically overestimates the statistical error.