Measuring the local Dark Matter density in the laboratory

Despite strong evidence for the existence of large amounts of dark matter (DM) in our Universe, there is no direct indication of its presence in our own solar system. All estimates of the local DM density rely on extrapolating results on much larger scales. We demonstrate for the first time the possibility of simultaneously measuring the local DM density and interaction cross-section with a direct detection experiment. It relies on the assumption that incoming DM particles frequently scatter on terrestrial nuclei prior to detection, inducing an additional time-dependence of the signal. We show that for sub-GeV DM, with a large spin-independent DM-proton cross section, future direct detection experiments should be able to reconstruct the local DM density with smaller than 50% uncertainty.


I. INTRODUCTION
A self-gravitating fluid that does not emit or absorb radiation at any observable wavelength, Dark Matter (DM) is the only coherent explanation for a number of otherwise anomalous phenomena [1,2]. These range from stellar motions in nearby dwarf spheroidal galaxies [3] to anisotropies in the cosmic microwave background radiation [4]. There is also strong evidence for the presence of DM in the Milky Way (MW), as inferred from kinematic measurements of stellar populations [5], microlensing events [6] and the dynamics of satellite galaxies [7].
While the evidence for DM in the Universe and in our own Galaxy is compelling, there is no direct indication of the presence of DM within about one parsec of the Sun [8]. The only available estimates divide into two classes: 1) local methods based on the vertical motion of stellar populations [9][10][11][12][13][14][15][16][17][18][19][20]; 2) global methods relying on mass models for the MW [21][22][23][24][25][26][27][28][29]. Each method comes with its own limitations as well as systematic and statistical errors [30][31][32][33][34][35][36][37]. However, no currently known astronomical tracer can directly probe the DM contribution to the MW gravitational potential with sub-parsec resolution [38,39]. Progress in understanding the particle properties of DM [40][41][42], the shape, composition and merger history of the MW [43,44] and, more broadly, the formation of galaxies is hampered by the lack of such ultra-local sub-parsec information about the DM density. In partic-ular, by combining global and ultra-local measurements of the DM density, we can constrain the shape of the MW halo [39]. This in turn may resolve a long-standing tension in standard ΛCDM cosmology: theory and simulations predict triaxial DM halos [45][46][47][48][49][50], while observations in the MW [51][52][53][54] point towards a roughly spherical halo. Although significant effort has recently been made to explain the observed halo shapes with hydrodynamic simulations [55][56][57][58][59], a direct measurement of the local DM density would provide a crucial, independent test of our understanding of Galaxy formation, with important implications for astronomy, astrophysics, cosmology and particle physics.
The lack of direct astronomical measurements of the DM density at the Earth's location also hinders the success of terrestrial 'direct detection' experiments [60][61][62]. These detectors search for DM-nucleus scattering events in underground laboratories, with an expected event rate depending on both the local DM density and the DM-nucleus scattering cross section.
Here, we explore a radically new approach to the problem of finding the local DM density at the Earth's location. We propose to exploit the diurnal variation of the DM flux after Earth-crossing to simultaneously measure the local DM density ρ χ and DM-nucleus scattering cross section σ with future direct detection experiments. This diurnal variation arises from distortions in the DM distribution, due to interactions of DM particles in the Earth before they reach the detector [63][64][65]. The amplitude of this modulation depends on the cross section [66,67], as we will demonstrate via Monte Carlo (MC) simulations, allowing us to break the degeneracy between ρ χ and σ. A similar method has recently been used to measure the high-energy neutrino-nucleon cross section with IceCube [68]. We show that using event timing information, combined with the energy spectrum of a hypothetical DM signal, can enable a measurement of the local DM density and cross-section with low-threshold experiments.
Throughout this paper, we emphasise the DM density measurement, given that our proposal is potentially the only method of directly pinning down ρ χ . 1 We find that the precision of this measurement depends on the detector's location and can be smaller than about 50% for DMproton scattering cross sections larger than 10 −32 cm 2 and a DM mass around 100 MeV, becoming much more precise for larger cross sections. Here, we focus on DMnucleus scattering, but, if extended to DM-electron interactions [70] or more exotic detection strategies (e.g. [71][72][73][74]), our method can be applied to DM candidates in the keV to sub-GeV range, covering a significant fraction of the parameter space of detectable DM candidates [75].
The associated simulation and statistics codes are publicly available at [76] and [77] respectively.

