Astrophysical uncertainties on the local dark matter distribution and direct detection experiments

The differential event rate in Weakly Interacting Massive Particle (WIMP) direct detection experiments depends on the local dark matter density and velocity distribution. Accurate modelling of the local dark matter distribution is therefore required to obtain reliable constraints on the WIMP particle physics properties. Data analyses typically use a simple Standard Halo Model which might not be a good approximation to the real Milky Way (MW) halo. We review observational determinations of the local dark matter density, circular speed and escape speed and also studies of the local dark matter distribution in simulated MW-like galaxies. We discuss the effects of the uncertainties in these quantities on the energy spectrum and its time and direction dependence. Finally we conclude with an overview of various methods for handling these astrophysical uncertainties.


I. INTRODUCTION
Direct detection experiments aim to detect cold dark matter in the form of Weakly Interacting Massive Particles (WIMPs) via their elastic scattering off target nuclei [1]. The signals expected in these experiments, and hence the resulting constraints on the WIMP mass and cross-sections, depend on the local WIMP density and velocity distribution. We first outline the expressions for the differential event rate, including its time and direction dependence, focussing on how it depends on these astrophysical quantities (Sec. II). In Secs. III and IV we review what observations and simulations tell us about the local dark matter distribution, before discussing the effect of astrophysical uncertainties on the signals expected in direct detection experiments in Sec. V. Finally, in Sec. VI, we discuss strategies for handling astrophysical uncertainties, including 'astrophysics independent' methods, parameterizing the WIMP speed distribution and building self-consistent models of the Milky Way (MW).

II. BACKGROUND
A. Differential event rate The differential event rate, per unit detector mass, as a function of recoil energy, E, and direction, Ω, under the standard assumption of elastic scattering due to contact interactions, is given by where ρ 0 is the local (i.e. at the Solar radius) dark matter density, m χ is the WIMP mass, µ χA the WIMP-nucleus reduced mass, σ SI, SD the spin independent/spin dependent WIMP-nucleus cross section at zero momentum transfer, F SI, SD (E) the spin independent/spin dependent form factor,r a unit vector in the direction of the recoiling nucleus and is the minimum WIMP speed that can cause a recoil of energy E. Finallyf lab (v min ,r) is the Radon transform of the lab frame WIMP velocity distribution, f lab (v), i.e. the integral of f lab (v) over the plane orthogonal to the directionr a distance v min from the origin. For the majority of direct detection experiments, which can measure the energy of the nuclear recoils but not their directions, Eq. (1) where g(v min ) is the velocity integral: Since f (v) is positive definite, g(v min ) must be a monotonically decreasing function of v min . The spin independent cross section can be written as where Z and A are the atomic and mass number of the nucleus and f p,n is the coupling between the WIMP and protons/neutrons. The cross section for scattering on a proton is where µ χp is the WIMP-proton reduced mass and in the case where f p = f n the total spin independent cross section can be written as For further details of the derivation of Eqs. (1) and (4) see e.g. Ref. [5]. The lab frame speed distribution f lab (v) is time dependent, due to the time dependent transformation from the Galactic rest frame to the lab frame: v → v + v e (t). The Earth's motion relative to the Galactic rest frame, v e (t), is made up of three components: the motion of the Local Standard of Rest (LSR) [151], v LSR = (0, v c , 0) [152] where v c is the local circular speed (see Sec. III B), the Sun's peculiar motion with respect to the LSR, v p ⊙ , (see Sec. III C) and the Earth's orbit about the Sun, v orb e [153]. The Earth's orbit leads to an annual modulation in the differential event rate [6]. The recoil rate is also strongly direction dependent. The recoils are tightly concentrated around the inverse of the direction of Solar motion, such that the recoil rate in the forward hemisphere is roughly an order of magnitude larger than that in the backwards one [7].
The interaction between the WIMP and the nucleus may also be inelastic (due to, for instance, the WIMP particle belonging to a multiplet of states) [8]. In this case the relationship between v min and E, and hence the dependence of the differential event rate on the velocity distribution, is modified. There has also recently been interest in using non-relativistic effective operators to consider a wider range of WIMP-nucleus interactions, and again the relationship between the velocity distribution and the event rate is changed [9,10].

