# Dark Matter Direct Detection with Non-Maxwellian Velocity Structure

Michael Kuhlena,b, Neal Weinerc, Jürg Diemandd, Piero Madaue, Ben Moored, Doug
Potterd, Joachim Stadeld, Marcel Zempf
a School of Natural Science, Institute for Advanced Study, Princeton, NJ 08540
b Theoretical Astrophysics Center, University of California Berkeley,
Berkeley, CA 94720
c Center for Cosmology and Particle Physics, Department of Physics, New York
University, New York, NY 10003
d Institute for Theoretical Physics, University of Zurich, 8057 Zurich,
Switzerland
e Department of Astronomy & Astrophysics, University of California Santa Cruz,
Santa Cruz, CA 95064
f Department of Astronomy, University of Michigan, Ann Arbor, Michigan 48109
E-mail  mqk<EMAIL_ADDRESS>
###### Abstract:

The velocity distribution function of dark matter particles is expected to
show significant departures from a Maxwell-Boltzmann distribution. This can
have profound effects on the predicted dark matter - nucleon scattering rates
in direct detection experiments, especially for dark matter models in which
the scattering is sensitive to the high velocity tail of the distribution,
such as inelastic dark matter (iDM) or light (few GeV) dark matter (LDM), and
for experiments that require high energy recoil events, such as many
directionally sensitive experiments. Here we determine the velocity
distribution functions from two of the highest resolution numerical
simulations of Galactic dark matter structure (Via Lactea II and GHALO), and
study the effects for these scenarios. For directional detection, we find that
the observed departures from Maxwell-Boltzmann increase the contrast of the
signal and change the typical direction of incoming DM particles. For iDM, the
expected signals at direct detection experiments are changed dramatically: the
annual modulation can be enhanced by more than a factor two, and the relative
rates of DAMA compared to CDMS can change by an order of magnitude, while
those compared to CRESST can change by a factor of two. The spectrum of the
signal can also change dramatically, with many features arising due to
substructure. For LDM the spectral effects are smaller, but changes do arise
that improve the compatibility with existing experiments. We find that the
phase of the modulation can depend upon energy, which would help discriminate
against background should it be found.

dark matter, direct detection, numerical simulations

##  1 Introduction

Direct detection experiments aim to detect the low energy nuclear recoil from
rare scattering events between dark matter (hereafter DM) particles, assumed
to be weakly interacting massive particles (WIMPs), and target nuclei. The
event rate and its energy spectrum depend on the properties of the DM
distribution at Earth’s location, about 8.5 kpc from the Galactic Center [1].
Typical calculations of the scattering rate assume a “standard halo model”
(SHM) consisting of a local DM density of $0.3\pm 0.1$ GeV cm-3 and a Maxwell-
Boltzmann (MB) velocity distribution function with a (three-dimensional)
dispersion of $270\pm 70$ km/s [2] truncated at an escape speed of $\sim 550$
km/s. Recent numerical simulations of the formation of Galactic-scale DM halos
have reached the necessary resolution to directly test these assumptions.

At the same time, absent a positive signal, a set of uniform halo assumptions
allows a simple means to compare different experiments. However, in light of
the recent results from DAMA/LIBRA [3], confirming earlier results from
DAMA/NaI [4], it is important to consider halo model uncertainties when
discussing exclusion limits from experiments with different targets, or energy
ranges. This is particularly important because proposals such as light dark
matter (LDM) [5, 6] and inelastic dark matter (iDM) [7, 8], which aim to
reconcile DAMA with null results from other experiments, sample the high
velocity component of the WIMPs preferentially, and it is especially here that
the Maxwell-Boltzmann distribution is expected to break down.

In the hierarchical structure formation paradigm of standard cold dark matter
cosmology the DM halo of a typical galaxy is built up through the merger of
many individual gravitationally bound progenitor halos, which themselves were
assembled in a hierarchical fashion. High resolution cosmological simulations,
such as Via Lactea [9, 10], GHALO [11], and Aquarius [12], have shown that
this merging process is in general incomplete, with the dense cores of many of
the merging halos surviving as subhalos orbiting within their respective host
halos. During pericenter passages tidal forces can strip off a large fraction
of a subhalo’s material, but the resulting cold tidal streams can readily be
identified as velocity space substructure. The resulting DM halos are not
perfectly phase mixed, and the assumption of a smooth halo is in general not a
good one, neither in configuration space nor in velocity space [13].

At $8.5$ kpc the Sun is located quite close to the Galactic Center, at least
when compared to the overall extent of the Milky Way’s DM halo ($r_{\rm
vir}\sim 200-300$ kpc). This central region is notoriously difficult for
cosmological numerical simulations to resolve, as very high particle numbers
and very short time steps are required to avoid the so-called “over-merging
problem” [14], which has until recently resulted in an artifically smooth
central halo devoid of any substructure. With the advent of
$\mathcal{O}(10^{9})$ particle simulations at the Galactic scale this problem
finally seems to have been overcome, with hundreds of subhalos identified at
$\lesssim 20$ kpc. Nevertheless the local phase space structure is far from
completely resolved, and likely never will be through direct numerical
simulation. Any estimation of the importance of local density or velocity
substructure based on cosmological simulation must thus rely on extrapolations
over many orders of magnitude below its resolution limit.

Local density variations due to the clumpiness of the DM halo are unlikely to
significantly affect the direct detection scattering rate. Based on the
Aquarius Project suite of numerical simulations, Vogelsberger et al. (2008)
[15] report that at more than 99.9% confidence the DM density at the Sun’s
location differs by less than 15% from the average over a constant density
ellipsoidal shell. Extrapolating from their numerical convergence study they
estimate a probability of $10^{-4}$ of the Sun residing in a bound subhalo of
any mass. Analytical work by Kamionkowski & Koushiappas (2008) [16] predicts a
positively skewed density distribution with local densities as low as one
tenth the mean value, but probably not much less than half.

The situation for velocity substructure is less clear. It is well established
that numerically simulated dark matter halos exhibit significant velocity
anisotropy and global departures from a Maxwell-Boltzmann distribution [17,
18, 19, 20]. The implications for direct detection experiments of these global
departures from the standard Maxwellian model have previously been
investigated in the context of a standard WIMP model [21, 22, 23, 24] and for
inelastic dark matter [25, 26], and they were found to result in appreciable
differences (factor of a few) in the total event rates and the annual
modulation signal. Most recently Vogelsberger et al. (2008) reported
significant structure (“wiggles”) in the velocity distribution function
measured in their high-resolution Aquarius simulations, which they attributed
to events in the halo’s mass assembly history. Their analysis concluded that
velocity substructure due to bound subhalos or unbound tidal streams, however,
does not influence the detector signals, since it makes up a highly sub-
dominant mass fraction locally.

The aim of this paper is take a closer look at this velocity space
substructure and to examine its impact on the direct detection signal for
models that are particularly sensitive to the high velocity tail, such as LDM
or iDM. In contrast to Vogelsberger et al. (2008), we find that both global
and local departures from the best-fit Maxwell-Boltzmann distribution can
significantly affect the total event rate, the annual modulation, and the
recoil energy spectrum. Parameter exclusion limits derived using a standard MB
halo model are likely to be overly restrictive.

This paper is organized as follows: in Section 2 we present velocity
distribution functions derived from the high-resolution numerical simulations
Via Lactea II and GHALO. In Section 3 we look at the implications of high
velocity substructure for direct detection experiments with directional
sensitivity. In Section 4 we consider iDM and LDM models and show how the
observed local and global departures from the MB model affect scattering event
rates and recoil spectra at several ongoing direct detection experiments, and
how this modifies parameter exclusion limits. A summary and discussion of our
results can be found in Section 5.

##  2 Results from Numerical Simulations

The nuclear scattering event rate depends on the size of the detector, the
type of target material, the scattering cross section, and the number density
and velocity distribution of the impinging DM particles. We defer calculations
of the expected event rate for various experimental setups and types of DM
models to section 4, and focus in this section on the particle velocity
distributions, which we determine directly from numerical simulations.

###  2.1 The Via Lactea and GHALO simulations

Our analysis is based on two of the currently highest resolution numerical
simulations of Galactic DM structure: Via Lactea II (VL2) [10] and GHALO [11].
Both are cosmological cold DM N-body simulations that follow the hierarchical
growth and evolution of a Milky-Way-scale halo and its substructure from
initial conditions in the linear regime ($z=104$ for VL2, $z=58$ for GHALO)
down to the present epoch. For details about the setup of the simulations we
refer the reader to the above references. The VL2 host halo is resolved with
$\sim 400$ million particles of mass $m_{p}=4,100\,\rm M_{\odot}$ within its
virial radius111$r_{\rm vir}$ is defined as the radius enclosing a density of
$\Delta_{\rm vir}=389$ times the background density. of $r_{\rm vir}=309$ kpc
and has a mass of $M_{\rm halo}=1.7\times 10^{12}\,\rm M_{\odot}$ and peak
circular velocity $V_{\rm max}=201.3$ km/s. The GHALO host is somewhat less
massive, $M_{\rm halo}=1.1\times 10^{12}\,\rm M_{\odot}$ and $V_{\rm
max}=152.7$ km/s, but even more highly resolved, with 1.1 billion particles of
mass $m_{p}=1,000\,\rm M_{\odot}$ within its $r_{\rm vir}=267$ kpc. For
reference we show the circular velocity of the two halos in Fig.1. In order to
facilitate a more direct comparison between the two halos, we have also scaled
GHALO to match VL2’s $V_{\rm max}$ by multiplying the simulation’s length and
velocity units by a factor $f=V_{\rm max}({\rm VL2})/V_{\rm max}({\rm
GHALO})=1.32$, and the mass unit by $f^{3}$. We refer to this model as GHALOs.
The circular velocity of these three halos at 8.5 kpc is 158.1, 121.7, and
148.9 km/s for VL2, GHALO, and GHALOs, respectively.

