Dark Matter Direct Detection Signals inferred from a Cosmological N-body Simulation with Baryons
PPreprint typeset in JHEP style - HYPER VERSION
ULB-TH/09-30
Dark Matter Direct Detection Signals inferred from aCosmological N-body Simulation with Baryons
F.-S. Ling , E. Nezri , E. Athanassoula , R. Teyssier , Service de Physique Th´eorique, Universit´e Libre de Bruxelles,CP225, Bld du Triomphe, 1050 Brussels, Belgium Laboratoire d’Astrophysique de Marseille,Observatoire Astronomique de Marseille ProvenceCNRS/Universit´e de Provence38 rue Joliot Curie, 13388 Marseille, France IRFU, CEA SaclayL’Orme des merisiers 91191 Gif-sur-Yvette, France Institute f¨ur Theoretische Physik, Universit¨at Z¨urich,Winterthurerstrasse 190, CH-8057 Z¨urich, SwitzerlandE-mails: [email protected] , [email protected] , [email protected] , [email protected] – 1 – a r X i v : . [ a s t r o - ph . GA ] M a r bstract: We extract at redshift z = 0 a Milky Way sized object including gas, starsand dark matter (DM) from a recent, high-resolution cosmological N-body simulation withbaryons. Its resolution is sufficient to witness the formation of a rotating disk and bulgeat the center of the halo potential, therefore providing a realistic description of the birthand the evolution of galactic structures in the ΛCDM cosmology paradigm.The phase-space structure of the central galaxy reveals that, throughout a thick region,the dark halo is co-rotating on average with the stellar disk. At the Earth’s location, therotating component, sometimes called dark disk in the literature, is characterized by a min-imum lag velocity v lag (cid:39)
75 km/s, in which case it contributes to around 25% of the totalDM local density, whose value is ρ DM (cid:39) .
37 GeV / cm . The velocity distributions alsoshow strong deviations from pure Gaussian and Maxwellian distributions, with a sharperdrop of the high velocity tail.We give a detailed study of the impact of these features on the predictions for DM signalsin direct detection experiments. In particular, the question of whether the modulationsignal observed by DAMA is or is not excluded by limits set by other experiments (CDMS,XENON and CRESST...) is re-analyzed and compared to the case of a standard Maxwellianhalo. We consider spin-independent interactions for both the elastic and the inelasticscattering scenarios. For the first time, we calculate the allowed regions for DAMA andthe exclusion limits of other null experiments directly from the velocity distributions foundin the simulation. We then compare these results with the predictions of various analyticaldistributions.We find that the compatibility between DAMA and the other experiments is improved.In the elastic scenario, the DAMA modulation signal is slightly enhanced in the so-calledchanneling region, as a result of several effects that include a departure from a Maxwelliandistribution and anisotropies in the velocity dispersions due to the dark disk. For theinelastic scenario, the improvement of the fit is mainly attributable to the departurefrom a Maxwellian distribution at high velocity. It is correctly modeled by a general-ized Maxwellian distribution with a parameter α (cid:39) .
95, or by a Tsallis distribution with q (cid:39) . Keywords: dark matter, N-body simulations, direct detection experiments. ontents
1. Introduction 12. A cosmological simulation with baryons 3
3. Direct detection experiments and signals 15
4. Summary & Perspectives 29
1. Introduction
Nowadays, the hypothesis of dark matter (DM) is considered as the most plausible expla-nation for various observations from galactic to cosmological scales. Flat rotation curves inspiral galaxies [1, 2], velocity dispersions in galaxy clusters [3, 4], gravitational lensing massreconstructions [5–7], hierarchical patterns in large scale structures [8, 9] and cosmic mi-crowave background anisotropies [10] all point towards the existence of a new neutral, nonbaryonic, cold ( i.e. with low thermal velocities) and weakly interacting massive particle(WIMP) in the universe.Theories beyond the Standard Model of particle physics [11–13] provide many candi-dates for DM particles (see for example Refs. [14–17]). However, its existence will not beproved and its properties will not be determined until a clear identification is made throughsignals in indirect or direct detection experiments or in particle accelerators. Recently, sev-eral astrophysical excesses have been reported [18–21], and associated with a DM signalin more or less exotic scenarios (see, among many others, Refs. [22–27]). However, up tonow, none of them requires DM as an irrefutable explanation. The INTEGRAL 511 keVline γ -ray signal [18] can be explained by the emission of low mass X-ray binaries [28].The EGRET diffuse γ -ray GeV anomaly [19] was most likely due to an error in the en-ergy calibration [29], and has since indeed been ruled out by the FERMI instrument [30].The positron and electron peak seen by ATIC [21] has been flattened by HESS [31] and– 1 –ERMI [32] measurements. The positron fraction rise observed by PAMELA [20], althoughstriking, can be accommodated in standard astrophysical scenarios with pulsars [33, 34], orwith an inhomogeneous distribution of supernovae remnants [35].Direct detection provides a unique opportunity to disentangle DM signals from otherastrophysical processes. The idea is to directly measure the energy deposited in a lownoise detector during the collision between a DM particle and a detector nucleus. Severalexperiments have been designed for this aim [36–40]. Up to now, all of them are compatiblewith null results, except the DAMA/LIBRA collaboration which claims to have observeda signal with an annual modulation characteristic of an DM signal [41]. This modulationis induced by the Earth rotation around the Sun, and is presented as a signature hard toreproduce by background events. The DAMA result is still quite controversial because itleads to tensions with other experiment measurements (for a recent status, see Refs. [42–45]).Predicting direct detection signals in any particle physics model requires astrophysicalassumptions in the form of a value for the local DM density as well as velocity distributionsin the solar neighborhood. Usually, the local DM density is taken as ρ DM = 0 . / cm =0 .
079 M
Sun / pc , a value which is supported by the Milky-Way rotation curve [46, 47], andstars spectrophotometric data [48]. A recent determination which takes numerous dynam-ical observables into account gives a larger value, ρ DM = 0 .
39 GeV / cm [49]. However, itshould be mentioned that determinations based on mass models only give average values.A local peak in the DM density, although improbable without a dynamical reason [50],is in principle not excluded by very local bounds based on planetary data [51, 52]. An-other common over-simplified assumption is to take the velocity distribution as isotropicand Maxwellian. Incomplete virialization, streams and interaction with the stellar diskcould contribute to create anisotropies or a non Maxwellian spectrum. Moreover, evenin a completely virialized structure at equilibrium, anisotropy is expected as a result ofthe inhomogeneous gravitational potential with a power law density profile [53, 54]. Theconsequences for direct detection of some of these departures from the standard halo casehave been analyzed by several authors [55–62].Numerical simulations have been used for decades as a powerful virtual laboratory toexplore complex dynamical systems. On cosmological scales, early studies enabled to assessthe role of DM in the hierarchical structure formation sheme [63,64]. Recent techniques [65,66] as well as improving computing capabilities have brought the possibility to extractgalactic DM halos with tremendous resolution. The so-called Via Lactea simulation [67]or the
Aquarius project [68] have resolved a Milky-Way sized galactic halo with morethan a billion particles. It appears that numerous sub-halos are present, anisotropies anddeviations from standard Maxwellian velocity distributions are also found and lead tosomewhat different signatures in direct detection experiments [43, 53, 54, 62, 69]. Despitetheir impressive resolution, these DM only simulations fail as a realistic description of agalactic halo for the simple reason that they totally neglect the baryonic components (starsand gas in the galactic disk and bulge), although these dominate are known to dominatethe dynamics near the galactic center. This shortcoming is being fixed by very recentcosmological simulations which include baryons [70, 71]. It appears that the DM spatial– 2 –istribution is neither spherical, or triaxial, but has a thick oblate shape around the galacticdisk, in what some authors have called a dark disk [72, 73], that is co-rotating with thegalactic stellar disk. These particular characteristics lead to a potential enhancement ofthe direct detection signal at low recoil energies [73], although there is an important spreadin the predictions, depending on the particular merger history for each galaxy.In this paper, a recent and advanced cosmological simulation with baryons is used as arealistic framework to extract detailed predictions for direct detection. The paper is orga-nized as follows. In Sec. 2, the simulation is presented and described. Velocity distributionsin the solar neighborhood are extracted, analyzed and compared with standard simplifiedassumptions. The presence of a co-rotating dark disk is also discussed. In Sec. 3, after abrief overview of the event rate formalism and the experimental status, the direct detectionpredictions are given for both the elastic and the inelastic scenarios. Finally, our resultsare summarized in Sec. 4.
2. A cosmological simulation with baryons
Although galaxy formation is far from being completely understood, both theory and ob-servation have made significant progresses in the recent years. Based on first principles, itis now possible to perform simulations of Milky-Way–like galaxies in reasonable agreementwith observed galaxies of similar circular velocities [74–76]. Some long-lasting problemsseem however to remain, especially for massive galaxies like our own: the simulated bulgeseems systematically too massive (sometimes referred to as the angular momentum prob-lem), and the total baryon fraction too large (the overcooling problem). Following previouswork initiated by Bruch et al. [73], we nevertheless would like to use a Milky–Way–likegalaxy simulation including both dark matter and baryons dynamics to infer the expectedsignal in dark matter direct detection experiments.
The simulation was performed using the cosmological Adaptive Mesh Refinement (AMR)code RAMSES [77]. We follow the evolution of 2 collisionless fluids, namely dark matterand stars, coupled through gravity to the dissipative baryonic component. Gas dynamics isdescribed using a second-order unsplit Godunov scheme. Gas is allowed to cool radiatively,based on a standard equilibrium photo-ionised mixture of Hydrogen and Helium [78], takeninto account the excess cooling due to metals [79]. Star formation is implemented using astandard Schmidt law, spawning stars as a Poisson random process [80]. Our rather stan-dard recipe for supernovae feedback and the associated chemical enrichment is describedin Ref. [81].We used the so-called “zoom-in” simulation technique, for which we identify a candi-date halo using a low resolution, dark matter only simulation, and then re-simulate it athigher resolution, degrading both spatial and mass resolutions as we move further awayfrom the central region. We checked that our final halo is not contaminated by high mass,low resolution dark matter particles. Our effective resolution in the initial conditions cor-responds to a Cartesian grid of 1024 elements covering our 20 Mpc/ h periodic box. The– 3 – igure 1: View of the simulation box (top), with a zoom on the central and better resolved region(bottom). The size of the simulation box is 20 Mpc /h . The central galaxy is embedded in a largescale filamentary structure. corresponding dark matter particle mass is therefore 7 . × solar masses, given we usedthe following cosmological model Ω m = 0 . , Ω Λ = 0 . , Ω b = 0 . , H = 70 km/s/Mpc.The initial Gaussian random field was generated according to the Eisenstein and Hu (1999)ΛCDM model [82] for the transfer function.At redshift z = 0, the complete simulation volume contains 8,376,152 dark matterparticles of various masses, 11,702,370 star particles and 19,677,225 AMR cells. We havefixed the spatial resolution to be 200 pc in physical units , which ensures that the verticalscale height of the stellar disc is barely resolved. The final halo Virial radius R vir isfound to be 264 kpc. It contains N DM = 842 ,
768 dark matter particles corresponding The Virial radius R vir is defined as a function of the cosmological matter density ρ m by M ( R Face-on and edge-on distributions of the baryonic material in our galaxy. The whitevertical line in the upper left panel has a length of 10 kpc. Note that both the stars (upper panels)and the gas (lower panels) have formed a thin disk, in good agreement with observations. A numberof satellites can also be seen, more numerous in the stars than in the gas, the gas having beenexpelled from a large fraction of the substructures. The filamentary structures in the gas edge-onview outlines past satellite trajectories. Reducing the bulge mass by a factor of 2 would result in a 10% decrease of the circularvelocity and in a 20% decrease of the dark matter density in the solar neighbourhood (seeRef. [84] for a quantitative estimate).The central dark halo, stars and gas are shown respectively on figure 2 and 3. Figure 2shows the dark matter distribution. As expected, it is centrally concentrated and has aconsiderable number of subhaloes. Using Adaptahop to extract the substructure informa-tion, we find that the central galaxy contains 92 subhaloes, of which 79 are inside the virialradius. The subhalo mass fraction is 4.8%, a rather small value due to the limited mass res-olution in the simulation. The closest clump to the galactic center is at a radial distance of12.7 kpc. Therefore, as no clump lies in the the solar neighborhood region (7 < R < − < z < In this section, we extract and analyze the velocity distributions given by our N-bodysimulation, and compare them with the usual Gaussian/Maxwellian assumption. As thegoal of this study is to calculate direct detection predictions with a realistic DM halo, theemphasis will be put on the local phase-space structure, in the solar neighborhood.The Sun being located at a distance R (cid:39) < R < N shell = 16 , σ ij . They are shown on Fig. 4.The Gaussian fit for each component and the Maxwellian fit for the velocity module areshown in red. A clear departure from these distributions is seen.In order to model and parameterize these deviations, one can introduce generalizedGaussian and Maxwellian distributions, which are convenient fitting functions. The gen-eralized Gaussian distribution with a mean velocity µ , a velocity dispersion parameter v ,and a Gaussianity parameter α is defined byf( v ) = 1 N ( v , α ) e − ( ( v − µ ) /v ) α , (2.1)with a normalization factor N ( v , α ) = 2 v Γ(1 + 1 / α ). For α = 1, the usual Gaussiandistribution is recovered. Similarly, the generalized Maxwellian distribution can be definedby f( (cid:126)v ) = 1 N ( v , α ) e − ( (cid:126)v /v ) α , (2.2)with a normalization factor N ( v , α ) = 4 πv Γ(1 + 3 / α ) / 3. For α = 1, we get the usualMaxwellian distribution. – 7 – (cid:61) (cid:45) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40000.51.1.52.2.5 v (cid:72) km (cid:144) s (cid:76) f (cid:72) v (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Μ (cid:61) (cid:45) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40000.51.1.52.2.5 v (cid:72) km (cid:144) s (cid:76) f (cid:72) v (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Μ (cid:61) (cid:45) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40000.51.1.52.2.5 v (cid:72) km (cid:144) s (cid:76) f (cid:72) v (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) d (cid:76) Μ (cid:61) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) v (cid:61) (cid:61) v (cid:61) Α (cid:61) (cid:61) 332 ; Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) Χ (cid:144) dof (cid:61) (cid:144) Χ (cid:144) dof (cid:61) (cid:144) Χ (cid:144) dof (cid:61) (cid:144) 480 100 200 300 400 500 60001.2.3.4. v (cid:72) km (cid:144) s (cid:76) f (cid:72) v (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Figure 4: Velocity distributions of dark matter particles ( N shell = 16 , ) in a spherical shell < R < kpc around the galactic center.a), b) and c) Velocity components along the principal axes of the velocity dispersion tensor, togetherwith the Gaussian (red) and a generalized Gaussian (green) distribution fits (cfr. Eq. (2.1)).d) Velocity module, with Maxwellian (red), Tsallis (green) and generalized Maxwellian (orange andpurple) fits (cfr. Eqs. (2.2,2.3)). µ , σ (both in km/s) and K stand for the mean, the standard deviation and the Kurtosis parameterof the distribution. The goodness of fit is indicated by the value of the χ vs. the number of degreesof freedom (dof ). For the velocity components, a fit with a generalized Gaussian distribution (in greenin Fig. 4, panels a), b) and c)) shows that the velocity distribution given by the simulationis systematically more flat than a Gaussian distribution. This property can be describedby the so-called Kurtosis parameter K , which compares the fourth-order moment with thesquare of the variance. A Kurtosis parameter K < i.e. more flat than a Gaussian distribution with the same standard deviation,– 8 – (cid:61) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40000.51.1.52.2.5 v r (cid:72) km (cid:144) s (cid:76) f (cid:72) v r (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Μ (cid:61) Σ (cid:61) (cid:61) Μ (cid:61) (cid:61) 250 , f (cid:61) Μ (cid:61) (cid:61) 120 , f (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40000.51.1.52.2.53. v Φ (cid:72) km (cid:144) s (cid:76) f (cid:72) v Φ (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Μ (cid:61) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40000.51.1.52.2.53. v z (cid:72) km (cid:144) s (cid:76) f (cid:72) v z (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) d (cid:76) Μ (cid:61) Σ (cid:61) (cid:61) Α (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) Χ (cid:144) dof (cid:61) (cid:144) 470 100 200 300 400 500 60000.51.1.52.2.53.3.54.4.5 v (cid:72) km (cid:144) s (cid:76) f (cid:72) v (cid:76) Figure 5: Velocity distributions of dark matter particles ( N ring = 2 , ) in a ring < R < kpc, | z | < kpc around the galactic plane.a) Radial velocity v r , with Gaussian (red) and generalized Gaussian (green) fits (cfr. Eq. (2.1)).b) Tangential velocity v φ , with a double Gaussian fit. f indicates the fraction of each component.c) Velocity across the galactic plane v z , with Gaussian (red) and generalized Gaussian (green) fits(cfr. Eq. (2.1)).d) Velocity module, with Maxwellian (red) and a generalized Maxwellian (green) fit (cfr. Eq. (2.2)). µ , σ (both in km/s) and K stand for the mean, the standard deviation and the Kurtosis parameterof the distribution. The goodness of fit is indicated by the value of the χ vs. the number of degreesof freedom (dof ). while K > i.e. more peaked around thecentral value compared to the Gaussian one. So velocity distributions extracted from thissimulation are systematically platykurtic, in contrast with results obtained in DM onlysimulations [69, 87].The goodness-of-fit for each theoretical distribution is calculated using the Pearson’s– 9 – test. The values of χ (calculated with a number of bins n bin = 50) printed in eachpanel of Fig. 4 show that the fit is considerably improved when using a generalized Gaus-sian distribution with α > v and α are estimated on the basis ofthe standard deviation σ and the Kurtosis parameter K extracted from the velocity distri-bution). We also notice that along the principal direction “3”, the velocity dispersion andthe deviation from the Gaussian distribution are smallest. It appears that this direction isactually very close to the galactic pole axis.For the velocity module, a clear deviation from a Maxwellian distribution (in red) isapparent, with a sharper drop of the high velocity tail, particularly crucial for the inelasticrecoil scenario (see Sec. 3.5). There is no evidence of bumps in the velocity distribution, asseen in Ref. [69]. A fit of the distribution with a generalized Maxwellian distribution with α > α enablesto model the simulation distribution in a satisfactory way. With two free parameters,only the mean value and the standard deviation can be accommodated, while the Kurtosisparameter is derived. It appears that the N-body velocity distribution is “fatter” than thebest-fit generalized Maxwellian ( v = 301 . α = 1 . K = 2 . 39 compared to K = 2 . 71 for the fit. In panel d) of Fig. 4, we showtwo possible fits with a generalized Maxwellian distribution, corresponding to α = 1 . α = 1 . 95. The former gives a better global fit for the whole distribution and for the lowvelocity bins, the latter gives a better fit of the high velocity tail.In fact, the sharp drop of the high velocity tail has been recognized as a universalbehavior of relaxed collisionless structures [54]. The presence of long-range gravitationalforces in dark matter structures indicates that these systems should be described by non-extensive statistical mechanics. The generalization of the usual Boltzmann-Gibbs entropyto non-extensive systems by Tsallis [88] leads to distribution functions of the formf( (cid:126)v ) = 1 N ( v , q ) (cid:18) − (1 − q ) (cid:126)v v (cid:19) q/ (1 − q ) , (2.3)where N ( v , q ) is a normalization constant. The Maxwell-Botzmann distribution is recov-ered by taking the limit q → 1. Equilibrated self-gravitating collisionless structures havebeen shown to exhibit Tsallis distributions both analytically [89], and numerically [54, 90].For the particles in the spherical shell around 8 kpc in this simulation, the velocity moduledistribution is indeed very well fit by a Tsallis distribution with v = 267 . q = 0 . K = 2 . β = 1 − σ t /σ r . Also, the departure from a Maxwellian distribution,parameterized by 1 − q for Tsallis distributions could show the same anisotropy. In thissimulation, for the DM particles in the spherical shell 7 < R < β (cid:39) . 12, together with a slightly smaller Kurtosis parameter in the radialdirection K = 2 . 41, in general agreement with the conclusions drawn from the universal– 10 – (cid:61) Σ (cid:61) (cid:61) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:61) Α (cid:61) Χ (cid:144) dof (cid:61) (cid:144) (cid:76) (cid:45) (cid:45) 200 0 200 40001.2.3.4.5.6.7.8. v r (cid:72) km (cid:144) s (cid:76) f (cid:72) v r (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Μ (cid:61) Σ (cid:61) (cid:61) Μ (cid:61) (cid:61) 30 , f (cid:61) Μ (cid:61) (cid:61) 72 , f (cid:61) (cid:76) (cid:45) (cid:45) 200 0 200 40001.2.3.4.5.6.7.8.9.10. v Φ (cid:72) km (cid:144) s (cid:76) f (cid:72) v Φ (cid:76)(cid:72) (cid:180) (cid:45) (cid:76) Figure 6: Velocity distributions of star particles ( N star = 143 , ) in a ring < R < kpc, | z | < kpc around the galactic plane.a) Radial velocity v r with Gaussian (red) and generalized Gaussian (green) fits (cfr. Eq. (2.1)).b) Tangential velocity v φ with a double Gaussian fit, with a thin (red) and a thick (green) diskcomponents. f indicates the fraction for each component. µ , σ (both in km/s) and K stand for the mean, the standard deviation and the Kurtosis parameterof the distribution. The goodness of fit is indicated by the value of the χ vs. the number of degreesof freedom (dof ). relation between the density slope and the velocity anisotropy presented in Ref. [53]. Tocalculate the velocity anisotropy, the radial velocity dispersion σ r is directly extracted fromthe radial velocity distribution, while the tangential velocity dispersion σ t can be deducedfrom the mean µ and the standard deviation σ of the distribution of the velocity moduleas σ r + 2 σ t = µ + σ .In dark matter-only simulations, one has a 4 π uncertainty on the Sun position. Thisfreedom is lifted, at least partially, in simulations with baryons where the galactic planecan be identified. To determine it, we select star particles with 3 < R < 10 kpc to avoidthe galactic bulge. We have N star = 1 , , 104 in this selection. The galactic plane is thenrotated onto the xy plane by diagonalizing the position tensor. In this new reference frame,denoted by { (cid:126) x ,(cid:126) y ,(cid:126) z } , we select DM particles with 7 < R < | z | < N ring = 2 , { r, φ, z } . Anisotropy is manifest. There is a significant rotationalong the galactic disk, with (cid:104) v φ (cid:105) = 35 . z velocity components,the mean of the distribution is compatible with zero. The velocity dispersions and Kurtosisparameters also exhibit some anisotropy in the z direction compared to the galactic plane.For the distribution of the velocity module, a strong deviation from a Maxwellian distri-bution is again visible, with a cut-off of the high velocity tail. However, despite a coarser– 11 – (cid:60) R (cid:60) (cid:45) z (cid:72) kpc (cid:76) (cid:60) Ρ D M (cid:62) (cid:72) G e V (cid:144) c m (cid:76) xyyzxz R (cid:72) kpc (cid:76) Ρ D M d i s k (cid:144) Ρ D M s h e ll Figure 7: (Left) Average density of DM particles with < R < kpc as a function of the heightfrom the galactic disk z (R is the spherical radius to the galactic center). The dashed line givesthe average value for the entire spherical shell. To select particles in z slices, we used a thickness δz = 2 kpc.(Right) Ratio of ring to shell densities as a function of distance from the galactic center for differentplanes. The ratio fluctuates around 1.2 for the galactic plane (blue), while it drops to a value ∼ . for other planes (green, magenta). For the plane yz , the sudden peak at R (cid:39) kpc is due to thepresence of a satellite halo, visible on Fig. 8.b. resolution, we can observe that the distribution for ring particles differs significantly fromthe Tsallis form that is expected for equilibrated systems. Therefore, we can already stressthat predictions of direct detection signals that use equilibrium distributions should betaken with some caution. As pointed out in Ref. [72], the rotation of the DM halo observed in Fig. 5 is expected inΛCDM cosmology with hierarchical structure formation. At relatively low redshift, z < dark disk . The verticaldistribution of the stars of the host galaxy thickens due to the impact, but all new starsform in a thin disc [91–93]. Observational evidence for a thick stellar disk in the MilkyWay has been first pointed out by Gilmore and Reid (1983) [94]. Although the detailedmorphology and some quantitative properties of the final galaxy depend on its mergerhistory, the existence of a thick stellar disk and a co-rotating dark disk in the ΛCDMparadigm seems quite robust. Using three cosmological hydrodynamics simulations, Read et al. [71] conclude that the accreted dark matter can contribute ∼ . − . v φ distribution is best-fitwith a double Gaussian, one for the non-rotating halo component, and one for the rotating– 12 – (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) y (cid:72) kpc (cid:76) z (cid:72) k p c (cid:76) Figure 8: Density maps of the dark matter halo in the planes a) xy (galactic plane), b) yz .Contours correspond to ρ DM = { . , . , . , . } GeV / cm . Sun0 20 40 60 80 1000102030405060 R (cid:72) kpc (cid:76) (cid:60) v Φ (cid:62) (cid:72) k m (cid:144) s (cid:76) (cid:45) (cid:45) (cid:45) 10 0 10 20 30 (cid:45) z (cid:72) kpc (cid:76) (cid:60) v Φ (cid:62) (cid:72) k m (cid:144) s (cid:76) Figure 9: Mean tangential velocity (cid:104) v φ (cid:105) as a function of R for z = 0 (left panel), and as a functionof z for r = 8 kpc (right), where r = (cid:112) x + y is the polar radius. accreted DM. The rotation lag of the rotating component is in the range 0 − 150 km/s. Alarge value of the rotation lag, corresponding to a small halo circular speed, is only foundfor galaxies which had no significant merger after a redshift z = 2. The importance ofthe dark disk for the prediction of DM direct detection rates has been acknowledged byseveral authors [71, 73]. It could lead to an increase by up to a factor of 3 in the 5 − 20 keVrecoil energy range. The signal modulation can also be boosted and the modulation phaseis shifted. However, in Ref. [95], the authors use high-resolution simulations of accretionevents in order to bracket the range of co-rotating accreted dark matter. They find that– 13 –he Milky Way’s merger history must have been unusually quiescent compared to medianΛCDM expectations. As a result, the fraction of accreted dark matter near the Sun isless than 20% of the host halo density, and does not change the likelihood of detectionsignificantly.To investigate this topic further, we select star particles with 7 < R < | z | < N star = 143 , v r and v φ are shown onFig. 6. We observe that the dark matter and the star particles are indeed co-rotatingin the solar neighborhood. The mean tangential velocity is (cid:104) v φ (cid:105) = 201 km/s but tendstowards (cid:104) v φ (cid:105) = 225 km/s for stars closer to the galactic plane, which is consistent withMilky Way rotation curve data [47]. Moreover, the v φ distribution is clearly bimodal, adouble Gaussian gives a reasonable fit. The low velocity component ( µ = 175 km/s) hasa larger velocity dispersion parameter ( v = 72 km/s) than the high velocity component( µ = 270 km/s, v = 30 km/s). To check whether these correspond to the thick and thindisks populations, we take subsets of the previous selected star particles with 0 < v φ < 150 km/s and v φ > 275 km/s respectively, and calculate the disc scale height as well as thetotal velocity dispersion. We obtain z thickd (cid:39) z thind (cid:39) . 67 kpc, σ thicktot (cid:39) . σ thintot (cid:39) . z thickMW (cid:39) z thinMW (cid:39) . 35 kpc [96]; for the velocitydispersion, σ thickMW (cid:39) 85 km/s [97], σ thinMW (cid:39) 40 km/s [98]. Although the thick and thinstellar disks could be revealed by applying severe cuts on v φ , the inferred values of the diskscale height and the total velocity dispersion are not reliable, they probably overestimatethe correct values by up to a factor 2. Hence, we stay inconclusive regarding the consistencyof stellar populations in this simulation with Milky Way observations.The existence of two DM components in the solar neighborhood cannot be assessedundoubtfully in this simulation due to the limited local resolution. Nevertheless, as shownon panel b of Fig. 5, a fit with a double Gaussian gives a better description of the v φ distribution than with a single Gaussian. From this fit, it appears that the mean velocityof the rotating component is comparable for DM and for stars. The dark disk particles alsohave a smaller velocity dispersion. also, it appears that the dark disk constitutes around25% of the total local density in the solar neighborhood. The value of this fraction canbe further checked by calculating average densities from the simulation. For DM particlesin the spherical shell 7 < R < (cid:104) ρ DM (cid:105) = 0 . 31 GeV / cm , in agreementwith what is found in DM only simulations and the commonly adopted value in this field, ρ DM = 0 . / cm . However, when we restrict DM particles to a ring of thickness δz around the galactic disk, the average density is higher. We find (cid:104) ρ DM (cid:105) = 0 . 37 GeV / cm for δz/ (cid:104) ρ DM (cid:105) (cid:39) . 39 GeV / cm for δz/ ≤ . < R < δz = 2 kpc parallel to theplane xy (galactic plane), and the average density is shown as a function of z . If the halowas close to spherical, there would be no bump around z = 0. On the right panel, weshow the disk to shell density ratio for different planes. For the galactic plane, the densityincrease due to the dark disk is typically between 15 and 30% for R ≤ 25 kpc. For other– 14 –lanes, the ratio oscillates around 0.9, except when satellite halos are encountered, this isthe case for the plane yz and R (cid:39) 13 kpc, as illustrated by the density maps of Fig. 8.With the dark disk component, the entire DM halo appears with an oblate shape. Thelocal density increase due to the dark disk is relatively mild in this simulation, comparedto more extreme cases [71]. Therefore, following the results of Ref. [73], we do not expecta strong increase of the direct detection signal in this simulation compared to a standardMaxwellian halo.Finally, we checked that the DM rotation is not limited to a thin galactic plane but ispresent in a rather thick dark disk. In particular, Fig. 9 shows the variation of the meantangential velocity vφ as a function of R for z = 0 (left panel), and as a function of z for r = 8 kpc (right), where r = (cid:112) x + y is the polar radius. A significant rotation isfound in the entire galactic plane, up to R = 100 kpc. For r = 8 kpc, the rotation remainsimportant up to heights of ∼ 15 kpc from the galactic plane. 3. Direct detection experiments and signals Direct detection experiments aim to measure the energy deposited during scatterings ofWeakly Interacting Massive Particles (WIMPs) with the detector material. Given the weakinteraction rate and the low expected signal, radiopurity of the material and control of otherbackgrounds are essential although experimentally challenging. In order to attenuate thesedifficulties, one possibility is to look for a typical signature of the signal that is hard tomimic with a background noise. The DAMA/Libra experiment [41] is the only one thathas observed a positive signal with such a signature, namely the annual modulation of theevent rate. However, as the signal could not be confirmed by any other experiment, neitherwith a heavier or a lighter target, the compatibility of the DAMA result and the exclusionlimits set by these other experiments becomes more and more controversial [99, 100].Ways of reconciliation have been sought in the detector characteristics (quenchingfactors uncertainties and channelling effects [101, 102]), in particle physics uncertainties(nuclear form factors [103], nucleon effective couplings [104]) or scenarios (elastic scat-tering [105–109], inelastic scattering [62, 110–114]), and in the Dark Matter halo velocitydistribution [43, 44, 62, 115–117].In this section, we give all the relevant steps necessary to compute the spin independentscattering signal of Dark Matter on nuclei in direct detection experiments. Reference valuesfor all the useful parameters and quantities will be given, and discussed in light of theDAMA controversy. In this paper, the accent will be put on the impact of a realistic velocitydistribution (as given by our cosmological simulation with gas and stars) on the directdetection signal and its modulation. In particular, we will focus on how the fit to DAMAand the other exclusion limits change compared to a standard Maxwellian halo. Both theelastic and the inelastic cases will be covered. For the influence of other parameters, thereader will be referred to the existing literature. We consider collisions between dark matter particles from the Galactic halo and the nuclei– 15 –f a given low background detector. The relevant characteristics of the halo are the localdensity of dark matter, taken to have the fiducial value ρ DM = 0 . at the Sun’slocation, and the local distribution of velocities f( (cid:126)v ) with respect to the Earth.The differential event rate of nuclear recoils as a function of the recoil energy E R isconveniently factored as d R dE R = ρ DM M DM dσdE R η ( E R , t ) (3.1)where M DM is the WIMP mass, dσ/dE R encodes all the particle and nuclear physicsfactors, and η ( E R , t ) is the mean inverse velocity of incoming DM particles that can deposita recoil energy E R . The time dependence of the velocity distribution is induced by themotion of the Earth around the Sun, which leads to a seasonal modulation of the eventrate [118, 119].The total recoil rate per unit detector mass in a given energy bin [ E , E ] is obtainedby integrating Eq. (3.1), R ( t ) = (cid:90) E E dE R (cid:15) ( E R ) (cid:18) d R dE R (cid:63) G ( E R , σ ( E R )) (cid:19) , (3.2)where (cid:15) is the efficiency of the detector and the finite energy resolution of the experimentis taken into account by convoluting the differential rate with a Gaussian distribution withspread σ ( E R ). For detectors made of several elements, the total rate is the average of therates R i ( t ) for each component i , weighted by its mass fraction f i R ( t ) = (cid:88) i f i R i ( t ) (3.3)Finally the expected number of observed events per unit time is the product of the totalrate times the detector mass M det .The particle physics is enclosed in the term dσ/dE R , which is generally parameterizedas dσdE R = M N σ n µ n (cid:16) Zf p + ( A − Z ) f n (cid:17) f n F ( E R ) , (3.4)where M N is the nucleus mass, µ n is the reduced WIMP/neutron mass, σ n is the zeromomentum WIMP-neutron effective cross-section, Z and A are respectively the number ofprotons and the atomic number of the element, f p ( f n ) are the WIMP effective coherentcoupling to the proton (resp. neutron), and F ( E R ) is the nuclear form factor. In thesequel, we will consider the case of scalar interactions, most relevant for the elastic scat-tering scenario, for which f p (cid:39) f n and the case of weak vector interactions through a Z boson, most relevant for the inelastic scattering scenario, for which f p = 4 sin θ W − f n = 1. The form factor F ( E R ) characterizes the loss of coherence for non zero momentumtransfer. We use the simple parameterization given by Helm [103, 120], defined as F ( E R ) = 3 e − q s / J ( qr ) qr , (3.5)– 16 –ith q = √ M N E R , s = 1 . r = { (1 . m N /m n ) / − . + 0 . π − s } / fm theeffective nuclear radius, J is a spherical Bessel function of the first kind. This form factoris optimal for scattering on Iodine [103]. For simplicity we use the same form factor forall targets (more accurate form factor are off by at most O (20%), for some targets and atlarge recoil energies [121]).Finally, the velocity distribution appears in the quantity η ( E R , t ) = (cid:90) v min d (cid:126)v f( (cid:126)v ( t )) v , (3.6)where (cid:126)v the WIMP velocity wrt the Earth and v min is the minimum velocity needed toprovoke a recoil inside the detector. Two cases need to be considered. In the elasticscattering scenario, a dark matter particle is simply scattered off a nucleus. In the inelasticscattering scenario, a dark matter particle DM is supposed to scatter into a slightly heavierstate DM , with a mass splitting δ = M DM − M DM ∼ 100 keV. In this scenario, whichhas been first proposed in Ref. [110] and confronted to the recent data in Refs. [62,111–114],a much broader range of dark matter candidates may both fit DAMA and be consistentwith the other experiments. The threshold velocity is given by v min = (cid:114) M N E R (cid:16) M N E R µ + δ (cid:17) , (3.7)with µ the WIMP-nucleus reduced mass. This formula encompasses both the elastic ( δ =0) and the inelastic ( δ (cid:54) = 0) scattering scenarios.When the velocity distribution in the galactic frame f gal ( v ) is isotropic, η is given by η = 2 πv ⊕ (cid:90) v + v − (cid:0) F( v esc ) − F( v ) (cid:1) dv , (3.8)with v ± = min { v esc , v min ± v ⊕ } , v esc is the galactic escape velocity, (cid:126)v ⊕ ( t ) is the Earth’svelocity in the galactic frame, which varies with the time of year t , and F( v ) = (cid:82) v f gal ( v ) dv .In deriving this formula, we have used the realistic assumption v esc > v ⊕ at any time t .Notice that F( v ) is an even function of v , because f gal only depends on the modulus | (cid:126)v | ofthe velocity. For a standard Maxwellian distribution with a mean velocity v ,f gal ( (cid:126)v ) = 1 N ( v esc ) π / v e − v /v for v < v esc , (3.9)where N ( v esc ) is a normalization factor, Eq. (3.8) reads η = 12 N ( v esc ) v ⊕ (cid:26) Erf (cid:16) v + v (cid:17) − Erf (cid:16) v − v (cid:17) − v + − v − ) √ π v e − v esc /v (cid:27) , (3.10)which agrees with the result of Ref. [122]. It might seem curious that Eq. (3.8) is expressedas a definite integral between v − and v + while particles with any velocity above v min should contribute to η . However, the function F ( v ) is itself an integral on the velocitydistribution. Also the dependence in v esc is now explicit. In the rest of this paper, we willassume an escape velocity v esc = 600 km/s, somewhat on the upper part of the expected– 17 – 00 150 200 250 300146148150152 v c (cid:72) km (cid:144) s (cid:76) t (cid:72) da y s (cid:76) Figure 10: Variation of the phase t corresponding to a maximal Earth velocity and event rate asa function of the galactic circular velocity v c . range, 498 km/s < v esc < 608 km/s [117]. Modifications of the velocity distributions andtheir impact on the fits, have been discussed in various works, see for instance Refs. [43,44, 62, 113].The seasonal modulation of the event rate is caused by the variation of the Earth’svelocity throughout the year. Following Refs. [57, 115], in the galactic reference frame suchthat (cid:126) X is pointing towards the galactic center, (cid:126) Y is along the disk rotation, and (cid:126) Z ispointing towards the galactic north pole, we write (cid:126)v ⊕ ( t ) = (cid:126)v (cid:12) + (cid:126)v EO ( t ) ,(cid:126)v (cid:12) = (cid:126)v c + (cid:126)v (cid:12) ,pec = (0 , , 0) + (10 , , 7) km / s ,(cid:126)v EO ( t ) = v EO ( (cid:126)e sin ω ( t − t ) − (cid:126)e cos ω ( t − t )) , (3.11)where the velocity (cid:126)v (cid:12) of the Sun with respect to the halo [123] is decomposed as the sumof the galactic disk circular velocity (cid:126)v c and a peculiar velocity (cid:126)v (cid:12) ,pec , the Earth meanorbital velocity is v EO = 29 . ω = 2 π/T and a period T = 1 year, t =79 . 62 days is the time of the vernal equinox, and (cid:126)e = ( − . , . , − . (cid:126)e =( − . , − . , . v EO (cid:28) v (cid:12) , we have v ⊕ ( t ) (cid:39) v (cid:12) + v EO sin γ cos ω ( t − t ) , (3.12)with γ (cid:39) . ◦ is the angle between the ecliptic plane perpendicular axis and (cid:126)v (cid:12) . TheEarth velocity is maximal for t = t , with t = t + π ω + 1 ω arctan (cid:126)v (cid:12) · (cid:126)e (cid:126)v (cid:12) · (cid:126)e (cid:39) . . (3.13)With these values, the Earth velocity modulation amplitude amounts to 6 . v c = 220 km / s– 18 –sually considered. Also, if the DM halo is co-rotating with the galactic disk, the relativevelocity relevant for calculating the direct detection signal is effectively reduced. In Fig. 10,the circular velocity is varied between 100 km/s and 300 km/s, the corresponding phase t increases with v c , from t = 144 to t = 153 days.For an isotropic velocity distribution, the differential rate also follows a temporal cosinefunction at first order in v EO /v (cid:12) . Indeed, the only time dependent factor in Eq. (3.1) is η ,for which we get, in the limit v esc > v + η ( E R , t ) (cid:39) η ( E R ) + η ( E R ) v EO v (cid:12) sin γ cos ω ( t − t ) , (3.14)with η ( E R ) = η ( E R , v EO = 0) ,η ( E R ) = 2 πv (cid:12) (cid:90) v v − F ( v ) dv − π (cid:16) F ( v ) + F ( v − ) (cid:17) , (3.15)with v ± = v min ( E R ) ± v (cid:12) . It should be stressed that the time of year t = t does notnecessarily correspond to a maximum of the event rate. For low energy events, Eq. (3.15)shows that it is a minimum instead, as η becomes negative.In this paper, we also need to consider discrete velocity distributions as the output ofa simulation gives a list of particles with their position and velocity. If we denote by N part the number of particles, and by (cid:126)v i the velocity of the particle i in the galactic frame, thedistribution can be written asf gal ( (cid:126)v ) = 1 N part (cid:88) i δ ( (cid:126)v − (cid:126)v i ) , (3.16)where δ is the Dirac delta distribution in 3 dimensions. Therefore, for η , we get in thiscase η = 1 N part (cid:88) i w i H ( w i − v min ) , (3.17)with w i = | (cid:126)v i − (cid:126)v ⊕ | and H ( x ) is the Heaviside step function. To compute the total rate, orthe modulation amplitude in some energy bin, we need to integrate the differential eventrate. For discrete distributions, we take advantage of the step-like shape of the differentialevent rate for a given flow (cid:126)v i to get accurate results. The former DAMA/NaI [36] and its follow-up DAMA/LIBRA [41] are experiments runat the LNGS at Gran Sasso, Italy using NaI(Tl) crystals as target. They are designed todetect the dark matter recoil off nuclei through the model independent annual modulationsignature which is due to the motion of the Earth around the Sun. The experimental resultsobtained by DAMA/LIBRA, with an exposure of 0.53 ton × yr collected over 4 annual cycles,combined with the ones of DAMA/NaI, for an exposure of 0.29 ton × yr collected over 7annual cycles, corresponding to a total exposure of 0.82 ton × yr, show a modulated signalwith a confidence level of 8.2 σ [41]. – 19 – (cid:45) Eee (cid:72) keV (cid:76) S m , k (cid:72) c pd (cid:144) k g (cid:144) ke V (cid:76) Figure 11: DAMA modulation amplitude data (in counts per day per kg per keV) as a function ofthe recoil energy, expressed in electron-equivalent keV units. The 36 bins, from 2 to 20 keV, showa modulation signal in the low energy part, from 2 to 6 keV. The upper part is compatible with zeromodulation, as shown in the two bins version of the data (from 2 to 6 keV and from 6 to 14 keV). The DAMA modulation data, reproduced on Fig. 11, are obtained by adjusting (witha likelihood function) the observed event rate in each energy bin k to a temporal cosinefunction S k = S ,k + S m,k cos ω ( t − t ) , (3.18)where S ,k is the constant part of the signal, S m,k is the modulation amplitude. The DAMAcollaboration uses the value t = 152 . t depends in general upon the velocity of the Sun, and upon the DM velocitydistribution. However, for an isotropic distribution, Eqs. (3.13-3.14) show that the phase t is completely determined once the velocity of the Sun in the galaxy is specified. For thelow energy region where the signal modulation is present, the experimental value of t hasbeen determined by the DAMA collaboration. When the phase is allowed to vary freely,and fitting again the event rate in the energy interval 2 ≤ Eee ≤ t is t = 144 . ± . σ ) . (3.19)The DAMA spectrum is given in keVee (electron-equivalent keV). The observed energyreleased in scintillation light is related to the nucleus recoil energy through a so-calledquenching factor q , E scint. = q · E recoil. This expresses the fact that a recoiling nucleusmay loose energy by collisions with other nuclei, hence in the form of heat, or throughcollisions with electrons, which create scintillation light. The reference values for Iodine andSodium are respectively q I = 0 . 09 and q Na = 0 . 3. However it has been pointed that the so-called channelled events may play a role [101,102]. This refers to events in which a recoilednucleus moves along the axis of the NaI crystal, losing most of its energy by collisions withelectrons, in which case the quenching factor may be larger, up to q ≈ 1. Once channelling– 20 –s taken into account, collisions of light WIMPs with Iodine become relevant, while recoilson Sodium become negligible in all scenarios [42, 43, 102, 122, 125, 126]. We use the fraction f of channelled events given in Ref. [43], f Na ( E R ) = e − E R / . E R , f I = e − E R / . E R , (3.20)where the recoil energy E R is in keV. We have verified that the other parameterizationsfound in the literature [42, 112, 122] give identical results. Finally, we use the energyresolution of the detector given by the collaboration [41] σ ( E ) E = 0 . √ E + 0 . (cid:15) = 1 [41, 122].To fit the observed energy spectrum of Fig. 11 with a theoretical model, we use agoodness-of-fit (GOF) method with a χ metric, for the 12 first bins of data, with width0 . . . χ = n (cid:88) i =1 ( S i − S obsi ) σ i (3.22)where S i is the predicted value of the modulation amplitude in the bin number i , S obsi isthe value reported by DAMA and σ i is the corresponding experimental uncertainty. Asdata above 8 . χ fit enlarges the allowed regions in the parameter space. For explicit numericalvalues of S obsi , we refer the reader to Table III of Ref. [122]. As emphasized earlier, thesevalues were obtained by the DAMA collaboration with the hypothesis t = 152 . 5. To beconsistent with their analysis and to facilitate comparisons with previous analyses, we willalso use this fixed value of t for the determination of the allowed regions in the parameterspace, although the true phase corresponding to the maximal event rate might be differentfor some energy bins. For two parameters ( M DM and σ in the elastic scenario, and δ and M DM in the inelastic scenario), the χ metric has 10 degrees of freedom, therefore the 90%(99% and 99 . χ < . χ < . χ < . So far all the other direct detection experiments searching for dark matter are compatiblewith null results. In this section we briefly describe the experiments that lead to the mostconstraining limits on both the elastic and the inelastic scenarios.In the elastic scenario, light nuclei, like Aluminum ( A = 27) and Silicon ( A = 29), aremore sensitive to light WIMPs scattering M DM ∼ multi-GeV, and provide the strongestupper bounds on the allowed parameter space favored by DAMA. In the inelastic scenario,which involves heavier candidates, experiments made of heavy nuclei, like Iodine ( A = 127),Xenon ( A = 129) and Tungsten ( A = 184) are the most constraining ones. Germanium( A = 73) made detectors fall in between. – 21 –n computing the rate of Eq. (3.2) for each experiment, we uniformly assume that thesmall number of events seen, if any, are signals from dark matter and we use the Poissonstatistics to find the parameter space excluded at a given confidence level (CL). Integrat-ing over the exposure time, the upper bounds are obtained by requiring that the totalnumber of events N tot is compatible with the number of observed events at 99 . 9% of CL,while Ref. [113] uses 99% CL. This means that we take less severe constraints, leadingto slightly smaller excluded regions. We have checked that the exclusion curves that weobtain reproduce fairly well the published 90% CL upper bounds for each experiment. Inthe following we describe the main features of the experiments considered, their particularcharacteristics, and the χ function adopted.CDMS: The Cryogenic CDMS experiment at Soudan Underground Laboratory op-erates Ge and Si made solid-state detectors. For a heavy WIMP, the Ge data are moreconstraining: we have considered the ensemble of the three released runs, with respec-tively an exposure of 19.4 kg-day [127], 34 kg-day [128] after cuts and a total exposureof 397.8 kg-days before cuts for the third run [37]. The sensitivity to nuclear recoils isin the energy range between 10-100 keV while the efficiencies and the energy resolutionof the detectors are given in [106]. For the third run, the efficiency is parameterizedas (cid:15) ( E R ) = 0 . 25 + 0 . E R − / < E R < 15 keV and (cid:15) ( E R ) = 0 . E R ≥ 15 keV. The searches on Si are more sensitive to light DM candidates. Forthe run released in Ref. [129], two 100 g Si ZIP detectors were used, and the data setcorresponds to 65.8 live days. The energy-dependent efficiency has been modeled inRef. [122] as (cid:15) ( E R ) = 0 . × . . 10 + 0 . E R − / 15) for 5 keV < E R < 20 keVand (cid:15) ( E R ) = 0 . × . . 40 + 0 . E R − / 80) for 20 keV < E R < 100 keV. Two can-didate events with energies E R (cid:39) 55 keV and E R (cid:39) 95 keV were observed, consistent withzero event once expected neutron background is taken into account. For the run releasedin Ref. [128], the total exposure is 12 kg-day after cuts, and the energy range has a higherminimum threshold 7 keV < E R < 100 keV.For zero observed event and Poisson distributed data, the χ function reduces to χ = 2 N tot (3.23)where N tot is the total predicted number of events in the entire range of energy, for givenvalues ( p , p , . . . ) of parameters ( P , P , . . . ) of the tested model. Therefore, for a 99 . N tot is lessthan 6 . χ for Poisson distributed data χ = 2 (cid:88) i (cid:40) N predi + B i − D i + D i log (cid:32) D i N predi + B i (cid:33)(cid:41) , (3.24)– 22 –ith N predi is the predicted number of events in the i th bin, B i = (0 . , . , . , . , . , . , . D i = (1 , , , , , , 4) is the number of detected events, andthe log term is zero if D i = 0. For a 99 . 9% CL limit and two parameters, the value ofthe χ function must be less than 13.8. As observed in Ref. [43], using a constant nuclearrecoil scintillation efficiency L eff = 0 . 19 gives a slightly underestimated energy thresholdof 4.5 keV. Correcting for this bias leads to a slightly less severe exclusion limit. In thispaper however, we will keep the energy bins as given by the XENON collaboration. Finally,we can notice that this experiment is also sensitive to inelastic dark matter in a similarway as the DAMA experiment, due to the close proximity of the target nucleus masses.CRESST: The CRESST-I and CRESST-II experiments at LNGS are mostly sensitiveto light and heavy WIMPs respectively. In the first data release (CRESST-I) [130], theprototype detector module uses Al O crystals as target, with a total exposure of 1.51kg-days covering the energy range between 0.6 - 20 keV. The energy resolution is given by σ ( E ) = (cid:112) (0 . keV ) + (0 . E ) and the collaboration reports eleven observed eventsafter cuts. Therefore, for Poisson distributed data, the number of predicted events cannotexceed 23 at 99.9% CL. For the second commissioned run of 2007 (CRESST-II) [40, 131],the target is changed to CaWO . The presence of heavy nuclei of Tungsten enhances thesensitivity to spin-independent inelastic scatterings. The energy range for the CRESST-IIZora and Verena detectors is 12-100 keV, with a total exposure of 47.9 kg-days before cutsand an acceptance of 0.9 on Tungsten recoil. The collaboration measured seven events,leading to a limit N tot < 16 at 99.9% CL.Other experiments are potentially relevant, like ZEPLIN, CoGeNT and TEXONO.Also, the total rate observed by DAMA provides a constraint that excludes part of theparameter space. However, to avoid the cluttering of the figures, we do not include theselimits in this paper, as they are weaker than the constraints presented. In the case of spin-independent elastic scatterings, the effective coherent couplings f p and f n of the WIMP to the proton and the neutron have similar values. Therefore, the onlyrelevant parameters are the WIMP mass M DM and the spin independent cross section onnucleon, here simply denoted as σ (see Eq. (3.4)). The main result for the elastic scenariois shown on the four panels of Fig. 12, where the regions in the plane M DM - σ favored byDAMA are plotted against the exclusion limits set by CDMS and XENON10 (the excludedregions are on the right and above the curves). From top left, the predictions for differenthalos: a) standard Maxwellian, b) this simulation, c) generalized Maxwellian with α = 1 . q = 0 . (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) M DM (cid:72) GeV (cid:76) Σ (cid:72) c m (cid:76) XENON10CDMS (cid:45) SiCDMS (cid:45) Ge b (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) M DM (cid:72) GeV (cid:76) Σ (cid:72) c m (cid:76) XENON10CDMS (cid:45) SiCDMS (cid:45) Ge c (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) M DM (cid:72) GeV (cid:76) Σ (cid:72) c m (cid:76) XENON10CDMS (cid:45) SiCDMS (cid:45) Ge d (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) M DM (cid:72) GeV (cid:76) Σ (cid:72) c m (cid:76) XENON10CDMS (cid:45) SiCDMS (cid:45) Ge Figure 12: Elastic Scenario, allowed regions in the plane M DM − σ (DM mass vs. spin independentcross-section on nucleon at zero momentum transfer) a) Standard Maxwellian halo ( ρ DM = 0 . / cm , v c = 220 km/s, v = 220 km/s, v esc =600 km/s)b) N-body simulation with baryons ( ρ DM = 0 . 37 GeV / cm , v c = 220 km/s)c) Generalized Maxwellian distribution ( ρ DM = 0 . 37 GeV / cm , v c = 190 km/s, v = 301 km/s, α = 1 . )d) Tsallis distribution ( ρ DM = 0 . 37 GeV / cm , v c = 190 km/s, v = 267 . km/s, q = 0 . )For c) and d), the circular velocity v c has been reduced compared to the standard value in orderto take into account a halo rotation. DAMA contours correspond to , and . % CL. Starsindicate local best-fit points. All other exclusion curves are at the . % CL. limit, so that the elastic solution becomes very marginal. Here we have chosen to fit theDAMA data on 12 bins, which yields a smaller allowed region than with 36 bins. For astandard Maxwellian halo, with v = 220 km/s, it appears that the 99.9% CL region ofDAMA is totally excluded by XENON10 at 99.9% CL.– 24 – DM (cid:61) Σ(cid:61) (cid:45) cm M DM (cid:61) Σ(cid:61) (cid:45) cm (cid:45) Eee (cid:72) keV (cid:76) S m , k (cid:72) c pd (cid:144) k g (cid:144) ke V (cid:76) M DM (cid:61) Σ(cid:61) (cid:45) cm M DM (cid:61) Σ(cid:61) (cid:45) cm (cid:45) (cid:45) t (cid:72) days (cid:76) R t o t (cid:72) c pd (cid:144) k g (cid:76) Figure 13: Elastic Scenario Modulation amplitude as a function of recoil energy (left) or time(right) for the best-fit points of Fig. 12. In the case of the time variation, events with a recoilenergy ≤ E R ≤ keV (in electron-equivalent units) were considered. For the simulation halo, the direct detection prediction is made by analyzing the phase-space of the central and best resolved Milky-Way sized galactic DM halo. The local struc-ture is obtained by selecting the particles ( N ring = 2 , ≤ R ≤ | z | ≤ R = 8 kpc in the galactic plane. Despite the small number offlows in this selection, which causes some numerical noise and artifacts in the χ plots,a few observations can be made. The DAMA regions are slightly smaller than in theMaxwellian case, and the best-fit points move to the left. Also, the region excluded byCDMS-Si slightly enlarges towards smaller cross sections. Moreover, the channeling regionextends outside the XENON exclusion limit, so that the compatibility between DAMA andthe other experiments is improved with the realistic halo from the simulation. While thebest-fit point of the channeling region is still excluded, the 99.9% and 99% CL regions arenot totally excluded anymore.Several factors contribute to this improvement. The departure from a strict Maxwelliandistribution has some impact, although moderate, as shown in panels c and d of Fig. 12.We considered a generalized Maxwellian halo with the two values α = 1 . α = 1 . q = 0 . ρ DM = 0 . v = 239 km / s,and smaller circular velocity v c = 190 km / s to reproduce the global average halo rotationvelocity), the 99.9% CL DAMA region is still totally excluded by XENON. Other possiblefactors include the velocity dispersion anisotropy, shown on Fig. 5, and a general anisotropyin the velocity field, induced by the co-rotating DM component, which has a small lagvelocity v lag (cid:39) 75 km/s compared to the galactic disk. A detailed modeling and discussionof the impact of these features will be presented in a subsequent paper [132].– 25 –he panels of Fig. 13 show the differential rate modulation for 2 ≤ Eee ≤ t (cid:39) 150 as expected. The number of flows is however too small to reliably indicateany deviation from t = 152 . The dominance of inelastic interactions over elastic ones is achieved in a natural way whenthe scattering is mediated by a massive gauge boson [113]. In the simplest scenario, themediator is the weak SU (2) L boson Z , so that no new gauge boson needs to be introducedbeyond the Standard Model. The free parameters are therefore the mass M DM of the DMcandidate, and its mass splitting δ with the next-to-lightest DM state. In this paper, wewill only consider this situation. More involved models, where the cross section differs fromthe weak cross section σ Z can be found in Ref. [113].Results for the inelastic scenario are summarized on Fig. 14. The DAMA favoredregions are shown together with the exclusion limits from CDMS and CRESST-II (theexcluded region extends on the left of each curve). As in the elastic scenario, the fourpanels are the predictions for different halos: a) standard Maxwellian, b) this simulation,c) generalized Maxwellian with α = 1 . 95, d) Tsallis with q = 0 . ρ DM = 0 . / cm , v =220 km / s, v esc = 600 km / s), there are three 99% CL regions compatible with DAMA,which are all mainly due to quenched scatterings on Iodine. The region with δ (cid:39) 30 keVis excluded by the constraints from both CDMS and CRESST. It is actually excludedby DAMA itself, as it leads to a year averaged rate higher than observed. The regionwith δ > 50 keV and M DM > Z , the standard out-of-equilibrium freeze-outpicture leads to a DM relic abundance higher than the value measured by WMAP [133],because the (co)annihilation rates are suppressed for M DM (cid:29) M Z . For scalar candidates,the adjustable mass of the charged SU (2) partner of the DM candidate enables to obtainthe correct relic density for a whole range of mass at the TeV scale [134]. Finally, theDAMA region with M DM < 100 GeV is not excluded by other experiments. For thesecandidates, due to the proximity of the Z pole, the relic abundance is lower than neededfor WMAP unless some asymmetry in the dark sector prevents DM density from collapsingduring the freeze-out [45].The impact of the numerical noise from the finite number of flows in the simulationbecomes more severe in the case of the inelastic scenario. Indeed, for inelastic interactions,only particles with a velocity high enough can produce a recoil. From Eq. (3.7), it followsthat the minimum velocity v ∗ corresponding to the energy E ∗ = µδ/M N is given by v ∗ =– 26 – (cid:76) ∆ (cid:72) keV (cid:76) M D M (cid:72) G e V (cid:76) CRESST IICDMS (cid:45) Ge b (cid:76) ∆ (cid:72) keV (cid:76) M D M (cid:72) G e V (cid:76) CRESST IICDMS (cid:45) Ge c (cid:76) ∆ (cid:72) keV (cid:76) M D M (cid:72) G e V (cid:76) CRESST IICDMS (cid:45) Ge d (cid:76) ∆ (cid:72) keV (cid:76) M D M (cid:72) G e V (cid:76) CRESST IICDMS (cid:45) Ge Figure 14: Inelastic Scenario, allowed regions in the plane δ − M DM (DM mass splitting vs. DMmass) a) Standard Maxwellian halo ( ρ DM = 0 . / cm , v c = 220 km/s, v = 220 km/s, v esc =600 km/s)b) N-body simulation with baryons ( ρ DM = 0 . 37 GeV / cm , v c = 220 km/s)c) Generalized Maxwellian distribution ( ρ DM = 0 . 37 GeV / cm , v c = 190 km/s, v = 332 km/s, α = 1 . )d) Tsallis distribution ( ρ DM = 0 . 37 GeV / cm , v c = 190 km/s, v = 267 . km/s, q = 0 . )For c) and d), the circular velocity v c has been reduced compared to the standard value in orderto take into account a halo rotation. DAMA contours correspond to , and . % CL. Starsindicate local best-fit points. All other exclusion curves are at the . % CL. (cid:112) δ/µ . As a consequence, in the inelastic case, the signal is due to the high velocity tailof the distribution. As panel b of Fig. 14 shows, the DAMA region with δ (cid:39) 30 keV isstrongly affected by numerical noise. In this region, only a small number of flows with avelocity around v ∗ (cid:39) 235 km / s contribute to the time dependent modulation. The region– 27 – DM (cid:61) ∆ (cid:61) 131 keVM DM (cid:61) ∆ (cid:61) DM (cid:61) ∆ (cid:61) (cid:45) Eee (cid:72) keV (cid:76) S m , k (cid:72) c pd (cid:144) k g (cid:144) ke V (cid:76) M DM (cid:61) ∆ (cid:61) 131 keVM DM (cid:61) ∆ (cid:61) DM (cid:61) ∆ (cid:61) (cid:45) t (cid:72) days (cid:76) R t o t (cid:72) c pd (cid:144) k g (cid:76) Figure 15: Inelastic Scenario Modulation amplitude as a function of recoil energy (left) or time(right) for the best-fit points of Fig. 14. In the case of the time variation, events with a recoil energy ≤ E R ≤ keV (in electron-equivalent units) were considered. with M DM < 100 GeV is also reduced, as it requires flows with a velocity higher than v ∗ (cid:39) 700 km / s. The third region is probably the most reliable. Its deformed shape, comparedto a Maxwellian halo, as well as its relative position to the exclusion curves are featuresthat show that the velocity distribution of the simulation halo is not Maxwellian and notisotropic. In particular, the compatibility between DAMA and the other experiments isstrongly improved. DAMA solutions at 90% CL are found, while only solutions at 99%CL were available in the case of the standard Maxwellian halo. This improvement of thecompatibility can also be attested by the different position of the best-fit point for eachcase.Results for a generalized Maxwellian or a Tsallis halo are shown on panels c and dof Fig. 14. The deviations from Maxwellianity cause the DAMA regions to shrink. Mostimportantly, in the region δ (cid:39) 120 keV and M DM (cid:39) δ , and can be seen to be in better agreement with the simulation.As in the elastic case, the differential rate modulation for 2 ≤ Eee ≤ . Summary & Perspectives In this paper, we have analyzed the phase-space structure of a galactic halo extracted froman advanced cosmological N-body simulation which contains stars, gas, and DM, and usedthis information to make direct detection predictions.Usually, such predictions are done with simplified astrophysical assumptions. The localphase-space structure in the solar neighborhood is taken as a smooth, isotropic Maxwelliandistribution. However, it is known for quite some time that deviations should be expectedin any realistic DM halo. The action of long-range gravitational forces, as well as thepresence of a non-virialized halo component in the form of cold streams from the contin-uous secondary galactic infall [115, 116] already contribute to distort the simple pictureof an isotropic Maxwellian halo. Cosmological numerical simulations done in the ΛCDMparadigm have demonstrated the hierarchical character of structure formation, leadingto galactic structures with numerous sub-halos [135, 136]. Recent N-body simulationswith very refined resolution, such as Aquarius [68] or Via Lactea [67] have even shownthat non-Gaussianity and anisotropy (the velocities are rather in the radial direction) arepresent [43, 54, 62, 137]. However, these simulations contain only DM particles, and cannottherefore provide a realistic description of a galaxy in a satisfactory way, despite the largenumber of particles. Dynamical interactions between the baryonic and the dark compo-nents of a galaxy are indeed expected to cause non negligible distortions on the phase-spacestructure of both components [86]. Recently, simulations with baryons have shown that thepresence of a galactic disk drags merging satellites in a preferential direction [71], whichleads to the formation of a thick stellar disc and a so-called dark disk.In this simulation, it appears that the DM halo contains a dark disk component thatis co-rotating with the galactic disk. The local DM density in the solar neighborhood has avalue ρ DM = 0 . 37 to 0 . 39 GeV / cm , consistent with the recent determination of Ref. [49],but higher than the usually quoted value of 0 . / cm . Despite large fluctuations, theaverage circular velocity differs significantly from zero throughout a thick region of thehalo. At the Sun’s location, a double Gaussian provides a reasonable fit of the distributionof the tangential velocity v φ . We infer that the rotating DM component has a mean lagvelocity v lag (cid:39) 75 km/s compared to the stellar disc, and contributes to around 25% of thetotal local density. This rather small value, compared to results of other ΛCDM simula-tions [71], rather points towards a quiescent merger history for the galactic halo consideredin our simulation. For the other velocity components and for the velocity module, we alsoobserve strong deviations from Gaussian and Maxwellian distributions, which we parame-terized by introducing different generalizations of Gaussian and Maxwellian distributions,see Eqs. (2.1–2.3). For the velocity module of DM particles around R = 8 kpc, we findthat a Tsallis distribution with a parameter q = 0 . 773 gives an excellent fit.The direct detection predictions are given as fits for the DAMA modulation signal andexclusion limits for null experiments. Both the elastic and the inelastic scenarios have beenconsidered with spin-independent cross-sections. As opposed to previous works [43,62], herethe regions in the parameter space are obtained directly from the velocity distributionsextracted from the simulation data, without modeling. They are then compared to the– 29 –redictions of various analytical models.As the co-rotating dark disk contribution to the local density is rather small, the ab-solute detection rates are not significantly enhanced compared to the standard Maxwelliancase. Nevertheless, the compatibility of DAMA vs. null experiments shows some improve-ment. For the elastic scenario, the rotating DM component causes the channeling regionto enlarge outside of the XENON constraint. However, the overall goodness of fit of aglobal solution is still poor as most of the channeling region remains excluded. For theinelastic scenario, the situation is much brighter, as two solutions with δ (cid:39) 130 keV arepermitted. The solution with a heavy candidate ( M DM ∼ TeV) becomes much moreacceptable with the realistic halo provided by our simulation. Moreover, such candidatecan also achieve naturally the relic abundance needed for WMAP [134], and is very littleconstrained by accelerator data [138]. If the DAMA modulation signal is taken seriously,all these hints concur to some interesting new physics beyond the Standard Model at theTeV scale. The solution with a mass around M DM (cid:39) 80 GeV is also acceptable from boththe direct detection data and the particle physics point of view. A correct relic abundancehowever requires some asymmetry in the dark sector [45] or some non thermal productionmechanism. The accelerator constraints could be more stringent, but have to be analyzedwithin a particular model [139]. These questions are beyond the scope of the present paper. Acknowledgments We are grateful towards Jean-Charles Lambert for his help in handling GLNemo to producethe pictures of the galaxy for this paper. We would also like to thank Chiara Arina andArman Khalatyan for helpful comments and discussions. RT would like to thank St´ephanieCourty for her help in selecting the simulated halo and generating the initial conditions.FSL work is supported by a belgian FNRS grant and the IAP. This work was grantedaccess to the HPC resources of CINES under the allocation 2009-SAP2191 made by GENCI(Grand Equipement National de Calcul Intensif). This work was partly supported by grantANR-06-BLAN-0172. References [1] Yoshiaki Sofue and Vera Rubin. Rotation Curves of Spiral Galaxies. Ann. Rev. Astron.Astrophys. , 39:137–174, 2001.[2] Albert Bosma. Dark matter in galaxies: Observational overview. Dark Matter in Galaxies,IAU symposium series, , 220, 2004.[3] F. Zwicky. Spectral displacement of extra galactic nebulae. Helv. Phys. Acta , 6:110–127,1933.[4] Aaron D. Lewis, David A. Buote, and John T. Stocke. Chandra Observations of Abell 2029:The Dark Matter Profile at ¡0.01 Rvir in an Unusually Relaxed Cluster. Astrophys. J. ,586:135–142, 2003.[5] Alexandre Refregier. Weak Gravitational Lensing by Large-Scale Structure. Ann. Rev.Astron. Astrophys. , 41:645–668, 2003. – 30 – 6] Douglas Clowe, Anthony Gonzalez, and Maxim Markevitch. Weak lensing massreconstruction of the interacting cluster 1e0657-558: Direct evidence for the existence ofdark matter. Astrophys. J. , 604:596–603, 2004.[7] Raphael Gavazzi et al. The Sloan Lens ACS Survey. VI: Discovery and analysis of a doubleEinstein ring. 2008.[8] Shaun Cole et al. The 2dF Galaxy Redshift Survey: Power-spectrum analysis of the finaldataset and cosmological implications. Mon. Not. Roy. Astron. Soc. , 362:505–534, 2005.[9] Uros Seljak et al. Cosmological parameter analysis including SDSS Ly-alpha forest andgalaxy bias: Constraints on the primordial spectrum of fluctuations, neutrino mass, anddark energy. Phys. Rev. , D71:103515, 2005.[10] D. N. Spergel et al. Wilkinson microwave anisotropy probe (wmap) three year results:Implications for cosmology. 0300.[11] C. Amsler et al. Review of particle physics. Phys. Lett. , B667:1, 2008.[12] Jean Iliopoulos. Physics Beyond the Standard Model. 2008.[13] Howard Baer. Physics Beyond the Standard Model. 2009.[14] R. D. Peccei and Helen R. Quinn. Constraints Imposed by CP Conservation in the Presenceof Instantons. Phys. Rev. , D16:1791–1797, 1977.[15] John R. Ellis, J. S. Hagelin, Dimitri V. Nanopoulos, Keith A. Olive, and M. Srednicki.Supersymmetric relics from the big bang. Nucl. Phys. , B238:453–476, 1984.[16] Hsin-Chia Cheng, Jonathan L. Feng, and Konstantin T. Matchev. Kaluza-Klein darkmatter. Phys. Rev. Lett. , 89:211301, 2002.[17] Marco Cirelli, Nicolao Fornengo, and Alessandro Strumia. Minimal dark matter. Nucl.Phys. , B753:178–194, 2006.[18] J. Knodlseder et al. The all-sky distribution of 511-keV electron positron annihilationemission. Astron. Astrophys. , 441:513–532, 2005.[19] S. D. Hunter et al. EGRET observations of the diffuse gamma-ray emission from thegalactic plane. Astrophys. J. , 481:205–240, 1997.[20] Oscar Adriani et al. An anomalous positron abundance in cosmic rays with energies 1.5.100GeV. Nature , 458:607–609, 2009.[21] J. Chang et al. An excess of cosmic ray electrons at energies of 300.800 GeV. Nature ,456:362–365, 2008.[22] Celine Boehm, Dan Hooper, Joseph Silk, Michel Casse, and Jacques Paul. MeV darkmatter: Has it been detected? Phys. Rev. Lett. , 92:101301, 2004.[23] W. de Boer et al. Excess of EGRET galactic gamma ray data interpreted as dark matterannihilation. 2004.[24] J. M. Frere et al. MeV right-handed neutrinos and dark matter. Phys. Rev. , D75:085017,2007.[25] Patrick J. Fox and Erich Poppitz. Leptophilic Dark Matter. Phys. Rev. , D79:083528, 2009. – 31 – 26] Marco Cirelli, Mario Kadastik, Martti Raidal, and Alessandro Strumia. Model-independentimplications of the e+, e-, anti-proton cosmic ray spectra on properties of Dark Matter. Nucl. Phys. , B813:1–21, 2009.[27] Alejandro Ibarra, David Tran, and Christoph Weniger. Decaying Dark Matter in Light ofthe PAMELA and Fermi LAT Data. 2009.[28] Georg Weidenspointner et al. An asymmetric distribution of positrons in the Galactic diskrevealed by γ -rays. Nature , 451:159–162, 2008.[29] F. W. Stecker, S. D. Hunter, and D. A. Kniffen. The Likely Cause of the EGRET GeVAnomaly and its Implications. Astropart. Phys. , 29:25–29, 2008.[30] T. A. Porter and for the Fermi LAT Collaboration. Fermi LAT Measurements of the DiffuseGamma-Ray Emission at Intermediate Galactic Latitudes. 2009.[31] H. E. S. S. Collaboration: F. Aharonian. Probing the ATIC peak in the cosmic-ray electronspectrum with H.E.S.S. 2009.[32] Aous A. Abdo et al. Measurement of the Cosmic Ray e+ plus e- spectrum from 20 GeV to1 TeV with the Fermi Large Area Telescope. Phys. Rev. Lett. , 102:181101, 2009.[33] Dan Hooper, Pasquale Blasi, and Pasquale Dario Serpico. Pulsars as the Sources of HighEnergy Cosmic Ray Positrons. JCAP , 0901:025, 2009.[34] Dmitry Malyshev, Ilias Cholis, and Joseph Gelfand. Pulsars versus Dark MatterInterpretation of ATIC/PAMELA. 2009.[35] Tsvi Piran, Nir J. Shaviv, and Ehud Nakar. Inhomogeneity in the Supernova Remnants as aNatural Explanation of the PAMELA/ATIC Observations. 2009.[36] 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:23–31, 2000.[37] Z. Ahmed et al. Search for Weakly Interacting Massive Particles with the First Five-TowerData from the Cryogenic Dark Matter Search at the Soudan Underground Laboratory. Phys. Rev. Lett. , 102:011301, 2009.[38] J. Angle et al. First Results from the XENON10 Dark Matter Experiment at the GranSasso National Laboratory. Phys. Rev. Lett. , 100:021303, 2008.[39] J. Angle et al. Limits on spin-dependent WIMP-nucleon cross-sections from the XENON10experiment. Phys. Rev. Lett. , 101:091301, 2008.[40] G. Angloher et al. Commissioning Run of the CRESST-II Dark Matter Search. 2008.[41] R. Bernabei et al. First results from DAMA/LIBRA and the combined results withDAMA/NaI. Eur. Phys. J. , C56:333–355, 2008.[42] Frank Petriello and Kathryn M. Zurek. DAMA and WIMP dark matter. JHEP , 09:047,2008.[43] Malcolm Fairbairn and Thomas Schwetz. Spin-independent elastic WIMP scattering andthe DAMA annual modulation signal. JCAP , 0901:037, 2009.[44] Christopher Savage, Katherine Freese, Paolo Gondolo, and Douglas Spolyar. Compatibilityof DAMA/LIBRA dark matter detection with other searches in light of new Galacticrotation velocity measurements. 2009. – 32 – 45] Chiara Arina, Fu-Sin Ling, and Michel H. G. Tytgat. The Inert Doublet Model andInelastic Dark Matter. 2009.[46] Ricardo A. Flores. DYNAMICAL ESTIMATES AND BOUNDS FOR THE LOCALDENSITY OF DARK MATTER. Phys. Lett. , B215:73, 1988.[47] Y. Sofue, M. Honma, and T. Omodaka. Unified Rotation Curve of the Galaxy –Decomposition into de Vaucouleurs Bulge, Disk, Dark Halo, and the 9-kpc Rotation Dip –.2008.[48] Alfred Bing-Chih Chen, Phillip K. Lu, Rene A. Mendez, and William F. van Altena. Darkmatter: Local volume density versus total surface density. Astron. J. , 126:762, 2003.[49] Riccardo Catena and Piero Ullio. A novel determination of the local dark matter density.2009.[50] J. Lavalle, Q. Yuan, D. Maurin, and X. J. Bi. Full Calculation of Clumpiness Boost factorsfor Antimatter Cosmic Rays in the light of LambdaCDM N-body simulation results. 2007.[51] Mauro Sereno and Ph. Jetzer. Dark matter vs. modifications of the gravitationalinverse-square law. Results from planetary motion in the solar system. Mon. Not. Roy.Astron. Soc. , 371:626–632, 2006.[52] J. M. Frere, Fu-Sin Ling, and G. Vertongen. Bound on the Dark Matter Density in theSolar System from Planetary Motions. Phys. Rev. , D77:083005, 2008.[53] Steen H. Hansen and Ben Moore. A universal density slope - velocity anisotropy relation forrelaxed structures. New Astron. , 11:333, 2006.[54] Steen H. Hansen, Ben Moore, Marcel Zemp, and Joachim Stadel. A universal velocitydistribution of relaxed collisionless structures. JCAP , 0601:014, 2006.[55] Anne M. Green. The WIMP annual modulation signal and non-standard halo models. Phys.Rev. , D63:043005, 2001.[56] Anne M. Green. Effect of halo modelling on WIMP exclusion limits. Phys. Rev. ,D66:083003, 2002.[57] Anne M Green. Effect of realistic astrophysical inputs on the phase and shape of the WIMPannual modulation signal. Phys. Rev. , D68:023004, 2003.[58] Christopher Savage, Katherine Freese, and Paolo Gondolo. Annual Modulation of DarkMatter in the Presence of Streams. Phys. Rev. , D74:043531, 2006.[59] J. D. Vergados, S. H. Hansen, and O. Host. The impact of going beyond the Maxwelldistribution in direct dark matter detection rates. Phys. Rev. , D77:023509, 2008.[60] J. D. Vergados. The NFW Dark Matter Density Profile Leads to Axially SymmetricVelocity Distribution-Imlications on Direct Dark Matter Searche. 2008.[61] Marc Kamionkowski and Savvas M. Koushiappas. Galactic substructure and directdetection of dark matter. arXiv:0801.3269 [astro-ph].[62] John March-Russell, Christopher McCabe, and Matthew McCullough. Inelastic DarkMatter, Non-Standard Halos and the DAMA/LIBRA Results. JHEP , 05:071, 2009.[63] Simon D. M. White, Carlos S. Frenk, Marc Davis, and George Efstathiou. Clusters,filaments, and voids in a universe dominated by cold dark matter. Astrophys. J. ,313:505–516, 1987. – 33 – 64] Giuseppe Tormen, Francois R. Bouchet, and S. D. M. White. The Structure and DynamicalEvolution of Dark Matter Halos. Mon. Not. Roy. Astron. Soc. , 286:865–884, 1997.[65] Romain Teyssier. Cosmological Hydrodynamics with Adaptive Mesh Refinement: a newhigh resolution code called RAMSES. 2001.[66] Volker Springel. The cosmological simulation code GADGET-2. Mon. Not. Roy. Astron.Soc. , 364:1105–1134, 2005.[67] J. Diemand et al. Clumps and streams in the local dark matter distribution. 2008.[68] Volker Springel et al. The Aquarius Project: the subhalos of galactic halos. Mon. Not. Roy.Astron. Soc. , 391:1685–1711, 2008.[69] Mark Vogelsberger et al. Phase-space structure in the local dark matter distribution and itssignature in direct detection experiments. 2008.[70] Brad K. Gibson et al. Hydrodynamical Adaptive Mesh Refinement Simulations of DiskGalaxies. 2008.[71] J. I. Read, L. Mayer, A. M. Brooks, F. Governato, and G. Lake. A dark matter disc in threecosmological simulations of Milky Way mass galaxies. 2009.[72] J. I. Read, G. Lake, O. Agertz, and Victor P. Debattista. Thin, thick and dark discs inLCDM. 2008.[73] T. Bruch, J. Read, L. Baudis, and G. Lake. Detecting the Milky Way’s Dark Disk. Astrophys. J. , 696:920–923, 2009.[74] Lucio Mayer, Fabio Governato, and Tobias Kaufmann. The formation of disk galaxies incomputer simulations. 2008.[75] Fabio Governato et al. Forming Disk Galaxies in Lambda CDM Simulations. Mon. Not.Roy. Astron. Soc. , 374:1479–1494, 2007.[76] Fabio Governato et al. Forming a Large Disk Galaxy from a z¡1 Major Merger. 2008.[77] R. Teyssier. Cosmological hydrodynamics with adaptive mesh refinement. a new highresolution code called ramses. A&A , 385:337–364, 2002.[78] Neal Katz, Lars Hernquist, and David H. Weinberg. Galaxies and gas in a cold dark matteruniverse. Astrophys. J. , 399:L109, 1992.[79] Ralph S. Sutherland and Michael A. Dopita. Cooling functions for low - densityastrophysical plasmas. Astrophys. J. Suppl. , 88:253, 1993.[80] Yann Rasera and Romain Teyssier. The History of the Baryon Budget: Cosmic Logistics ina Hierarchical Universe. 2005.[81] Y. Dubois and R. Teyssier. Cosmological MHD simulation of a cooling flow cluster. 2008.[82] Wayne Hu and Daniel J. Eisenstein. The Structure of structure formation theories. Phys.Rev. , D59:083509, 1999.[83] George R. Blumenthal, S. M. Faber, Ricardo Flores, and Joel R. Primack. Contraction ofDark Matter Galactic Halos Due to Baryonic Infall. Astrophys. J. , 301:27, 1986.[84] Oleg Y. Gnedin, Andrey V. Kravtsov, Anatoly A. Klypin, and Daisuke Nagai. Response ofdark matter halos to condensation of baryons: cosmological simulations and improvedadiabatic contraction model. Astrophys. J. , 616:16–26, 2004. – 34 – 85] Mario G. Abadi, Julio F. Navarro, Mark Fardal, Arif Babul, and Matthias Steinmetz.Galaxy-Induced Transformation of Dark Matter Halos. 2009.[86] E. Athanassoula. On the nature of bulges in general and of box/peanut bulges in particular.Input from N -body simulations. Mon. Not. Roy. Astron. Soc. , 358:1477–1488, 2005.[87] Jurg Diemand, Michael Kuhlen, and Piero Madau. Dark matter substructure andgamma-ray annihilation in the Milky Way halo. Astrophys. J. , 657:262, 2007.[88] C. Tsallis. Possible generalization of boltzmanngibbs entropy. J. Stat. Phys. , 52:479, 1988.[89] J. A. S. Lima, R. Silva, and A. R. Plastino. Nonextensive thermostatistics and theh-theorem. Physical Review Letters , 86:2938, 2001.[90] Steen H. Hansen, Daniel Egli, Lukas Hollenstein, and Christoph Salzmann. Dark matterdistribution function from non-extensive statistical mechanics. New Astron. , 10:379, 2005.[91] Ian R. Walker, J. Christopher Mihos, and Lars Hernquist. Quantifying the Fragility ofGalactic Disks in Minor Mergers. Astrophys. J. , 460:121, 1996.[92] Jorge Penarrubia, Alan McConnachie, and Arif Babul. On the formation of extendedgalactic disks by tidally disrupted dwarf galaxies. Astrophys. J. , 650:L33–L36, 2006.[93] Stelios Kazantzidis, James S. Bullock, Andrew R. Zentner, Andrey V. Kravtsov, andLeonidas A. Moustakas. Cold Dark Matter Substructure and Galactic Disks I:Morphological Signatures of Hierarchical Satellite Accretion. 2007.[94] G. Gilmore and N. Reid. New light on faint stars. III - Galactic structure towards the SouthPole and the Galactic thick disc. Mon. Not. Roy. Astron. Soc. , 202:1025–1047, 1983.[95] Chris W. Purcell, James S. Bullock, and Manoj Kaplinghat. The Dark Disk of the MilkyWay. 2009.[96] Mario Juric et al. The Milky Way Tomography with SDSS. 1. Stellar Number DensityDistribution. Astrophys. J. , 673:864–914, 2008.[97] C. Soubiran, O. Bienayme, and A. Siebert. Vertical distribution of Galactic disk stars I -Kinematics and metallicity. Astron. Astrophys. , 398:141–152, 2003.[98] B. Nordstrom et al. The Geneva-Copenhagen survey of the Solar neighbourhood: Ages,metallicities, and kinematic properties of 14,000 F and G dwarfs. Astron. Astrophys. ,418:989, 2004.[99] Graciela B. Gelmini. Search for Dark Matter. Int. J. Mod. Phys. , A23:4273–4288, 2008.[100] Dan Hooper. TASI 2008 Lectures on Dark Matter. 2009.[101] R. Bernabei et al. Possible implications of the channeling effect in NaI(Tl) crystals. Eur.Phys. J. , C53:205–213, 2008.[102] A. Bottino, F. Donato, N. Fornengo, and S. Scopel. Zooming in on light relic neutralinos bydirect detection and measurements of galactic antimatter. Phys. Rev. , D77:015002, 2008.[103] J. D. Lewin and P. F. Smith. Review of mathematics, numerical factors, and corrections fordark matter experiments based on elastic nuclear recoil. Astropart. Phys. , 6:87–112, 1996.[104] Sarah Andreas, Thomas Hambye, and Michel H. G. Tytgat. WIMP dark matter, Higgsexchange and DAMA. JCAP , 0810:034, 2008. – 35 – Phys. Rev. , D71:123520, 2005.[107] A. Bottino, N. Fornengo, and S. Scopel. Light relic neutralinos. Phys. Rev. , D67:063519,2003.[108] A. Bottino, F. Donato, N. Fornengo, and S. Scopel. Lower bound on the neutralino massfrom new data on CMB and implications for relic neutralinos. Phys. Rev. , D68:043506, 2003.[109] A. Bottino, F. Donato, N. Fornengo, and S. Scopel. Light neutralinos and WIMP directsearches. Phys. Rev. , D69:037302, 2004.[110] David Tucker-Smith and Neal Weiner. Inelastic dark matter. Phys. Rev. , D64:043502, 2001.[111] David Tucker-Smith and Neal Weiner. The status of inelastic dark matter. Phys. Rev. ,D72:063509, 2005.[112] Spencer Chang, Graham D. Kribs, David Tucker-Smith, and Neal Weiner. Inelastic DarkMatter in Light of DAMA/LIBRA. Phys. Rev. , D79:043513, 2009.[113] Yanou Cui, David E. Morrissey, David Poland, and Lisa Randall. Candidates for InelasticDark Matter. JHEP , 05:076, 2009.[114] Douglas P. Finkbeiner, Tongyan Lin, and Neal Weiner. Inelastic Dark Matter andDAMA/LIBRA: An Experimentum Crucis. 2009.[115] Graciela Gelmini and Paolo Gondolo. WIMP annual modulation with opposite phase inlate-infall halo models. Phys. Rev. , D64:023504, 2001.[116] Fu-Sin Ling, Pierre Sikivie, and Stuart Wick. Diurnal and Annual Modulation of Cold DarkMatter Signals. Phys. Rev. , D70:123503, 2004.[117] Martin C. Smith et al. The RAVE Survey: Constraining the Local Galactic Escape Speed. Mon. Not. Roy. Astron. Soc. , 379:755–772, 2007.[118] A. K. Drukier, K. Freese, and D. N. Spergel. Detecting Cold Dark Matter Candidates. Phys. Rev. , D33:3495–3508, 1986.[119] Katherine Freese, Joshua A. Frieman, and Andrew Gould. Signal Modulation in Cold DarkMatter Detection. Phys. Rev. , D37:3388, 1988.[120] Richard H. Helm. Inelastic and Elastic Scattering of 187-Mev Electrons from SelectedEven-Even Nuclei. Phys. Rev. , 104:1466–1475, 1956.[121] Gintaras Duda, Ann Kemper, and Paolo Gondolo. Model independent form factors for spinindependent neutralino nucleon scattering from elastic electron scattering data. JCAP ,0704:012, 2007.[122] Christopher Savage, Graciela Gelmini, Paolo Gondolo, and Katherine Freese. Compatibilityof DAMA/LIBRA dark matter detection with other searches. JCAP , 0904:010, 2009.[123] Walter Dehnen and James Binney. Local stellar kinematics from Hipparcos data. Mon. Not.Roy. Astron. Soc. , 298:387–394, 1998.[124] M. J. Reid et al. Trigonometric Parallaxes of Massive Star Forming Regions: VI. GalacticStructure, Fundamental Parameters and Non- Circular Motions. Astrophys. J. , 700:137–148,2009. – 36 – Phys. Rev. , D78:083520, 2008.[126] Spencer Chang, Aaron Pierce, and Neal Weiner. Using the Energy Spectrum atDAMA/LIBRA to Probe Light Dark Matter. 2008.[127] D. S. Akerib et al. First results from the cryogenic dark matter search in the SoudanUnderground Lab. Phys. Rev. Lett. , 93:211301, 2004.[128] D. S. Akerib et al. Limits on spin-independent WIMP nucleon interactions from thetwo-tower run of the Cryogenic Dark Matter Search. Phys. Rev. Lett. , 96:011302, 2006.[129] D. S. Akerib et al. New results from the Cryogenic Dark Matter Search experiment. Phys.Rev. , D68:082002, 2003.[130] G. Angloher et al. Limits on WIMP dark matter using sapphire cryogenic detectors. Astropart. Phys. , 18:43–55, 2002.[131] G. Angloher et al. Limits on WIMP dark matter using scintillating CaWO-4 cryogenicdetectors with active background suppression. Astropart. Phys. , 23:325–339, 2005.[132] Fu-Sin Ling. Is the Dark Disc contribution to Dark Matter Signals important. Submitted toPhys. Rev. D , 2009.[133] G. Hinshaw et al. Five-Year Wilkinson Microwave Anisotropy Probe (WMAP)Observations:Data Processing, Sky Maps, & Basic Results. 2008.[134] T. Hambye, F. S. Ling, L. Lopez Honorez, and J. Rocher. Scalar Multiplet Dark Matter.2009.[135] E. Athanassoula, F. S. Ling, E. Nezri, and R. Teyssier. Gamma ray fluxes from acosmological dark matter simulation. Astropart. Phys. , 31:37–45, 2009.[136] Jurg Diemand and Ben Moore. The structure and evolution of cold dark matter halos. 2009.[137] Ben Moore et al. Dark matter in Draco and the Local Group: Implications for directdetection experiments. Phys. Rev. , D64:063508, 2001.[138] Riccardo Barbieri, Lawrence J. Hall, and Vyacheslav S. Rychkov. Improved naturalnesswith a heavy Higgs: An alternative road to LHC physics. Phys. Rev. , D74:015007, 2006.[139] Erik Lundstrom, Michael Gustafsson, and Joakim Edsjo. The Inert Doublet Model andLEP II Limits. 2008., D74:015007, 2006.[139] Erik Lundstrom, Michael Gustafsson, and Joakim Edsjo. The Inert Doublet Model andLEP II Limits. 2008.