B. Phase space distribution function
The steady-state phase space distribution, f (x, v), of a collection of collisionless particles is given by the solution of the collisionless Boltzmann equation: df /dt = 0 (see e.g. Ref. [11]). The standard halo model (SHM) is an isotropic, isothermal sphere with density profile ρ(r) ∝ r −2 . In this case the solution to the collisionless Boltzmann equation is a Maxwellian velocity distribution where N is a normalization constant and σ v is the one-dimensional velocity dispersion, which is related to the circular speed, v c , (see Sec. III B) by σ v = v c / √ 2. Formally this distribution extends to infinity, therefore it is usually truncated by hand at the escape speed, v esc (Sec. III D).
For spherical isotropic systems the Eddington equation (see e.g. Ref. [11]) provides a unique relationship between the density profile and the distribution function. However for triaxial and anisotropic systems there is no unique relationship between the density profile and the velocity distribution. For an anisotropic system the phase space distribution function is a function of both energy and the magnitude of the angular momentum. In this case solutions can be found using specific ansatzes for the form of the anisotropy [12][13][14][15]. Another approach, used in Ref. [16] for triaxial systems, is to use the Jeans equations (which are derived by taking moments of the collisionless Boltzmnn equation) to find the components of the velocity dispersion and approximate the velocity distribution as a multi-variate gaussian.
The formation of dark matter halos is hierarchical, rather than monolithic, and therefore additional non-virialised contributions to the phase space distribution function, such as tidal streams, are possible. We discuss these features in Sec. IV.

A. Density
In the context of WIMP direct detection experiments, the standard value of the local dark matter density is ρ 0 = 0.3 GeV cm −3 (or equivalently 0.008M ⊙ pc −3 ). Observational probes of the dark matter density at the Solar radius can be divided into two classes: local and global. Local measures use the vertical motion of tracer stars in the Solar neighbourhood, while global measures use an ensemble of data sets to constrain a parameterized model of the MW. The local dark matter density could in principle be larger than the global (spherically averaged) value if the MW halo is flattened or there is a dark disc. For a detailed review see Ref. [17].
There is a long history of local measurements of the local dark matter density, dating back to Kapteyn [18], Jeans [19] and Oort [20] in the 1920s and 30s. In recent years there have been a number of studies, using different data sets and different methodologies (and hence different assumptions). These studies agree with each other within their quoted uncertainties, finding best-fit values in the range (0.22 − 0.33) GeV cm −3 [21][22][23][24][25][26]. The errors depend on the method used and lie, roughly, in the range (0.05 − 0.5) GeV cm −3 . See Fig. 2, Table 4 and Sec. 5.2 of Ref. [17] for further details. Recent determinations of the global dark matter density typically lie in the range (0.2−0.6) GeV cm −3 [27][28][29][30][31][32][33][34]. The size of the statistical errors (which in some cases are as small as 10%) depend on the assumptions made. See Table 4 and Sec. 5.3 of Ref. [17] for further details.
As discussed in Ref. [17], significant improvement in the accuracy (and precision) of determinations of the local dark matter density is expected in the near future using data from the Gaia satellite.

B. Circular speed
The standard value of the local circular speed is v c = 220 km s −1 . As discussed in Sec. 5.2 of Ref. [2], the local circular speed can be measured via two different general methods: by measuring the Sun's velocity with respect to an object assumed to be at rest with respect to the Galactic centre or by directly measuring the local radial force (via the Oort constants or observations of tidal streams).
The measured proper motion of Sgr A ⋆ [35], combined with a Solar radius of R 0 ≈ 8 kpc, gives a total Solar velocity of v φ,⊙ = 242 km s −1 . The φ component of the Sun's motion with respect to the Rotational Standard of Rest (RSR) then has to be subtracted to find v c . Ref. [2] measures the Milky Way rotation curve using stellar tracers. They find v c = (218 ± 6) km s −1 and v φ,⊙ = 242 +10 −3 km s −1 , with a value for the φ component of the Sun's velocity relative to the RSR of ∼ 24 km s −1 which is somewhat larger than the φ component of the Sun's velocity with respect to the LSR [36]. Ref. [37] finds v c = 220 ± 18 km s −1 from the orbit of the GD-1 stellar stream. A larger value of v c = (254 ± 16) km s −1 was claimed using the kinematics of masers in the Galactic disk [38], however subsequent reanalyses, taking into account the non-random orbital phases of the masers, have found a lower central value and significantly larger uncertainties [39,40].