Figure 1: Circular velocity profiles of the VL2, GHALO, and GHALOs host halos.

###  2.2 Velocity Modulus Distributions

The DM-nucleon scattering event rate is directly proportional to

$g(v_{\rm min})=\int_{v_{\rm min}}^{\infty}\frac{f(v)}{v}dv,$ (1)

where $f(v)$ is the DM velocity distribution function in the Earth’s rest
frame and $v_{\rm min}(E_{R})$ is the minimum velocity that can result in a
scattering with a given nuclear recoil energy $E_{R}$. For a target with
nuclear mass $m_{N}$ and a WIMP/nucleon reduced mass
$\mu=m_{N}m_{\chi}/(m_{N}+m_{\chi})$, $v_{\rm min}(E_{R})$ is given by

$\left(\frac{v_{\rm
min}}{c}\right)^{2}=\frac{1}{2}\frac{m_{N}E_{R}}{\mu^{2}}\left(1+\frac{\mu}{m_{N}E_{R}}\delta\right)^{2}.$
(2)

The $\delta$ refers to the possible mass splitting between the incoming and
outgoing DM particle, which would be 0 for standard and light DM and
$\mathcal{O}$(100 keV) for inelastic DM.

We determine $f(v)$ in the halo rest frame directly from the particle
velocities in our numerical simulations, and in the Earth’s rest frame by
first applying a Galilean velocity boost by $v_{\oplus}(t)$. The Earth’s
velocity with respect to the Galactic center is the sum of the local standard
of rest (LSR) circular velocity around the Galactic center, the Sun’s peculiar
motion with respect to the LSR, and the Earth’s orbital velocity with respect
to the Sun,

$\vec{v}_{\oplus}(t)=\vec{v}_{\rm LSR}+\vec{v}_{\rm pec}+\vec{v}_{\rm
orbit}(t).$ (3)

We follow the prescription given in Chang et al. (2008) [27] and set
$\vec{v}_{\rm LSR}=(0,220,0)$ km/s, $\vec{v}_{\rm pec}=(10.00,5.23,7.17)$ km/s
[28], and $\vec{v}_{\rm orbit}(t)$ as specified in reference [29]. The
velocities are given in the conventional $(U,V,W)$ coordinate system where $U$
refers to motion radially inwards towards the Galactic center, $V$ in the
direction of Galactic rotation, and $W$ vertically upwards out of the plane of
the disk. We associate these three velocity coordinates with the
$(v_{r},v_{\theta},v_{\phi})$ coordinates of the simulation particles.

The Earth’s orbital motion around the Sun results in the well-known annual
modulation of the scattering rate, which the DAMA collaboration claims to have
detected at very high statistical significance [3]. For the SHM the peak of
this modulation occurs around June 2nd, when the Earth’s relative motion with
respect to the Galactic DM halo is maximized.

Figure 2: Velocity distribution functions: the left panels are in the host
halo’s restframe, the right panels in the restframe of the Earth on June 2nd,
the peak of the Earth’s velocity relative to Galactic DM halo. The solid red
line is the distribution for all particles in a 1 kpc wide shell centered at
8.5 kpc, the light and dark green shaded regions denote the 68% scatter around
the median and the minimum and maximum values over the 100 sample spheres, and
the dotted line represents the best-fitting Maxwell-Boltzmann distribution.

We have measured the DM velocity distribution from all particles in a 1 kpc
wide spherical shell (8 kpc $<r<$ 9 kpc), containing 2.1, 5.4, and 3.6 million
particles in VL2, GHALO, and GHALOs, respectively. The large particle numbers
in these measurements result in a very small statistical uncertainty, but fail
to capture any local variations. To address this we have also determined
$f(v)$ from the particles in 100 randomly distributed sample spheres centered
at 8.5 kpc. These sample spheres have radii of 1.5 kpc for VL2 and 1 kpc for
GHALO and GHALOs, and contain a median of 31,281, 21,740, and 14,437 particles
in the three simulations.222Tables of $g(v_{\rm min})$ determined from the
spherical shell and the 100 sample spheres, and tracing the annual modulation
over 12 evenly spaced output times, are available for download at
http://astro.berkeley.edu/$\sim$mqk/dmdd/.

The resulting distributions, both in the halo rest frame and translated into
Earth’s rest frame, are shown in Fig. 2. The shell averaged distribution is
plotted with a solid line, while the light and dark green shaded regions
indicate the 68% scatter around the median and the absolute minimum and
maximum values of the distribution over the 100 sample spheres. For comparison
we have also overplotted the best-fitting Maxwell-Boltzmann (hereafter MB)
distributions, with 1D velocity dispersion of $\sigma_{\rm 1D}=$ 130, 100, and
130 km/s. These clearly underpredict both the low and high velocity tails of
the actual distribution. This is not a new result and has previously been
found in cosmological numerical simulations [17, 18, 19, 15]. Actually there
is no reason to assume that a self-gravitating, dissipationless system would
have a locally Maxwellian velocity distribution, and in fact it has been shown
that self-consistent, stable models of cuspy DM structures require just such
non-Gaussianity [30, 31].

In addition to its overall non-Maxwellian nature, we notice several broad
bumps present in both the shell averaged and, at very similar speeds, in the
sub-sample $f(v)$. Similar features were reported by Vogelsberger et al.
(2008) [15] for the host halos of their completely independent Aquarius
simulations, and thus appear to be robust predictions of hierarchically formed
collisionless objects. Vogelsberger et al. also showed that the broad bumps
are independent of location and persistent in time and hence reflect the
detailed assembly history of the host halo, rather than individual streams or
subhalos. The extrema of the sub-sample distributions, however, exhibit
numerous distinctive narrow spikes at certain velocities, and these are due to
just such discrete structures. Note that although only a small fraction of
sample spheres exhibits such spikes, they are clearly present in some spheres
in all three simulations. The Galilean transform into the Earth’s rest frame
washes out most of the broad bumps, but the spikes remain visible, especially
in the high velocity tails, where they can profoundly affect the scattering
rates for inelastic and light DM models (see Section 4).

In order to assess the dependence of these features on the sample sphere size,
we also considered for the VL2 simulation sphere radii of 1 and 2 kpc,
containing a median of 9,200 and 74,398 particles, respectively. The coarse
features of the distributions persist, but of course the prevalence of the
spikes increases with sample sphere size, as more of the substructure is
probed. It is difficult to assess with our simulations the true likelihood of
significant local velocity substructure, as it depends on the abundance and
physical extent of subhalos and tidal streams many orders of magnitude below
the length scales that we can accurately resolve. Higher resolution numerical
simulations, as well as analytical models [32], perhaps in conjunction with
simulations [33], will be necessary to settle this question.

###  2.3 The effects of neglected baryonic physics

The values of $\sigma_{\rm 1D}$ we report here may appear surprisingly low to
a reader familiar with the standard isothermal MB halo assumption of $\langle
v^{2}\rangle=3\,\sigma_{\rm 1D}^{2}=3/2\,v_{0}^{2}$, where $v_{0}$, the peak
of the MB distribution, i.e. the most probable speed, is assumed to be equal
to the rotation velocity of the Sun around the galaxy, $v_{0}\simeq 220$ km/s.
In this standard model $\sigma_{\rm 1D}$ would be 156 km/s, considerably
higher than our values of 130 km/s and 100 km/s, respectively. In fact, the
local circular velocity $v_{c}$ and velocity dispersion $\sigma$ are only
indirectly related and not necessarily equal. While $v_{c}$ is set by the
local radial gradient of the potential, $\sigma$ depends on the shape of the
potential at exterior radii. For a non-rotating spherical system the relation
between $v_{c}(r)$ and the radial velocity dispersion $\sigma_{r}(r)$ is given
by

$v_{c}^{2}=-\sigma_{r}^{2}\left(\frac{d\ln\rho}{d\ln
r}+\frac{d\ln\sigma_{r}^{2}}{d\ln r}+2\beta\right),$ (4)

where $\frac{d\ln\rho}{d\ln r}\equiv\gamma(r)$ is the logarithmic slope of the
density profile and $\beta(r)\equiv 1-\sigma_{\theta}^{2}/\sigma_{r}^{2}$ is
the velocity anisotropy. In a singular isothermal sphere ($\gamma=-2$,
$\frac{d\ln\sigma_{r}^{2}}{d\ln r}=0$, $\beta=0$) we have $v_{c}=v_{0}$, but
for an NFW profile $v_{c}/v_{0}\approx 0.88$ at $r=r_{s}/2$. In the VL2 and
GHALO host halos the relation is $v_{c}/v_{0}=0.85$ and $0.86$, respectively.

A central baryonic condensation in the form of a Galactic disk and bulge will
deepen the central potential, raise the local circular velocity to $\sim 220$
km/s, and increase the velocity dispersion of DM particles at the Sun’s
location. Since the VL2 and GHALO simulations do not include baryonic physics,
it is not surprising that the values of $v_{c}$ and $v_{0}$ at 8.5 kpc are
lower than appropriate for our Milky Way galaxy. The main focus of our work
here is to investigate the effects of global and local variations from the MB
assumption, and therefore we compare our results to the best-fitting MB model
($v_{0}=184$ km/s for VL2 and GHALOs and 141 km/s for GHALO) instead of the
standard MB halo model ($v_{0}=220$ km/s). In principle we could have scaled
just the velocities up to give $v_{0}=220$ km/s, but this would remove many of
the effects of substructure by pushing it above the escape velocity. Thus, we
use VL2 and GHALOs simulations for the parameter exclusions plots in Section
4, but it should be recognized, that the velocity dispersion is somewhat lower
than in conventional MB halo parameterizations.

