Observationally inferred dark matter phase-space distribution and direct detection experiments

We present a detailed analysis of the effect of an observationally determined dark matter (DM) velocity distribution function (VDF) of the Milky Way (MW) on DM direct detection rates. We go beyond local kinematic tracers and use rotation curve data up to 200 kpc to construct a MW mass model and self-consistently determine the local phase-space distribution of DM. This approach mitigates any incomplete understanding of local dark matter-visible matter degeneracies that can affect the determination of the VDF. Comparing with the oft used Standard Halo Model (SHM), which assumes an isothermal VDF, we look at how the tail of the empirically determined VDF alters our interpretation of the present direct detection WIMP DM cross section exclusion limits. While previous studies have suggested a very large difference (of more than an order of magnitude) in the bounds at low DM masses, we show that accounting for the detector response at low threshold energies, the difference is still significant although less extreme. The change in the number of signal events, when using the empirically determined DM VDF in contrast to the SHM VDF, is most prominent for low DM masses for which the shape of the recoil energy spectrum depends sensitively on the detector threshold energy as well as detector response near the threshold. We demonstrate that these trends carry over to the respective DM exclusion limits, modulo detailed understanding of the experimental backgrounds. With the unprecedented precision of astrometric data in the GAIA era, use of observationally determined DM phase-space will become a critical and necessary ingredient for DM searches. We provide an accurate fit to the current best observationally determined DM VDF (and self-consistent local DM density) for use in analyzing current DM direct detection data by the experimental community.


INTRODUCTION
One of the most popular candidates for the dark matter (DM) particle is the hypothetical weakly interacting massive particle (WIMP) [1][2][3][4][5][6]. The determination of the density of these particles in the halo of the Milky Way (MW) galaxy in general and the solar neighborhood in particular is crucial for many direct detection experiments which attempt to measure the rate of nuclear recoil events caused by the WIMPs scattering off of the detector target nuclei. The expected scattering behavior is strongly dependent on the local astrophysical properties of DM. In particular, the scattering rate is directly proportional to the local DM density. In addition, the velocity distribution function (VDF) of the DM crucially affects the shape of the nuclear recoil energy spectrum. Thus, it is imperative to have precise knowledge of the local phase-space distribution of DM from observations in order to set precise bounds on the DM particle physics parameter space [7]. DM direct detection experiments usually assume the simplest possible 'Standard Halo Model' (SHM) for the DM halo, in which the velocity distribution is Maxwellian. This model assumes the halo to be an isotropic, isothermal sphere -hypotheses that are unlikely to be valid in reality. Moreover, N-body simulations produce halos with velocity distributions which de-viate systematically from a Maxwellian [8][9][10][11][12][13]. One can also construct more realistic, analytical VDFs that differ from the predictions of the SHM [14][15][16][17]. However, a self consistent connection to observations is generally missing -for example, simulations obtain a VDF of a Milky Way 'like' halo identified inside a simulation box using user defined criteria (such as halo mass or circular velocity); on the other hand, analytic models typically cannot incorporate the effect of the baryonic mass component of the MW on the DM VDF.
In this work, we advocate for a more realistic, observationally-driven approach which follows a three step procedure -(i) first, visible matter tracers in the Milky Way are used to map out the gravitational potential, Φ(r), in the region of interest, (ii) second, a multicomponent mass density model, having both DM and visible matter (VM) in various configurations is obtained consistent with the observed gravitational potential, and (iii) finally, since the density is an integral of the total phase-space distribution function, one can invert the equation connecting the two to obtain the VDF. A convenient way to get the inversion is to use the Eddington formalism [15]. The full procedure results in a selfconsistent determination of both the DM density and its velocity distribution. This was first done by Bhattacharjee et al. [18]; this work is based on similar, but more detailed, MW DM phase-space analysis using a larger dataset [19].
Our approach, of using rotation curve (RC) data up to ∼ 200 kpc as the visible matter tracer, makes us sensitive to the DM and VM distributions throughout the halo. In this approach it is easier to separate the contributions of the different components, modulo the analytic form of the DM distribution that is assumed [67]. Estimates of the local DM phase-space based on local dynamics rely on the accuracy of separating the VM contributions and are prone to contaminations, viz. without the leverage of a large number of galactic radial bins, a Bayesian analysis marginalizing over the VM contribution using only local data has large degeneracies. Our method of using data at a large number of radial points leads to substantially better breaking of the degeneracy between the VM and DM distributions.
For simplicity, experimental DM direct detection exclusion curves in the literature are almost always calculated by assuming a naïve SHM expectation for the MW DM halo. There have been a number of theory papers that have tried to go beyond the SHM and assess the impact of astrophysical uncertainties on DM direct detection experiments (for example, see [20,21]). Some of these attempts have examined the impact of an uncertain local dark matter density (such as [22]), which results in a trivial re-scaling of the exclusion curves. Others have attempted to look at the effect of the DM velocities via changes in the local DM escape speed and the local DM velocity dispersion (see [23]). Studies that tried to incorporate the full local DM VDF, other than SHM, have been mainly restricted to ansatzes or VDFs extracted from simulations [20,23,24].
We claim that the right approach should be to use an observationally inferred determination, along with the related uncertainties, of the local DM phase-space in a selfconsistent manner (also see [68]). In this work, we take a detailed look at this fundamentally important element of DM detection results. The main result of our work is the first re-estimation the DM exclusion curves, for some of the major DM direct detection experiments, using observationally determined local DM phase-space.