II. DIRECT DETECTION FORMALISM
The differential recoil rate for a DM particle of mass m χ with a nucleus A of mass m A can be written [78,79] with local DM density ρ χ and local DM velocity distribution in the laboratory f (v). Neglecting the effect of Earth scatterings, the usual choice for f (v) in the context of direct detection is the standard halo model (SHM) [78,80,81], a Maxwell-Boltzmann distribution in the galactic frame, truncated at the local galactic escape speed v esc ≈ 544 km s −1 [82,83]. We integrate over v > v min , the minimum speed kinematically required to produce a nuclear recoil of energy E R . It is a crucial feature of this work that for large enough cross sections both ρ χ and f (v) are modified by underground scatterings, thereby modifying the rate in Eq. (1). The true recoil energy E R does not directly correspond to the detected energy deposit E D . We account for a finite energy resolution by transforming the theoretical recoil spectrum of Eq. (1) into the observed spectrum Here, we model the detector response function as a Gaussian with mean E R and standard deviation or energy resolution σ E . For a given energy threshold E th , a finite energy resolution means that a nuclear recoil below the threshold, E R < E th might fluctuate above the threshold and be detectable. However, the approximation of a Gaussian breaks down for energies far below the threshold [84], which is why we set E min R = E th − 2σ E and thereby only include up-fluctuations of 2σ to avoid unphysical signal rates.
We also assume standard spin-independent (SI) interactions for the differential scattering cross section, Here, σ SI p is the DM-proton cross section at zero momentum transfer and A the nucleus' mass number. We consider light DM, m χ m A , so we set the nuclear form factor F 2 (E R ) = 1.

III. EARTH SCATTERING
Above a certain DM-proton cross section σ SI p 10 −37 cm 2 , the probability for a DM particle to scatter on a terrestrial target becomes non-negligible. In this regime, underground scatterings prior to passing through the detector decelerate and deflect the incoming DM particles and thereby change the local DM density and distribution. These distortions grow with the cross section, and the signal thus depends non-linearly on σ SI p . In the single-scattering regime of moderate cross sections, the impact of Earth scatterings on the local DM properties can be quantified analytically [66]. However, the precise contributions of multiple scatterings require the use of MC simulations of underground DM particle trajectories [64,65,67], where we use the numerical tool DaMaSCUS [76] 2 . The simulation details are described extensively in [67,109], and we briefly review the essentials here.
The shape of a DM particle's trajectory is primarily determined by the local mean free path, where n i (x) is the local number density of isotope i, and σ SI i is the total DM-nucleus scattering cross section for that nucleus. The number densities depend on the Earth's mass density profile ρ ⊕ (r), taken from the Preliminary Reference Earth Model (PREM) [110], and the relative nuclear abundances [111]. Furthermore, the distribution of the DM-nucleus scattering angle θ arises from the differential cross section in Eq. (3) and the relation between θ and the recoil energy, E R ∝ (1 − cos θ).
The simulated system features an axial symmetry around the direction of the Earth's velocity v ⊕ . This symmetry allows us to define the isodetection angle Θ [64,65], the polar angle from the symmetry axis as illustrated in Fig. 1. The time-dependent local isodetection angle of a terrestrial observer at x obs reads where r ⊕ ≈ 6370 km is the Earth's radius, and d ∼ 1 km is the underground depth of the observer. It varies over a sidereal day, as described e.g. in App. A of [67].
To extract local estimates based on the MC simulations, we define isodetection rings of finite size ∆Θ = 5 • . By counting the particles passing through each isodetection ring, we obtain an MC estimation of the local DM densityρ χ . By recording their speeds, we obtain a (weighted) histogram estimate of the local speed distributionf (v, Θ) [67]. Finally, these estimates are used to determine the local nuclear recoil spectrum expected for a given value of Θ via Eq. (1). We performed a grid of 45 MC simulations and evaluated the recoil spectra for DM parameters in the ranges m χ ∈ [0.058, 0.5] GeV and σ SI p = [10 −38 , 10 −30 ] cm 2 , accounting for the crucial impact of Earth scatterings. Below m χ ≈ 0.058 GeV, the experimental setups we consider begin to rapidly lose sensitivity, due to the exponential suppression of events above the energy threshold.