###  2.4 Local Escape Speed

The escape speed from the Sun’s location in the Galactic halo is another
factor that can strongly affect scattering rates, especially for inelastic and
light DM models. The RAVE survey’s sample of high-velocity stars constrains
the Galactic escape velocity to lie between 498 and 608 km/s at 90%
confidence, with a median likelihood of 544 km/s [34]. This is in good
agreement with the highest halo rest frame speed of any particle in our 8.5
kpc spherical shells, namely 550 km/s in VL2 and 586 km/s in GHALOs. The lower
$V_{\rm max}$ of the GHALO host is reflected in a significantly lower escape
speed, only 433 km/s. The corresponding maximum speeds in the Earth rest frame
are 735-761, 773-802, and 634-660 km/s, where the range refers to the
modulation introduced by the Earth’s orbit around the Sun.

###  2.5 Radial and Tangential Distributions

Figure 3: Radial and tangential velocity distribution functions. The solid
lines show the shell average, the shaded range the 68% scatter around the
median, and the dotted curve shows the best-fit MB model.

For comparison with past work [25, 24] we also present separately the
distribution functions of the radial $v_{r}$ and tangential
$v_{t}=\sqrt{v_{\theta}^{2}+v_{\phi}^{2}}$ velocity components in Fig. 3. We
fit these distributions with functions like the ones in Fairbairn & Schwetz
(2009) [24], except that we don’t normalize the velocities by the square root
of the gravitational potential at the particles’ location and instead provide
separate fits for each simulation. We also include an estimate of the variance
due to local velocity substructure by separately fitting the distribution in
each of the sample spheres. The fitting functions are

$\displaystyle f(v_{r})$ $\displaystyle=$
$\displaystyle\frac{1}{N_{r}}\exp\left[-\left(\frac{v_{r}^{2}}{\bar{v}_{r}^{2}}\right)^{\alpha_{r}}\right]$
(5) $\displaystyle f(v_{t})$ $\displaystyle=$
$\displaystyle\frac{v_{t}}{N_{t}}\exp\left[-\left(\frac{v_{t}^{2}}{\bar{v}_{t}^{2}}\right)^{\alpha_{t}}\right],$
(6)

and the parameters of these fits are listed in Table 1. The normalizations
$N_{r}$ and $N_{t}$ can readily be obtained numerically for a given set of
parameters by ensuring that the distributions integrate to unity.

|  | radial | tangential
---|---|---|---
|  | shell | median | $16^{\rm th}$ | $84^{\rm th}$ | shell | median | $16^{\rm th}$ | $84^{\rm th}$
VL2 |  $\bar{v}_{r,t}$ [km/s] | 202.4 | 199.9 | 185.5 | 212.7 | 128.9 | 135.1 | 124.2 | 148.9
$\alpha_{r,t}$ | 0.934 | 0.941 | 0.877 | 0.985 | 0.642 | 0.657 | 0.638 | 0.674
GHALO |  $\bar{v}_{r,t}$ [km/s] | 167.9 | 163.6 | 156.4 | 173.0 | 103.1 | 114.3 | 93.21 | 137.0
$\alpha_{r,t}$ | 1.12 | 1.11 | 1.02 | 1.20 | 0.685 | 0.719 | 0.666 | 0.819
GHALOs |  $\bar{v}_{r,t}$ [km/s] | 217.9 | 213.8 | 202.3 | 226.6 | 138.2 | 162.2 | 125.1 | 183.1
$\alpha_{r,t}$ | 1.11 | 1.11 | 1.01 | 1.18 | 0.687 | 0.759 | 0.664 | 0.842
Table 1: Radial and tangential velocity distribution fit parameters (see Eqs.
5 and 6). The columns labeled “shell” refer to the fit for all particles in
the 1 kpc wide shell centered on 8.5 kpc, whereas the following three columns
give the median as well as the $16^{\rm th}$ and $84^{\rm th}$ percentile of
the distribution of fits over the 100 sample spheres (see text for more
detail).

###  2.6 Modulation Amplitude and Peak Day

The annual modulation of the scattering rate $g(v_{\rm min})$ (Eq. 1) grows
with $v_{\rm min}$, since a reduction in the number of particles able to
scatter makes the summer-to-winter difference relatively more important. At
sufficiently high velocities the modulation amplitude can even reach unity,
when during the winter there simply aren’t any particles with velocities above
$v_{\rm min}$.

The expected modulation amplitude for a given experiment is then determined by
the relation between the measured recoil energy $E_{R}$ and $v_{\rm min}$, see
Eq. 2. There are important qualitative difference between standard DM models
($\delta=0$) and inelastic models ($\delta\sim 100$ keV). Inelastic DM models
require a higher $v_{\rm min}$ for a given $E_{R}$, resulting in a lower total
scattering rate and a more pronounced modulation. Furthermore, while for
standard DM $v_{\rm min}$ grows with $E_{R}$, in inelastic models $v_{\rm
min}$ typically falls with $E_{R}$ (for $\delta>E_{R}m_{N}/\mu$).

Figure 4: The amplitude (left panels) and peak day (right) as a function of
$v_{\rm min}$. The solid red line indicates the shell averaged quantities, and
the dotted line the best-fit MB distribution. The light and green shaded
regions cover the central 68% region around the median and the minimum and
maximum values of the distribution over the 100 sample spheres. The thin black
line shows the behaviour of one example sample sphere.

In the left panels of Fig. 4 we show the $v_{\rm min}$ dependence of the
annual modulation amplitude, $({\rm max}(g(v_{\rm min}))-{\rm min})/({\rm
max}+{\rm min})$ as measured in our simulations. We focus on the high $v_{\rm
min}$ region that is relevant for inelastic and light dark matter models, as
well as directional detection experiments. In VL2 and GHALOs the shell
averaged modulation amplitude (solid red line) rises from about 20% at $v_{\rm
min}=400$ km/s to unity at $\sim 750$ km/s. The GHALO amplitudes are shifted
to lower $v_{\rm min}$, growing from 40% at 400 km/s to 100% already at $\sim
600$ km/s. The strong high velocity tails of $f(v)$ (Fig. 2) result in
somewhat lower modulation amplitudes compared to the best-fit MB distributions
(dotted line). At the very highest velocities $f(v)$ drops below the
Maxwellian distribution in VL2 and GHALO, and this leads to the rise in
amplitudes above the Maxwellian case for $v_{\rm min}>600$ and 550 km/s,
respectively. In GHALOs the distribution more closely follows the MB fit, and
only barely rises above it at $v_{\rm min}>670$ km/s. Interestingly, all three
simulations exhibit a pronounced dip in the modulation amplitude at close to
the highest $v_{\rm min}$. These correspond to bumps in $f(v)$ discussed in
Section 2.2.

As before, the light and dark shaded green regions in Fig. 4 cover the 68%
region around the median and the extrema of the distribution over the 100
sample spheres. Over most of the range of $v_{\rm min}$, the typical
modulation amplitude in a sample sphere differs from the spherical shell
average by less than 10%. The extrema of the distribution, however, can differ
much more significantly, especially at higher velocities. If the Sun happens
to be passing through a fast moving subhalo or tidal stream, the typical
velocity of impinging DM particles can greatly differ from the smooth halo
expectation, leading to an increase (or decrease) in the overall scattering
rate and modulation amplitude. For a conventional massive ($>10$ GeV) DM
particle, scattering events are dominated by the peak of the velocity
distribution function and the effects from streams or subhalos are washed out
and typically negligible [15]. However, whenever scattering events are
dominated by particles on the high velocity tail of the distribution, either
because of a velocity threshold (inelastic DM), a particularly low particle
mass (light DM), or for detectors intrinsically sensitive to only high recoil
energies (e.g. many directionally sensitive detectors), the presence of
velocity substructure can have a profound impact on the event rates.

Another potentially interesting signature of velocity space structure is the
shift in the peak day of the annual modulation, as was explored in [22]. For
an isotropic velocity distribution (e.g. MB) the peak day is independent of
$v_{\rm min}$ and occurs around the beginning of June. Any kind of departure
from isotropy, however, would leave its signature as a change in the phase of
the modulation. In the most extreme case of a very massive DM stream moving in
exactly the same direction as the Sun, the phase of modulation could flip
completely, with the maximum scattering amplitude occurring in the winter.

In the right panel of Fig. 4 we show the $v_{\rm min}$ dependence of the peak
day as measured in our simulations. The spherical shell average is consistent
with a phase shift of zero, except at the very highest velocities where a few
discrete velocity structures dominate the shell average and introduce small
velocity dependent phase shifts. The peak days in the sample spheres, however,
often differ by $\sim 20$ days from the fiducial value. As an example we have
plotted a curve for one of the sample spheres. The peak day determination
appears to be quite noisy here, but since the particle numbers are still
fairly large ($N(>v_{\rm min})=9114,802,105$ for $v_{\rm min}=400,600,670$
km/s in the VL2 sample sphere shown), we believe that the variations in peak
day arise from actual velocity structure rather than discreteness noise. The
68% scatter of the phase shifts over the 100 sample spheres remains remarkable
constant at $\pm 10-20$ days throught most of the $v_{\rm min}$ range. At
higher velocities ($v_{\rm min}>600-700$ km/s) the typical phase shifts grow
to $\sim 50$ days, which could be due to the increasing relative importance of
individual streams, but may also be due to particle discreteness noise. Based
on our measurements we would expect some amount of variation in the peak day
as a function recoil energy. As it is hard to imagine any background
contamination to exhibit such a phase shift, this raises the possibility of
such a measurement confirming a DM origin of the DAMA modulation signal.