D. Escape speed
The escape speed is the speed above which stars and particles are not gravitationally bound to the MW: v esc (r) = 2|Φ(r)| where Φ(r) is the potential. The local escape speed, v esc ≡ v esc (r = R 0 ), is estimated from the speeds of local high velocity halo stars.
The most recent estimate, v esc = 533 +54 −41 km s −1 at 90% confidence, comes from Ref. [41] using data from the Radial Velocity Experiment (RAVE). They define the escape speed as the minimum speed required to reach three virial radii [154] and parametrise the tail of the stellar velocity distribution as [42], which Ref. [41] found to be a good fit to cosmological simulations with k in the range 2.3 < k < 3.7. Note that, as emphasised by Ref. [43], the escape speed is correlated with other astrophysical parameters, in particular the circular speed, v c .

IV. SIMULATIONS
Since the local dark matter distribution can not be directly 'observed' (apart from by lab based dark matter detection experiments!) numerical simulations are a useful source of information about this quantity. For a detailed review of the state of the art (as of 2012) and future prospects see Ref. [44].
The resolution of simulations is many orders of magnitude larger than the scales probed by Earth based experiments (the Earth moves at ∼ 200 km s −1 ∼ 0.1 mpc yr −1 with respect to the Galactic rest frame). This raises the question of whether the dark matter distribution on these scales contains ultra-fine structure that can not be resolved by simulations. Vogelsberger and White [45] addressed this issue by devising a technique which probes the evolution of the fine-grained phase space distribution within a numerical simulation. They found that the typical density of the resulting streams is of order 10 −7 times the local halo density, and hence the local dark matter distribution consists of a huge number of overlapping streams, producing a close to smooth distribution. Similar conclusions have been reached by Schneider et al. [46] by considering the evolution of the first WIMP microhalos [47][48][49][50] that form.

A. N-body
In the 2000s various high resolution, dark matter-only simulations of MW-like halos in a cosmological context were carried out (e.g. Aquarius [51], GHALO [52] and Via Lactea [53]). The local velocity distributions of these halos were found to deviate systematically from a multi-variate Gaussian, having more low speed particles and a lower, flatter peak [54][55][56][57]. These distributions are well fit by a phenomenological form proposed by Mao et al. [58]: The parameter (v 0 /v esc ) varies with radius in a single halo and p lies in the range [0, 3] with significant scatter from halo-to-halo [58,59]. In additional to the dominant smooth component, the dark matter distribution contains stochastic components, in particular at high speed. At some locations there are tidal streams, giving narrow spikes in the velocity distribution [57]. There are also broad features, which vary from halo to halo, but do not vary with position within a single halo. These features, which have been dubbed 'debris flow', are believed to reflect the formation history of the halo [60,61].