SELF-CONSISTENT DETERMINATION OF THE DM PHASE-SPACE DISTRIBUTION
The rate of nuclear recoil events, in direct detection searches, depends crucially on the local (i.e. solar neighborhood) density and velocity distribution of the WIMPs in the Galaxy, which are a priori unknown. In contrast to the density, not much knowledge directly based on observational data is available on the likely form of the velocity distribution function (VDF) of the WIMPs in the Galaxy. The standard practice is to use what is often referred to as the Standard Halo Model (SHM), in which the DM halo of the Galaxy is described as a single-component isothermal sphere, for which the VDF is assumed to be isotropic and of the Maxwell-Boltzmann form. High resolution cosmological simulations of DM halos give strong indications of significant departure of the VDF from the Maxwellian. On the other hand, these cosmological simulations do not yet satisfactorily include the gravitational effects of the visible matter components of the real galaxy, namely, the central bulge and the disk.
One approach to determining the local density of DM is to use the rotation curve data to find the likelihood of the parameters characterizing the density distributions of the various mass components of the galaxy. In general, the visible matter (VM) parameters are fixed (from observational data) and the dark matter (DM) parameters are obtained from a likelihood maximization. A full likelihood analysis (DM and VM) was first done [18] by taking the Navarro-Frenk-White (NFW) profile [25] for the DM halo, a spheroidal bulge and an axisymmetric disk [26][27][28][29][30]. The NFW DM density profile is given where r is the distance from the galactic center, ρ DM, is the local DM density, r s the scale radius of the halo and R is the distance of the sun from the galactic center. The bulge and disk density profiles are given respectively by ρ b (r) = where ρ b0 is the normalization of the bulge density, Σ is the local disk surface density, and r b and R d are the scale radii of the bulge and the disk respectively. The parameter z d is the scale height of the disk. The visible matter parameterizations are based on fits to local kinematical data. This fiducial model of the MW consisting of a dark matter halo, visible matter bulge and a single disk is a minimal model of the mass distribution in the MW [69]. For a given choice of the density profiles of both DM and VM, we can use the Poisson equation to obtain the total gravitational potential Φ(R) at a given radius on the galactic plane, and thus we get the circular rotation speed at R from the relation, Next, a Markov Chain Monte Carlo (MCMC) analysis is carried out to determine the most likely values of the density parameters (along with their 1σ uncertainties). In this work, we adopt the Python framework CosmoHammer [31] which embeds the emcee package by Foreman-Mackey et. al. [32] that is based on an improved MCMC algorithm by Goodman and Weare [33]. The χ 2 test statistic used within the MCMC is where v obs,i is the observed circular velocity value, v theo,i is the one theoretically calculated, σ i the error in the observed velocity value and N is the total number of binned data points at different distances from the galactic center. The best fit density parameters are used to estimate the full spatial density of the DM particles ρ DM (r) and the total gravitational potential Φ(r).
Under isotropic conditions , the phase space distribution function F of the DM component, at a position r, depends only on the total specific energy E = 1 2 v 2 +Φ(r), with v = |v|, r = |r|. This function F can be uniquely determined using the Eddington formula [34], 2 v 2 , and ρ(r) is the total density. The VDF at a radius r can be obtained as, Also, for E > 0, F > 0, and for E < 0, F = 0; this ensures that the VDF truncates naturally at the escape velocity v esc = 2|Ψ(r)|. We work with the normalized, A key assumption in using Eddington's method is that of isotropy (i.e, the net potential has spherical symmetry). Due to the axisymmetric nature of the VM disk, the total potential in our MW mass model is non-spherical. In order to use the Eddington approximation, we use a spherical approximation [35] for the VM potential given by Φ V M (r) The DM density ρ DM (r) is obtained by integrating the single-particle phase-space distribution function (which in the simplest case is a function of total energy per unit mass of a DM particle) over the velocities of the DM particle. Under the spherical approximation, ρ DM (r) = F[E = (1/2)v 2 +Φ(r)]4πv 2 dv where v is the DM velocity and F is the phase-space distribution function, which is inverted give the DM VDF. Note that ρ DM is also implicitly present on the right hand side of the above equation in Φ(r), and hence the equation demands a selfconsistent solution. For a given DM and VM distribution that is consistent with the particular Φ(r) 'at all r', a unique solution for the VDF f (r, v) ≡ F(E)/ρ DM (r) 'at all r' exists in the spherically symmetric case.
The VDF (in eq. (3)) determined using the technique outlined in the preceding discussion is completely empirical, without any reference to any simulations; it is connected to the local DM density in a self-consistent manner through the Eddington formula and we only need to model the potential contribution from the disk as arising due to a spherical mass distribution. We find that the spherical approximation induces corrections of the order of ≈ 10% in the value of v 2 (R )).
We stress that the arguments presented in this section apply to virialized DM halos. Kinematic outliers associated with local DM substructure such as streams or debris flows [36][37][38][39][40][41][42][43] could also impact the interpretation of the experimental results. Because the origin of these outliers is unknown, however, we choose to not include them. Additionally, recent mergers of satellite galaxies could lead to spatial or kinematic substructure. We also ignore the possibility of velocity spikes in the VDF due to local substructure [8].