We conclude this section by noting that our simulations reveal both global
(shell averaged) and local (sample spheres) departures from the standard halo
model in velocity space. These can have a significant effects on the overall
scattering rate, as well as on the amplitude and peak day of the annual
modulation thereof. In the next section we go on to explore how these effects
translate into predicted detection rates for actual direct detection
experiments, and how they modify parameter exclusions plots based on existing
null-detections.

##  3 Velocity Substructure and Directional Detection

Up to this point we have focused on the integrated features of the halo, i.e.,
$f(v)$ or $g(v_{\rm min})$, but there are many more structures that appear
that are not pronounced in these measures. In particular, the presence of
clumps and streams, while contributing to the bumps and wiggles of these
functions, can be more pronounced when the direction of the particles’
motions, rather than just the amplitude, is considered. These features are
especially important for directional WIMP detectors, such as DRIFT [35],
NEWAGE [36] and DMTPC [37], where the presence of such structures could show
up as a dramatic signal. Although we do not discuss it here, such directional
experiments have been argued to be especially important for testing inelastic
models [38, 39].

Figure 5: $f(v)$ and HealPix skymaps of the fraction of particles above
$v_{\rm min}$ coming from a given direction, for VL2 sample sphere #03 which
contains a fast-moving subhalo. The top row is in the halo rest frame ($v_{\rm
min}=400$ km/s), and the middle translated into the Earth rest frame ($v_{\rm
min}=500$ km/s). For comparison the bottom row shows the Maxwell-Boltzmann
halo case without substructure.

To quantify this, we begin by searching for “hotspots” on the sky. To do this,
we make a map of the sky in HealPix, and consider the flux of WIMPs from each
direction in the sky. We divide the sky into 192 equal regions of 215 square
degrees, and determine $p_{i}$, the fraction of particles above $v_{\rm min}$
with a velocity vector pointing towards bin $i$. We take the same 100 sample
spheres as before, but note that beyond determining sample sphere membership
we do not consider the location of a particle: the assignment to a given sky
pixel is based solely on the direction of its velocity vector. All of the
structure in a given sample sphere is considered local, i.e. able to influence
the signal at an Earth-bound direct detection experiment.

As an example of the kind of effects high velocity substructure can produce,
we show in Fig. 5 the speed distribution $f(v)$ and skymaps for VL2 sample
sphere #03. In the halo rest frame (top row, $v_{\rm min}=400$ km/s) a very
pronounced feature is visible due to the presence of a subhalo moving with
galacto-centric velocity modulus of $\sim 440$ km/s. This feature persists in
the Earth rest frame (center row, $v_{\rm min}=500$ km/s), where the direction
of the subhalo’s motion is “hotter” than the “DM headwind” hotspot in the
direction of Earth’s motion. When translating to the Earth’s rest frame, there
is an additional degree of freedom arising from the unspecified plane of the
Galactic disk. We associate radial motion towards (away from) the Galactic
Center with the center (anti-center) of the map, but the hotspot arising from
Earth’s motion is free to be rotated around this axis. In this example we have
chosen this rotation to maximize the angle between the subhalo hotspot and
Earth’s motion. In the bottom row we show for comparison the Earth rest frame
map for the MB-case without any substructure.

Figure 6: The distribution of HR$(v_{\rm min})$ as a function of $v_{\rm min}$
(left) and of $\psi$ for $v_{\rm min}=0$ and 500 km/s (right) for the VL2
simulation.

This is just one strong example, and we would like to understand what sorts of
hotspots we might expect. For this purpose we define the hotspot ratio

${\rm HR}(v_{\rm min})=\frac{{\rm max}\left\\{p_{i}(v_{\rm
min})\right\\}}{{\rm max}\left\\{p^{\rm MB}_{i}(v_{\rm min})\right\\}},$ (7)

the ratio of the hottest pixel in the sphere sample skymap above $v_{\rm min}$
to the hottest pixel in the corresponding MB-case, and the hotspot angle
$\psi$ between these two pixels. We calculate HR and $\psi$ for all 100 sample
spheres, and in each sphere for a full $2\pi$ rotation (in one degree
increments) of the direction of Earth’s motion. We show in Fig. 6 for VL2 the
distribution of HR as a function of $v_{\rm min}$ and the distribution of
$\psi$ for the case without a velocity threshold and for $v_{\rm min}=500$
km/s. For small velocities the mean of HR is unity and the r.m.s. variation is
only 10%. As $v_{\rm min}$ increases, HR grows: at $v_{\rm min}=600$ km/s, the
mean HR is $1.3\pm 0.35$, and at $v_{\rm min}=700$ km/s it’s $3.1\pm 1.4$. The
downturn at the very highest velocities is caused by running out of particles
in the sample spheres. There are also marked changes in the direction of the
hottest pixel. Even without a velocity threshold ($v_{\rm min}=0$ km/s) in 38%
of all cases the direction of the hottest spot on the sky is more than 10∘
removed from the direction of Earth’s motion (i.e. the MB halo expectation).
At $v_{\rm min}=500$ km/s the bulk of DM particles are more likely to be
coming from a local hotspot with $\psi>10^{\circ}$ than from the direction of
Earth’s motion. In this case there is only an 18% chance of having
$\psi<10^{\circ}$. Not all of this structure is due to individual subhalos;
some of it may arise from tidal streams, and some from the anisotropic
velocity distribution of the host halo particles.

##  4 Implications for DAMA

In light of the DAMA/LIBRA results [3], two scenarios which have garnered a
great deal of attention of late are light dark matter [5, 6], and inelastic
dark matter [7, 8]. These scenarios were proposed some time ago to address the
conflicts between DAMA [4] and CDMS [40] as well as other experiments at the
time.

Recent studies [27, 25, 41, 42] have shown that iDM remains a viable
explanation of the DAMA data, consistent with recent results from CDMS [43],
XENON10 [44], KIMS [45], ZEPLIN II [46], ZEPLIN III [47] and CRESST [48]. Such
models have some tension with CRESST (see the discussion in [27, 49]), which
observed 7 events in the tungsten band, while approximately 13 would be
expected [27]. While lighter masses have no tension with CDMS, higher mass
($\gtrsim 250$ GeV) WIMPs do.

Light dark matter no longer works in its original incarnation [5, 6], but
instead relies upon “channeling,” [50, 51], a difficult to quantify effect
whereby some fraction of nuclear recoils have essentially all of the energy
converted into scintillation light, rather than just a fraction. Even
including channeling, such scenarios seem to have strong tension with the data
[52, 24, 53], both in the spectrum of the modulation, and constraints from the
unmodulated event rate at low (1-2 keVee) energies333keVee stands for
“electron equivalent keV”, a unit of energy used for scintillation light,
which is produced by interactions of the recoiling nucleus with electrons. It
is related to the full nuclear recoil energy (in keVr) through the quenching
factor $q$: $E({\rm scintillation})/{\rm keVee}=q\;E({\rm recoil})/{\rm
keVr}$. $q$ is a material-dependent quantity that must be experimentally
determined.. If one disregards the lowest (2-2.5 keVee) bin from DAMA,
however, the fit improves [52, 54], but there is still tension for the
lightest particles from the presence of modulation above 4 keVee and an
overprediction of an unmodulated signal at lower energies. It is important to
note that the uncertainties in $L_{\rm eff}$, the scintillation efficiency of
liquid Xenon (see [55] for a discussion), are not completely included in these
analyses, and the most recent measurements of $L_{\rm eff}$ [55] suggest that
the low energy analysis threshold may be somewhat higher than the 4.5keVr used
with $L_{\rm eff}=0.19$, weakening the limits. On the other hand, the most
recent analysis from XENON10 [56] using a more advanced rejection of double-
scatter events shows no events at all in the signal region all the way down to
the S2 threshold, strengthening the limits. In light of these uncertainties,
and in view of the importance of the result, we believe it is best to maintain
an open mind about the viability of the light WIMP explanation.

For our purposes, the crucial feature is that these scenarios sample the high
($v\sim 600$ km/s) component of the velocity distribution, and so are
especially sensitive to the departures from a simple Maxwell-Boltzmann
distribution. Consequently, it is important to investigate whether these
changes affect the tensions described above.

###  4.1 Inelastic Dark Matter

The inelastic Dark Matter scenario [7] was proposed to explain the origin of
the DAMA annual modulation signal. The scenario relies upon three basic
elements: first, the presence of an excited state $\chi^{*}$ of the dark
matter $\chi$, with a splitting $\delta=m_{\chi^{*}}-m_{\chi}\sim
m_{\chi}v^{2}$ comparable to the kinetic energy of the WIMP. Second, the
absence of, or at most a suppressed elastic scattering cross section off of
nuclei, i.e., the process $\chi N\rightarrow\chi N$ should be small. Third, an
allowed cross section for inelastic scatterings, i.e., $\chi
N\rightarrow\chi^{*}N$, with a size set roughly by the weak scale.

Although these properties may seem odd at first blush, they are in fact
perfectly natural if the scattering occurs through a gauge interaction [7, 8],
where the splitting is between the two Majorana components of a massive
pseudo-Dirac fermion, or between the real and imaginary components of a
complex scalar. Simple models can be constructed where the mediating
interaction is the Z-boson [7, 8, 57, 41]. Models with new vector interactions
to explain the PAMELA positron excess also naturally provide models of iDM
[58], and often where all scales arise naturally from radiative corrections
[58, 59, 60, 61, 62]. Composite models provide a simple origin for the excited
states as well [63, 64].

The principle change is a kinematical requirement on the scattering,

$\displaystyle\beta_{\rm min}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{2m_{N}E_{R}}}\left(\frac{m_{N}E_{R}}{\mu}+\delta\right),$
(8)