B. Hydrodynamical
The contribution of baryons and dark matter to the potential in the Solar neighbourhood is comparable. Therefore the baryons are likely to have a non-negligible effect on the local dark matter distribution. Simulating baryonic physics requires formulating prescriptions for physical processes which occur on scales that are smaller than the resolution of the simulation. This is extremely challenging, however there have been significant advances in recent years in producing simulated galaxies which resemble the Milky Way.
Refs. [62][63][64] found velocity distributions in hydrodynamical simulations of MW-like galaxies which deviated significantly, in different ways, from a Maxwellian distribution. In early 2016 several studies of the local dark matter distribution in MW-like galaxies were published [65][66][67]. Borzognia et al. [65] selected galaxies from the EAGLE [68] and APOSTLE [69] projects with masses in the range (10 12 − 10 13 M ⊙ ) that satisfied observational constraints on the MW rotation curve and stellar mass. Kelso et al. [66] studied two galaxies with global properties similar to the MW from the MaGICC [70] simulations. Sloane et al. [67] studied four MW-like galaxies, with masses in the range 0.7 − 0.9 × 10 12 M ⊙ with a variety of different merger histories (i.e. quiescent or with one or more relatively recent major mergers). The resolution of their simulations is lower than those studied in Ref. [65,66].
In general the baryons deepen the potential and increase the average speed of the dark matter particles. Refs. [65,66] found that a Maxwellian speed distribution is a good fit to the speed distribution in their hydrodynamical simulations. This may be because adiabatic contraction makes the logarithmic slope of the density profile in the Solar neighbourhood closer to that of an isothermal sphere [66]. Even though the Maxwellian speed distribution is a good fit, there is still significant halo to halo variation in the peak speed/velocity dispersion. Ref. [67] also found that a Maxwellian speed distribution is a better fit to their hydrodynamical simulations than to their dark matter only simulations, however they still found a deficit of high speed particles relative to the Maxwellian. They also found that the Earth frame speed distribution in the hydrodynamical simulations contained features and the velocity distributions varies more between halos in the hydrodynamical simulations than in the dark matter only ones.
In some baryonic simulations a co-rotating dark matter disc is formed, due to late merging sub-halos being preferentially dragged towards the stellar disc by dynamical friction [62,71,72]. Ref. [73] argued that to be consistent with the observed properties of the MW's thick disc, the MW's merger history must be quiescent compared with typical ΛCDM merger histories, and hence the density of its dark disc must be relatively small. Most of the MW-like galaxies formed in recent hydrodynamical simulations do not have a significant dark disc [65,66].
Ref. [74] studied the ongoing disruption of the Sagittarius dwarf galaxy, and in particular the question of whether one of its tidal streams passes through the Solar neighbourhood (as suggested, in the context of direct detection, by Ref. [75]). They found that the dark matter in the stream is more extended, and not coaxial with, the stellar stream and can hence have a non-negligible contribution to the local dark matter distribution even if the stellar stream does not pass through the Solar neighbourhood. Ref. [76] studied the dynamical effects of a galactic bar on the dark matter distribution, finding that the DM density at the Solar radius varies depending on the Earth's location relative to the stellar bar and a quadrupole wake in the DM that lags the stellar bar forms.
In summary the simulations discussed above use different prescriptions for sub-grid physics as well as different criteria for selecting MW-like galaxies. Recent studies [65][66][67] suggest that the speed distribution does not deviate hugely from a Maxwellian, however there may still be non-negligible deviations and a firm consensus about the exact form of the local dark matter velocity distribution has still to emerge. In particular there may be significant variation from halo to halo.

V. CONSEQUENCES
The normalization of the differential event rate is proportional to the product of the density and cross-section, so the accuracy of constraints on the cross-section is limited by the accuracy with which the local density is known. The differential event rate is proportional to an integral over the velocity distribution, Eq. (4), and hence it depends fairly weakly on the detailed shape of the velocity distribution [77,78]. Therefore the resulting uncertainty in exclusion limits [79][80][81][82] and determinations of the WIMP mass [83,84] is usually fairly small. The main exception to this is the case of light WIMPs. If the WIMP is much lighter than the target nuclei, or the experimental energy threshold is large, then the experiment is only sensitive to the high speed tail of the speed distribution. In this case the precise value of the escape speed, and also the shape of the speed distribution just below it, can have a significant effect on the energy spectrum, and hence exclusion limits or determinations of the WIMP parameters [80,85]. Astrophysical uncertainties also affect the level of the 'neutrino floor' [86], the cross-section below which neutrino backgrounds strongly limit the sensitivity of non-directional experiments [87].
The annual modulation arises from the small shift in the speed distribution in the lab frame which occurs due to the Earth's orbit. The properties of the annual modulation are therefore far more sensitive to the speed distribution than the time averaged differential event rate. The amplitude and phase of the modulation, and hence the regions of parameter space compatible with an observed signal, can vary significantly [88][89][90][91][92][93][94][95][96][97][98][99]. The detailed angular dependence of the directional event rate is sensitive to the exact shape of the velocity distribution, however the 'forwards-backwards' anisotropy is robust [100][101][102].
A sufficiently high density, high speed tidal stream would produce a step in the differential event rate [103], the position and height of which vary annually [104]. In a directional experiment a stream would produce recoils which are strongly peaked in the opposite direction [105], allowing the stream properties to be measured [106]. A dark disc with a sufficiently high density would increase the differential event rate at low energies and also change the phase and amplitude of the annual modulation signal [107].