IV. EXTRACTING THE LOCAL DM DENSITY FROM DATA
We express the sensitivity of direct detection experiments to the local DM density in terms of contours of constant p-value. We can then reject a point θ = (σ SI p , ρ χ ) on these contours in favour of the alternative, benchmark point θ = (σ SI p , ρ χ ) with a statistical significance of Φ −1 (1 − p), where Φ is the standard normal distribution. For the local DM density, we assume ρ χ = 0.4 GeV cm −3 .
We calculate such p-value contours by using t θ = −2 ln λ(θ) as a test statistic, where λ(θ) is the profile likelihood ratio, defined in Eq. (7) of [112]. We account for the unknown DM mass by maximising the likelihood (at fixed θ) with respect to m χ ∈ [0.058, 0.5] GeV. The p-value calculation requires the probability density function (pdf) of t θ under the assumption that the true model parameters are θ or θ . We denote these pdfs by f (t θ |θ) and f (t θ |θ ), respectively. Following [112], we approximate f (t θ |θ) as a chi-square distribution with k = 2 degrees of freedom and f (t θ |θ ) as a non-central chi-square distribution with the same number of degrees of freedom [112] and non-centrality parameter Λ = −2 ln λ(θ). Here, we restrict ourselves to "Asimov data", defined as the hypothetical dataset such that the maximum likelihood estimator,θ, and benchmark point, θ , coincide. The p-value for rejecting the hypothesised point θ in favour of θ is then given by where t med is the median of f (t θ |θ ). The profile likelihood ratio, λ(θ), depends on the expected number of nuclear recoils from DM signal and background events in the i-th energy bin and in the j-th time bin, s ij and b ij , respectively (see Eq. (7) in [112]). We calculate s ij , i = 1, . . . , N , j = 1, . . . , M , by integrating Eq. (1) over N = 12 (M = 12) energy (time) bins of equal size. We consider two experimental setups. Motivated by existing experiments [113,114], the first detector we consider is a germanium detector, with a Gaussian energy resolution of 18 eV. The energy bins in the analysis cover the energy interval from the assumed threshold, 60 eV, to a maximum energy of 500 eV. The second detector we consider is a cryogenic calorimeter with a sapphire target (Al 2 O 3 ), inspired by the ν-cleus experiment [115,116]. For the energy threshold and resolution, we assume E th = 10 eV and σ E = 3 eV respectively, which should be achievable for sapphire targets with some improvements in detector performance [115]. We assume perfect detection efficiency for both detectors.
The exposure spans a total of 30 days, starting from January 1st 2020, folded onto a single sidereal day, which is then divided into M = 12 time bins. For both detectors, we assume a target mass of 35 g, leading to a total exposure of 1 kg day. 3 We calculate b ij assuming a timeindependent background consisting of a flat component and an exponentially falling component, as observed by EDELWEISS-Surf [113]. We assume that both detectors are operated at a depth d = 1000 m underground.
We consider two benchmark masses for the DM particle. The first benchmark is m χ = 400 MeV, for which the DM-proton scattering cross section is constrained to be σ SI p 10 −37 cm 2 by current direct searches [117]. We consider searches for this particle with the germanium detector (E th = 60 eV). The second benchmark is m χ = 100 MeV, which is significantly less constrained: cross sections of σ SI p 5 × 10 −31 cm 2 are still allowed by current constraints [118]. A very low threshold is required for sensitivity to such light DM and we therefore consider searches for this particle with the sapphire detector (E th = 10 eV).