where $m_{N}$ is the mass of the target nucleus and $\mu$ is the reduced mass
of the WIMP/target nucleus system. For $\delta\sim\mu\beta^{2}/2$, the
consequences can be significant. This requirement has three principal effects:
first, the kinematical constraint depends on the target nucleus mass, and is
more stringent for lighter nuclei. If we are sampling dominantly the tail of
the velocity distribution, the relative effect between heavy and light targets
(e.g., iodine versus germanium) can be significant. Second, again because the
signal is sampling the tail of the velocity distribution, the modulated
amplitude can be significantly higher than the few percent expected for a
standard WIMP. In fact, in the cases where there are particles kinematically
scattering in the summer, but not the winter, the modulation can reach 100%.
Third, the inelasticity suppresses or eliminates the event rate at lower
energies. Because standard WIMPs have exponentially more events at lower
energies, most experiments have focused on controlling background in this
region and lowering the threshold. In particular, the XENON10 and ZEPLINIII
experiments have few events at low energies, but a significant background at
higher energies. Consequently, their limits for standard WIMPs are quite
strong, but significantly weaker for inelastic WIMPs. The combination of these
three elements allows an explanation of the DAMA result that is consistent
with all other current experiments [27, 25, 41].

####  4.1.1 Fitting the DAMA signal

DAMA reports an annual modulation in the range of $2-6$ keVee range. This can
be interpreted as arising from either Na or I scattering events. In the former
case (Na), it has been shown that the “light inelastic” region (an
approximately 15 GeV WIMP with $\delta\approx 30{\rm~{}keV}$) can open up
significant parameter space [52, 49]. In the latter case (I), we find the
“heavy inelastic” region (a 100+ GeV WIMP with
$\delta\mathop{}_{\textstyle\sim}^{\textstyle>}100{\rm~{}keV}$) opens
significant parameter space. Since the constraints are stronger on the heavy
case, we focus on this region.

When determining the precise values of parameters that might agree with DAMA,
we must convert from keVee to keVr, which already introduces significant
uncertainties, a point which has been recently discussed by [25]. The
quenching factor for iodine has been found to have a range of different
values, including $q=0.09\pm 0.01$ [65], $q=0.05\pm 0.02$ [66], $q=0.08\pm
0.02$ [67] and $q=0.08\pm 0.01$ [68] 444The value quoted by [68] is generally
$q=0.086\pm 0.007$, but this value averages over a wide range of energies. The
two measured values for the DAMA energy range specifically are $q=0.08\pm
0.01$ and $q=0.08\pm 0.02$.. The first two measurements of the quenching
factor are somewhat indirect, fitting event distributions at low energies. The
last two are more direct. Of the last two, we should note that the first
includes non-linearity in its stated uncertainty, while the second explicitly
fixes to a signal electron energy, and thus this effect is not included. We
adopt $q=0.08\pm 0.02$ as the value for quenching, which corresponds to a
range of $25-75$ keVr for DAMA.

Figure 7: The ratio of DAMA signal in a given simulation to that in a MB
distribution. The solid red line is the spherical shell average, the dashed
line the median of the distribution over the sample spheres.

####  4.1.2 Constraining the iDM interpretation of DAMA

While iDM allows DAMA to exist consistently with the other experiments, the
reach of other experiments is at or near the predicted DAMA level. In
particular, as we see in Fig. 14, within the context of a MB halo, at high
masses, CDMS excludes the entire range of the parameter space that fits DAMA.
At the same time, the CRESST results (with a tungsten target) seem very tense
with the data over the whole mass range. These two experiments provide the
greatest present tension with the iDM interpretation of DAMA. An immediate
question we can ask is then whether the velocity-space distribution of
particles in the simulations makes the constraints more or less significant
than a MB distribution.

This is not an easy question to answer. Since there can be significant spatial
variation in the velocity distribution, we would like to quantify this without
making full exclusion plots for every data point. The three principle
constraints on the iDM interpretation of DAMA come from a) Xe experiments
(XENON10, ZEPLIN-II and ZEPLIN-III), b) CDMS (Germanium) and c) CRESST
(Tungsten). While varying the halo model can have significant effects on each
experiment, including rates and spectra, the variation of velocity structures
in the different halos affects these limits differently, and it is difficult
to quantify in aggregate what the effects are.

Before proceeding into a detailed analysis of the DAMA signal, we can already
study a preliminary question: how does the cross section needed to explain
DAMA compare between a MB halo and a simulation? We consider the ratio of the
modulation between MB and the simulations in Fig. 7. Overall, we see that the
signal at DAMA is increased in the simulation, as much as a factor of a few.

Of course, increasing the DAMA signal does not change the constraints if the
signal in any of the other experiments is changed, as well. We proceed by next
considering a ratio of ratios (RoR). Since Germanium (CDMS) has a higher
velocity threshold than Iodine (DAMA), while Tungsten (CRESST) has a lower
velocity threshold, a simple test is to look at the variation of the signal at
different experiments as a fraction of the DAMA modulation amplitude. We thus
calculate for each of the simulation samples (spherical shell and 100
spheres), as well as for the best-fitting MB model, the ratios of the CDMS and
CRESST signals, integrated from 10-100 keV, to the total DAMA modulation in
either of the high-$q$ and low-$q$ benchmark regions described above. We then
divide the simulation ratio by the MB ratio: (CDMS/[DAMA])sim/(CDMS/[DAMA])MB
and (CRESST/[DAMA])sim/(CRESST/[DAMA])MB. If these RoR’s are smaller than one,
it suggests that using the simulation’s velocity distribution instead of the
best-fitting MB distribution weakens the limits relative to DAMA, while values
larger than one imply stronger limits. The effects on Xenon limits are harder
to quantify, because of the added uncertainties in the conversion from keVee
to keVr at those experiments, and where, precisely, the backgrounds lie. Thus,
we consider first only CDMS and CRESST limits.

Figure 8: Scatter plots of the ratios of ratios (RoR, see text for details)
for the CRESST and CDMS experiments. RoR’s less than one indicate that using
the simulation’s velocity distribution instead of the best-fitting MB model
weakens the CRESST or CDMS limits relative to the DAMA (high-$q$) signal.
Figure 9: Left: The differential recoil energy distribution (in
dru$=$events/kg/day/keV) at DAMA for selected spheres from VL2, GHALO and
GHALOs, normalized to unity in the DAMA signal range 2-6 keVee Right: The
distribution at CRESST for the same spheres, normalized to unity from 0-100
keV. The models are $m_{\chi}=100{\rm~{}GeV}$ and $\delta=130\,({\rm
top}),150\,({\rm middle}),170\,({\rm bottom}){\rm~{}keV}$. Figure 10: As in
Fig. 9, but with $m_{\chi}=700{\rm~{}GeV}$.

In Fig. 8 we show scatter plots of the CRESST RoR against the CDMS RoR. Large
filled symbols indicate the spherical shell sample and small symbols are used
for the sample spheres. Note that in some cases the CDMS RoR is zero,
indicating that the simulation velocity distribution resulted in no CDMS
signal at all. A few things are immediately obvious from this plot. First, the
limits from CDMS can vary wildly between simulations, and even between
different spheres within a single simulation. This simply represents how
dramatically the velocity distribution can change at the highest velocities.
Second, we see that the CRESST rate is much more weakly affected, with
typically suppressions of (0.6-1), (0.7-1.1), and (0.7-1.3) for Via Lactea II,
GHALO and GHALOs respectively. Thus, from the perspective of Poisson limits,
these results would suggest that those derived from Maxwellian halos are
possibly excessively aggressive by almost a factor of two. However, many
limits are placed using one of Yellin’s techniques [69], where not only the
overall rate, but also the distribution of signal versus background is
important. Here we find that the halo uncertainties can be at their largest.
We show in Fig. 9 the spectra at DAMA and CRESST for a 100 GeV WIMP with
$\delta=130,150,170{\rm~{}keV}$. We employ energy smearing at DAMA by assuming
that the smearing reported for the one of 25 targets of DAMA/LIBRA [70] is
characteristic of all of them. We assume a smearing of 1 keV at CRESST.

One can see that the peak positions and properties can change by quite a large
amount. At DAMA, the effect of this is principally to shift the peak. Such an
effect is largely degenerate with the quenching uncertainty, which can
reasonably range from $q=0.06$ to $q=0.09$. (For a lower quenching value, for
instance, the peak will shift to lower energy in keVee.) The effects would be
similar at XENON10, where the energy smearing will eliminate most interesting
structures, and the shift in the location of the peak will be comparable to
the uncertainties induced from $L_{\rm eff}$. In contrast, the spectrum for
CRESST can vary dramatically, with the peak location moving from below
$30{\rm~{}keV}$ to above $60{\rm~{}keV}$. The implication of this is that
techniques such as optimum interval, maximum gap and ${\rm p_{max}}$ are
likely all overly aggressive in that a peak in the spectrum might arise at a
specific location in the real Milky Way halo, but that would not be reproduced
at the appropriate position, for instance in a Maxwellian halo. In general,
any technique that relies on knowing precisely the predicted spectrum is
unreliable when considering these variations.

Figure 11: Allowed parameter space for DAMA at 90% (purple) and 99% (blue)
with $m_{\chi}=70{\rm~{}GeV}$. For comparison the 90% and 99% regions for
DAMA/LIBRA rates only (as opposed to DAMA/NaI + DAMA/LIBRA) are shown in red
and green lines. Constraints are CDMS (solid), ZEPLIN-III (long dashed, thin),
ZEPLIN-II (long-dashed, thick), CRESST (medium-dashed) and XENON10 (short-
dashed). The regions are MB with $v_{0}=220{\rm~{}km/s}$ and $v_{\rm
esc}=550{\rm~{}km/s}$ (top left), VL2 (top right), and a sample sphere from
VL2 (bottom, left) and GHALOs (bottom, right). Figure 12: As in 11, but with
$m_{\chi}=150{\rm~{}GeV}$. Figure 13: As in 11, but with
$m_{\chi}=300{\rm~{}GeV}$. Figure 14: As in 11, but with
$m_{\chi}=700{\rm~{}GeV}$.