VI. HOW TO HANDLE THE UNCERTAINTIES?
In light of the uncertainties in the local WIMP density and velocity distribution discussed in Sec. III and IV above, it is crucial to develop methods for analysing and comparing direct detection data that are independent of these uncertainties. Conventionally direct detection experimental results are presented in terms of constraints on the WIMP mass m χ and cross-sections, σ SI and σ SD , however this requires assumptions to be made about the local density and velocity distribution. The energy dependence of the differential event rate depends on both the WIMP mass and velocity distribution, so that it is impossible to measure the WIMP mass with a single experiment without making assumptions about the velocity distribution [108,109]. In Sec. VI A we discuss astrophysics independent methods, in Sec. VI B we explore parameterizing the WIMP speed distribution and finally, in Sec. VI C, we discuss using astrophysical data to constrain mass models of the MW.

A. Astrophysics independent methods
In principle, with data from two experiments using different targets, the WIMP mass can be be determined, without any assumptions about f (v) by taking moments of the energy spectra [108,109]. However the WIMP mass is systematically underestimated by this method if it is comparable with or heavier than the mass of the target nuclei.
Fox et al. [110] showed that an astrophysics independent comparison of different experiments can be made using the rescaled velocity integral,g(v min ) = (ρ 0 σ SI p /m χ )g(v min ), where g(v min ) is defined in Eq. (5). Different experiments probe different sections of the speed distribution, and it is possible to assess whether they overlap and whether or not the measured rates are consistent. Because the relationship between recoil energy E and v min depends on the (unknown) WIMP mass, the comparison has to be done for a selection of fixed values of the WIMP mass. Alternatively the comparison can be done in terms of the nuclear recoil momentum which, unlike v min , is independent of m χ [111].
The original method has been generalized to include annual modulation [112][113][114][115], experiments containing multiple isotopes [113], energy-dependent experimental efficiency [115] as well as 'non-standard' particle interactions [116][117][118][119][120]. This method is an extremely powerful tool for examining whether signals and limits from different experiments are consistent with each other. In particular it has been used to show that the excesses over backgrounds seen by CRESST [121], CoGeNT [122] and CDMS Si [123], and the annual modulations seen by DAMA/LIBRA [124] and CoGeNT [125] are not compatible with each other and also the exclusion limits from LUX [126], PandaX [127] and CDMSlite [128] in the absence of fine-tuned non-standard interactions (see e.g. Ref. [129]). This method has been developed so that it can be applied to non-binned data [130,131] and the degree of agreement between different experiments quantified [131][132][133][134]. This involves writing the rescaled velocity integral as a series of monotonically decreasing steps (which is physically equivalent to the velocity distribution being composed of a series of streams with different densities and velocities) and finding the optimum step properties [130][131][132]135]. It is also possible to make astrophysics independent comparisons of direct detection experiments and neutrino telescopes [135,136].

