Direct Detection of Dark Matter Bound to the Earth
DDirect Detection of Dark Matter Bound to the Earth
Riccardo Catena ∗ Chalmers University of Technology, Department of Physics, SE-412 96 G¨oteborg, Sweden
Chris Kouvaris † CP -Origins, University of Southern Denmark, Campusvej 55, DK-5230 Odense, Denmark We study the properties and direct detection prospects of an as of yet neglected population of darkmatter (DM) particles moving in orbits gravitationally bound to the Earth. This DM population isexpected to form via scattering by nuclei in the Earth’s interior. We compute fluxes and nuclear recoilenergy spectra expected at direct detection experiments for the new DM population consideringdetectors with and without directional sensitivity, and different types of target materials and DM-nucleon interactions. DM particles bound to the Earth manifest as a prominent rise in the low-energypart of the observed nuclear recoil energy spectrum. Ultra-low threshold energies of about 1 eV areneeded to resolve this effect. Its shape is independent of the DM-nucleus scattering cross-sectionnormalisation.
Preprint: CP3-Origins-2016-036 DNRF90
I. INTRODUCTION
The detection of Milky Way dark matter (DM) parti-cles is one of the most pressing research questions in As-troparticle Physics. The experimental technique knownas direct detection will play a crucial role in this con-text in the coming years [1]. It searches for nuclear recoilevents induced by the non-relativistic scattering of MilkyWay DM particles in low-background detectors [2]. Thegoal is to disentangle the expected DM signal, i.e. afew nuclear recoil events per ton per year, from back-ground events induced by environmental radioactivity,muon-induced neutrons or solar and atmospheric neutri-nos [3]. In order to achieve this goal, different experimen-tal read-out strategies are currently under investigation,including the detection of scintillation light, athermalphonons, ionisation charge, and bubble nucleation [4]. Analternative to background discrimination is the detec-tion of an annual modulation in the observed rate of nu-clear recoil events, which would allow to identify the DMorigin of the observed signal unambiguously [5–7]. Thefirst ton-scale detectors for DM direct detection exploit-ing liquid Xenon or Argon are currently in a construc-tion or commissioning stage [8]. The first data release ofXENON1T is for instance expected in 2017, with greatexpectations for groundbreaking discoveries [9]. At thesame time, detectors with directional sensitivity, i.e. de-signed to measure anisotropies in the distribution of nu-clear recoil events, are currently in a research and de-velopment stage, and some first encouraging results havealready been achieved [10].Low-threshold detectors are a priority in the designof DM direct detection experiments. A first motivationfor low-threshold detectors arises from models of lightDM [11]. A DM particle of mass m χ moving at a speed ∗ [email protected] † [email protected] of 10 − in natural units can deposit at most an energy2 × − m χ m N / ( m χ + m N ) in the scattering by nuclei ofmass m N . Therefore, it is required a threshold energy ofabout 1 keV (1 eV) to detect a 1 GeV (1 MeV) DM par-ticle in DM-nucleus elastic collisions. This can be some-what improved by looking at inelastic channels [12]. Cur-rently none of the operating direct detection experimentshas reached threshold energies of 1 eV yet. However, var-ious strategies are under consideration, ranging from theinitial proposal of Drukier and Stodolsky for the detec-tion of neutrinos via neutral-current interactions [13] tomore recent studies by where DM detection is achievedvia excitations in superfluid helium [14] or semiconduc-tors [15].We have recently argued that low-threshold direct de-tection experiments are crucial for a second importantreason [16]. They would allow for the detection of an as ofyet neglected population of DM particles gravitationallybound to the Earth, for which we have calculated the ex-pected flux and induced event rate at detector. This newpopulation of DM particles would manifest in a directdetection experiment as a prominent spectral feature inthe low-energy part of the observed nuclear recoil energyspectrum. Such a population of bound DM particles canform if DM interacts with the nuclei in the Earth andscatters to orbits gravitationally bound to the planet,where it accumulates over the whole history of the solarsystem until the present time, when it is eventually de-tected. The velocity distribution of this new population ofDM particles peaks just below the Earth’s escape veloc-ity, and the induced nuclear recoil spectrum at detectoris maximum for values of the DM particle mass close tothe mass of abundant elements in the Earth, since in thismass range the probability of scattering to bound orbitsis larger.The literature on the capture of DM particles in orbitsbound the solar system is considerable. Most of thesestudies focus on the capture of DM particles by the Sun,and on the subsequent accumulation and annihilation of a r X i v : . [ a s t r o - ph . C O ] A ug such particles at the Sun’s centre, resulting in energeticneutrinos observable on Earth, e.g. [17–20]. The directdetection of DM particles from orbits bound to the Sunis studied in [21, 22]. It is found that the expected rate ofnuclear recoils is small due to the large Earth to Sun dis-tance. The capture of DM particles in orbits bound to theEarth is investigated in, e.g. [23–26]. Most of the workson this topic focus on the neutrino signal produced byDM annihilation at the Earth’s centre. To the best of ourknowledge, the direct detection of DM particles bound tothe Earth is addressed in two articles only, besides ourrecent publication [16]. In the pioneering work by Gouldet al. [27], the direct detection of DM particles bound tothe Earth is studied assuming a modified isothermal ve-locity distribution for DM. This study carefully accountsfor various effects related to the Sun’s gravitational po-tential, but focuses on standard spin-independent darkmatter-nucleon interactions only. In a subsequent pub-lication [28], an explicit expression for the velocity dis-tribution at the Earth’s surface of DM particles in or-bits bound to the planet is found. Our work [16] extendsthese first investigations by considering a broader set ofdark matter-nucleon interactions, a refined chemical com-position for the Earth, and detectors with and withoutdirectional sensitivity. In the present study, we furtherextend the results presented in [16] by providing signifi-cantly more general expressions for fluxes and rates nowvalid for arbitrary dark matter-nucleon interactions, andconsidering different target materials for the assumed ter-restrial detectors.This paper is organised as follows. In Sec. II we will re-view and significantly extend the calculations presentedin Ref. [16], providing all details needed to compute theflux of DM particles bound to the Earth potentially ob-servable in a terrestrial detector. In Sec. III we will con-vert this flux into a rate of nuclear recoil events, consid-ering both non-directional and directional detectors, andexpressing all equations in terms of general DM-nucleusscattering cross-sections. In Sec. IV we will numericallyevaluate the main equations previously derived and dis-cuss how our conclusions depend on assumptions regard-ing the direct detection of DM particles bound to theEarth. Finally, we will conclude in Sec. V. II. DM CAPTURE BY THE EARTH
The capture of DM particles by stellar objects and theEarth has been studied extensively in the past [17, 25]. Inparticular, the capture of DM in the Sun and its sub-sequent distribution in bound elliptical orbits has beenstudied both analytically [21, 22] and numerically [19,20]. Here we focus on DM capture by the Earth. Thekey point for the DM capture is that the particle shouldscatter underground to velocities that are below the es-cape velocity of that particular point of the Earth, thusleading to a gravitational bound orbit.Let us review the capture rate of halo DM particles to gravitationally bound orbits in the Earth after a scat-tering with a nucleus inside the Earth starting from firstprinciples. Let us assume that the DM particle densityinside the Earth at the scattering point just before thescattering takes place is dn χ = f ( (cid:126)x, (cid:126)v ) d xd v, (1)where f ( (cid:126)x, (cid:126)v ) is the DM distribution right before the col-lision (at position (cid:126)x with velocity (cid:126)v ). The number ofDM scatterings per time per center of mass solid angle d Ω that takes place within an infinitesimal volume d x inside the Earth with nuclei of the element A of density n A ( x ) is given by d ˙ N A = d xn A ( (cid:126)x ) d vf ( (cid:126)x, (cid:126)v ) v dσ A dE R dE R , (2)where dσ A /dE R is the differential cross section per recoilenergy E R . Not all scatterings lead to capture. The cap-ture condition for a scattering is for the particle to loseenergy larger than the kinetic energy it had asymptot-ically far away from the Earth. The energy before thecollision (i.e. kinetic plus potential one) is E before = 12 m χ ( v − v ( r )) = 12 m χ v ∞ , (3)where v esc ( r ) is the escape velocity from the Earth at aradius r from the center of the Earth (i.e. at the place ofthe scattering) and v ∞ is the velocity of the particle atan asymptotically far away distance from the Earth. Thetotal energy after the collision must be negative in orderfor the DM particle to remain in a bound orbit aroundthe Earth. Its value is E after = 12 m χ ( v (cid:48) − v ( r )) = − Gm χ M ⊕ a ≡ − m χ α, (4)where a is the major semi-axis of the elliptical orbit afterthe collision and α is defined as α ≡ GM ⊕ /a ( G being thegravitational constant and M ⊕ the mass of the Earth).Using Eqs. (3) and (4) we get the energy transfer E R E R = 12 m χ ( v − v ( r ) + α ) = 12 m χ ( v ∞ + α ) . (5)Eq. (5) gives dE R = (1 / m χ dα and Eq. (2) now reads d ˙ N A = 12 d xn A ( (cid:126)x ) d vf ( (cid:126)x, (cid:126)v ) v dσ A dE R m χ dα Θ α , (6)where Θ α represents a step function that enforcesthe kinematic constraint E R ≤ β A + E kb where E kb =(1 / m χ v is the kinetic energy before the collision. Wedefine β A ± = 4 m χ m A ( m χ ± m A ) . (7)Using Eq. (5) the above condition can be written as2 E R /m χ = v ∞ + α ≤ β A + v = β A + ( v ∞ + v ( r )). Since1 /β A + − /β A − = 1 the above constraint can be rewrittenas Θ α ≡ Θ (cid:20) β A − (cid:18) v ( r ) − αβ A + (cid:19) − v ∞ (cid:21) , (8)where it is understood that the step function Θ( x ) = 1if x ≥ f ∞ ( (cid:126)v ∞ ) = n χ π / v exp (cid:18) − ( (cid:126)v ∞ + (cid:126)v e ) v (cid:19) , (9)where v = 220 km s − is the local standard of rest, v e = 232 km s − is the Earth velocity in the galac-tic rest frame, and n χ the DM number density in theEarth’s neighborhood. Liouville’s theorem states thatthe distribution function remains constant along the tra-jectory of a particle, i.e. f ( (cid:126)x, (cid:126)v ) = f ∞ [ v ∞ ( (cid:126)x, (cid:126)v )] where f ∞ is the DM distribution far away from the Earth and v ∞ = v − v ( r ). Taking the angular average of f ( (cid:126)x, (cid:126)v )defined as (cid:82) f ( (cid:126)x, (cid:126)v ) d v = 4 π (cid:82) v ¯ f ( r, v ) dv we get¯ f ( v ) dv = n χ π / v E v (cid:112) v − v (cid:32) e − v − v − e − v v (cid:33) dv , (10)where v ± = (cid:112) v − v ± v e . Note that we have droppedthe variable r from ¯ f . The escape velocity of theEarth varies from 15 km s − at the Earth’s centre to11.2 km s − at the Earth’s surface. Since the variationis small, we simplified our calculation, by setting the es-cape velocity to its surface value v =11.2 km s − . Thismakes ¯ f ( r, v ) independent of r (leading to Eq. (10)).Upon making the isotropic approximation, we can sim-plify further Eq. (6). The specific angular momentum ofthe particle after the collision is J = rv sin θ where r isthe distance from the center of the Earth, v the velocityafter the collision and θ the angle subtended by (cid:126)r and (cid:126)v . Since we assume that cos θ is uniformly distributed,and J = J (1 − cos θ ), the distribution rewrittenin terms of J is d cos θ = dJ / (2 J (cid:112) − J /J )where J max = rv is the maximum possible specific angu-lar momentum after the collision. Within this approxi-mation we can now rewrite Eq. (6) as d ˙ N A = πd xn A ( r ) v dv ¯ f ( v ) dσ A dE R m χ (cid:32) J (cid:115) − J J (cid:33) − (11) × dαdJ Θ α Θ J , where Θ J = Θ( J max − J ) is a step function enforcing J ≤ J max . With the use of Eq. (4), J max = r ( v ( r ) − α ) / .One can easily check that in the case of spin-independentinteractions where dσ A dE R = m A σ A µ A v F A ( E R ) , (12)where µ A is the DM-nucleus reduced mass, Eq. (11) be-comes the one derived in [22] d ˙ N A = 2 πσ A v ¯ f ( r, v ) n A ( r ) J β A + (cid:18) − J J (cid:19) − / F A ( E R ) × Θ α Θ J (cid:0) d x dv (cid:1) dα dJ . (13)The form factor F A ( E R ) accounts for the loss of coher-ence and it is usually approximated by F A ( E R ) = exp( − E R /Q A ) , (14)where E R is the energy transferred during the collisionand Q A = 3 / (2 m A R A ), m A being the nucleus mass and R A = 10 − cm (cid:104) . . (cid:0) m A GeV (cid:1) / (cid:105) the radius of thenucleus. In this paper since we will present results fordifferent types of DM-nuclei interactions, we will useEq. (11) which can be used for any generic interactionand form factor.Eq. (11) can be written in a more convenient form interms of new more useful variables for the purposes ofthis study. Instead of using J and α , we will use theperihelion (minimum distance of the elliptical orbit tothe center of the Earth) r m and the ellipticity of the orbit e . Recall that the semi-major axis for an ellipse is a = r m / (1 − e ) and consequently α = GM ⊕ (1 − e ) /r m . Notealso that J = r m ( v − α ). From these two expressionswe can calculate the Jacobian and get dJ dα = 2 GM ⊕ (cid:18) v − GM ⊕ (1 − e ) r m (cid:19) dedr m . (15)Eq. (11) can be written in terms of the new variables r m and e as d ˙ N A = 2 πGM ⊕ d xn A ( r ) v dv ¯ f ( v ) dσ A dE R m χ (16) × (cid:32) r (cid:114) − r m r (cid:33) − Θ r m Θ e dr m de. The condition J = r m ( v − α ) / ≤ J max imposed by Θ J becomes Θ r m ≡ Θ( r − r m ) and Θ e is Θ α having subsi-tuted α = GM ⊕ (1 − e ) /r m . Recall that the semi-majoraxis a = r m / (1 − e ). For the typical spin-independentDM-nucleus cross section of Eq. (12), Eq. (16) takes theform provided in [16] d ˙ N A = 4 πGM ⊕ σ A vf ( v ) n A ( r ) r β A + (cid:18) − r m r (cid:19) − / F A ( E R ) × Θ r m Θ e (cid:0) d x dv (cid:1) de dr m . (17)Since we consider generic DM-nuclei interactions, we aregoing to use the more generic form of Eq. (16).Eq. (16) should be summed over all elements abun-dant in the Earth. In practice we take into account themost abundant elements, i.e. O, Si, Mg, Fe, Ca, Na, S, Ni, and Al assuming the standard com-position and density profile of chemical elements in theEarth n A ( r ) provided in [35]. Integrating Eq. (16) over d x dv and summing over elements gives d ˙ N = 8 π GM ⊕ m χ (cid:88) A K A ( r m , e ) × (cid:90) R ⊕ r m dr n A ( r ) (cid:18) − r m r (cid:19) − / de d r m ≡ g ( r m , e ) de dr m . (18)Eq. (18) gives the rate of accumulation of trapped DMparticles into bound elliptical orbits of ellipticity within[ e, e + d e ], and perihelion within [ r m , r m + d r m ]. In thederivation of Eq. (18), we have assumed spherical sym-metry, i.e. d x = 4 πr d r . K A ( r m , e ) is defined as K A ( r m , e ) ≡ (cid:90) v v dv v ¯ f ( v ) dσ A dE R . (19)The upper limit v comes from the step function Θ e andit given by v = (cid:115) (1 + β A − ) v − GM ⊕ r m (1 − e ) β A − β A + . (20)The lower limit of intergration is obviously the escape ve-locity v since a DM particle with zero speed at asymp-totic far distances from the Earth, will acquire v onceit reaches the Earth. dσ A /dE R depends generally on E R (either explicitly or via the form factor F A ( E R ). In sucha case E R = (1 / m χ (cid:18) v − v + GM ⊕ (1 − e ) r m (cid:19) (21)is the energy loss in the collision that must be used inthe evaluation of K A ( r m , e ). III. RECOIL ENERGY SPECTRUM OF BOUNDDARK MATTER
In order to estimate the rate of events of bound DMparticles scattering off a detector, we need to estimate theprobability of DM particles that follow a specific ellipticorbit to scatter off the detector as well as the numberof bound DM particles per specific elliptical orbit. Tosimplify our estimate, we are going to consider DM par-ticles that have scattered in the Earth once in order toget captured and a second time in the detector creatinga recoil signal. Multiple scatterings that take place un-derground diminish further the kinetic energy of the DM particle leading to recoil energies that are practically be-low any experimental threshold. Therefore within thisapproximation, we estimate the number of DM particlesthat can accumulate in different orbits and have scat-tered only once. We can now estimate the number ofperiods N required for a bound DM particle to scatterfor a second time N = (cid:32)(cid:88) A (cid:90) θ n A ( r ) σ A ξ ( r m , e ) dθ (cid:33) − , (22)where ξ ( r m , e ) dθ is an infinitesimal path along the el-liptic trajectory of the orbit. The length of the path thata DM particle travels underground is (cid:90) d(cid:96) = 2 (cid:90) θ dθ (cid:115)(cid:18) drdθ (cid:19) + r ≡ (cid:90) θ ξ ( r m , e ) dθ. (23)Using the parametric equation for the elliptic orbit Pr = 1 + e cos θ, (24)where P is a constant, e the ellipticity of the orbit and θ the angle subtented from a point of the orbit with dis-tance r from the center and the perihelion, it is easilyfound that ξ ( r m , e ) = 2 r m (1 + e ) (cid:112) e + 2 e cos θ/ (1 + e cos θ ) . (25)The limit of integration θ is given bycos θ = r m R ⊕ (1 + e ) e − e (26)and corresponds to the angle subtended by the perihelionand the point where the orbit crosses the Earth ( r = R ⊕ ) from the Earth’s center. It can be found by setting r = R ⊕ and solve for θ in Eq. (24) The condition − < cos θ < − e e ≤ r m R ⊕ ≤ . (27)For a given orbit, the time T ( r m , e ) a DM particle canspend without scattering for a second time until today ison average T ( r m , e ) ≡ min[ N × τ ( r m , e ) , τ ⊕ ] , (28)where τ ⊕ (cid:39) . × years is the age of the Earth and τ ( r m , e ) = (cid:115) π GM ⊕ r m (1 − e ) (29)is the period of the elliptical bound orbit. We will referto T as accumulation time. A. Non-Directional Detectors
The differential event rate in a non-directional detectorfor a given orbit characterized by r m and e is dR r m ,e dE R = N T dσ N dE R F = N T dσ N dE R d ˙ N πl c T ( r m , e ) τ ( r m , e ) , (30)where N T is the number of target nuclei in the detec-tor. F is the flux of bound DM particles in orbits ofperihelion r m and ellipticity e crossing the detector. Theflux is equal to the rate d ˙ N with which a particular or-bit is populated (see Eq. (18)) multiplied by the time T ( r m , e ) this orbit can accumulate DM particles dividedby τ ( r m , e ) / π(cid:96) c ( (cid:96) c being the distance between the detector and the centerof the Earth). We have assumed that the elliptical orbitscross the surface of the Earth isotropically, i.e. there areno bound DM particles crossing a particular patch of theEarth’s surface with a higher rate than another patch.This gives the factor 4 π(cid:96) c . Since generically dσ N /dE R depends on the DM particle velocity, it is needed to knowthe velocity before the scattering with the detector. Itis completely determined by r m and e and can be easilyshown to be v = (cid:115) GM ⊕ (cid:18) r − − e r m (cid:19) . (31)with r = (cid:96) c . Note that dσ N /dE R refers to DM scatteringoff a detector nucleus and it should not be confused with dσ A /dE R that was the scattering that lead to the captureof DM by a random underground nucleus.Combining Eqs. (18), (28) and (30) we obtain the dif-ferential rate of events dRdE R = N T π(cid:96) c (cid:90) (cid:90) R ⊕ − e e R ⊕ dedr m g ( r m , e ) dσ N dE R T ( r m , e ) τ ( r m , e ) dr m de. (32)We stress again that in general dσ N /dE R depends on v ,and v should be evaluated at the value given by Eq. (31).Eq. (32) represents the main equation that gives the eventrate in non-directional detectors. If one assumes spin-independent interactions (Eq. (12)), the spectrum recoilbecomes dRdE R = κ (cid:90) (cid:90) R ⊕ − e e R ⊕ g ( r m , e ) v T ( r m , e ) τ ( r m , e ) dr m de, (33) where κ = N T m N σ n A N F ( E R ) / (4 π(cid:96) c µ N ).Eq. (32) must be contrasted to the recoil events comingfrom direct halo DM scatterings off nuclei targets in thedetectors. The rate is as usually given by dRdE R = N T n χ (cid:90) v esc + v e v min dσ N dE R f ( v ) vd v, (34)where n χ is the local DM density in the Earth, and v min = (cid:113) m N E R / (2 µ N ) (35)is the minimum velocity that can produce nuclear recoilof energy E R . For f ( v ) we use the usual Maxwell-Boltzmann of Eq. (9) with v esc and v e being the escapevelocity of the Galaxy and the velocity of the Earth inthe rest frame of the Galaxy respectively. B. Directional Detectors
We also study the spectrum of bound DM scatteringoff directional detectors. By choosing an appropriate re-coil direction, directional detectors have the advantageof minimizing the rate of events coming from the haloDM particles. Pointing the cone of detection along withthe DM wind, one looks at particles that have velocities (cid:126)v − (cid:126)v e . This leads to overall smaller particle fluxes andconsequently to smaller rate of events. On the contrarythis choice does not affect the rate of events of boundDM particles. In particular we will consider the spec-trum of recoils coming from a direction perpendicular tothe vector that connects the center of the Earth with thedetector. We have found that such horizontal directionscan give an enhancement in the bound/halo ratio of DMevents in the detector. Generically the directional ratefor energy recoil E R and recoil direction within the solidangle d Ω q is dRdE R d Ω q = N T (cid:90) dσdE R d Ω q d Φ , (36)where d Φ is the flux of particles arriving at the detector.For a generic DM-nucleus interaction, the cross sectionper nuclear recoil energy per recoil solid angle is dσ N dE R d Ω q = dσ N dE R π δ (cid:16) cos θ q − v min v (cid:17) , (37)where θ q is the angle between the nuclear recoil and theinitial DM velocity and v min is given by Eq. (35). Eq. (36)can be rewritten with the help of (37) as dRdE R d Ω q = N T πδ(cid:96) c (cid:90) dσ N dE R g ( r m , e ) T ( r m , e ) τ ( r m , e ) δ (cid:16) cos θ q − v min v (cid:17) dr m de d cos θdφ π dω π . (38)Eq. (38) requires some explanation. The flux of boundDM particles is proportional to g ( r m , e ) T ( r m , e ) τ ( r m , e )as in the case of non-directional detectors divided by theeffective area of the detector δ(cid:96) c . Eventually we will showthat the result will be independent of δ(cid:96) c . In the case ofnon-directional detectors we were interested in the to-tal flux of particles passing through the detector withoutcaring about the direction. Therefore once we knew thedensity of bound particles per orbit, we had to integrateover all possible orbits (i.e. r m and e ) in order to esti-mate the total rate. In the case of directional detection,not only do we care about the total number of eventsper time, but we need to know from what direction DMparticles come from. Since we care about detecting par-ticles that scatter off nuclei in the detector creating anuclear recoil to a particular direction, r m and e are notthe only variables we need to achieve that. In addition tothe characteristics of the elliptical orbit, we need to knowwhat is the location of the perihelion of the orbit com-pared to the detector location.Therefore we parametrizethe orbits by r m , e , the polar angles θ and φ that definethe location of the perihelion with respect to the detec-tor (i.e. the detector is along the z -axis) and the angle ω between the plane of the orbit and the plane definedby the perihelion the center of the Earth and the detec-tor. We expect an isotropic distribution of the perihelionaround the Earth and a uniform distribution for ω . Thisis why we divide the corresponding quantites by 4 π and2 π respectively in Eq. (38). The δ function enforces therecoil angle θ q to be the one that kinematics dictates. Wenow need to find the orbits that pass from the detector’slocation and can create a nuclear recoil to a particular horizontal direction. Eq. (24) evaluated at θ = 0 givesthe perihelion r = r m . Therefore trading P for r m andusing r = (cid:96) c (the distance of the detector from the centerof the Earth) we rewrite Eq. (24) as r m = (cid:96) c e cos θ e . (39)For a given orbit where the perihelion forms an angle θ with the center of the Earth and the detector, r m mustbe given by the above equation in order for the particleto pass from the detector’s location. Varying the valueof the perihelion while keeping e and θ fixed leads to δr m = δ(cid:96) c e cos θ e . (40)The integration over dr m can be substituted approxi-mately by δr m which is related to the size of the detector.On the other hand in order for the orbit to pass throughthe detector (of dimension δ(cid:96) c ), (cid:96) c sin θδω = δ(cid:96) c ⇒ δω = δ(cid:96) c (cid:96) c sin θ . (41)Since θ takes values from 0 to π , it is always positive.Since δ(cid:96) c << (cid:96) c δω is extremely small unless one consid-ers very small values of θ (practically locating the peri-helion inside the detector). If we ignore this tiny patchof surface for the perihelion, we can substitute the inte-gration over dω by δω . Using Eqs (40) and (41) we canwrite (38) as dRdE R d Ω q = N T π (cid:96) c (cid:90) dσ N dE R g ( r m , e ) T ( r m , e ) τ ( r m , e ) δ (cid:16) cos θ q − v min v (cid:17) ez e √ − z dzdφde, (42)where r m is given by Eq. (39). We defined z ≡ cos θ .Note that the rate does not depend anymore on the char-acteristic size of the detector δ(cid:96) c . We will eventually usethe delta function to perform the integral over z . Beforewe do this, we need to find the relation of θ q with thevariables of the problem i.e. e , φ and z . Let us considerfor the moment an orbit with φ = 0 and an angle θ sub-tended by the detector, the center of the Earth and theperihelion of the orbit. If we use cartesian coordinateswith the perihelion being along the x -axis, a point in theorbit has coordinates x = a e + r cos θ, y = r sin θ (43)with a being the focal point. Let us choose a horizontaldirection at the location of the detectorˆ θ = − sin θ ˆ x + cos θ ˆ y. (44)A bound DM particle that follows a particular ellipti-cal orbit reaches the detector with a velocity that has a direction ˆ (cid:96) = dxd(cid:96) ˆ x + dyd(cid:96) ˆ y. (45)With the help of Eq. (43) dxd(cid:96) = dr cos θ − r sin θdθ (cid:112) dx + dy = dθ (cid:0) drdθ cos θ − r sin θ (cid:1) dθ (cid:113)(cid:0) drdθ (cid:1) + r . (46)Canceling the dθ from numerator and denominator andcalculating dr/dθ from Eq. (24) we get the final result dxd(cid:96) = − sin θ √ e + 2 e cos θ . (47)Similarly dyd(cid:96) = e + cos θ √ e + 2 e cos θ . (48)Using Eqs. (44), (45), (47) and (48) we getcos θ q = ˆ θ · ˆ (cid:96) = ± e cos θ √ e + 2 e cos θ . (49)Recall that ˆ θ is the recoil direction and ˆ (cid:96) the directionof the velocity of the bound DM particle. The ± refersto the two possibilities that the particle is orbiting theellipse (counter)clockwise. It is not difficult to show thatfor a nonzero value of φ cos θ q = ˆ θ · ˆ (cid:96) = ± e cos θ √ e + 2 e cos θ cos φ. (50)We show the details of the derivation in the case ofnonzero φ in the appendix. We assume that there isequally probable to have clockwise or counterclockwiseorbits. Let us consider first the orbits with a plus signin Eq. (50). We will multiply the corresponding rate bya factor of 1/2 since there is 50% probability. In orderto evaluate the dz integration using the delta function inEq. (42), we will use the well known property δ [ h ( z )] = δ ( z − z ) | h (cid:48) ( z ) | , (51)where h ( z ) is a function of z , z is the solution of theequation h ( z ) = 0 and h (cid:48) ( z ) is the derivatize of h ( z )with respect to z evaluated at z . In our particular case h ( z ) = cos θ q − v min v = (1 + ez ) cos φ √ e + 2 ez − v min (cid:113) GM ⊕ (cid:96) c (cid:113) − − e ez ) , (52)where v is given by Eq. (31) with r m given by Eq. (39).Recall that z = cos θ . The equation h ( z ) = 0 has the solution z = − cos φ + γe cos φ , (53)where γ = v (cid:96) c GM ⊕ . (54)It is also easy to show that | h (cid:48) ( z ) | = e cos φ (cid:112) γ − (1 − e ) cos φ . (55)The constraint − < z < (cid:114) γ e < cos φ < (cid:114) γ − e . (56)From Eq. (52) it is clear that 0 < cos φ < (cid:112) γ/ (1 + e ) <
1. This last condition can be rewritten as e > γ − . (57)Recall that 0 < e < γ − < ⇒ γ <
2. Using the definition of γ (Eq. (54)) and v min from Eq. (35), the constraint γ < E R < µ N GM ⊕ m N (cid:96) c . (58)This condition in fact sets the upper limit in the recoil en-ergy spectrum that bound DM particles can contribute.We can now rewrite Eq. (42) performing the integra-tion over z by using the delta function as we prescribedabove dRdE R d Ω q = N T π (cid:96) c (cid:90) e (cid:90) φ φ dσ N dE R g ( r m , e ) T ( r m , e ) τ ( r m , e ) 1 + ez e (cid:112) − z (cid:112) γ − (1 − e ) cos φ cos φ Θ(2 − γ ) dφde. (59) e = Max[ γ − ,
0] is derived from the constraint ofEq. (57) and the fact that e >
0. The step functionΘ(2 − γ ) ensures that γ < φφ = cos − (cid:20) Min (cid:18) , (cid:114) γ − e (cid:19)(cid:21) φ = cos − (cid:114) γ e . (60)Note that r m is evaluated at the value r m = (cid:96) c ez e , (61) (see Eq. (39)). Eq. (59) is our final result for the recoilspectrum in directional detectors. Comparing the overallcoefficient of Eq. (59) with respect to that of Eq. (42),the former is larger by a factor of 4 (there is a factor of1 / /
16 respectively). A factor of 2 comes fromthe integration of φ . Note that Eq. (56) is satisfied in tworegions i.e. one with positive and one with negative valueof φ . Since cos φ always appears as a square, we integrateonly over positive φ and multiply by 2. The second factorof 2 comes from the fact that the orbits with the oppositedirection (i.e. with a minus sign in Eq. (50)) give exactlythe same contribution as the orbits with the plus sign.This is easy to show: The solution of Eq. (52) is still m χ [GeV]20 40 60 80 100 K A [ ( s / k m ) c m - ] -39 -38 -37 -36 -35 -34 -33 Oxygen / e=0.6Oxygen / e=0.95Iron / e=0.6Iron / e=0.95 m χ [GeV]20 40 60 80 100 K A [ ( s / k m ) c m - ] -44 -43 -42 -41 -40 -39 -38 -37 Oxygen / e=0.6Oxygen / e=0.95Iron / e=0.6Iron / e=0.95
FIG. 1. K A as a function of the DM particle mass m χ for two elements in the Earth, namely Oxygen and Iron. K A isproportional to the probability of scattering towards a bound orbit of given ellipticity e and perihelion r m . In the figure wevary e as reported in the legends, and fix r m to R ⊕ /
2. The left panel refers to the interaction O with c = 2 /m V and c = 0,whereas the right panel to the interaction O with c = 2 /m V and c = 0. The parameter m V = 246 . given by (53) even for the orbits with a minus sign in(50). The only difference is that in this case cos φ < φ →− cos φ . However since cos φ appears always as cos φ inEq. (59), one can change variable φ (cid:48) ≡ π − φ keeping inmind that cos φ (cid:48) = cos φ . The constraint on φ (cid:48) is thesame of Eq. (56) with φ → φ (cid:48) since cos φ (cid:48) = − cos φ . Thevalue of | h (cid:48) ( z ) | is the same as before and therefore theoverall contribution of the “negative sign” orbits is thesame as the ones with positive sign.Eq. (59) can be used for any generic form of DM-nucleon interactions. For the spin-independent interac-tion of Eq. (12), (59) becomes dRdE R d Ω q = κ d (cid:90) e de (cid:90) φ b φ a dφ v g ( r m , e ) τ ( r m , e ) T ( r m , e ) × ez e (cid:112) − z (cid:112) γ − (1 − e ) cos φe cos φ , (62)where κ d = N T m N σ n A N F ( E R ) / (8 π µ N (cid:96) c ).Eq. (59) describes the recoil spectrum of bound DMscattering off an underground detector. This spectrummust be contrasted to the usual directional spectrum ofhalo DM. Using Eq. (37) we get dRdE R d Ω q = N T n χ π (cid:90) dσ N dE R δ (cid:16) ˆ v · ˆ q − v min v (cid:17) f ( v ) vd v. (63) In the case of spin indeppendent interactions (seeEq. (12)), it takes the form dRdE R d Ω q = κ h ˆ f ( v min , ˆ q ) , (64)where κ h = N T n χ m N σ n A N F N ( E R ) / (4 πµ N ) andˆ f ( v min , ˆ q ) is the so-called Radon transfromation of f ( v )defined as [29]ˆ f ( v min , ˆ q ) = (cid:90) δ ( (cid:126)v · ˆ q − v min ) f ( v ) d v. (65) IV. RESULTS
The main equations derived in the previous sectionsare Eqs. (32) and (59). They describe the rate of nuclearrecoil events expected in non-directional and directionaldetectors, respectively. Now we numerically evaluate andinterpret these expressions under different assumptionsregarding the cross-sections dσ A /dE R (for scattering inthe Earth) and dσ N /dE R (for scattering in a terrestrialdetector). We will also investigate the dependence ofour results on the type of target nuclei composing thedetector in analysis. r m [m]10 T [ y ea r s ] -8 -7 -6 -5 -4 -3 -2 -1 m χ =15 GeV / e=0.6m χ =15 GeV / e=0.95m χ =50 GeV / e=0.6m χ =50 GeV / e=0.95 r m [m]10 T [ y ea r s ] -5 -4 -3 -2 -1 m χ =15 GeV / e=0.6m χ =15 GeV / e=0.95m χ =50 GeV / e=0.6m χ =50 GeV / e=0.95 FIG. 2. Accumulation time T as a function of the perihelion r m for two reference values of the DM particle mass m χ and of theellipticity e . The left panel refers to the interaction O with c = 2 /m V and c = 0, whereas the right panel to the interaction O with c = 2 /m V and c = 0. Coupling constants are expressed in terms of the electroweak scale m V = 246 . O is significantly larger than that of O . Overall, T × K A for O is larger than thecorresponding of O . A. General considerations
The rate of nuclear recoil events in Eq. (32) dependson the cross-section dσ A /dE R through the functions K A and T (defined in Eqs. (19) and (28), respectively). Italso depends on the differential cross-section dσ N /dE R ,which appears in Eq. (32) directly. We can therefore char-acterise each single scattering event at detector as the re-sult of a complex three stage physical process. Each stageexplicitly depends on how DM interacts with nuclei andis briefly described below:1. Capture of the DM particle χ by the Earth. The el-ement A contributes with probability proportionalto K A .2. Motion of the DM particle χ along the bound orbitcharacterised by r m and e . This motion lasts onaverage for a time T , i.e. the accumulation timedefined in Eq. (28).3. Scattering of the particle χ at detector (with cross-section given by dσ N /dE R ).In all numerical applications, we will assume the cross-section dσ A dE R = 2 m A (2 j A + 1) v (cid:88) τ =0 , (cid:88) τ (cid:48) =0 , (cid:104) c τ c τ (cid:48) W ττ (cid:48) M ( E R )+ 2 j χ ( j χ + 1)3 m A E R m n c τ c τ (cid:48) W ττ (cid:48) Φ (cid:48)(cid:48) ( E R ) (cid:105) , (66) and an analogous expression for dσ N /dE R . The isotope-dependent nuclear response functions W ττ (cid:48) M and W ττ (cid:48) Φ (cid:48)(cid:48) inEq. (66) are quadratic in nuclear matrix elements andare defined in Ref. [30]. They have been calculated forthe 16 most abundant elements in the Sun, including O, Si, Mg, Fe, Ca, Na, S, Ni, and Al,in Ref. [31] and for various isotopes of Xe and Ge, andfor Na in Ref. [32]. The labels 1 and 11 in Eq. (66)refer to the non-relativistic effective operators O and O introduced in Ref. [30]. The former corresponds tothe familiar spin-independent interaction operator, thelatter to the momentum-dependent interaction operator O = ( q /m n ) · S χ , where m n is the nucleon mass, and q and S χ are the momentum transfer and DM particlespin operators, respectively. They are explicitly definedin Ref. [31]. A comparison of Eqs. (66) and (12) allows toexpress σ A and F A in terms of the coupling constants andresponse functions in Eq. (66). For the isoscalar couplingconstants, c and c , we assume the reference values2 /m V , with m V = 246 . O or O . At the same time, we set the isovectorcoupling constants to zero: c = c = 0. Finally, j A and j χ are the A element and DM particle spins, respectively.Knowledge of the Earth’s chemical composition isneeded in order to evaluate K A and T . In this study, weconsider the nine elements: O, Si, Mg, Fe, Ca, Na, S, Ni, and Al, with mass fractions as givenin Ref. [33], and the radial density given in Ref. [34] andimplemented in Ref. [35]. We have verified numerically0 D r u k i e r D A M I C CD M S li t e E R [keV]10 -3 -2 -1 R a t e [ k g - da y - k e V - ] m χ =15 GeV / halom χ =15 GeV / totalm χ =25 GeV / halom χ =25 GeV / totalm χ =50 GeV / halom χ =50 GeV / total D r u k i e r D A M I C CD M S li t e E R [keV]10 -3 -2 -1 R a t e [ k g - da y - k e V - ] Ge / haloGe / totalNa / haloNa / totalXe / haloXe / total
FIG. 3. Rate of nuclear recoil events dR/dE R as a function of E R . We assume dark matter-nucleon interactions of type O and c = 2 /m V , c = c = c = 0 ( m V = 246 . m χ = 50 GeV, and considerdifferent target materials for dσ N /dE R , namely, Xenon, Germanium and Sodium. In both panels, solid lines correspond to thetotal rates, including the contribution from halo and bound DM particles. Dashed lines represent the contribution to dR/dE R from halo DM particles. Vertical lines show illustrative energy thresholds of running or proposed dark matter direct detectionexperiments. D r u k i e r D A M I C CD M S li t e E R [keV]10 -3 -2 -1 R a t e [ k g - da y - k e V - ] -4 -3 -2 -1 m χ =15 GeV / halom χ =15 GeV / totalm χ =25 GeV / halom χ =25 GeV / totalm χ =50 GeV / halom χ =50 GeV / total D r u k i e r D A M I C CD M S li t e E R [keV]10 -3 -2 -1 R a t e [ k g - da y - k e V - ] -4 -3 -2 -1 Ge / haloGe / totalNa / haloNa / totalXe / haloXe / total
FIG. 4. Same as for Fig. 3 but now for the interaction O . that changes in the mass fraction of single elements inthe Earth have a negligible impact on the scattering rateevaluation.Fig. 1 shows K A as a function of the DM particle mass m χ for two elements in the Earth, namely Oxygen and Iron, and for two reference values of the ellipticity e . Wefind that K A increases for m χ → m A since in this massrange the upper limit v in Eq. (19) tends to infinity, i.e.maximum momentum transfer in the scattering. Fig. 1also shows that for large values of e the range of masses1 m χ [GeV]10 R a t e r a t i o -4 -3 -2 -1 d σ /dE ∝ d σ /dE ∝ d σ /dE ∝ q /v m χ [GeV]10 R a t e r a t i o -6 -5 -4 -3 -2 -1 TotalOxygenSiliconMagnesiumIronCalciumSodiumSolfurNickelAluminium
FIG. 5. Left: Ratio of Eqs. (32) and (34) as a function of m χ for three different dark matter-nucleon interactions. From thetop to the bottom in the legend: O , a modified version of O obtained by replacing c with c /v (i.e. resonant scattering),and O . In all cases we set c and c to zero, and assume Germanium as a target material. In the figure, we introduce thesymbol dσ/dE to characterise the scaling of dσ A /dE R and dσ N /dE R as a function of the dark matter-nucleus relative velocity v and of the momentum transferred q . Right: Contribution of O, Si, Mg, Fe, Ca, Na, S, Ni, and Al to theratio of Eqs. (32) and (34) as a function of the DM mass m χ . We assume O as dark matter-nucleon interaction, c = 2 /m V ( m V = 246 . c = 0. In both panels we assume a nuclear recoil energy of 1 eV in the evaluation of the scatteringrates. where K A (cid:54) = 0 is broader than for e (cid:39)
0. The reasonis that for a given r m , the upper limit v (Eq. 20) inthe integral defining K A grows with e and the integrandin Eq. (19) is proportional to v which also grows with e . Fig. 1 has been obtained by setting r m = R ⊕ / T as a function of r m for two reference val-ues of m χ and e . As expected, T grows when r m → R ⊕ and e → r m → R ⊕ and e →
1, i.e. forthe orbits contributing the most to the rate of nuclear re-coil events presented in what follows. At the same time,Fig. 2 is quantitatively reliable in the limit r m → R ⊕ and e → B. Non-directional detectors
In this section we focus on the rate of nuclear recoilevents dR/dE R in Eq. (32). Fig. 3 shows dR/dE R as a function of E R assuming dark matter-nucleon inter-actions of type O . We have obtained this figure underthe assumption c = 2 /m V , c = 0. The left panel re-ports results obtained for three different values of theDM particle mass and assuming Germanium as a targetmaterial, whereas in the right panel we consider differ-ent target materials for dσ N /dE R , namely Xenon, Ger-manium and Sodium, and fix the DM particle mass at m χ = 50 GeV. In both panels, solid lines are the totalrates, including the contribution from halo and boundDM particles. Dashed lines represent the contributionto dR/dE R from halo DM particles. The case of darkmatter-nucleon interactions of type O is discussed inFig. 4, where in the left (right) panel we have reportedresults obtained for different DM particle masses (tar-get materials). In both figures, vertical lines correspondto the threshold energies of present (CDMSlite [36] andDAMIC [37]) or proposed (Drukier [13]) direct detectionexperiments.From Figs. 3 and 4 we conclude that DM particles inorbits bound to the Earth can be revealed in future directdetection experiments as pronounced features in the low-energy part of the induced nuclear recoil spectrum. As inthe case of halo DM, DM particles bound to the Earthcan produce a larger number of nuclear recoil events atlow-energies if they are light, and in detectors composedof heavy nuclei.2In order to assess the significance of the predicted spec-tral features, we evaluate the ratio of Eqs. (32) and (34)as a function of m χ . The result of this calculation is re-ported in Fig. 5. The left panel shows the rate ratio forthree dark matter-nucleon interaction types: O , a mod-ified version of O obtained by replacing c with c /v inthe equations above (e.g. resonant scattering [38]), andfinally O .The rate ratio can be as large as 0.1 for the interaction O , and 0.4 for its resonant analogous. Notably, for theinteraction O the value can be up to ∼
200 at the Ironresonance (i.e. m χ ∼
50 GeV). The large value found forthe operator O is related to the large accumulation time T that DM particles interacting with nuclei via O canspend on bound orbits before a second scattering occurs(see right panel in Fig. 2). We have verified numericallythat the ratio of Eqs. (32) and (34) is independent of thecoupling constants c and c when a single interactionat the time is considered.Finally, the right panel in Fig. 5 shows the contributionof O, Si, Mg, Fe, Ca, Na, S, Ni, and Alto the ratio of Eqs. (32) and (34) as a function of the DMparticle mass, and assuming O as dark matter-nucleoninteraction. The overall shape of the rate ratio reflects theresonant form of the function K A , and contributions fromdistinct elements in the Earth can easily be identified inthe figure. C. Directional detectors
We conclude this section with a quantitative analysisof Eq. (59), which describes the double differential rateof nuclear recoil events induced by the scattering of DMparticles bound to the Earth in directional detection ex-periments [39–41].In Fig. 6, the left panel shows the ratio of the doubledifferential rates in Eqs. (59) and (63) as a function of theDM particle mass, assuming O as dark matter-nucleoninteraction. The solid blue line refers to a hypotheticaldetector composed of Fluorine, whereas the dashed redline corresponds to a second hypothetical detector whichuses He as a target material. The right panel in Fig. 6shows the same ratio, now evaluated for the operators O and O , and assuming Fluorine as a target material. Asfor the case of non-directional detectors, the predictedspectral feature is more pronounced for the operator O than for O by roughly three orders of magnitude.For directional detectors the size of the effect is of theorder of 0.1 (10 ) for O ( O ), and it is generically largerthan that of the non-directional detectors (e.g. compareFig. 5 with Fig. 6). Furthermore, Fig. 6 shows that thepredicted spectral feature is slightly more pronounced fora He based detector than for a detector adopting F as atarget material.
V. CONCLUSION
We have studied the properties and detection prospectsof DM particles bound to the Earth. The new DM pop-ulation forms via scattering of Milky Way DM parti-cles by nuclei in our planet’s interior. We have derivedfluxes and nuclear recoil event rates at directional andnon-directional detectors expected for the new popula-tion of DM particles. The equations presented in thiswork are valid for arbitrary dark matter-nucleon inter-actions, and extend those found in Ref. [16]. We havenumerically evaluated such expressions under differentassumptions regarding the scattering of DM in the Earthand at detector, carefully modelling the Earth internalcomposition, and considering different target materialsfor the assumed directional and non-directional DM di-rect detection experiments.We have found that future DM direct detectionexperiments with an ultra-low energy threshold ofabout 1 eV (and equally low energy resolution) havethe potential to reveal the population of DM particlesstudied in this paper with the same exposure needed todetect the associated Milky Way DM component. DMparticles bound to the Earth manifest as a prominentfeature in the low-energy part of the observed nuclearrecoil energy spectrum. In particular we have found thatDM-nucleus operators like O can give rates in recoilevents of bound DM in detectors up to a few hundredtimes higher than the corresponding Milky Way DM inlow energies. The existence and the shape of this featureare independent of the dark matter-nucleus scatteringcross-section normalisation. This work provides anadditional important motivation to invest in the designand development of a new class of ultra-low thresholdenergy detectors. Acknowledgments.
CK is partially funded by theDanish National Research Foundation, grant numberDNRF90 and by the Danish Council for Independent Re-search, grant number DFF 4181-00055.
VI. APPENDIX
We show what is the value of cos θ q in the generic caseof a nonzero φ . Recall that cos θ q = ˆ θ · ˆ (cid:96) with ˆ θ and ˆ (cid:96) aregiven in Eqs. (44) and (45) (see also Eqs. (47) and (48)).Note the coordinate system convention we have used: theperihelion, the center of the Earth and the detector lieon the x − y plane with the perihelion being on the x -axis. In order to find cos θ q for a nonzero value of φ , weneed to go to a reference system that is rotated by anangle φ around the axis that connects the center of theEarth and the detector. Practically this can be achievedby the following coordinate transformations: i) We rotatearound the z -axis by an angle θ . This will make the x -axispass through the detector. ii) We rotate around the new x -axis (the axis passing from the detector and the center3 m χ [GeV]10 R a t e r a t i o -2 -1 Fluorine He m χ [GeV]10 R a t e r a t i o -2 -1 O O FIG. 6. Left: Ratio of the double differential rates in Eqs. (59) and (63) as a function of the DM particle mass. We assume O as dark matter-nucleon interaction. The solid blue line refers to a hypothetical detector composed of Fluorine, whereas thedashed red line corresponds to a second hypothetical detector which uses He as a target material. Right: Same ratio as inthe left panel now evaluated for the operators O (blue solid line) and O (red dashed line) for comparison. In both cases weassume Fluorine as a target material. In the two panels we assume a nuclear recoil energy of 1 eV. of the Earth) by an angle φ . iii) We rotate around thenew z -axis by an angle − θ . The new coordinate system will be given in terms of the old one asˆ x (cid:48) = C · C · C · ˆ x , (67)where ˆ x (cid:48) = (ˆ x (cid:48) , ˆ y (cid:48) , ˆ z (cid:48) ), ˆ x = (ˆ x, ˆ y, ˆ z ) and C , , are the3 × x (cid:48) = (cos θ + sin θ cos φ )ˆ x + sin θ cos θ (1 − cos φ )ˆ y − sin θ sin φ ˆ z ˆ y (cid:48) = sin θ cos θ (1 − cos φ )ˆ x + (sin θ + cos θ cos φ )ˆ y + cos θ sin φ ˆ z ˆ z (cid:48) = sin θ sin φ ˆ x − sin φ cos θ ˆ y + cos φ ˆ z. (68)The velocity of the particle that follows an elliptic orbitwhere the perihelion is rotated around the detector axisby φ should be given by Eq. (45) with ˆ x and ˆ y substituted by ˆ x (cid:48) and ˆ y (cid:48) respectivelyˆ (cid:96) = dxd(cid:96) ˆ x (cid:48) + dyd(cid:96) ˆ y (cid:48) , (69)with dx/d(cid:96) and dy/d(cid:96) given from Eqs. (47) and (48).Using Eq. (68) we calculate cos θ q = ˆ θ · ˆ (cid:96) and obtainEq. (50). [1] L. Baudis, Phys. Dark Univ. (2012) 94doi:10.1016/j.dark.2012.10.006 [arXiv:1211.7222 [astro- ph.IM]]. [2] M. W. Goodman and E. Witten, Phys. Rev. D (1985)3059. doi:10.1103/PhysRevD.31.3059[3] J. D. Lewin and P. F. Smith, Astropart. Phys. (1996)87. doi:10.1016/S0927-6505(96)00047-3[4] T. Marrodn Undagoitia and L. Rauch, J. Phys. G (2016) no.1, 013001 doi:10.1088/0954-3899/43/1/013001[arXiv:1509.08767 [physics.ins-det]].[5] A. K. Drukier, K. Freese and D. N. Spergel, Phys. Rev.D , 3495 (1986). doi:10.1103/PhysRevD.33.3495[6] K. Freese, J. A. Frieman and A. Gould, Phys. Rev. D (1988) 3388. doi:10.1103/PhysRevD.37.3388[7] Bernabei, R., Belli, P., Cappella, F., et al. 2013, Euro-pean Physical Journal C, 73, 2648[8] L. Baudis, Phys. Dark Univ. (2014) 50doi:10.1016/j.dark.2014.07.001 [arXiv:1408.4371 [astro-ph.IM]].[9] E. Aprile et al. [XENON Collaboration], JCAP (2016) no.04, 027 doi:10.1088/1475-7516/2016/04/027[arXiv:1512.07501 [physics.ins-det]].[10] F. Mayet et al. , Phys. Rept. (2016) 1doi:10.1016/j.physrep.2016.02.007 [arXiv:1602.03781[astro-ph.CO]].[11] R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto,T. Volansky and T. T. Yu, JHEP (2016) 046doi:10.1007/JHEP05(2016)046 [arXiv:1509.01598 [hep-ph]].[12] C. Kouvaris and J. Pradler, arXiv:1607.01789 [hep-ph].[13] A. Drukier and L. Stodolsky, Phys. Rev. D (1984)2295.[14] K. Schutz and K. M. Zurek, arXiv:1604.08206 [hep-ph].[15] Y. Hochberg, T. Lin and K. M. Zurek, arXiv:1608.01994[hep-ph].[16] R. Catena and C. Kouvaris, Phys. Rev. D ,no. 2, 023527 (2016) doi:10.1103/PhysRevD.94.023527[arXiv:1602.00006 [astro-ph.CO]].[17] W. H. Press and D. N. Spergel, Astrophys. J. , 679(1985). doi:10.1086/163485[18] Krauss, L. M., Freese, K., Spergel, D. N., & Press, W. H.1985, Astrophys. J., 299, 1001[19] A. H. G. Peter, Phys. Rev. D , 103531 (2009)doi:10.1103/PhysRevD.79.103531 [arXiv:0902.1344[astro-ph.HE]].[20] A. H. G. Peter, Phys. Rev. D , 103533 (2009)doi:10.1103/PhysRevD.79.103533 [arXiv:0902.1348[astro-ph.HE]].[21] T. Damour and L. M. Krauss, Phys. Rev. Lett. , 5726 (1998) doi:10.1103/PhysRevLett.81.5726 [astro-ph/9806165].[22] T. Damour and L. M. Krauss, Phys. Rev. D ,063509 (1999) doi:10.1103/PhysRevD.59.063509 [astro-ph/9807099]. [23] K. Freese, Phys. Lett. B (1986) 295.doi:10.1016/0370-2693(86)90349-7[24] L. M. Krauss, M. Srednicki and F. Wilczek, Phys. Rev.D (1986) 2079. doi:10.1103/PhysRevD.33.2079[25] A. Gould, Astrophys. J. , 571 (1987).doi:10.1086/165653[26] J. Lundberg and J. Edsjo, Phys. Rev. D (2004) 123505doi:10.1103/PhysRevD.69.123505 [astro-ph/0401113].[27] A. Gould, J. A. Frieman and K. Freese, Phys. Rev. D (1989) 1029. doi:10.1103/PhysRevD.39.1029[28] J. I. Collar, Phys. Rev. D (1999) 063514doi:10.1103/PhysRevD.59.063514 [astro-ph/9808058].[29] P. Gondolo, Phys. Rev. D (2002) 103513doi:10.1103/PhysRevD.66.103513 [hep-ph/0209110].[30] A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbersand Y. Xu, JCAP (2013) 004 doi:10.1088/1475-7516/2013/02/004 [arXiv:1203.3542 [hep-ph]].[31] R. Catena and B. Schwabe, JCAP , no.04, 042 (2015) doi:10.1088/1475-7516/2015/04/042[arXiv:1501.03729 [hep-ph]].[32] N. Anand, A. L. Fitzpatrick and W. C. Hax-ton, Phys. Rev. C (2014) no.6, 065501doi:10.1103/PhysRevC.89.065501 [arXiv:1308.6288[hep-ph]].[33] W.F. Mcdonough, Treatise on Geochemistry, Vol 2, Else-vier, 2003. (The values for the Earth composition are veryclose to those in The Encyclopedia of Geochemistry, Eds.Marshall and Fairbridge, Klower Acadmic Publ., 1998.)[34] The Earth: its properties, composition, and structure ,Britannica CD, Version 99 c (cid:13) (2004) 008doi:10.1088/1475-7516/2004/07/008 [astro-ph/0406204].[36] R. Agnese et al. [SuperCDMS Collabora-tion], Phys. Rev. Lett. (2014) 4, 041302doi:10.1103/PhysRevLett.112.041302 [arXiv:1309.3259[physics.ins-det]].[37] A. E. Chavarria et al. , Phys. Procedia (2015)21 doi:10.1016/j.phpro.2014.12.006 [arXiv:1407.0347[physics.ins-det]].[38] Y. Bai and P. J. Fox, JHEP (2009) 052doi:10.1088/1126-6708/2009/11/052 [arXiv:0909.2900[hep-ph]].[39] R. Catena, JCAP (2015) 07, 026 doi:10.1088/1475-7516/2015/07/026 [arXiv:1505.06441 [hep-ph]].[40] B. J. Kavanagh, Phys. Rev. D (2015) 2, 023513doi:10.1103/PhysRevD.92.023513 [arXiv:1505.07406[hep-ph]].[41] C. Kouvaris, Phys. Rev. D93