Although we are still limited by our sample of simulations, we can study the
effects on various limits by looking at the allowed parameter space and limits
in a variety of halos. To do so, we largely follow [27] in calculating the
allowed parameter space and limits, with a few important differences. As
pointed out by [25], the uncertainty in quenching factor for iodine can be
extremely important. We thus consider the range $q=0.06$ to $q=0.09$ and take
the 90% allowed region to be the union of 90% allowed regions, and similarly
for the 99% region. We smear the signal as described above.

For CDMS we employ the maximum gap method with the data set specified in [27].

For limits arising from XENON10, we use the data presented in the recent
reanalysis of [56], and calculate the maximum gap limit (as in ZEPLIN-II and
ZEPLIN-III) for both [71] and [55] values of $L_{\rm eff}$, taking at every
point the more conservative of the two. For limits, we employ the maximum gap
method. Since the data are smeared and there are great uncertainties in
$L_{\rm eff}$, the detailed spectral information arising from structure is
lost.

For ZEPLIN-II, we take the data described in [46], but (conservatively) take
the energy values to be shifted by a factor of $0.24/0.19$ for a maximum gap
analysis, which employs the value of $L_{\rm eff}$ from [71] at higher
energies.555We could parameterize $L_{\rm eff}$ as a function of energy, and
use different shifts at different energies, but the limits are essentially the
same, being dominated by the events at the highest energy.

For ZEPLIN-III, we utilize the data within the published acceptance box up to
15 keVee for the maximum gap analysis. We additionally employ data outside the
blind acceptance box up to 30 keVee, the maximum to which efficiencies were
provided in [47]. We use the delineated $1\sigma$ region for the data in this
high energy range since the $3\sigma$ region is not specified. To convert from
keVee to keVr, we use

$E_{r}=(0.142E_{\rm ee}+0.005)\,\exp(-0.305\,E_{\rm ee}^{0.564})$ (9)

below 10 keVee [25, 42], and a constant ratio of 0.55 above that. (We use a
larger value than [42] at high energies to account for the energy dependent
value of $L_{\rm eff}$ suggested by [71], which tends to weaken limits
slightly.) Note that we do not employ a background subtraction using the
science data as in [72], nor do we do so for ZEPLIN-II as the same basic
technique seems to have led to a significant overprediction of background for
ZEPLIN-III.

For all Xe based experiments, we take a fixed smearing of $\sqrt{30}$ keVr,
which is typical of the smearing of the experiments in the range of the peak
signal.

For the CRESST experiment, we consider the data from the commissioning run
[48], as well as the additional events at $\sim$ 22 keV, 33 keV and 88 keV
presented for Run 31 by [73], and take a smearing of 1 keV. As described
above, Yellin-style techniques are ironically inappropriate for experiments
with high energy resolution. Because of the significant variation in the
possible spectrum when varying halo models with the high energy resolution, we
use a binned Poisson analysis, rather than a maximum-gap technique. We divide
up the signal region into a low region ($E_{R}<35{\rm~{}keV}$), a middle
region ($35{\rm~{}keV}<E_{R}<50{\rm~{}keV}$) and a high region
$(50{\rm~{}keV}<E_{R}<100{\rm~{}keV})$. In this way, we do not overweight the
events near the zero of the form factor (in the mid region), but still include
some broad spectral information. We require that the rate be lower than 95%
Poisson upper limits for each bin, which is a 95% limit when the signal is
exclusively in the low bin and $\sim$ 90% when there is signal in both the low
and high bins, and $\sim$ 85% in the (rarer) instances when there is signal in
all three bins. Although stronger limits can be set using for instance,
optimum interval, because of the uncertainties in the halo distribution, this
seems inappropriate. In other words, if some of the events at CRESST are real,
then the distribution of events may well be telling us the properties of the
halo.

We show the results of these limits in Figs. 11-14. Even when accounting for
smearing, and using a binned Poisson limit instead of optimum interval, it is
remarkable how significantly the allowed parameter space can change from
simulation to simulation, and between spheres in the same simulation. We see
that CRESST remains the most constraining, as found previously, but we see
important differences. First, agreement is generally improved in these
spheres. This is due to the presence of structure at high velocities which
increases the modulated fraction at DAMA. As a consequence, small regions of
allowed parameter space exist at high mass, the size of which depends
significantly on the particular halo. Larger regions exist at smaller mass
($\sim 150{\rm~{}GeV}$). Moreover, the allowed region of parameter space can
shift dramatically in $\sigma$ and $\delta$, with pockets appearing at large
$\delta$ from high velocity particles. Should iDM be relevant for nature, the
detailed nature of the halo is clearly significant in determining the
properties of the direct detection signals.

Figure 15: Spectra for a sample of spheres best fit to the DAMA data (shown)
within the different simulations for LDM models with $m_{\chi}=3{\rm~{}GeV}$
(top, left), $7{\rm~{}GeV}$ (top, right), and $13{\rm~{}GeV}$ (bottom). Figure
16: The allowed parameter space for LDM models. Top, left: MB with
$v_{0}=220{\rm~{}km/s}$ and $v_{\rm esc}=550{\rm~{}km/s}$. Top, right: VL2.
Middle, left: GHALO. Middle, right GHALOs. Bottom two sample spheres taken
from the simulations.

###  4.2 Light dark matter

Light ($\sim$ GeV) dark matter [5] has been proposed as a means to reconcile
DAMA with the other existing experiments [6, 74, 51]. The current iteration of
LDM involves “channelling” [75, 50] whereby for some small fraction of the
time the entirety of the nuclear recoil energy is converted to scintillation
light. Since the observed energy is lower, it allows lighter WIMPs to scatter
at the DAMA experiment. Such light WIMPs may be incapable of depositing energy
above the threshold at XENON10 or CDMS-Si, for instance, allowing DAMA to
evade these bounds.

Such an explanation is not without difficulty, however. These light WIMPs have
a rapidly falling event spectrum, due to the decreasing probability of
channeling at higher energies, and the increasingly suppressed number of
particles with higher kinetic energy. As a consequence the spectrum of LDM
seems generally to fall too rapidly at light masses [52, 24, 53]. At higher
masses, the spectrum is acceptable, but is excluded by other experiments. We
should note again here that this is under the assumption that the errors in
the data are statistical only. Should there be, e.g., significant changes in
the efficiency, it is possible that LDM could give a better fit. As noted
earlier, the uncertainties in $L_{\rm eff}$ should weaken these limits, while
the most recent analysis from XENON10 (with no events to the S2 threshold)
would strengthen the limits. Thus, while the models do not seem to describe
the data well, there are adequate uncertainties that this scenario remains an
interesting possibility as an explanation of the DAMA signal.

Since these models are also extremely sensitive to the presence of high
velocity particles, one might wonder whether, just as in the iDM case, the
properties of halos at high velocities might modify the spectrum and the
relative signal strength at other experiments.

We begin again by addressing the question of the effects on the spectrum. We
show these for three different masses in Fig. 15. One can see that, while
there are noticeable changes in the spectrum, its overall shape remains
basically unchanged. The difference between LDM and iDM here is simple: in
iDM, the “threshold” (i.e., minimum velocity) scattering is actually at large
energy, and to go to lower energies one requires higher velocities. Competing
against this are the significant form factor effects. The competition between
these two sensitive quantities leads to dramatic changes and features. For
LDM, the threshold scattering is at zero, and all effects (channeling, form
factor, velocity distribution) contribute to the falling (unmodulated)
spectrum. Thus, it can fall somewhat more or less rapidly, depending on the
number of particles at high velocity, but without very significant changes.

To determine the spectrum, we follow the analysis of [52]. The most proper
analysis would be to include the updated measurement of $L_{\rm eff}$ (which
will weaken the limits) and the reanalysis to remove multi-hit events from
[56] (which strengthens the limits). However, incorporating the S2 threshold,
resulting in contribution to an energy-dependent acceptance at low energies is
complicated. Lacking detailed information on this, we will restrict ourselves
to the previous XENON10 limits. Since our point is the relative change from
halo model to halo model, this is reasonable, so long as we recognize that the
limits should be taken with a grain of salt.

We show in Fig. 16 the allowed parameter space for the LDM scenario for the
canonical MB model, for the (averaged) simulations VL2, GHALO and GHALO-
scaled, as well as two sample spheres from the halos. We see that the small
spectral effects translate into small effects in the allowed parameter space.
There can be some improvement, however, and although the exclusion curves
technically cover the allowed regions, uncertainties in the low-energy
scintillation of Xe, for instance, and the ambiguity of employing statistical
errors without a covariance matrix suggest that a small region may yet be
allowed near 10 GeV. We should again emphasize that there is no sense in which
the sampling we have done is statistically complete, and our failure to find
an allowed region does not preclude the possibility that one exists, or might
arise in another simulation. A proper reanalysis using the data from [56] is
warranted.

##  5 Discussion

Direct detection of dark matter intimately ties up questions of astrophysics
and particle physics. In recent years, the range of particle physics models
has exploded, many with new interactions and properties. Some of these, in
particular light dark matter and inelastic dark matter, are particularly
sensitive to the high velocity tail of the WIMP velocity distribution. The
same is true, even for standard WIMPs, for many directionally sensitive
detection experiments (e.g. DRIFT, NEWAGE, DMTPC), which rely on high recoil
energies in order to measure the direction of the scattering DM particle. It
is especially at these high velocities that deviations from the standard
isothermal Maxwell-Boltzmann assumption are expected to be most significant.
Velocity substructure, arising from nearby subhalos and tidal streams, as well
as from the host halo’s velocity anisotropy, can have profound effects on the
expected event rates, the recoil energy spectrum, and the typical direction of
scattering DM particles.

