{'fraction_non_alphanumeric': 0.030387468464058442, 'fraction_numerical': 0.009827917309937116, 'mean_word_length': 4.6117274167987325, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The local dark matter phase-space density and impact on WIMP direct detection

We present a new determination of the local dark matter phase-space density. This result is obtained implementing, in the limit of isotropic velocity distribution and spherical symmetry, Eddington's inversion formula, which links univocally the dark matter distribution function to the density profile, and applying, within a Bayesian framework, a Markov Chain Monte Carlo algorithm to sample mass models for the Milky Way against a broad and variegated sample of dynamical constraints. We consider three possible choices for the dark matter density profile, namely the Einasto, NFW and Burkert profiles, finding that the velocity dispersion, which characterizes the width in the distribution, tends to be larger for the Burkert case, while the escape velocity depends very weakly on the profile, with the mean value we obtain being in very good agreement with estimates from stellar kinematics. The derived dark matter phase-space densities differ significantly--most dramatically in the high velocity tails--from the model usually taken as a reference in dark matter detection studies, a Maxwell-Boltzmann distribution with velocity dispersion fixed in terms of the local circular velocity and with a sharp truncation at a given value of the escape velocity. We discuss the impact of astrophysical uncertainties on dark matter scattering rates and direct detection exclusion limits, considering a few sample cases and showing that the most sensitive ones are those for light dark matter particles and for particles scattering inelastically. As a general trend, when adopting a self-consistent phase-space density, we find that rates are larger, and hence exclusion limits stronger, than with the standard Maxwell-Boltzmann approximation. Tools for applying our result on the local dark matter phase-space density to other dark matter candidates or experimental setups are provided.


Introduction
The direct detection technique has played a major role in the quest for the identification of the dark matter (DM) component of the Universe. The goal is measure the recoil energy due to the scatterings on a target material of the particles forming the Milky Way DM halo. The prime focus is on weakly interacting massive particles (WIMPs), non-relativistic thermal relics from the Universe with mass in the range, say, between few GeV and few TeV, which would scatter elastically on a nucleus, depositing energies of the order of few (tens of) keV (for a recent review on the DM problem and direct detection see, e.g., [1]). In recent years several experiments have reached the sensitivity to start probing the level of scattering cross sections expected for WIMPs. In particular, three collaborations have published results compatible with a positive signal: DAMA and DAMA/LIBRA [2] detected an annual modulation in the total event rate consistent with the effect expected from WIMPs scatterings because of the orbit of the Earth around the Sun; CoGeNT has recently confirmed [3] the detection of a low-energy exponential tail in their count rate consistent with the shape predicted for the signal from a WIMP slightly lighter than 10 GeV, as already found previously [4], showing in addition a 2.8 σ indication in favor of an annual modulation effect; finally, CRESST-II [5] has just reported the detection of a number events in the acceptance region of low energy nuclear recoils, rejecting the hypothesis that they are due to background only, and showing the indication for a component of mass either around 25 or 11 GeV. On the other hand, other experiments such as CDMS [6] and Xenon100 [7] have not found any evidence for a DM signal and produced exclusion plots in the plane WIMP mass versus scattering cross section which seem to disfavor the region of the parameter space preferred by DAMA, CoGeNT or CRESST-II (in this respect, see also the recent reanalysis of the data taken at an early stage by CDMS at a shallow site [8], and the result of Xenon10 [9]).
The interpretation of a count rate or an annual modulation effect in a given experiment in terms of properties of DM particles (and hence obviously the comparison of results obtained by different collaborations, possibly, with different techniques and/or different target materials) depend on a number of assumptions implemented in the analysis of the data. Besides issues regarding understanding the target material, such as nuclear form factors or whether channeling occurs, and the performance of the detector, such as determining the energy threshold, the quenching factors and background rejection/contamination, there are uncertainties related to the DM signal itself. First of all there is an uncertainty in the normalization of the incident DM particle flux, which scales with the local halo density, often quoted to be unknown within a factor of 2 or so. Large uncertainties also affect the energy spectrum of the DM particles in the detector frame, in turn connected to their velocity distribution in the Galactic frame and to the proper motion of the Sun/Earth system. Concerning these, the vast majority of the analyses adopt a standard paradigm in which the velocity distribution is assumed to be Maxwell-Boltzmann with velocity dispersion scaled up of a factor of 3/2 with respect to the Sun circular velocity (usually taken as the IAU standard value of 220 km s −1 ) and truncated to the value assumed for the escape velocity.
The Maxwellian distribution is the configuration maximizing the entropy for a selfgravitating collisionless system (complete wash out of the initial conditions after gravitational collapse) and is associated to the spherical isothermal sphere density profile, which declines as r −2 at large radii and hence supports a flat rotation curve. It is a well-motivated form but unlikely a fair description of the Milky Way DM halo. In particular cosmological N-body simulations find that DM halos have density profiles falling more rapidly at large radii, as r −3 , and velocity distributions showing significant departures from the Maxwell-Boltzmann shape, see, e.g., the results from the high-resolution simulations Via Lactea [10] and Aquarius [11]. A few analyses have discussed the impact on direct detection results when using DM phasespace distribution function as directly read out from the numerical simulations or from other specialized approaches devised to describe in detail the fine-grained structure of the DM velocity distribution, see, e.g., [12][13][14][15][16][17][18][19][20][21]. The main shortcoming of these approaches is that, since it is highly non-trivial to include baryons and baryonic feedback in the simulations, these analyses treat the Galaxy as if made of DM only (in some cases, normalizing the circular velocity obtained in the model for the DM particles to the locally measured circular velocity). In reality, the stellar and gas components dominate the potential well in the inner Galaxy and out to, at least, our Galactocentric distance, with, e.g., the DM contribution to local surface mass density standing essentially within the error bar from star population counts [22].
Taking an observationally driven point of view, in Ref. [23] (hereafter Paper I) we reconsidered the problem of constructing mass models for the Milky Way, and showed that, when considering an updated set of dynamical tracers and an efficient way of scanning a multidimensional parameter space, a high precision determination of the local DM halo density can be obtained, namely around 0.39 GeV cm −3 with a 1-σ error bar of about 7%. This result relies partially on the assumption of spherical symmetry we imposed on the DM component (and on the same prescription to compute the gravitational potential described in section 2 of the present paper), and would be revised if one argues, for instance, for some level of flattening of the DM halo, see, e.g., [24][25][26][27]. At the same time, it is a rather solid result since there is no tension between the simplified picture we proposed and any of the very var-iegated sets of complementary observables implemented in the analysis (this, together with the effectiveness of Bayesian inference applied to our mass model for the Galaxy, is reflected in the tiny error we found).
In case one assumes that the DM distribution is also isotropic, there is a one-to-one correspondence between the spherically symmetric density profile and the underlying distribution function, with the latter that can be computed from the former through Eddington's formula [28]. Although this requires some heavy numerical integrals, it is nowadays possible to perform this inversion on very large samples of trial cases. Relying on the same mass models introduced in Paper I and an analogous Bayesian approach with a Markov Chain Monte Carlo scan of the underlying parameter space, we study here the family of isotropic distribution functions one can associate to the spherical DM density profiles. Following this approach, we can address here for the first time the theoretical uncertainty on the direct detection signal within a framework in which the value of the local halo density, the shape of the distribution function and its truncation at the escape velocity, and the circular velocity of the Sun are taken self-consistently and in agreement with the available dynamical constraints. Such method is way more powerful than deriving an overall theoretical error assuming that the local halo density, the velocity dispersion, the escape velocity and the Sun circular velocity have given uncertainties to be propagated as uncorrelated errors, as actually done in most analyses in the literature.
Results presented in this paper are valid in the limit of spherical symmetry 1 and isotropy, indeed rather strong assumptions, and should be regarded as a first step towards a study allowing for, at least, axisymmetric configurations (the case with axisymmetric models can be in principle treated in a specular way, but it is computationally much more demanding; such case is subject of ongoing work). On the other hand, the inner regions of a galaxy, such as at own position within the Milky Way DM halo, are those for which the simulations find weaker evidence for departure from spherical symmetry, and, favored also by the presence of large amount of baryonic matter, those with the largest chance for gravitational relaxation of the collisionless DM system, and hence where the distribution function is expected to be close to isotropic.
The paper is organized as follows: In Section 2 we introduce the Galactic model and the procedure adopted to compute the local distribution function. In Section 3 we review the computation of direct detection rates, emphasizing the connection to the distribution function, and detail the method implemented to compute the exclusion limit in one sample case, the one from the recent data release from the Xenon Collaboration. Section 4 illustrates the statistical method implemented, while Section 5 contains our results. Section 6 concludes.

Galactic model and DM phase-space distribution function
Assuming that the distribution of DM particles in the Galaxy is spherically symmetric and isotropic, and in the limit of spherical symmetry for the underlying gravitational potential for the Galaxy Φ(r), Eddington's formula [28] gives an one-to-one correspondence between the DM halo density profile ρ h (r) and its phase-space distribution function F h : where we defined the relative potential Ψ(r) = −Φ(r) + Φ(r = R vir ), with the virial radius R vir being the radius at which the DM halo as a virialized object is truncated. F h depends on particle velocity and position in the Galaxy only through the value of the relative energy E = −E + Φ(r = R vir ) = −E kin + Ψ(r), with E and E kin , respectively, the total and kinetic energy. From the numerical point of view, it is actually simpler to implement Eq. (2.1) by changing the integration variable from Ψ to the radius of the spherical system r; to compute F h at the local Galactocentric distance, it is sufficient to specify the radial dependence of ρ h and Φ from the local position in the Galaxy out to R vir . We will then make an ansatz for a mass model for the Galaxy, including a parameterization for the contributions of the stellar bulge/bar and stellar disc, compare it against available dynamical constraints and implement a prescription to extrapolate Φ(r). The Galactic mass model we adopt is the same introduced in Paper I; we summarize it briefly here. The DM halo component takes the generic form: where f (x) is the function defining the shape of the DM density profile. We will consider three separate choices, namely: the profile originally proposed by Navarro, Frenk and White [29] as the universal profile describing DM halos in numerical N-body simulations, and extensively used in the literature: the Einasto profile [30,31], favored by the latest simulations: and, finally, the Burkert profile [32]: which has core, possibly reflecting the wash out of the central DM enhancement induced by the baryon infall. In Eq. (2.2), ρ ′ and a h fix, respectively, a reference normalization and length scale; most often these are given in terms of the virial mass M vir and the concentration parameter c vir , inverting the relations: where, in the first equation, which defines also the virial radius R vir , the virial overdensity ∆ vir is computed according to Ref. [33], whileρ 0 is the mean background density today. Ω DM and Ω b , the dark matter and baryon energy densities in units of the critical density, enter in this expression since we are assuming that only a fraction equal to Ω DM /(Ω DM + Ω b ) of the total virial mass consists of dark matter; their values are fixed according to the mean values from the fit of the 7-year WMAP data [34]. Finally, in the equation for c vir , r −2 is the radius at which the effective logarithmic slope of the profile is −2.
For what regards the luminous mass components, the stellar disk mass density profile is assumed to take the form: where Σ d is the central disk surface density, R d and z d are the length scales in the cylindrical coordinates R and z, while R dm is the truncation radius. Among these, only Σ d and R d will be treated as free parameters. R dm will be assumed to scale with the local Galactocentric distance R 0 according to the prescription R dm = 12[1 + 0.07(R 0 − 8 kpc)] kpc [35]; we will also fix the vertical height scale to the best fit value suggested in Ref. [35], z d = 0.340 kpc, since the dynamical constraints we consider in our analysis are insensitive to a slight variation around this value. The bulge/bar region is modeled through the mass density profile: (2.10) The first term in Eq. (2.9) represents an axisymmetric nucleus [36], while the second term describes a triaxial bar. When comparing against dynamical constraints, we will implement an axisymmetrized version of Eq. (2.9), and assume x b ≃ y b = 0.9 kpc · (8 kpc/R 0 ), z b = 0.4 kpc·(8 kpc/R 0 ) and q b = 0.6. As for the stellar disk, we fixed these parameters because of the lack of observables to discriminate among these values and small deviations around them. Such values are in agreement with the best fit obtained in Ref. [37] assuming R 0 = 8 kpc, and then scaled to an arbitrary R 0 . Rather than using the two mass normalization scales Σ d andρ bb , in Paper I the parameter space of our reference Galactic model was scanned varying the disc mass M d and the balge/bar mass M bb ; here we re-parameterize these in terms of two dimensionless quantities, namely, the fraction of collapsed baryons f b and the ratio between the bulge/bar and disk masses Γ, i.e., respectively: In Eq. (2.11) we also included the sub-leading contributions to the total virial mass associated with the atomic (H I ) and the molecular (H 2 ) Galactic gas layers, with profiles as given in [38], see also Paper I. Finally, in our mass model, the baryons which do not collapse in the disc are distributed with the same profile of the DM component. The mass density profiles can be used as source terms in a Poisson equation to compute the Galactic gravitational potential. A rigorous procedure would require the solution of partial differential equations in cylindrical coordinates, a method which would be numerically challenging to apply to an extensive scan of our parameter space; moreover, it would actually give more information that those required in Eq. (2.1) which assumes a spherically symmetric Φ. In this respect, since the goal is to use Eddington's formula at our position in the Galaxy, i.e. in the outskirts of the stellar Galactic plane, it is a fairly good approximation to infer the local deepness of the gravitational potential well, and profile of Φ at larger radii, by applying the following prescription 2 : First, varying the Galactocentric distanceR, we calculate the total mass profile M (R) obtained integrating the sum of all density components we introduced within a sphere of radiusR centered in the Galactic center. This quantity is then used in a Poisson equation in the spherically symmetric limit, whose solution Φ(r) reads: With this last step, one has all the ingredients to estimate the DM phase-space distribution function according to Eq. (2.1). In particular, since we wish to discuss DM direct detection, we will focus on the local DM velocity distribution: having factorized out the local DM density One quantity which will be useful to compute to compare with previous analyses is the local velocity dispersion: while one of the advantages of Eddington's approach is that the local escape velocity v esc is not a free parameter, being instead directly related to the gravitation potential through the relation: As already mentioned, the DM velocity profile most often used in DM direct detection studies is the Maxwell-Boltzmann distribution. This is obtained from the distribution function of the isothermal sphere, taking a form in which radial and velocity dependence can be factorized: and 0 otherwise. Neglecting the small correction due to the truncation at the escape velocity (which is assumed here as an external input), this distribution has a constant velocity dispersion, which, after properly normalizing, is found to be equal to σ v . The isothermal sphere profile has an asymptotically flat rotation curve, with the circular velocity at large radii tending to the value Θ ∞ = 2/3 σ v . Assuming a flat rotation curve for the Milky Way down to our position in the Galaxy, that is Θ 0 = Θ ∞ , one usually constrains the Maxwell-Boltzmann velocity distribution imposing: Taking into account the escape velocity in the normalization factor, the local DM velocity distribution in the isothermal sphere case reads: where N esc is given by: 3 Direct-detection rates and phase-space densities The differential direct detection rate for a DM particle χ in a given material (per unit detector mass) is given by the convolution of the incident flux of the χ particles and the differential cross section for their scattering off the target nucleus dσ/dQ [39,40]: where Q is the energy deposited in the detector, M N the nucleus mass and ρ χ · f χ the χ phase-space density at the detector and in its reference frame. The lower limit in the integral is the minimum velocity required for a χ particle to deposit the energy Q. Its expression depends from the kinematics of the scattering: In the case of an elastic scattering, such as for WIMP-nucleus interactions, the minimum velocity is given by the χ-nucleus reduced mass, being M χ the mass of the particle χ. If the scattering is instead inelastic, such as for inelastic DM, the minimum velocity also depends from the mass split δ between the DM particle and the final state; in this case one finds In this paper we will mainly focus on spin-independent WIMP interactions. Rewriting accordingly the differential scattering cross section in terms of the WIMP-nucleon scattering cross section at zero momentum transfer σ χn and the nucleus form factor F 2 is the WIMP-nucleon reduced mass and m p the proton mass, the differential rate takes the form: where A is the mass number of the target nucleus. In the prefactor there are quantities depending on the particle physics model, while the integrand depends on astrophysical quantities only, though the limit of integration contains again a dependence on M N and M χ .
The local χ phase-space density is computed assuming that the χ particles account for the entire DM component and is directly associated to F h (E): with v(t) = v obs (t) + u. In the last equation v obs takes into account the motion of the detector frame compared to the Galactic frame in which F h is computed. One needs to take into account the motion of the Local Standard of Rest (LSR), namely the rotation of the sun around the Galactic center v LSR , the motion of the Earth around the Sun V ⊕ and the Sun's peculiar velocity with, in Galactic coordinates (wherex is the direction to the Galactic Center,ŷ the direction of disk rotation, andẑ the North Galactic Pole) and using (partially) the notation of Ref. [41]: In the computations the motion of the Earth around the Sun has been modeled as in the DarkSUSY package [42]. Neglecting here for simplicity the ellipticity of the Earth's orbit, one where V ⊕ = 29.8 km/s and t 1 = 0.218 is the fraction of the year before the Spring equinox, whileε 1 andε 2 are the directions of the Earth's velocity at times t 1 and t 1 + 0.25 years. From the previous equations, one can see that astrophysical uncertainties on the event rate computation are all encoded in one function, which we will denote by g χ (u, t) and it is defined as follows This function is the key object in our analysis. For convenience, we will refer to it as to the DM phase-space density function, though it is actually an angular integral of the true phase-space density boosted in the detector rest frame and multiplied by |u|.

Differential rate and exclusion limit for the Xenon100 experiment
In a real experiment Eq. (3.5) has to be modified in order to account for experimental limitations related to the efficiency of the detector and its finite energy resolution. In discussing the impact of astrophysical uncertainties on the differential event rate, we will consider a sample toy-model case, focussing on specifications which would apply, as a first approximation, to the case of the Xenon detector, taking a constant efficiency ǫ = 0.3 and including finite energy resolution effects by assuming a Gaussian probability density ξ(E, Q) that an event depositing an amount of energy Q in the detector is detected with an associated energy in the interval [E, E + dE] (the same assumptions are implemented for the Xenon experiment, e.g. see [44] and references therein). When plotting below observed differential event rates in an actual experiment, we will refer hence to the quantity: 3 with the energy dependent dispersion estimated for the Xenon detector: We will also discuss how astrophysical uncertainties affect the estimation of the exclusion limits in the scattering cross section -particle mass plane for a direct detection experiment, focusing, as for the differential event rate, on the sample case applicable to the Xenon detector. To do that, however, we need to treat finite energy resolution effects in a more refined way than what just sketched above. We will follow the approach discussed in [45], and take into account that the quantity measured by the Xenon experiment is the number of photoelectrons associated with an event, and not the corresponding recoil energy. To translate such a signal in to a recoil energy, one has to use the conversion formula [46] which gives the number of expected photoelectrons (PE) as a function of the recoil energy.
There has been some discussion in the literature on which is the correct form to assume for the relative scintillation efficiency of the detector L ef f . In our analysis we will refer to its most recent determination [46] providing an estimate for such a quantity at energies as low as 3 KeV, while at lower energies we follow [7] assuming that L ef f goes logarithmically to zero. The differential event rate is then given by the following expression where the convolution with the Poisson distribution of mean ν(Q) accounts for the fact that a given recoil energy has a probability ν(Q) n /n! e −ν(Q) of generating n observable PE. The Gaussian weighted sum, instead, describes the finite resolution in the observed number of PE: the probability of observing S 1 PE when the actual number is n is assumed to be Gaussian, with dispersion equal to √ n σ PMT and σ PMT = 0.5 PE. Following [7], we define as signal region the interval in the observed number of photoelectrons between 4 PE and 30 PE. In this region the Xenon collaboration has reported three candidate DM events associated with recoil energies of 12.1 KeV, 30.2 KeV and 34.6 KeV. To extract an upper bound on the cross section σ χn from the three Xenon100 events, one can calculate, as a function of σ χn , the probability -denoted by C(σ χn ) -that a random experiment produces a number of events in the signal region equal or larger than four. Then, one can exclude cross sections such that, for instance, C(σ χn ) ≥ 0.9. The probability C(σ χn ) can be written as the sum of different Poisson distributions of mean µ: where µ(σ χn ) = µ S (σ χn )+µ B is the total number of events in the signal region, including any known background term µ B . In the case of the Xenon100 experiment, we include a constant background contribution µ B equal to 1.8 events. The signal contribution, instead, is given by: where η = 48 × 100.9 [kg days] is the exposure of the Xenon100 experiment. Introducing the function: we can define the Xenon100 upper limitσ(M χ ) as: 4 where in the last line we emphasized the dependence of the exclusion limitσ from the DM particle mass.

Bayesian analysis of the parameter space
Eddington's inversion formula Eq. (2.1) establishes -under the assumptions discussed -a direct relation between the local DM phase-space density and the parameter space of the underlying Galactic model. In turn, the local DM phase-space density, along with the coordinate transformation in Eqs. (3.7) and (3.8), defines the factor in the DM scattering rate depending on astrophysical assumptions. To address the astrophysical uncertainties on DM direct detection we need therefore to study the variability range of the Galactic parameters. The model we summarized in Section 2 has 7 or 8 free parameters: the local Galactocentric distance R 0 ; the DM halo parameters c vir and M vir , plus the exponent α E in case of the Einasto profile; the parameters defining mass profiles for the baryonic components R d , f b and Γ (recall Eqs. (2.11) and (2.12)); the halo star anisotropy parameter β (introduced to make predictions for the velocity dispersion of halo stars and then compare it to the data of [47], see Paper I for further details). To set constraints on this parameter space, we consider a broad and variegated sample of dynamical constraints for the Milky Way, including several recent results, the same set of observables implemented in the analysis on the local halo density of Paper I (to which we refer for a comprehensive description): the proper motion of stars in the outer Galaxy, the radial velocity dispersion of halo stars, terminal velocities for the rotation curve in the inner part of the Galaxy, the monitoring of stellar orbits around the Galactic Center, the peculiar motion of SgrA * , the Oort's constants, the total mean surface density within |z| < 1.1kpc, the local disk surface mass density and, finally, the total mass inside 50 kpc and 100 kpc. The scan of the parameter space to address variability regions of the parameters is performed implementing a Bayesian analysis, again specular to the one carried out in Paper I. We sample the Galactic model posterior probability density function (pdf) implementing a Markov Chain Monte Carlo scan of the underlying parameter space. By construction, the pdf reflects the change of our prejudice about the most probable values for the model parameters g after seeing the data d. It is related to the Likelihood L(d|g) and the prior distribution π(g) through Bayes' theorem, which states where p(d) is the Bayesian evidence and in the present discussion simply plays the róle of a normalization constant. For a detailed description of our choice of the Likelihood and of the prior distribution we refer the reader to Paper I. In the present work, we generated 16 chains for each DM profile, accumulating (after burn in) approximately 1.5 × 10 5 accepted points in the parameter space for every considered profile. Rather than to the posterior pdf of the Galactic model parameters, in this work we are interested in the marginal one or two-dimensional pdf of specific functions of the Galactic model parameters. We therefore used in our analysis the fact that the marginal posterior pdf for any quantity f depending on the model parameters -denoted by p(f, g|d) -is related to the posterior pdf p(g|d) by the relation which follows from the definition of conditional pdf. Expectation values and variances of any function f of the model parameters with respect to p(g|d) can be then calculated as follows Eqs. (4.2) and (4.3) apply to functions of the Galactic model parameters only. The phase-space density, the direct detection differential event rate and the Xenon100 exclusion limit are instead all functions of both the Galactic model parameters and a further continuous variable. Respectively: the detector rest frame velocity u 5 , the observed recoil energy E and the DM particle mass M χ . To keep the present discussion as general as possible, we will denote during this section such an additional continuous variable by x.
The problem is now that of determining the marginal posterior pdf of objects ζfunctions of the Galactic model parameters -with an explicit dependence from x, i.e. ζ : x −→ ζ(x). In our analysis we compute such pdf's discretizing the variable x (see also [48] for a different approach to this problem). We therefore use Eq. (4.2) to calculate the marginal posterior pdf for a set of N quantities ζ i = ζ(x i ), with i = 1, · · · N , where now each ζ i is a function of the Galactic model parameters only. If the steps of such discretization {x 1 · · · x N } are small enough, this procedure leads to a good determination of the function ζ(x). As a result, to any given point x k in the range of variability of x, we can associate instead of a single number, a marginal posterior pdf, namely the pdf related to ζ k through Eq. (4.2).
Finally, given a marginal posterior pdf for the quantities ζ i , one can calculate from Eq. (4.3) the means ζ i and construct the required γ% credibility intervals, namely the shortest intervals including the γ% of the total probability. In our analysis we considered γ% equal to either 68% or 95%.

Results
The main result of the present analysis is a Bayesian determination of the local DM phasespace density based on the Eddington's inversion formula and a broad sample of dynamical constraints for the Milky Way. As already mentioned, this quantity incorporates all the astrophysical uncertainties which affect the expected event rate at a DM direct detection  Table 1. Means, standard deviations and confidence intervals for the circular velocity Θ 0 , the escape velocity v esc and the dark matter velocity dispersion σ v . We also show our results for the parameter f b , which affects the shape of the local gravitational potential and therefore the local dark matter phase-space density. Other Galactic model parameters have been already discussed in Paper I. experiment. We therefore first focus on this quantity; then, in a few illustrative examples, we also show how astrophysical uncertainties translates in to error bars for the predicted differential event rate and exclusion limit in the scattering cross section -WIMP mass plane. Concerning the exclusion limit, we will concentrate on the last Xenon100 data release, however, our results for the phase-space density can be readily applied to any other case.

Phase-space density
To better understand our results for the DM phase-space density, it is instructive to analyze first a few related local quantities, namely the circular velocity Θ 0 , the escape velocity v esc and, finally, the DM velocity dispersion σ v . In Table 1 we show means and credibility intervals for these quantities, while their marginal posterior pdf are shown in Fig. 1, where each quantity has been plotted in three cases, corresponding to the three DM profiles considered in this work. As we explained in detail in the previous section, such pdf's have been determined applying Eq. (4.2) to the case of the posterior pdf of the Galactic model parameters sampled through our MCMC scan. The local circular velocity (Fig. 1, left panel) has been very well reconstructed in all the three cases. Indeed the three means of Θ 0 are all close to the corresponding experimental value which has been implemented in the Likelihood function, i.e. Θ 0 = (245 ± 10.4) km/s [49]. The minor differences between the three cases are related to the fact that different profiles have slightly different mean values for R 0 and also different total mass within R 0 .
For what concern the escape velocity (Fig. 1, central panel), we find a mean value close to 555 km/s for the Einasto and the Burkert profiles, and to 550 km/s in the NFW case. Interestingly enough, although we did not consider any direct constraint on the escape velocity, our determination of this quantity is in a very good agreement with the central value quoted by the RAVE collaboration [50] (median likelihood of 544 km/s and 90% confidence interval in the range 498-608 km/s). We also notice that in the Einasto case the v esc pdf has a broader shape. This result is related to the larger uncertainties -compared to the other cases -in the fraction of collapsed baryons that we find in a scan based on an Einasto profile (see Table 1).
The DM velocity dispersion is connected to the width of the phase-space density as a function of the velocity. From Fig. 1 (right panel), we can see that in the case of a Burkert profile the mean value of σ v is significantly larger that in the Einasto and NFW cases. This result can be understood as follows. A spherically symmetric gravitational potential generates a velocity dispersion for the local population of DM particles which can be estimated through the Jeans equation. One finds By numerical inspection we find that the first derivative of the gravitational potential is larger -in the relevant integration region -for a Burkert profile than in the other cases. This is related to the fact that, according to our posterior pdf, this profile is in a better agreement with observations when a large fraction of collapsed baryons is considered, as one can see in Table 1. A larger value of f b -with respect to the other profiles -is in turn responsible for the larger velocity dispersion which emerges from our marginal posterior pdf 6 . We can now focus on the DM phase-space density. In Fig. 2 (left panel) we show the marginal posterior pdf for the DM phase-space density evaluated at a reference velocity in the detector rest frame (240 km/s). In the notation of the previous section, we are therefore considering the case in which x = u and ζ(x) = g χ (u). The u variable has been discretized and the reference velocity corresponds to one of the points of such a N-point discretization, i.e. u k = 240 km/s with k ∈ {1 · · · , N }. Let us note that the Burkert pdf is systematically shifted towards larger velocities. This is a feature which is present at each point of the discretization. From this marginal posterior pdf (Fig. 2, left panel) we can determine the  mean value of g χ (u k ) and the corresponding 68% and 95% credibility intervals. By proceeding in this way for all the N points in the discretization we can plot a phase-space density "band", instead of the usual phase-space density curve. This has been done in Fig. 3 for the three profiles under analysis. The red band is associated with the 68% credibility interval, while the yellow one corresponds to the 95% credibility interval. This figure makes quantitative what we already mentioned in the previous section: given a value of the velocity, say u, we can now associate to u a pdf -and not just a number -which accounts for astrophysical uncertainties in the determination of the phase-space density at that point.
In Fig. 2 (central panel) we show the "mean" phase-space densities corresponding to the three profiles considered in this paper. In this figure to each point of the velocity discretization we associated the mean of the corresponding pdf. It should be now clear that, because of a larger velocity dispersion, the Burkert case is characterized by a broader phase-space density.
We now compare our MCMC predictions with the usual Maxwell-Boltzmann approximation. As already mentioned, a Maxwell-Boltzmann distribution is identified by two parameters: the escape velocity v esc and the velocity dispersion σ v . Then, as explained in section 2, the Maxwell-Boltzmann distribution is normalized imposing σ v = 3/2 Θ 0 , and this parameter, together with the sharp cutoff imposed at the assumed value for the escape velocity, sets the shape of the high velocity tail of the distribution. In Fig. 2 (right-panel) we show the ratio between the distribution function found for an isothermal sphere profile (Maxwell-Boltzmann distribution times the local density, see Eq. (3.6)) with Θ 0 = 220 km/s, v esc = 544 km/s, ρ χ = 0.3 GeV/cm 3 and σ v = 3/2 Θ 0 , and the tails of our MCMC phase-space densities. From this figure one can see that different profiles lead to distributions behaving differently at high velocities, crossing each other in a non trivial way. Interestingly enough, we also notice that, in the high velocity tail, our MCMC mean distributions are systematically above the ones obtained in the Maxwell-Boltzmann approximation considered in this example (and often adopted in the literature as a benchmark distribution), with a mismatch which is significantly larger than the mismatch in the numerical value for the local halo density (which is about 0.4 GeV/cm 3 for the three profiles considered).
In Fig. 4 we compare our mean MCMC phase-space densities with the distribution functions of different isothermal spheres, one for each profile, constructed setting the required Galactic parameters at their mean values according to our MCMC scan (see Table 1). From this figure we can appreciate the fact that the phase-space density associated with a Burkert profile has a shape which -compared to the other cases -is more difficult to approximate with the distribution function of an isothermal sphere. As already pointed out, this result is related to the higher velocity dispersion found in the Burkert case.
Finally, in Fig. 5 we show the two-dimensional marginal posterior pdf in the (Θ 0 , σ v ) plane and compare the shape of this 2D-distribution with the straight line defined by σ v = 3/2 Θ 0 . As one can see, the correlation pattern between these parameters suggested by the data is different from the one described by this linear relation. Indeed, the eigenvalues of the covariance matrix for the parameters (Θ 0 , σ v ) point towards directions which are different from the one identified by the relation σ v = 3/2 Θ 0 .

Differential rate
As a first application of our results on the DM phase-space density, we now discuss how uncertainties on this quantity affect the expectations for the differential event rate at a direct detection experiment. The differential event rate can be calculated using Eq. (3.11) and depends on the Galactic model parameters through the phase-space density. It is however Similarly to what we did for the phase-space density, we now discretize the variable E and study a set of quantities defined by ζ k = dR dE (E k ), where E k is a generic point in the observed energy discretization. Thus, to each energy E k we can now associate a marginal posterior pdf. The associated credibility intervals define the "signal bands" which incorporate the astrophysical uncertainties.
In this section we present the results that we obtained following this approach for three illustrative cases: a 100 GeV WIMP (elastic scattering), a 10 GeV WIMP (elastic scattering) and, finally, a 45 GeV WIMP (inelastic scattering). The corresponding results are shown for the case of an Einasto profile in Fig. 6. In these plots a given credibility interval is the region between two curves of the same type.
From these plots one can appreciate that for light candidates (Fig. 6, central panel), and in general in the case of an inelastic scattering (Fig. 6, right panel), astrophysical uncertainties can be extremely important and completely spoil the accuracy of the predicted event rate. This phenomenon is related to the fact that such candidates -light or inelastic -probe only the high velocity tail of the phase-space density. Therefore, a few per cent change in v esc , for instance, can modify the integral over the phase-space density in Eq. (3.5) by a factor O(1). Instead, heavier candidates scattering off the target nuclei elastically are less sensitive to astrophysical uncertainties because these are sensitive to a much larger portion of the phase-space density. Thus, in this case, the uncertainties on the differential ratewhich are due to the uncertainties on the total area subtended by the phase-space densityare much less pronounced and are of the order of ∼ 20%.
In Fig. 7 we compare the differential rates obtained assuming different DM profiles both for the elastic and for the inelastic scenario. In these plots all the curves have been normalized to the differential rate computed assuming a Maxwell-Boltzmann distribution with Θ 0 = 220 km/s, v esc = 544 km/s and ρ χ = 0.3 GeV/cm 3 . In the case of a 100 GeV WIMP the curves shown in the plot differ mostly because the three corresponding profiles have different local DM densities. In the light and inelastic WIMP scenario, instead, the distance between the three curves also depend on the different high velocity behavior of the phase-space densities corresponding to the considered profiles. In all cases the Maxwell-Boltzmann approximation leads to lower differential rates.

Xenon100 exclusion limit
Finally, we study how astrophysical uncertainties affect the estimation of the exclusion limit in the scattering cross section -mass plane. In this section we focus on the recent Xenon100 data, implementing a method based on the Poisson statistics, as explained in section 3. We expect, however, similar conclusions, i.e similar credibility intervals, for different experiments analyzed with more refined statistical tools.
In the case of the exclusion limit, we discretize the mass M χ of the DM candidate and study the quantities ζ k =σ(M k χ ), where M k χ is a point in the discretization of the WIMP mass, andσ is the upper bound defined in Eq. (3.19). From the marginal posterior pdf's of the quantitiesσ(M k χ ) we derived a mean exclusion limit and the corresponding "exclusion bands" calculated from the credibility intervals associated with the differentσ(M k χ ). The results are shown in Fig. 8 (left panel), where we focus on the case of an Einasto profile. Together with the 95% MCMC credibility intervals we also plot for comparison the exclusion limits which one finds in the Maxwell-Boltzmann approximation with Galactic parameters set at their upper and lower 95% C.I. values (see Table 1). In agreement with our results for the differential events rates, uncertainties at small masses are rather important while at larger masses they are less strong (in this respect the logarithmic scale in the plot might be misleading).
In the right panel of Fig.8 we show the ratios between the MCMC exclusion limits computed for the Einasto, NFW and Burkert profiles, and the one obtained in the Maxwell-Boltzmann approximation with Θ 0 = 220 km/s, v esc = 544 km/s and ρ χ = 0.3 GeV/cm 3 . As for the differential rate, the curves differ because the three profiles are charcterized by different local densities and velocity distributions, being the tails of such distributions particularly relevant for low mass WIMPs. Fig.8 also shows that in all cases the Maxwell-Boltzmann approximation leads to upper bounds for the scattering cross section less strong than in the MCMC computation, in particular at low masses where the Maxwell-Boltzmann distribution is artificially truncated at an escape velocity v esc = 544 km/s. We also studied how astrophysical uncertainties affect the exclusion limit in the scattering cross section -mass plane for the inelastic DM scenario. The results obtained assuming an Einasto profile are shown in Fig.9 (left panel). As for the case of DM elastically scattering off the detector nuclei, together with the 95% MCMC credibility intervals we also plot the results derived in the Maxwell-Boltzmann approximation with Galactic parameters set at their upper and lower 95% C.I. values (see Table 1). As expected from our results on the differential rates, in the inelastic DM scenario uncertainties on the exclusion limit are considerably larger than in the case in which the DM -nucleus interaction occurs elastically. Finally, in the right panel of Fig.9, we show the ratios between the MCMC inelastic exclu- sion limits obtained for the Einasto, NFW and Burkert profiles, and the one computed in the Maxwell-Boltzmann approximation with Θ 0 = 220 km/s, v esc = 544 km/s and ρ χ = 0.3 GeV/cm 3 . As already mentioned presenting the analogous results for the elastic DM scenario, the curves differ bacause of the different local densities and velocity distributions characterizing the three considerd DM profiles.

Discussion and Conclusions
We have presented a new determination of the local DM phase-space density based on the Eddington's inversion formula. In this approach, assuming an isotropic velocity distribution, and in the limit of spherically symmetric DM density profile ρ h (r) and gravitational potential for the Galaxy Φ(r), there is an one-to-one correspondence between ρ h and the DM phasespace density.
Implementing, within a Bayesian framework, a broad and variegated sample of dynamical constraints for the Milky Way, and sampling through a MCMC algorithm the posterior pdf of the Galactic models parameters, we computed the local DM phase-space density. We have then applied this result to study a few illustrative cases of direct detection rates, and have discussed the impact on the exclusion limit one can infer from the latest data release from the Xenon experiment; the extension of our analysis to other DM models or to other experimental configurations is straightforward, and similar pattern are expected. The advantage of our approach is that it allows to keep a close inspection on the correlation among the astrophysical quantities setting the phase-space density of DM particles (and hence in turn the direct detection signals) and propagate uncertainties in a consistent way.
We studied separately three cases corresponding to three possible choices for the DM density profile, namely the NFW, Burkert and Einasto profiles. We have found that, on average, the velocity dispersion, which characterize the width in the velocity distribution, tends to be larger for the Burkert profile than for the NFW or the Einasto. Studying two-dimensional marginal posterior probability in the plane velocity dispersion versus circular velocity, we have verified that for any of the three distributions these two quantities do not follow the linear correlation pattern suggested in the isothermal sphere model, and adopted in most analyses of direct detection rates. We have computed the distributions of the escape velocities in our samples and found mean values which are consistent with the mean value obtained from observations of high velocity stars by the RAVE Survey; on the other hand, we have found that the shape of the high velocity tails in our distributions is sensibly different with respect to what you find assuming a Maxwell-Boltzmann distribution truncated at a given value of the escape velocity.
Concerning the DM direct detection signals, astrophysical uncertainties have the largest impact on models for which the scattering rate depends mainly on the tails of the phasespace density, namely for light WIMPs or DM candidates inelastically scattering, given that in such cases the integral of the phase-space density can be considerably different for different Galactic models. We have also compared our result with those obtained with a Maxwell-Boltzmann defined by setting circular velocity, escape velocity and local halo density to the "standard" reference values adopted in most direct detection studies, respectively, Θ 0 = 220 km/s, v esc = 544 km/s and ρ χ = 0.3 GeV/cm 3 . Regardless of the assumed DM profile, when adopting a self-consistent MCMC DM phase-space density, we have found that the differential rates are higher and the exclusion limits stronger than in the standard Maxwell-Boltzmann approximation.
Results in this paper are derived under the hypotheses of spherical symmetry and isotropy of the velocity distribution, two rather strong assumptions which may be (slightly) inaccurate for describing the DM distribution in the local neighborhood (see the note referring to this aspect in section 2). Having demonstrated here the power of the Bayesian approach in this analysis, we remark that such approach is general enough to be implemented within more general Galactic geometries, where extensions of the Eddington's inversion are required in order to relate the phase-space density to the parameter space of the underlying Galactic model; a closer inspection on this point will be carried out in future work. good for ∼100 GeV WIMPs and is accurate at the per cent level for light WIMPs or inelastic DM candidates in the range of energies where the corresponding signals are significantly different from zero (see Fig. 10). Error bars on direct detection signals can be then estimated calculating an upper limit σ (+) R and a lower limit σ (−) R analogously: where σ gχ can be computed as in Eqs. (4.3). In Fig. 10 we explicitly checked the validity of this approximation by studying how the quantities 7 depend from the recoil energy. As already mentioned, when the signal is non-negligible, the relative errors ∆ and ∆ ± are of a few per cent or less. Eqs. (A.1) and (A.2) should be therefore useful for a quick estimation of DM signals which depend from the local DM phase-space density. Files including time averaged g χ , g χ ± σ gχ and g χ ± 2σ gχ for the three considered DM profiles are provided with the electronic version of the paper.

B Data files provided with the paper
With the electronic version of the paper we provide fifteen data files. They contain different tabulations of g χ as a function of the velocity u in the detector rest frame as calculated by using Eq. (3.10). There are five files for each profile considered in the paper, corresponding to g χ , g χ ± σ gχ and g χ ± 2σ gχ respectively. Each file has a name of the form: uDF − halotype − curvetype.dat with halotype=nfw,burkert,einasto and curvetype=mean, 1splus, 1sminus, 2splus, 2sminus. 1splus and 1sminus correspond to g χ + σ gχ and g χ − σ gχ respectively, while 2splus and 2sminus are instead associated with g χ + 2σ gχ and g χ − 2σ gχ . The velocity is in km/s while g χ is in (GeV/cm 3 )(s 2 /km 2 ) and has been time averaged. The format of the files is the one required by the DarkSUSY function dshmuDFnum, which means that the first two lines are two empty headers and the data just start from the third line.

C The role of the priors
In this appendix we will focus on the impact of the choice of the priors on our conclusions.
In the present analysis we sampled the posterior pdf using flat priors for the Galactic model parameters, however, this does not guarantee that the priors actually imposed on the derived quantities are also non-informative. To determine the effective priors used for the local circular velocity, escape velocity, velocity dispersion and the values of the function g χ , following [51], we sampled the posterior pdf for these quantities with a MCMC scan where no constraints were imposed on the underlying Galactic model parameters. The pdf found in this way provides an estimation of the effective priors used for the derived quantities. In Fig.  11 one can appreciate that, though the actual priors on these quantities are not "rigorously flat", the obtained results are not dominated by the priors as well, since a high probability is a priori assigned to a large range of possible values for the derived quantities.