B. Parameterizing the speed distribution
If the underlying shape of the speed distribution is known, it is possible to make unbiased reconstructions of the WIMP particle physics parameters by marginalising over the parameters of the speed distribution (for instance the peak speed) [137,138]. However Peter [139] showed that erroneous assumptions about the shape of the local speed distribution lead to biased determinations of the WIMP particle physics parameters. She suggested the alternative approach of using an agnostic empirical parameterization of the speed distribution, and jointly constraining this, along with the WIMP particle parameters, using data from multiple experiments. A simple parameterization, as a series of bins with constant values, performed better than assuming an incorrect functional form, however the reconstructed parameter values were still biased [139]. Kavanagh and Green [140] showed that this was because with fixed speed bins, the width of the bins in energy space can be reduced, so as to give a better fit to the measured energy spectrum, by artificially reducing the WIMP mass. They proposed parameterising the momentum distribution instead to avoid this problem. They subsequently found [141] that parameterizing the natural logarithm of the speed distribution in terms of Legendre polynomials allows an accurate, unbiased reconstruction of the WIMP mass to be made, even for speed distributions containing a high density stream or dark disc. It is also possible to reconstruct the speed distribution itself, however the uncertainties are quite large [141].
It is not possible to make unbiased measurements of the cross-sections using direct detection experiments alone, however, as we do not know what fraction of the WIMP have speeds which are too low to cause recoils with energies above the experimental threshold energy. Ref. [142] showed that this problem can be avoided by combining data from direct detection experiments with neutrino telescope data. This is because it is the low speed WIMPs, that direct detection experiments are not sensitive to, that are captured in the Sun and then annihilate. So by combining these data sets the entire speed distribution is probed. Unbiased measurements of both the WIMP mass and cross sections can then be made using either the binned or polynomial speed parameterisations. The polynomial parameterisation allows for a wider range of possible speed distributions, and hence leads to larger uncertainties in the particle physics parameters. These methods have recently been extended to parameterizing the full 3d velocity distribution as probed by directional experiments [143], see also Refs. [144,145]. See Ref. [146] for a detailed study of how the WIMP astrophysics and particle physics can be constrained using multiple direct detection experiments.

C. Modelling the Milky Way
In early work in this area Strigari and Trotta [147] constructed a mass model of the MW and used stellar kinematics and (simulated) direct detection data to simultaneously fit the WIMP particle physics parameters and the free parameters of the speed distribution, which they assumed to be an isotropic Maxwellian. Catena and Ullio [148] built mass models of the MW, including its luminous components and a range of possibilities for the density profile of the DM halo. They used a large ensemble of astronomical data sets to constrain the parameters of the mass model, and then used the Eddington formalism to self-consistently calculate the speed distribution (see also Ref. [149]). Ref. [15] extended this approach to include anisotropy of the velocity distribution. The resulting mean velocity distribution has a high-speed tail that is enhanced relative to the isotropic case and is similar to the empirical form, Eq. (10), found by Mao et al. to give a good fit to N-body simulations [58]. This approach can be extended to include information about the shape of the density profile from gamma-rays produced by WIMP annihilation in the MW halo [150].

VII. SUMMARY
The differential event rate in WIMP direct detection experiments depends on the local dark matter density and velocity distribution. Experimental data analyses usually assume a simple Standard Halo Model, an isothermal sphere with isotropic Maxwellian speed distribution. However this model may not be an accurate description of the actual dark matter distribution in the Milky Way. We first reviewed observational determinations of the relevant quantities: the local density, circular speed and escape speed. Typically the statistical errors can be fairly small, ∼ 10%, however the systematic errors may be substantially larger. We then turned our attention to numerical simulations. There has recently been significant progress in carrying out hydrodynamical simulations including baryons and it appears that the speed distribution is closer to a Maxwellian distribution than previously thought. However there may still be non-negligible deviations and significant variation from halo to halo We then discussed the consequences of these astrophysical uncertainties. The energy spectrum depends on an integral over the WIMP speed distribution, therefore exclusion limits are usually relatively insensitive to its exact shape. However the annual modulation and detailed direction dependence are far more dependent on the detailed form of the velocity distribution. Finally we concluded by discussing methods for handling astrophysical uncertainties. Astrophysics independent methods allow model independent comparisons of data sets from different experiments. Alternatively with an appropriate parameterization of the speed distribution it will be possible, using data from multiple experiments, to make an unbiased measurement of the WIMP mass without knowledge of, or assumptions about, the speed distribution. Finally another approach is to build a mass model of the Milky Way, constrain its parameters using observational data and calculate the speed distribution self-consistently. In the event of convincing signals being observed in direct detection experiments all three of these approaches will likely be useful in establishing their consistency and making reliable measurements of the WIMP particle physics parameters.