We have studied the effects of these deviations by employing data from N-body
simulations directly, rather than with analytic fits that miss interesting
structures. Our analysis is based on two of the highest resolution numerical
simulations of Galactic dark matter structure, Via Lactea II and GHALO. We
confirm previously reported global departures from the Maxwell-Boltzmann
distribution, with a noticeable excess of particles on the high velocity tail.
In many of our local samples we find additional discrete features due to
subhalos and tidal streams.

These global and local departures from the standard Maxwellian model have a
number of interesting consequences:

  * •

For directionally sensitive detection, we have found that the fraction of
particles coming from the “hottest” direction in the sky increases with the
velocity threshold. More importantly, the direction from which most DM
particles are coming can be affected significantly. At high velocity
thresholds ($\gtrsim 500$ km/s in the Earth’s frame) the typical direction is
more often than not shifted by $>10^{\circ}$ away from the direction of
Earth’s motion, and up to $80^{\circ}$ in extreme cases.

  * •

For inelastic DM the effects are dramatic. The relative strength of CDMS
limits versus DAMA can change by an order of magnitude, while the relative
strength of CRESST limits can change by almost a factor of two. For CRESST
most of this effect is due to an enhanced modulation, which can be a factor of
two larger than for a Maxwellian halo. Such effects would persist at Xenon
detectors as well, changing their sensitivity relative to DAMA by
$\mathcal{O}(1)$.

  * •

For light DM, we have seen that small changes in the spectrum are possible,
although not so much as to allow very light ($\sim$ 1 GeV) WIMPs to fit the
existing data, but only somewhat heavier
($\mathop{}_{\textstyle\sim}^{\textstyle>}$ 10 GeV).

  * •

In general, dramatic variations in the spectrum can appear for iDM models,
complicating attempts to employ Yellin-type analyses, which rely upon a
detailed knowledge of the predicted spectrum of the model. Instead we advocate
a binned Poisson analysis, which is less sensitive to assumption about the
shape of the signal spectrum.

  * •

A further important effect that we have found is the possible variation of the
phase of modulation as a function of energy. This has been noted before in the
context of streams [76]. Since most background modulations would be expected
to have the same phase as a function of energy, such a change could be a
strong piece of evidence that a modulated signal was arising from DM.

It is important to keep in mind that our simulations are nowhere close to
resolving the phase-space structure at the relevant sub-parsec scales, and we
rely on extrapolation from a coarser sampling (1-1.5 kpc radius spheres). If
velocity substructure does not persist on smaller scales, then our analysis
may overestimate the likelihood of these effects. On the other hand, we are
probably underestimating the significance of the effects, since at the moment
we are diluting the substructure signal by the host halo particles, while it
would completely dominate the background host halo, should the Earth be
passing through one of these substructures.

We have seen here how, should DM be conclusively discovered with direct
detection experiments, it may begin to help unlock questions about the
formation of the galaxy, precisely because of these dramatic sensitivities.
With the fantastic improvements in direct detection on the horizon it is more
important than ever to increase our awareness of the significant impact that
the detailed structure of our Galaxy’s dark matter halo may have.

## Acknowledgments

The authors thank S.Chang and A.Pierce for discussions and their collaboration
in the development of the DM direct detection analysis tools. MK gratefully
acknowledges support from William L. Loughlin at the Institute for Advanced
Study, and from the Theoretical Astrophysics Center at UC Berkeley. NW is
supported by DOE OJI grant # DE-FG02-06ER41417 and NSF CAREER grant
PHY-0449818.