VELOCITY PROFILES
To determine the RC of the Galaxy, one has to use kinematical and positional information for some tracer objects moving in the gravitational potential of the galaxy. In general, one measures line-of-sight (LOS) quantities (positional and kinematical data) and the RC is derived from these. To determine the RC for the disk region, we have to adopt a value for the the local standard of rest (LSR) which corresponds to the position (R 0 ) and velocity (v c0 ) of the sun with respect to the galactic center, and make the assumption that the tracer objects follow a circular orbit around the galactic center. From this, the positions and velocities about the galactic center can be obtained. The choice of LSR plays a crucial role in determining the parameters of the mass model of the MW -it affects the value of local DM density and our estimation of other MW properties, like the mass of the MW. Following [18,44], we pick two popular LSRs used in MW studies: (i) R 0 = 8 kpc and v c0 = 200 km/s and (ii) R 0 = 8.5 kpc and v c0 = 220 km/s. The various tracer objects used for constructing our RC include HI and HII regions (CO emissions from the latter), Cepheids, Planetary nebulae, etc. For regions extending beyond the visible disk of the galaxy, we look at tracers distributed in the halo of the Milky Way (like dwarf spheroidals, globular clusters, K-giant stars, etc.). The latter tracers are of the non-disk kind, and their motion around the galactic center is typically unsystematic. Under the assumption that these objects are isotropically distributed in the halo, one can define an effective circular velocity v c at a galactocentric distance r, and use the Jeans equation to relate it to the observed number densities and velocity dispersions [19,34]. One disadvantage of using these non-disk tracers is that they can have non-negligible velocity anisotropy. Currently, only the line-of-sight velocity dispersion is known to precision at large distances, the velocity anisotropy cannot be well determined [70]. In this work, we neglect the effect of any such anisotropy in constructing the RC. For the effect of velocity anisotropy on the RC at large distances, we refer the reader to [19,45,46].
After the raw RC is generated as above, it is suitably radially binned and averaged. The binning strategy used by Bhattacharjee et.al. [19] is two-fold. In the first step, they average over each individual data set, choosing different bin sizes at different radius ranges (smaller bins at low r where there is more data, etc). These bin sizes are manually optimized to best reflect the overall behavior of the raw data points. Once these individual datasets are binned and averaged, they are compiled into a larger dataset, and the above process is repeated once again to arrive at the final dataset.
Following the detailed procedure described above, we performed an MCMC analysis on the RC compiled by Bhattacharjee et al. [19] to obtain the best fitting DM plus VM density distributions of our galaxy. There is a wealth of knowledge available on the distribution of visible matter in the MW based on decades of astrophysical observations, and it is prudent to add some of this information in the form of VM priors in our MCMC analysis. In this Bayesian approach, the final best fit DM distribution depends on the imposed priors on the VM distribution.
Using the Eddington formula, we estimated the local DM VDF, which is self-consistently related to the DM and VM mass distribution, and in particular to the DM local density. Note, that the extracted DM VDF implicitly depends on the choice of VM prior and LSR. Although flat priors give the most 'unbiased' best fit parameters, we use the wealth of knowledge on the VM distribution to impose comparatively tighter local VM priors. For all the results quoted in the rest of the paper, we chose our priors based on observational constraints on the local VM density [28,[47][48][49], which can be expressed in terms of constraints on the disk parameters Σ and R d . For our analysis, we adopt the following gaussian priors on the disk parameters: Σ = 67 ± 8 M pc −2 and R d = 2.3 ± 0.6 kpc. Mandal and Majumdar [44] have studied the impact of the choice of local VM prior and LSR on the estimates of the MW DM distribution and, in this work, we adopt their method and results.
Together with our choice of VM prior and two sets of LSRs, we determine the corresponding local DM densities and VDFs: 1. B200-8.0-67: For the choice of LSR with R 0 = 8 kpc and v c0 = 200 km/s, we find the best fit local DM density to be ρ DM, = 0.18 ± 0.02 GeVcm −3 .
The corresponding VDF is shown in Fig. 1 with the high velocity cutoff corresponding to the local escape velocity v esc = 536.83 km s −1 .
In the rest of the paper, the observationally determined VDF B220-8.5-67 is selected as our fiducial VDF and we denote it as obs for the sake of brevity. From Fig. 1, we can see that the VDF B200-8.0-67 is very close to our fiducial VDF and the main difference between the expected DM signal rates for the two choices of LSR enters through a trivial rescaling of the rate due to the different best-fit local DM densities.
The fractional uncertainty on the local DM density (as listed above) is between 5 − 10% irrespective of the choice of LSR or any observational priors on the VM. In contrast, estimates of ρ DM, using local kinematics typically has an uncertainty of ∼ 30% [50]. We attribute this difference to our use of data at a larger number of radial bins which leads to the breaking of any DM-VM degeneracies. In the rest of the paper, we neglect this small uncertainty on ρ DM, since the VDF uncertainty becomes the dominant astrophysical uncertainty for DM direct detection searches.
Note, that the VDF is determined from the mean values of the determined model parameters, and the 1-σ errors on these parameters gives us an uncertainty band on the VDF. The best fit parameters of the MW mass model have correlated uncertainties and these correlations are naturally computed in the MCMC analysis. These correlations lead to correlations in errors on the VDF at different velocities [71]. In our analysis, we use the above RC datasets with the error bars on the velocities (a) having the reported values (which we call current), and (b) reduced by 1/3 (which we call third ). The reduction of the error bars on the RC unsurprisingly leads to the narrowing of the error bands on the VDF. Our analysis with the reduced error bars models upcoming data from the GAIA satellite [51]. Compared to the ∼ 50 RC data points [19] that we have used in this work, GAIA data will potentially have an order of magnitude larger number of RC data points within 200 kpc which will reduce the relative error on the velocity by a factor of σv v ∼ 1 √ N ∼ 1 3 , where N = 10 is a benchmark increase in the number of tracers that we expect, thus giving tighter constraints on the VDF and the local DM density.
As a benchmark for comparison, we introduce the SHM which is conventionally used in all direct detection DM experimental analyses. The VDF is assumed to be given by a Maxwellian distribution. Note, that the corresponding self-consistent DM density distribution assumes that the DM halo is a singular isothermal sphere. Moreover, the SHM assumes a single component halo mass model, an assumption which breaks down in the presence of VM. The isothermal velocity profile of dark matter as a function of the dark matter velocity in the galactic rest frame v is given by, where v d is the velocity dispersion taken to be 220 km/s, and k iso is a numerically determined normalization constant found by integrating the velocity profile over d 3 v, with v < v esc . The canonical value of the galactic escape velocity v esc is taken to be 544 km/s. We denote the SHM VDF with these parameters as iso. The local DM density is conventionally chosen to be ρ DM, = 0.3 GeVcm −3 ; however, this value is not self-consistent with the assumed VDF parameters.

Analytic fits to the observational VDFs
Finally, for the sake of ease-of-use of our VDFs estimated from MW observations, we fit the VDFs to the following functional form where Here k , p, v and v esc are fit parameters and B is an appropriate normalization factor. The best fit parameters are given in Table I ; the fit parameters for the error envelopes are given in Appendix I. We show these fits for the obs VDF in Fig 2. The analytic expression in eq. (5) fits the obs VDF to roughly 1% accuracy over most of the velocity range, and is a better fit to the DM VDFs over other fitting forms proposed in the literature by Bhattacharjee et al. [18] and Mao et al. [24]. A table of the VDFs along with the upper and lower uncertainty envelopes can be found at

DIRECT DETECTION RATE
We will now examine the effect of the difference between the observationally determined (obs) velocity profile and the canonical isothermal SHM profile that we have considered in the previous section on the dark matter direct detection rate. We will consider here only elastic scattering of DM with a target nucleus. Assuming isotropic scattering in the center-of-mass frame of the DM-nucleus system, the rate of direct detection signal events per unit recoil energy (E R ), per unit detector mass, is given by [52], The nominal rate R 0 is given by, where 1 tru is 1 count/kg/day. Here, the dark matter mass m D and target nucleus mass m T are expressed in GeV/c 2 , σ 0 is the DM-target-nucleus cross-section, ρ DM, is the local dark matter density and v 0 is a "typical relative velocity" parameter which is representative of the dark matter velocity in the detector rest frame. The factor of v 0 is introduced only for dimensional convenience and cancels out in the full expression for the rate.
An explanation of the other factors in eq. (6) is in order here. In accordance with [52], E 0 = 1 2 m D v 2 0 is the characteristic recoil energy of the nucleus (typically of the order of a few keV for a 100 GeV dark matter particle) and r = 4m D m T /(m D + m T ) 2 . F 2 (E R ) is the nuclear form factor which is target dependent, here we use the Helm form-factor from [52]. (E R ) is the detector efficiency as a function of recoil energy.
For spin-independent (and isospin independent) scattering the WIMP-nucleus cross-section σ 0 can be expressed in terms of the DM-nucleon cross-section σ n by, where A is the nucleon number of the target, µ D,N (µ D,n ) is the reduced mass of the DM particle and the nucleus (nucleon). We will only consider spin-independent interactions in this work. The case of spin-dependent scattering is a trivial extension.
In eq. (6), the factor I(E R ) is a dimensionless velocity averaged integral given by, Here, the integrand is evaluated over the relative velocity of a DM particle and the detector v r = v gal + v e , where v gal is the DM velocity in the galactic rest frame and v e is the earth's velocity. Since we will only focus on the time integrated recoil signal in this work [72], we take the earth's velocity to be a constant with magnitude v e = 240 km/s. Demanding that the dark matter particle should have a large enough relative velocity to be able to cause a recoil energy E R in the detector gives us the lower bound on the relative velocity, v min = (2E R /(rm D )) 1/2 . This introduces an explicit dependence of I on the recoil energy as well as the DM mass. The velocity integral fully captures the way that the VDF of the dark matter affects the detection rate. Our VDFs by definition are cut-off above the galactic escape velocity. Thus, when v min exceeds v esc + v e , no recoils should occur. For a fixed DM mass, this implies that there is a maximum nuclear recoil energy given by E max R = 1 2 rm D (v esc + v e ) 2 above which no signal is expected in the detector. We note that since the obs VDF has a lower escape velocity compared to the isothermal VDF, we expect lower recoil energy cut-offs for a fixed DM mass when using the obs velocity profile.
Conversely, for a fixed recoil energy, if m D is sufficiently large, then the velocity integral receives contributions from all possible relative velocities. In the absence of the factor v 0 /v r in the integrand, the integral would evaluate to unity, independent of the velocity profile. Thus, for large DM masses, the distinction between the various VDFs arises mainly from the spread of DM relative velocities. However, for low DM masses, the support of the velocity integral shrinks and the rate becomes highly sensitive to the differences between the high velocity tails of the VDFs.
We can see from Fig. 1 that the obs and isothermal VDFs have similar widths but the obs velocity profile has a more suppressed high velocity tail. Thus, for a fixed recoil energy, we expect that the difference between the direct detection rates for the obs and isothermal VDFs will be most dramatic for low dark matter masses.
We define ζ ≡ I obs /I iso as the ratio of I factors for the fiducial obs velocity profile and the isothermal profile. In Fig. 3, we plot ζ as a function of the dark matter mass, for several different target nuclei, for a fixed recoil energy E R = 5 keV. For low DM masses compared to the target mass, v min → (E R m T /2(m D ) 2 ) 1/2 . Thus, for a given target and a fixed recoil energy, the I factor is sensitive to the tail of the VDF. The tail of the obs VDF falls below that of the isothermal profile at a galactic velocity v dev 330 km s −1 (see Fig. 1). Hence, the I ratio drops below unity for low DM masses, m D → In Fig. 4, we show how ζ varies as a function of dark matter mass for a Si target at various recoil energies. We can see from the figure that at higher recoil energies the deviation of ζ from unity occurs for correspondingly higher values of the dark matter mass. The error bands in For a given recoil energy it thus seems that the Ifactors can vary by several orders of magnitude for low dark matter masses. We would thus expect a large change in the total direct detection signal rates and consequently the exclusion bounds for low dark matter masses if we used the obs VDF rather than the canonical isothermal profile. However, as we shall see next this expectation is tempered by the fact that at low recoil energies we actually get a small contribution to the overall rate due to the low detector efficiency at these energies.

RESULTS
In this section we would like to see whether our expectation of the strong sensitivity of the recoil rate to the VDF continues to hold upon including detector efficiency effects. In Fig. 6 we plot the detector efficiency for several existing experiments CRESSTII [53], LUX [54], and PICO [55]. We also show the efficiencies for the proposed SuperCDMS experiment using published projections for the Silicon and Germanium High Voltage (HV) detectors [56][57][58]. Each efficiency function can be characterized by a threshold energy, (E T ) which can be defined as the recoil energy for which the efficiency drops to 50% of maximum efficiency, and the width of the detector response near the threshold (conventionally defined by parameterizing the response as an error-function, see discussion on detector specifications in Appendix II).
Using these efficiencies, we can now compute the re-  [54], PICO [55], CRESSTII [53] and the proposed SuperCDMS Si and Ge High Voltage detectors [56-58]. coil energy spectra for different detectors for the obs and SHM VDFs using eq. (6). Note that we need to choose the appropriate value of the local DM density corresponding to the VDF that we are using, which changes the overall normalization of the rate for different VDFs. It is the strength of our present approach, as detailed in previous sections, that we have self-consistent pairs of local DM density and VDFs. In Fig. 6 we plot the recoil spectra due to DM-nucleus interactions for the SuperCDMS Silicon detector and the LUX detector for DM masses of 6 GeV and 10 GeV, assuming a WIMP-nucleon crosssection σ n = 10 −40 cm 2 .
We note a few interesting features of the recoil spectra. The high energy tail of the recoil spectra is sensitive to the tail of the VDFs. However, at low recoil energies the shape of the recoil spectrum is determined by the threshold energy and efficiency of the detector near the threshold. We note that each detector has a minimum dark matter mass below which no events are seen in the detector. This minimum mass is given by m min where v esc is the local escape velocity and E min R is the minimum deposited energy that can be detected. Since CRESST and Super-CDMS are low threshold experiments, they are sensitive to much lower recoil energies and hence to much lower dark matter masses.
In Fig. 7 we show the differential recoil rate dR/dv gal as a function of the magnitude of the dark matter velocity in the galactic rest frame (v gal = v r + v e ) for DM particles scattering in LUX and SuperCDMS with the same benchmark cross-section. Comparing this plot with the VDFs in Fig. 1 allows us to visualize the contribution to the total recoil signal for different ranges of DM galactic velocities. From the figure, we can clearly see the recoils of low mass DM probe the high velocity tail of the VDFs.
Finally, we can examine the difference between the   7: Plot of the contribution to total direct detection rate (solid lines) for a 6 GeV and 10 GeV dark matter particle scattering in LUX and SuperCDMS (Si) as a function of v gal (the speed of DM in the galactic rest frame). We have assumed a WIMP-nucleon scattering cross-section σn = 10 −40 cm 2 . The thick and thin envelopes indicate the uncertainty on the recoil rate estimated by a propagation of the current and third velocity envelopes respectively, of the obs VDF. The dashed lines show the corresponding rate assuming the iso VDF.
fully integrated signal rates for the obs and isothermal VDFs. In Fig. 8, we plot as a function of dark matter mass the relative difference (in percent) of the expected number of signal events using the observationally determined VDFs and the expected number obtained by using the SHM. To derive a realistic upper bound on the DM-nucleon cross-section (σ n ) from experiment requires an understanding of the detector background recoil spectra shape and corresponding systematic uncertainties, understanding of the detector efficiency and threshold as well as knowledge of the detector exposure. A detailed study of these effects is typically possible only within each experimental collaboration. Here, we have estimated a simple bound on σ n for the LUX, PICO and CRESSTII [53] experiments. We have also computed an expected exclusion limit for the SuperCDMS Silicon and Germanium high-voltage detectors.
Our procedure was as follows: We first assumed detector exposures and efficiencies for several detectors based on their published results where applicable. These exposures and efficiencies are summarized in the appendix, along with the references from which they were obtained. For all experiments we assumed a background rate of 1 dru ≡ 1 count/keV/kg/day which is constant over recoil energies from 0-200 keV. We then calculated a simple median 90% CL estimated bound on the dark matter nucleon cross-section σ n by doing a simple counting pseudoexperiment over the entire recoil energy range from 0-200 keV, assuming that only background is observed.
The upper end of the recoil energy range of 200 keV is an arbitrary choice. This value is much higher than the maximum recoil energy expected for low DM masses. We are unable to perform a profile likelihood analysis [59] which would take into account shape differences between the signal and background without a detailed understanding of the detector backgrounds, but a detailed experimental analysis would optimize the choice of this upper cut-off for every candidate DM mass. We assumed a simple flat background spectrum, and therefore our choice of the upper recoil energy cut-off of 200 keV added a fixed amount of background that could realistically be reduced. We can see that the difference is most significant for masses close to the detector threshold mass m min D . This difference is significant when VDF errors improve to third errors (shown); with current errors (not shown) the difference cannot be resolved.
We plot our estimated upper-bound on σ n for the obs and SHM velocity profiles for the selected experiments in the upper panel of Fig. 9. Note that it is important to use appropriate choice of local DM density corresponding to the selected VDF when computing the exclusion curves. We also show the uncertainty bands around the exclusion curves for the obs VDF corresponding to a straight forward propagation of the third VDF uncertainty envelope from Fig. 1 . The lower panel of Fig. 9 shows the percentage difference between the exclusion bounds set by assuming the obs and SHM local DM phase-space.
There are several interesting features in Fig. 9 . We note that there is a few percent normalization difference in the bounds for the different velocity profiles, which is apparent at high DM masses, due to the difference in the local DM density for each profile choice. In addition to the density effect, the estimated bounds on σ n differ in shape, at low dark matter masses, due to impact of the different velocity profiles . This difference is most stark at the threshold mass m min D of each detector where m min D (E T m T /2) 1/2 /(v esc + v e ) and E T is the detector threshold energy. The exact value of m min D , as well as the shape of the exclusion curve near m min D , depend on the tails of the VDFs. Additionally, they also depend on the detector sensitivity near the detector threshold, which is different for each detector.
For example, taking the case for LUX, which has m min D 4 GeV, the exclusion limit obtained by assuming the obs VDF profile as opposed to the SHM VDF yields > 200% difference near threshold, i.e. the constraints are weakened by a factor of three. Similarly, for PICO, m min D 3 GeV, and the exclusion limit differs by > 150% near threshold. For the lower threshold experiments, such as CRESSTII and the near future Super-CDMS experiments, m min D < 1 GeV, and the difference in the exclusions are 50%.
In order to judge whether these differences in the exclusion curves when using the obs versus iso VDFs are significant, we need to consider the uncertainty on the obs VDF profile. Given the current errors, the difference between the central values of the exclusion curves is well within the VDF uncertainty. However, assuming the same central values of the VDFs but with errors reduced (by future astrophysical measurements) to third errors, we find that for LUX and PICO, this difference could be ∼ 5σ significant.
In a direct detection experiment, there are four main sources of systematic errors : (i) astrophysical uncertainties on the local DM density and the VDF (ii) detector response uncertainty, (iii) uncertainty of the nuclear form factors and (iv) uncertainty on the detector background. In this work, we have discussed the uncertainty on the local DM density and VDF and argued that the VDF uncertainty is the dominant astrophysical unknown. To assess the impact of other uncertainties on the DM exclusion bound requires a careful understanding of the detector and is beyond the scope of this work.
To get an idea of the relative importance of the astrophysical vs detector uncertainties, we use the published uncertainties on the LUX exclusion [54] as a benchmark. An examination of their bounds indicate ∼ 50% detector related uncertainties, at all candidate DM masses. Note, that they do not include astrophysical uncertainties in their exclusion limits. In contrast, expected uncertainties in the mean DM exclusion curve due to third errors on the obs VDFs are ∼ 30% and, hence, it is expected to be subdominant to detector systematics once precision astrophysics results are used to determine the DM VDF.
Given that the application of the central obs VDF results in a systematic deviation of up to 200% in the mean DM exclusion curve, along with an estimate of the combined expected uncertainty from astrophysics as well the currently known experimental systematics of ∼ 60%, it is clearly important to use the best available observationally determined VDF when presenting the results of DM direct detection experiments.

DISCUSSION AND CONCLUSIONS
We have constructed an observationally driven determination of the local DM density and velocity distribution and used this to interpret the null results of DM direct detection experiments. Milky Way astrophysical data is poised for unprecedented precision measurements with the release of GAIA data and can be used to get precise estimates of the MW (and in particular the local) DM phase-space distribution. This approach goes beyond the simplistic, and incorrect, isothermal VDF in the so-called SHM of the MW. Previous attempts to go beyond the SHM have mainly relied on simulations of a MW-like DM halo. These simulations, although state-ofthe-art, are not guaranteed to describe exactly our MW. Despite the promise shown by these simulation results, it is prudent to use observationally inferred DM phasespace distributions. Our work is unique among other attempts in modelling the local DM phase-space in that it goes beyond just local kinematic data (prone to large DM-VM degeneracies) and uses rotation curve data up to ∼ 200 kpc for a full Bayesian reconstruction of the MW mass model and its corresponding self-consistent local DM density and VDF.
The mean observational VDF that we have determined differs from the SHM isothermal VDF in that it has a lower escape velocity and is significantly non-Maxwellian especially at the tails. Current RC data has large error bars for the majority of radial bins and results in a large current error band around the mean VDF. However, even with this large error bar, the obs VDF clearly differs from the isothermal VDF, in particular, at the high velocity tail. A similar determination of the DM VDF from the GAIA astrometric observations [73] are expected to reduce the error on VDF significantly, thus potentially differentiating the obs VDF from the isothermal VDF over most of the DM velocity range.
Low mass (∼ few GeV) DM has received considerable interest due to claims of a long-standing detection of an annual modulation signal by DAMA [60], as well as claims of excess events seen in other direct detection experiments such as CDMS-II [61], COGENT [62] and CRESST-II [63]. In addition, the observed excess of gamma rays from the galactic center could also be explained by the annihilation of a low mass DM species [64,65]. Although the DM origin of these anomalies is far from certain, these results are indicative of the need to precisely interpret the results of direct detection experiments for low DM masses.
In this work, we have shown that the difference between the observationally determined VDFs and the conventionally used isothermal VDF can yield very different interpretations of direct detection experiments for low candidate DM masses. Using the right DM VDF becomes especially pertinent when using VDFs inferred from future measurements from the GAIA telescope (where for simplicity we assume in this work that the mean obs VDF remains the same whereas the uncertainty around the mean reduces). For example, in such a scenario, for DM experiments like LUX or PICO, the DM exclusion limit using the obs DM VDF is expected to deviate by up to 5σ from the limit inferred from the SHM VDF at the detector threshold DM mass sensitivity. For future low threshold experiments, like SuperCDMS, accurate knowledge of the shape of the detector response is crucial to compute the impact of using the observationally determined DM VDF.
We emphasize that it is imperative that DM experiments use the best observationally estimated DM VDF when setting exclusion limits. For this purpose we have provided an accurate analytic fit to the obs VDF given in eq. (5) and a github link to the tables of the actual VDFs.
We provide the best fit parameters for the upper and lower edges of the error envelope for the VDFs shown in Fig. 1 in Tables II & III. Note that we have used the same parameterizations as eq. 5 although the fitting function is ideal for the mean VDF.   For experiments other than LUX, the detector response function was taken to have the parametric form, in terms of the error-function. Here E T is the detector threshold energy at which the efficiency drops to 50% and σ is the width of the efficiency near threshold. For the LUX experiment we used the efficiency curve given in Fig.  1 of ref. [54] with a hard low energy cut-off of 1.1 keV. We present for completeness in Table IV a listing of the detector properties for the experiments used in the results presented in the main text. We also list the references from which these specifications were extracted.