V. RESULTS
The projected p = 0.05 contours in the (σ SI p , ρ χ )-plane are shown in Fig. 2. The upper (lower) panel shows reconstructions for a hypothetical direct detection experiment in the Northern (Southern) Hemisphere, located at Laboratori Nazionali del Gran Sasso (LNGS, 46 • N) and Stawell Underground Physics Laboratory (SUPL, 37 • S) respectively [119,120]. In the left (right) panels, we show results for 400 MeV DM particle and germanium detector (100 MeV DM particle and sapphire detector).
While we have performed the analysis over a wider range of cross sections, we do not find closed contours in the (σ SI p , ρ χ ) plane for cross sections smaller than σ SI p . This indicates that for our heavier benchmark mass of 400 MeV (left panels), there remains no unconstrained region of parameter space for SI interactions, for which a substantial modulation effect due to Earth-scattering can be observed. We therefore focus in the remainder of this work on the light benchmark of mass 100 MeV (right panels).
We note also that it is also not possible to obtain closed contours in the (σ SI p , ρ χ ) plane when taking only the recoil energies of the observed signal events into account, assuming no knowledge of their timing information. This is because of the strong degeneracy between σ SI p and ρ χ , which can not be broken by the energy data alone.
In contrast, keeping track of the signals' timing and accounting for their modulation signature improves the situation drastically, as seen in the colored contours of Fig. 2, for the scenario of a light DM particle (right panels). In the case of the Northern experiment and cross sections above about 10 −33 cm 2 , the degeneracy between DM density and scattering cross section starts to become weaker. For higher benchmark cross sections, the true density as well as the cross section itself can be reconstructed with increasing precision. For example, for σ SI p = 10 −32 cm 2 (σ SI p = 10 −31 cm 2 ) we could determine the local DM density to be ρ χ = 0.40 +0. 26 −0.03 GeV cm −3 (ρ χ = 0.40 +0.01 −0.02 GeV cm −3 ) at 95% CL. In these cases, the cross section would be constrained to σ SI p = 1.00 +0.07 −0.33 × 10 −32 cm 2 (σ SI p = 1.00 +0.03 −0.01 × 10 −31 cm 2 ). Projected contours for an experiment at SUPL (lower panel) show a similar evolution. However, the reconstruction of the local DM density and cross section is generally less precise than in the Northern Hemisphere. In the case of a benchmark cross section of σ SI p = 10 −32 cm 2 , a closed contour is obtained though the constraints on ρ χ are very wide (extending over the entire range of our analysis, from 0.01 GeV/cm 3 to 1.0 GeV/cm 3 ). Instead, for a benchmark cross section of σ SI p = 10 −32 cm 2 , we find ρ χ = 0.40 +0. 23 −0.19 GeV cm −3 and σ SI p = 1.00 +3.77 −0.21 × 10 −31 cm 2 . These results indicate that with future ultra-low threshold detectors, it should be possible to reconstruct both the local DM density and cross section, if the DMproton cross section lies within a few orders of magnitude of current constraints, for a DM mass of 100 MeV.