## References

  * [1] F. J. Kerr and D. Lynden-Bell, Review of galactic constants, Mon. Not. Roy. Astron. Soc. 221 (1986) 1023.
  * [2] M. Kamionkowski and A. Kinkhabwala, Galactic halo models and particle dark matter detection, Phys. Rev. D57 (1998) 3256–3263, [hep-ph/9710337].
  * [3] DAMA Collaboration, R. Bernabei et al., First results from DAMA/LIBRA and the combined results with DAMA/NaI, Eur. Phys. J. C56 (2008) 333–355, [arXiv:0804.2741].
  * [4] DAMA Collaboration, R. Bernabei et al., Search for WIMP annual modulation signature: Results from DAMA / NaI-3 and DAMA / NaI-4 and the global combined analysis, Phys. Lett. B480 (2000) 23–31.
  * [5] A. Bottino, F. Donato, N. Fornengo, and S. Scopel, Light neutralinos and WIMP direct searches, Phys. Rev. D69 (2004) 037302, [hep-ph/0307303].
  * [6] P. Gondolo and G. Gelmini, Compatibility of DAMA dark matter detection with other searches, Phys. Rev. D71 (2005) 123520, [hep-ph/0504010].
  * [7] D. R. Smith and N. Weiner, Inelastic dark matter, Phys. Rev. D64 (2001) 043502, [hep-ph/0101138].
  * [8] D. Tucker-Smith and N. Weiner, The status of inelastic dark matter, Phys. Rev. D72 (2005) 063509, [hep-ph/0402065].
  * [9] J. Diemand, M. Kuhlen, and P. Madau, Dark matter substructure and gamma-ray annihilation in the Milky Way halo, Astrophys. J. 657 (2007) 262, [astro-ph/0611370].
  * [10] J. Diemand et al., Clumps and streams in the local dark matter distribution, Nature 454 (2008) 735–738, [arXiv:0805.1244].
  * [11] J. Stadel, D. Potter, B. Moore, J. Diemand, P. Madau, M. Zemp, M. Kuhlen, and V. Quilis, Quantifying the heart of darkness with GHALO - a multibillion particle simulation of a galactic halo, MNRAS 398 (Sept., 2009) L21–L25, [arXiv:0808.2981].
  * [12] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, and S. D. M. White, The Aquarius Project: the subhaloes of galactic haloes, MNRAS 391 (Dec., 2008) 1685–1711, [arXiv:0809.0898].
  * [13] M. Zemp, J. Diemand, M. Kuhlen, P. Madau, B. Moore, D. Potter, J. Stadel, and L. Widrow, The graininess of dark matter haloes, MNRAS 394 (Apr., 2009) 641–659, [arXiv:0812.2033].
  * [14] A. A. Klypin, S. Gottlober, and A. V. Kravtsov, Galaxies in N-body Simulations: overcoming the overmerging problem, Astrophys. J. 516 (1999) 530–551, [astro-ph/9708191].
  * [15] M. Vogelsberger, A. Helmi, V. Springel, S. D. M. White, J. Wang, C. S. Frenk, A. Jenkins, A. Ludlow, and J. F. Navarro, Phase-space structure in the local dark matter distribution and its signature in direct detection experiments, MNRAS 395 (May, 2009) 797–811, [arXiv:0812.0362].
  * [16] M. Kamionkowski and S. M. Koushiappas, Galactic Substructure and Direct Detection of Dark Matter, Phys. Rev. D77 (2008) 103509, [arXiv:0801.3269].
  * [17] J. Diemand, B. Moore, and J. Stadel, Velocity and spatial biases in CDM subhalo distributions, Mon. Not. Roy. Astron. Soc. 352 (2004) 535, [astro-ph/0402160].
  * [18] R. Wojtak, E. L. Lokas, S. Gottloeber, and G. A. Mamon, Radial velocity moments of dark matter haloes, Mon. Not. Roy. Astron. Soc. 361 (2005) L1, [astro-ph/0503391].
  * [19] S. H. Hansen, B. Moore, M. Zemp, and J. Stadel, A universal velocity distribution of relaxed collisionless structures, JCAP 0601 (2006) 014, [astro-ph/0505420].
  * [20] S. H. Hansen and B. Moore, A universal density slope Velocity anisotropy relation for relaxed structures, New Astronomy 11 (Mar., 2006) 333–338, [arXiv:astro-ph/0411473].
  * [21] B. Moore, C. Calcáneo-Roldán, J. Stadel, T. Quinn, G. Lake, S. Ghigna, and F. Governato, Dark matter in Draco and the Local Group: Implications for direct detection experiments, Phys. Rev. D 64 (Sept., 2001) 063508–+, [arXiv:astro-ph/0106271].
  * [22] C. Savage, K. Freese, and P. Gondolo, Annual Modulation of Dark Matter in the Presence of Streams, Phys. Rev. D74 (2006) 043531, [astro-ph/0607121].
  * [23] J. D. Vergados, S. H. Hansen, and O. Host, Impact of going beyond the Maxwell distribution in direct dark matter detection rates, Phys. Rev. D 77 (Jan., 2008) 023509–+, [arXiv:0711.4895].
  * [24] M. Fairbairn and T. Schwetz, Spin-independent elastic WIMP scattering and the DAMA annual modulation signal, JCAP 0901 (2009) 037, [arXiv:0808.0704].
  * [25] J. March-Russell, C. McCabe, and M. McCullough, Inelastic dark matter, non-standard halos and the DAMA/LIBRA results, Journal of High Energy Physics 5 (May, 2009) 71–+, [arXiv:0812.1931].
  * [26] F. Ling, E. Nezri, E. Athanassoula, and R. Teyssier, Dark Matter Direct Detection Signals inferred from a Cosmological N-body Simulation with Baryons, ArXiv e-prints (Sept., 2009) [arXiv:0909.2028].
  * [27] S. Chang, G. D. Kribs, D. Tucker-Smith, and N. Weiner, Inelastic Dark Matter in Light of DAMA/LIBRA, Phys. Rev. D79 (2009) 043513, [arXiv:0807.2250].
  * [28] W. Dehnen and J. Binney, Local stellar kinematics from Hipparcos data, Mon. Not. Roy. Astron. Soc. 298 (1998) 387–394, [astro-ph/9710077].
  * [29] J. D. Lewin and P. F. Smith, Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil, Astropart. Phys. 6 (1996) 87–112.
  * [30] L. Hernquist, N-body realizations of compound galaxies, Astrophys. J. Suppl. 86 (1993) 389–400.
  * [31] S. Kazantzidis, J. Magorrian, and B. Moore, Generating Equilibrium Dark Matter Halos: Inadequacies of the Local Maxwellian Approximation, Astrophys. J. 601 (2004) 37–46, [astro-ph/0309517].
  * [32] N. Afshordi, R. Mohayaee, and E. Bertschinger, Hierarchy in the Phase Space and Dark Matter Astronomy, ArXiv e-prints (Nov., 2009) [arXiv:0911.0414].
  * [33] M. Vogelsberger, S. D. M. White, A. Helmi, and V. Springel, The fine-grained phase-space structure of cold dark matter haloes, MNRAS 385 (Mar., 2008) 236–254, [arXiv:0711.1105].
  * [34] M. C. Smith et al., The RAVE Survey: Constraining the Local Galactic Escape Speed, Mon. Not. Roy. Astron. Soc. 379 (2007) 755–772, [astro-ph/0611671].
  * [35] D. P. Snowden-Ifft, T. Ohnuki, E. S. Rykoff, and C. J. Martoff, Neutron recoils in the DRIFT detector, Nuclear Instruments and Methods in Physics Research A 498 (Feb., 2003) 155–164.
  * [36] H. Nishimura, K. Miuchi, K. Hattori, S. Iwaki, C. Ida, S. Kabuki, H. Kubo, S. Kurosawa, H. Sekiya, A. Takada, M. Takahashi, A. Takeda, T. Tanimori, K. Tsuchiya, and K. Ueno, Tracking performance of the NEWAGE-0.3a direction-sensitive dark matter detector, Astroparticle Physics 31 (Apr., 2009) 185–191.
  * [37] G. Sciolla and the DMTPC Collaboration, The DMTPC project, Journal of Physics Conference Series 179 (July, 2009) 012009–+, [arXiv:0903.3895].
  * [38] D. P. Finkbeiner, T. Lin, and N. Weiner, Inelastic Dark Matter and DAMA/LIBRA: An Experimentum Crucis, arXiv:0906.0002.
  * [39] M. Lisanti and J. G. Wacker, Disentangling Dark Matter Dynamics with Directional Detection, arXiv:0911.1997.
  * [40] CDMS Collaboration, R. Abusaidi et al., Exclusion limits on the WIMP nucleon cross-section from the Cryogenic Dark Matter Search, Phys. Rev. Lett. 84 (2000) 5699–5703, [astro-ph/0002471].
  * [41] Y. Cui, D. E. Morrissey, D. Poland, and L. Randall, Candidates for Inelastic Dark Matter, JHEP 05 (2009) 076, [arXiv:0901.0557].
  * [42] C. Arina, F.-S. Ling, and M. H. G. Tytgat, The Inert Doublet Model and Inelastic Dark Matter, arXiv:0907.0430.
  * [43] CDMS Collaboration, Z. Ahmed et al., Search for Weakly Interacting Massive Particles with the First Five-Tower Data from the Cryogenic Dark Matter Search at the Soudan Underground Laboratory, Phys. Rev. Lett. 102 (2009) 011301, [arXiv:0802.3530].
  * [44] XENON Collaboration, J. Angle et al., First Results from the XENON10 Dark Matter Experiment at the Gran Sasso National Laboratory, Phys. Rev. Lett. 100 (2008) 021303, [arXiv:0706.0039].
  * [45] KIMS Collaboration, H. S. Lee. et al., Limits on WIMP-nucleon cross section with CsI(Tl) crystal detectors, Phys. Rev. Lett. 99 (2007) 091301, [arXiv:0704.0423].
  * [46] G. J. Alner et al., First limits on WIMP nuclear recoil signals in ZEPLIN-II: A two phase xenon detector for dark matter detection, Astropart. Phys. 28 (2007) 287–302, [astro-ph/0701858].
  * [47] V. N. Lebedenko et al., Result from the First Science Run of the ZEPLIN-III Dark Matter Search Experiment, arXiv:0812.1150.
  * [48] G. Angloher et al., Commissioning Run of the CRESST-II Dark Matter Search, arXiv:0809.1829.
  * [49] K. Schmidt-Hoberg and M. W. Winkler, Improved Constraints on Inelastic Dark Matter, arXiv:0907.3940.
  * [50] R. Bernabei et al., Possible implications of the channeling effect in NaI(Tl) crystals, Eur. Phys. J. C53 (2008) 205–213, [arXiv:0710.0288].
  * [51] F. Petriello and K. M. Zurek, DAMA and WIMP dark matter, JHEP 09 (2008) 047, [arXiv:0806.3989].
  * [52] S. Chang, A. Pierce, and N. Weiner, Using the Energy Spectrum at DAMA/LIBRA to Probe Light Dark Matter, Phys. Rev. D79 (2009) 115011, [arXiv:0808.0196].
  * [53] C. Savage, G. Gelmini, P. Gondolo, and K. Freese, Compatibility of DAMA/LIBRA dark matter detection with other searches, JCAP 0904 (2009) 010, [arXiv:0808.3607].
  * [54] D. Hooper, F. Petriello, K. M. Zurek, and M. Kamionkowski, The New DAMA Dark-Matter Window and Energetic-Neutrino Searches, Phys. Rev. D79 (2009) 015010, [arXiv:0808.2464].
  * [55] A. Manzur et al., Scintillation efficiency and ionization yield of liquid xenon for mono-energetic nuclear recoils down to 4 keV, arXiv:0909.1063.
  * [56] X. Collaboration, Constraints on inelastic dark matter from XENON10, arXiv:0910.3698.
  * [57] C. Arina and N. Fornengo, Sneutrino cold dark matter, a new analysis: Relic abundance and detection rates, JHEP 11 (2007) 029, [arXiv:0709.4477].
  * [58] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, A Theory of Dark Matter, Phys. Rev. D79 (2009) 015014, [arXiv:0810.0713].
  * [59] N. Arkani-Hamed and N. Weiner, LHC Signals for a SuperUnified Theory of Dark Matter, JHEP 12 (2008) 104, [arXiv:0810.0714].
  * [60] M. Baumgart, C. Cheung, J. T. Ruderman, L.-T. Wang, and I. Yavin, Non-Abelian Dark Sectors and Their Collider Signatures, JHEP 04 (2009) 014, [arXiv:0901.0283].
  * [61] C. Cheung, J. T. Ruderman, L.-T. Wang, and I. Yavin, Kinetic Mixing as the Origin of Light Dark Scales, arXiv:0902.3246.
  * [62] A. Katz and R. Sundrum, Breaking the Dark Force, JHEP 06 (2009) 003, [arXiv:0902.3271].
  * [63] D. S. M. Alves, S. R. Behbahani, P. Schuster, and J. G. Wacker, Composite Inelastic Dark Matter, arXiv:0903.3945.
  * [64] M. Lisanti and J. G. Wacker, Parity Violation in Composite Inelastic Dark Matter Models, arXiv:0911.4483.
  * [65] R. Bernabei et al., New limits on WIMP search with large-mass low- radioactivity NaI(Tl) set-up at Gran Sasso, Phys. Lett. B389 (1996) 757–766.
  * [66] K. Fushimi et al., Application of a large volume NaI scintillator to search for dark matter, Phys. Rev. C47 (1993) 425–428.
  * [67] G. Gerbier et al., Pulse shape discrimination and dark matter search with NaI(Tl) scintillator, Astropart. Phys. 11 (1999) 287–302.
  * [68] D. R. Tovey et al., Measurement of scintillation efficiencies and pulse-shapes for nuclear recoils in NaI(Tl) and CaF-2(Eu) at low energies for dark matter experiments, Phys. Lett. B433 (1998) 150–155.
  * [69] S. Yellin, Finding an upper limit in the presence of unknown background, Phys. Rev. D66 (2002) 032005, [physics/0203002].
  * [70] DAMA Collaboration, R. Bernabei et al., The DAMA/LIBRA apparatus, arXiv:0804.2738.
  * [71] P. Sorensen et al., The scintillation and ionization yield of liquid xenon for nuclear recoils, Nucl. Instrum. Meth. A601 (2009) 339–346, [arXiv:0807.0459].
  * [72] D. B. Cline, W. Ooi, and H. Wang, Constraints on Inelastic Dark Matter Signal using ZEPLIN- II Results, arXiv:0906.4119.
  * [73] CRESST Collaboration, F. Probst, , Talk given at the 12th Marcel Grossman Meeting, July 12-18, 2009, Paris, France.
  * [74] A. Bottino, F. Donato, N. Fornengo, and S. Scopel, Interpreting the recent results on direct search for dark matter particles in terms of relic neutralino, Phys. Rev. D78 (2008) 083520, [arXiv:0806.4099].
  * [75] A. Bottino, F. Donato, N. Fornengo, and S. Scopel, Zooming in on light relic neutralinos by direct detection and measurements of galactic antimatter, Phys. Rev. D77 (2008) 015002, [arXiv:0710.0553].
  * [76] R. Bernabei et al., Dark matter particles in the galactic halo: Results and implications from DAMA/NaI, Int. J. Mod. Phys. D13 (2004) 2127–2160, [astro-ph/0501412].