VI. DISCUSSION
From these results, the necessity of timing information is obvious. Without time-tagging, it is always possible to reabsorb a change of the local DM density into a rescaling of the cross section such that ρ χ × σ SI p remains constant. The time-dependence of the local DM distribution in the laboratory, caused by underground scatterings, introduces an additional dependence of the signal on σ SI p . Since this dependence manifests itself through diurnal modulation, knowledge of event timings is the key to disentangling the local DM density and the cross section.
Contrary to our intuition 4 , we find that experiments in the Northern Hemisphere are generally better suited to measuring the local DM density. For an experiment at LNGS, Θ(t) varies in the range [4 • , 84 • ], so the bulk of the incoming DM flux reaches the laboratory directly from above at a certain time of day, while it has to pass through a small fraction of Earth's mantle 12 hours later. Therefore, the experiment switches continuously between being totally exposed to and partially shielded from incoming DM particles, as illustrated in Fig. 1. Increasing the cross section increases the modulation amplitude and steadily improves the reconstruction of ρ χ .
However, in order to reach a direct detection experiment in the Southern Hemisphere, most DM particles need to traverse the planet's bulk mass throughout the day. An experiment at SUPL is always partially shielded, with Θ(t) ∈ [86 • , 167 • ]. For cross sections between 10 −34 − 10 −32 cm 2 , the Earth's stopping power renders the majority of the DM wind undetectable. The slower sub-component, which arrives from the opposite direction passing the atmosphere and overburden only, is not yet affected (see again Fig. 1). In this regime, the modulation amplitude depends only weakly on the cross section, and estimates of ρ χ and σ SI p are less precise than in the Northern Hemisphere. This demonstrates that the determining factor for reconstructing the DM density is not so much the diurnal modulation's amplitude, but rather its sensitivity to changes in the cross section. To better understand the contour shapes in Fig. 2, we focus on the benchmark point m χ = 100 MeV, σ SI p ≈ 3 × 10 −33 cm 2 and a sapphire experiment in the Southern Hemisphere. The left panel of Fig. 3 shows the projected contour and best fit mass at each point in parameter space (σ SI p , ρ χ ), while the right panel shows the loglikelihood across three fixed cross section slices (A-C). In each slice, the curved region where the log-likelihood peaks corresponds roughly to a region where the total number of signal events is constant. The dominant effect is that increasing the DM mass from the benchmark value exponentially increases the number of events above threshold. This is because the typical recoil energy deposited by a 100 MeV particle in a sapphire detector O(8 eV), is so close to our assumed threshold of 10 eV. This exponential increase in the number of signal events must be compensated by a decrease in ρ χ , as seen in the right panel of Fig. 3.
Focusing on slice C, overestimating the scattering cross section would mean predicting more events and a greater modulation amplitude. The former is compensated by a lower best-fit value of ρ χ , while the latter is compensated by an overestimate of m χ (red triangle). Though increasing m χ increases the Earth's stopping power 5 , it also increases the typical recoil energy which can be deposited by a DM particle. For signals close to threshold, this latter effect wins out, meaning that a larger DM mass reduces the overall attenuation of the signal and lowers the modulation amplitude due to Earth scattering. Thus, in slice C the log-likelihood peaks in a region of higher mass, over a restricted range of values for ρ χ . In general, the slope of the log-likelihood contours around the peak then determines the uncertainty on ρ χ and thereby the contour shape in the left panel of Fig. 3.

VII. CONCLUSIONS
Provided that DM-matter interactions are sufficiently strong for underground scatterings to occur frequently, signals in direct detection experiments should show a diurnal modulation, which can be exploited to break the degeneracy between ρ χ and σ SI p . We have explored this possibility for a number of benchmarks using MC simulations.
For the case of a DM mass of 400 MeV (left panel of Fig. 2), we find that both ρ χ and σ SI p can be reconstructed only for cross sections close to 10 −30 cm 2 , substantially larger than allowed by current constraints. Lighter DM of mass 100 MeV (right panel of Fig. 2) should be accessible to near-future low-threshold nuclear recoil searches. In this case, it should be possible to disentangle the local DM density and the scattering cross section for models just below current constraints, in the range σ SI p ∈ [10 −33 , 10 −30 ] cm 2 . For this, detectors in the Northerm Hemisphere offer the best prospects, with a O(50%) (O(5%)) reconstruction of ρ χ possible for σ SI p = 10 −32 cm 2 (σ SI p = 10 −31 cm 2 ). This is the first demonstration that it is possible to measure the local Dark Matter density directly in the laboratory and further motivates the search for light, strongly-interacting DM with low-threshold detectors.