A self-consistent phase-space distribution function for the anisotropic Dark Matter halo of the Milky Way

Dark Matter (DM) direct detection experiments usually assume the simplest possible 'Standard Halo Model' for the Milky Way (MW) halo in which the velocity distribution is Maxwellian. This model assumes that the MW halo is an isotropic, isothermal sphere, hypotheses that are unlikely to be valid in reality. An alternative approach is to derive a self-consistent solution for a particular mass model of the MW (i.e. obtained from its gravitational potential) using the Eddington formalism, which assumes isotropy. In this paper we extend this approach to incorporate an anisotropic phase-space distribution function. We perform Bayesian scans over the parameters defining the mass model of the MW and parameterising the phase-space density, implementing constraints from a wide range of astronomical observations. The scans allow us to estimate the precision reached in the reconstruction of the velocity distribution (for different DM halo profiles). As expected, allowing for an anisotropic velocity tensor increases the uncertainty in the reconstruction of f(v) but the distribution can still be determined with a precision of a factor of 4-5. The mean velocity distribution resembles the isotropic case, however the amplitude of the high-velocity tail is up to a factor of 2 larger. Our results agree with the phenomenological parametrization proposed in Mao et al. (2013) as a good fit to N-body simulations (with or without baryons), since their velocity distribution is contained in our 68% credible interval.

The direct detection energy spectrum and its annual modulation, due to the Earth's orbit [24], depend on the local velocity distribution f (v) = F(x ⊙ , v)/ρ(x ⊙ ), where x ⊙ denotes the position of the Sun and ρ(x) is the DM density (see Refs. [25][26][27][28][29][30], among others). Direct detection experimental data are usually analysed assuming the so-called Standard Halo Model (SHM), which describes the MW as an isotropic isothermal sphere with local density ρ 0 = 0.3 GeV cm −3 and a Maxwellian-Boltzmann velocity distribution with velocity dispersion σ = Θ 0 / √ 2 and Θ 0 = 220 km s −1 . This model is unlikely to be a good approximation to the real MW DM halo. N-body simulations produce halos with velocity distributions which deviate systematically from a Maxwellian [26,31].
Finding an appropriate phase-space distribution for the DM halo of the MW when you know its gravitational potential (i.e. obtaining a complete dynamical model for the Galaxy from its mass model) can be done under certain simplifying assump-tions. For instance the phase-space distribution of a spherically symmetric system, with an isotropic velocity tensor, can be written as function of the energy E alone [7]. In this case one can solve for F(E), using the Eddington equation [32]. The solution will be self-consistent, in that F(E) and the gravitational potential of the system, Φ(x), satisfies Boltzmann's equation (e.g. Refs. [13,33,34]).
For the more general case of a spherical system with an anisotropic velocity tensor, the phase-space distribution function also depends on the modulus of the angular momentum, L. Often a parametric form is considered for F(E, L) and one can still find the set of parameters that corresponds to a self-consistent solution. Ref. [35] assumed that the phase-space distribution is separable in the two variables (i.e. F(E, L) = F E (E)F L (L)) and proposed a particularly convenient expression for F L (L) that depends on three parameters and can fit the radial dependence of the velocity anisotropy parameter β(r) in the case of halos formed in N-body simulations, where σ t and σ r are the tangential and radial velocity dispersions. In this paper, we apply the formalism developed in Ref. [35] in the context of the DM halos of galaxy clusters to the MW DM halo (see also recent work in Ref. [20] for an alternative approach to anisotropy). The first step is to build a mass model of the MW, c.f., e.g., Refs. [11,13,15,18], using a wide range of astronomical observations to constrain the gravitational potential of the MW and, therefore, the DM density profile. Our inferred knowledge of Φ(x) will then be used to derive self-consistent solutions for F(E, L), using the threeparameter form of F L (L) introduced in Ref. [35]. This approach can be thought of as a generalization of the Eddington equation to the case of a system with an anisotropic velocity tensor. Moreover, it extends the approach used when constructing a mass model of the MW, where the density profiles of the different components of the MW are written as functions of a number of free parameters which are constrained using astronomical observations. In this case, we also parametrize F L (L) and use our knowledge of the gravitational potential to derive self-consistent solutions for F (E, L), and therefore f (v). This is a different approach to that which has previously been used for anisotropic halos (e.g. Refs. [25,30,36]) where the components of the DM velocity dispersion tensor have been found by solving the Jeans equation [7]. The velocity distribution is reconstructed with a remarkable precision, but the resulting solutions are not necessary self-consistent. In our approach, the Jeans equation is automatically satisfied (thanks to the Jeans theorem) without having to impose it explicitly.
The paper is structured as follows. In Sec. II we introduce our mass model for the MW, listing the free parameters of the model in Sec. II A and the observations we use to constrain the parameters in Sec. II B, while Sec. II C describes the statistical techniques employed in the scans. Sec. III is devoted to the discussion of the resulting constraints on the mass model parameters. In Sec. IV we present our technique for obtaining self-consistent anisotropic phase-space distributions, F(E, L), and we apply it to the MW model found in the previous sections. In Sec. V we discuss our results and in Sec. VI we summarize our conclusions.

II. MASS MODELS OF THE MILKY WAY
In this section we discuss how we obtain a viable mass model for the MW. Our general approach follows previous work e.g., Refs. [11,13,15,18], with some differences in the details of how observations are implemented and in the modelling of the mass components.
The basic idea is to model the dynamically important components of the MW with physically motivated parametrisations, and then constrain the free parameters using a range of observations. We use a nested sampling algorithm to search the parameter space, and find the Bayesian probability distribution functions for the free parameters. If the observational data used are informative enough, the final results will be a precise MW model in agreement with observations, as well as estimates of the residual uncertainties in the model parameters.
In the follow subsection (Sec. II A) we describe the components of our MW mass model, including the free parameters. In Sec. II B we list and discuss the different observations and their implementation. Finally in Sec. II C we give some details about the sampling technique.

A. Milky-Way mass contributors
Our mass model of the MW follows closely that proposed in Ref. [11], and has four components: • stellar disk: Following Ref. [11], we model the stellar disk with a single component thin disc with density profile (in cylindrical coordinates): which is in agreement with the fit to the COBE Diffuse Infrared Background Experiment data discussed in Ref. [37]. The scale length in the z-direction is fixed at z d = 0.34 kpc, since the dynamical constraints considered here are insensitive to small variations in its value, while the normalization, Σ d , and the radial scale length, R d , are left as free parameters.
Ref. [15] considered a mass model with two disks, a thin and a thick one. Since the stellar components are not the focus of our investigation, we consider a model with a single disk which has fewer free parameters (see also Sec. II B 2).
The gravitational potential produced by Eq. (3) is axisymmetric (see, e.g. Ref. [19]), however, for simplicity, we work under the assumption of spherical symmetry, leaving the investigation of non-spherical Galactic models to future work. Thus, the disk gravitational potential at a certain distance r from the center of the MW can be well approximated by G ∞ r dr M d (< r)/r, where M d (< r) is the disk mass enclosed in a sphere of radius r. The deviation of the spherical gravitational potential from the true axisymmetric one (for the best-fit point for a Navarro-Frenk-White DM halo, see later) is maximal near the Galactic Center and is less than 10% for distances larger than 2.2 kpc 1 .
• bulge/bar: As in Ref. [11], we consider an axisymmetric version of the model proposed in Ref. [38]: with and The two terms represent the bar and the bulge, respectively. Their parameters are set to q b = 0.6, z b = 0.4 kpc (8 kpc/R 0 ) and x b = 0.9 kpc (8 kpc/R 0 ), rescaling to arbitrary R 0 , the distance of the Sun from the Galactic Center, as suggested by Ref. [38]. The normalization ρ bb (0) is left free. Another viable model for the bulge/bar system can be found in Ref. [39]. This component is always subdominant in our results and therefore we do not expect our results for the local DM distribution to be sensitive to the details of the bulge/bar density parametrization. As for the disk, the gravitational potential of the bulge/bar is assumed to be spherical and is obtained by computing the enclosed mass.
• interstellar medium: The model for the interstellar medium is kept fixed, without any free parameters. The mass density of molecular hydrogen H 2 , as well as the HI and HII components, is modelled as in Ref. [40], based on the observations presented in Ref. [41].
• DM halo: Insight into the density profile of the DM halo comes mainly from the study of synthetic halos formed in N-body simulations [42][43][44][45]. We consider three different spherically symmetric parametrizations for the DM density profile: a Navarro-Frenk-White (NFW) profile [46] ρ χ (r) = ρ s r r s an Einasto profile [47] and a Burkert profile [48] The scale radius, r s , is related to the radius at which the logarithmic derivative of the density profile is equal to −2, while ρ s fixes the normalization. The NFW and Burkert profiles only have two free parameters (r s and ρ s ) while the Einasto profile has an additional parameter, α, which controls the curvature of the profile. The NFW and Einasto profiles have inner cusps and provide good fits to the density profiles of halos formed in DM-only N-body simulations. Baryons are likely to play an important role in determining the DM distribution in the inner regions. However, simulating baryonic physics, and forming realistic galaxies, is a difficult problem (see, e.g., Refs. [49,50] for recent progress) and it is not yet clear how baryonic physics will affect the DM density profile. The Burkert profile has a central core, a possibility that seems to be preferred by observations of dwarf Spheroidal [51] or Low-Surface Brightness galaxies [52].
To compare with other mass models and other DM halo constraints present in the literature, we will also calculate the concentration parameter c = r vir /r s , where r vir is the virial radius, i.e. the radius within which the average density of the halo is ∆, the virial overdensity, times the critical density of the Universe. In a flat Λ Cold DM Universe with matter density parameter Ω m = 0.227 [53] at z = 0, ∆ = 94 [54] 2 .
Unlike Ref. [13] we do not include uncollapsed baryons in our MW model, as their distribution is highly uncertain and the majority are thought to be in the warm-hot intergalactic medium (see e.g. Ref. [55]). See Sec. III for further discussion.
There are five additional quantities that will be needed when implementing the observational constraints (see Sec. II B). These are R 0 , β * (the velocity anisotropy of stellar halo tracers, see Sec. II B 6) and the three components of the velocity of the Sun with respect to the Rotation Standard of Rest (see Sec.
, in a system of coordinates where the first axis points towards the Galactic Center, the second along the direction of Galactic rotation and the third is perpedicular to the Galactic plane.

B. Experimental constraints
As emphasised by Ref. [15] the fact that different experimental constraints use different underlying assumptions is an issue when constructing a MW mass model. In principle the best approach would be to use the raw data, rather than values of derived quantities, however in practice this is not possible. Still, where possible, we avoid using constraints which make specific assumptions, e.g. a fixed value of the solar radius, R 0 .

Local circular velocity
Measurements of the local circular velocity, Θ 0 , can be divided into two categories: those that measure the rotation velocity of the Sun, V φ,⊙ , by observing the proper motion of an object (or a population of objects) at rest at the Galactic Center, and those that measure the difference between the two Oort constant, A − B = Θ 0 /R 0 , from the proper motions of tracers.
Ref. [56] measured the proper motion of Sgr A ⋆ with an extremely good accuracy of approximately 0.4%: µ l = (−6.379±0.026) mas yr −1 . The local circular velocity can then be calculated using R 0 , and the velocity of the Sun with respect to the so-called Rotation Standard of Rest (RSR), V RSR ⊙ , i.e. the rotation velocity of a circular orbit in the axisymmetric approximation of the gravitation potential [57,58]. On the other hand, Ref. [59] measured A − B with 3% accuracy using the motion of 220 Cepheids detected by the Hipparcos satellite: The two techniques appear to lead to inconsistent values of Θ 0 , depending on the value of V RSR ⊙ assumed. Traditionally it has been assumed that the Local Standard of Rest (LSR, 2 Using the more recent value of Ω m from the Planck collaboration would not affect our results significantly. i.e. the orbit of local stars with "zero velocity dispersion", obtained by extrapolating to σ R = 0 the definition of the asymmetric drift, [58,60]) moves on a circular orbit. If this is true then the rotational component of the Sun's velocity with respect to the LSR, V LSR ⊙ , coincides with V RSR ⊙ , and can be used to estimate Θ 0 from the measurement of V φ,⊙ . Using V LSR ⊙ = 5 km s −1 , from the analysis in Ref. [5] of Hipparcos data and R 0 = 8.0 kpc from Ref. [61], Ref. [56] find Θ 0 /R 0 = (29.4 ± 0.2) km s −1 kpc −1 , significantly larger than the value quoted in Ref. [59] of (27.2 ± 0.9) km s −1 kpc −1 . However, using line-of-sight velocities of more than 3000 stars observed by the APOGEE survey, Ref. [58] found V RSR ⊙ = V φ,⊙ − Θ 0 = 23.9 +5.1 −0.5 km s −1 (assuming a flat rotation curve). This is significantly larger than even the revised value of V LSR ⊙ of 13 km s −1 [62], found taking into account the radial metallicity gradient. Ref. [58] discussed two possible reasons for this discrepancy. For instance, the LSR would differ from the RSR if the orbit of the LSR is not circular (due, for instance, to large-scale non-axisymmetric streaming motions). Alternatively, V LSR ⊙ could be significantly larger than previously thought.
Using the value of V RSR ⊙ = 23.9 km s −1 quoted in Ref. [58], the value of Θ 0 /R 0 found from the measurement of the proper motion of Sgr A ⋆ in Ref. [56] drops to 27.1 km s −1 kpc −1 , consistent with both the value in Ref. [59] value and the measurement of Θ 0 from Ref. [58].
In light of this, we constrain the circular velocity by imposing the measurement of the proper motion of Sgr A ⋆ by Ref. [56], assuming the value of V RSR ⊙ quoted before from Ref. [58] (see also Sec. II B 7). We also follow Ref. [11] by using A + B = ∂Θ(r)/∂r| r=R 0 = (0.18 ± 0.47) km s −1 kpc −1 .

Local surface density
The total integrated local surface density, within a vertical distance z of the Galactic plane, is defined as where ρ tot is the total mass density of the MW. A demonstration that this quantity continues increasing for z larger than a few times the disk scale length, z d , would be very strong evidence for the presence of DM at the Solar radius. The local surface density can be estimated by means of the Poisson equation, once the vertical force, K z , is determined. We use the values, derived from the kinematics of stellar tracers, from Ref. [63] and Ref. [64] of Σ * (R 0 ) = (48 ± 8)M ⊙ pc −2 and Σ(R 0 , 1.1 kpc) = (71 ± 6)M ⊙ pc −2 for the visible component and the total mass within 1.1 kpc, respectively. These values are consistent with more recent analyses, e.g. Refs. [65][66][67][68][69].
Ref. [70] used the data from Refs. [71][72][73] to derive lower limits on the local surface density up to 4 kpc. We do not use these results since they are, strictly speaking, only lower limits and also because of the large residuals in the fit to UW (see Fig. 2 of Ref. [71]).

Terminal velocities
The inner rotation curve of the MW, i.e. inside the Solar radius, can be constrained by measurements of the so-called terminal velocities: along each line-of-sight towards Galactic longitude l there is a point at which the distance from the Galactic Centre is smallest. Under the assumption of circular motion, the modulus of the line-of-sight velocity is largest for objects at this minimum distance that are moving on an orbit that is tangential to the line-of-sight. This maximal velocity is normally referred to as the terminal velocity and can be used to directly constrain the rotation curve of the MW, at that specific minimal distance.
Measurements of terminal velocities are obtained from the observation of the spectral line of atomic hydrogen HI (see, e.g. Ref. [74]) or of CO [75]. We consider the data set in Ref.
[74], excluding all the points with | sin l| < 0.35 • , where the assumption of circular motion is not valid due to the presence of the Galactic bar. A constant experimental error of 7 km s −1 is assumed for each of the remaining data points (following Ref. [11]).

Microlensing
Microlensing observations constrain the gravitational potential of the MW since they provide us with a probe of the mass density in compact objects in the direction(s) of observation. The impact of microlensing data on the reconstruction of the MW potential has been studied in Ref. [14]. We consider the same 10 measurements of the optical depth τ discussed in that paper, coming from the MACHO [76], OGLE-II [77] and EROS [78] collaborations. The distribution of the gravitational lenses is assumed to follow the matter density of the disk and the bulge/bar, while the distribution of sources depends on the particular microlensing events (see Ref. [14] for more details).

Proper motion of masers in high-mass star-forming regions
Accurate measurements of the proper motion (µ l and µ b ) and of the line-of-sight velocity with respect to the LSR (v LSR ) are available in the literature for a number of masers. Their positions can also be determined from their Celestial coordinates (which we consider to be exact) and their parallaxes, π. The masers are found to be in almost circular orbits and Ref. [79] showed how, once one fixes the peculiar (i.e. noncircular) velocity (U s , V s , W s ) of the maser, it is possible to predict its proper motion. The transformation from peculiar to proper motion requires knowledge of the position of the maser, the circular velocity at the position of the maser and the velocity of the Sun with respect to the RSR. The latter is included in our nuisance parameters (see Sec. II B 7) so that, for each choice of (U s , V s , W s ), the predicted proper motion can be compared to the measured values for µ l , µ b and v LSR , providing an indirect constraint on the rotation curve at the position of the maser.
Ref. [79] found that the 18 masers they analyzed were orbiting around the Galactic Center with (on average) V s ∼ −15 km s −1 . However Ref. [80] argue that it is more likely that V LSR ⊙ ≈ 11 km s −1 , as advocated by Ref. [62], rather than the value of V LSR ⊙ ≈ 5 km s −1 [5] used by Ref. [79]. Note that the value of V LSR ⊙ proposed in Ref. [62] is still smaller than the value in Ref. [58] of V RSR ⊙ = 23.9 km s −1 which we adopt, see Sec. II B 1.
We implement the information from the motions of masers, following Ref. [80], by marginalizing over the peculiar velocities and parallaxes of the masers. We assume that the components of the peculiar velocities have a Gaussian probability distribution with zero mean and ∆v = 10 km s −1 . A Gaussian distribution is also assumed for the parallax, with the mean and dispersion coinciding with the measured value of π and the experimental error, respectively, for each maser.

Velocity dispersion of halo stars
Ref. [91] selected a sample of more than 2000 Blue Horizontal-Branch stars observed by the Sloan Digital Sky Survey and computed the dispersion of the line-of-sight velocity in 10 bins in distance from the Galactic Center, from 5.0 to 60.0 kpc. They compare these measurements with the results of two cosmological galaxy formation simulations of MW-like galaxies to infer the rotation curve of the MW up to 60.0 kpc.
We follow Ref. [11] and directly use the binned line-ofsight velocity dispersion data. We do not use the results of the simulations, however our analysis is based on the following three assumptions, that receive validation from the simulations used in Ref. [91]: i) the dispersion of the line-ofsight velocity can be used as an estimate of the radial velocity dispersion, ii) the stellar density ρ * is well fitted by a r −3.5 power law and iii) the stellar velocity anisotropy parameter β * is constant over the range of distances considered. Under these assumptions, one can solve the Jeans equations for the Blue Horizontal-Branch stars analyzed by Ref. [91] to obtain an expression for the radial velocity diversion of the stars [11]: The 10 velocity dispersion data points in Ref. [91] can thus be used to constrain the circular velocity at large radii. Note that the constant velocity anisotropy parameter of the stars, β * , is a free parameter in our scans.

Nuisance parameters
The Sun's distance from the Galactic Center, R 0 , as well as the three components of the velocity of the Sun with respect to the RSR are included in the scan as nuisance parameters. We assume the U and W components of the Sun's velocity with respect to the RSR and LSR are identical (see Sec. II B 1).
We also assume a Gaussian probability distribution for each nuisance parameter, with a mean value and dispersion corresponding to the measured value and experimental error, respectively: • R 0 = (8.33 ± 0.35) kpc, inferred from the observation of stellar orbits around the Galactic Center [92].

Other observations
Refs. [93,94] estimated the total mass of the MW from the kinematics of satellite galaxies, globular clusters and (in the case of Ref. [94]) individual stars, considered as tracers of the underlying gravitational potential. Their results have an accuracy of approximately a factor of 2, pointing towards a total DM halo mass of around ∼ 10 12 M ⊙ . However these estimates are obtained using assumptions which we do not make in our analysis (e.g. a fixed value of R 0 and of the Solar RSR velocity). Moreover, the results in Ref. [93] become very prior dependent as soon as the value of Θ 0 is left free, as it is in our scans since it depends on the mass model. Thus, we decide not to consider such results.
Ref. [95] found, using high velocity stars from the RAVE survey, that the local escape velocity lies between 492 and 587 km s −1 (at 90% confidence). No constraint on the escape velocity is included in our scan since Ref. [15] argued that assuming that these stars are in a steady state (as done in Ref. [95]) is probably unrealistic. Moreover, the range quoted could be even larger depending on the parameterisation of the high speed tail of the stellar velocity distribution.
Finally, the mass model in Ref. [15] used two additional constraints from simulations. Namely the concentration of MW-like halos from Ref. [96] and the following relation for the ratio of stellar mass, M ⋆ , to DM mass, M χ ,: found by Ref. [97] to be a good fit to the Millennium-II simulation. We prefer to include only observational data and thus we do not consider the information coming from N-body simulations. However, in Sec. III we will show the predictions of our mass models for the observables discussed here (i.e. total mass of the MW within 50 kpc and within the virial radius, the escape velocity, the concentration of the MW DM halo and the DM/stellar mass ratio) but not included as data.

C. Sampling technique
The free parameters that we scan over are summarised in Tab. I. The parameters of our MW mass model were intro- are the components of the velocity of the Sun with respect to the RSR. The final section contains the parameters which appear in the part of the phase-space distribution function which depends on the angular momentum: β 0 and β ∞ are the velocity anisotropy of the DM halo at r = 0 and r = ∞ respectively and k L 0 is a proportionality constant between the parameter L 0 (entering in the definition of the phase-space density) and r s Θ(r s ) √ ln 2 − 1/2 (see Sec. IV). We do not indicate here the upper limit for β 0 since it depends on the form of the density profile (see Sec. IV). duced in Sec. II A (note that the quantity α is only applicable in the case of the Einasto profile). The distance of the Sun from the Galactic centre, R 0 , and the components of the Sun's motion with respect to the RSR, U RSR ⊙ , V RSR ⊙ and W RSR ⊙ , are considered as nuisance parameters (see Sec. II B 7). The last three parameters, β 0 , β ∞ and k L 0 , parameterise the part of the phase-space distribution function which depends on the angular momentum and will be introduced in detail in Sec. IV. β 0 and β ∞ are the velocity anisotropy parameters at the center of the MW and at infinity, respectively. The third parameter entering in the L-dependent part of the phase-space density is a characteristic scale L 0 . The parameter considered in the scan, however, is the ratio between L 0 and "scale angular momen- The second column indicates the range of values considered for each of the parameters and the third column indicates the shape of the prior probability distribution.
The scan is performed using the public code MultiNest [98], which uses a nested sampling algorithm to determine the Bayesian posterior probability distributions for the parameters in the scan and functions of these parameters. The core principle is Bayes' theorem: which shows how the so-called prior probability Pr(Θ|H) of the parameters Θ (describing our knowledge of the quantities in the context of model H before the implementation of the experimental data D) is updated once we consider the observational information encoded in the likelihood L(Θ) = Pr(D|Θ, H). The result is the posterior probability distribution function (pdf) Pr(Θ|D, H), which also depends on the socalled evidence Pr(D|H). We do not need to include the evidence as it only depends on the data and therefore acts as a normalization constant. Different choices of prior distributions can affect the final posterior pdf if the data considered are not constraining enough to overcome the shape of the prior. The third column in Tab. I indicates for each parameter whether we assumed a Gaussian prior distribution (denoted by Gaussian) or a uniform distribution on a linear scale (flat) or logarithmic scale (log). In the case of the NFW DM halo for some parameters we used both flat and log priors, in order to check that the posterior pdfs do not depend significantly on the choice of the prior. When not explicitly mentioned, all the results presented in the following sections refer to the choice of priors indicated in Tab. I.
We are interested in the probability distributions of particular subsets (normally 1-or 2-dimensional) of the parameters in Tab. I. Two different statistics can be considered to measure such distributions. The so-called posterior pdf of parameter Θ i is found by marginalizing over all the other parameters: where N is the total number of parameters in the scan. On the other hand, the profile likelihood (PL) of parameter Θ i is found by maximizing over the other parameters 3 : The PL is very sensitive to fine-tuned regions in the parameter space with a very large likelihood, while the pdf takes into account volume effects where large regions with a moderate likelihood are integrated over. It is normally useful to consider both quantities when studying the characteristics of a parameter space [99,100]. For each parameter we will compute the so-called x% credible interval, defined so that the each of the two tails outside the interval has a probability 0.5x%. We also determine the x% confidence levels as the region with a likelihood at most ∆χ 2 (x) lower than the likelihood of the best-fit point, where ∆χ 2 (x) is determined by solving where χ 2 n (y) is the χ 2 distribution with n degrees of freedom. We used the SuperBayes package 4 [101,102] to compute the credible intervals and confidence regions and produce the plots presented in the following sections.
We performed our scans using 2000 live-points and a tolerance of 10 −4 . Ref. [99] showed that such a set-up allows a good reconstruction of the PL in the context of SuperSymmetric models of Particle Physics. The pdf and PL distributions for our scans do not differ much from each other, which confirms that we have achieved a reliable evaluation of the PL.

III. PARAMETER ESTIMATION FOR THE MILKY WAY MASS MODELS
In this section we present the results of our scans, determining the parameters of our MW mass model. Tabs. II, III and IV summarize the results, including mean and best-fit values and 68% and 95% credible intervals, for the NFW, Einasto and Burkert DM halo profiles, respectively. Note that in order to facilitate comparisons with other studies we include derived values of quantities (marked with * ) that are not considered as experimental data in our scans.
The three mass models associated with the different DM halos share the following properties: in the (ρ s , r s ) plane, lines of constant M χ (< R 0 ), the DM mass within the solar radius, are approximately parallel to lines with constant Σ χ (R 0 , 1.1 kpc), as well as lines with constant ρ 0 . Thus, imposing the constraints on Σ(R 0 , 1.1 kpc) and Σ * (R 0 ) (fixing, indirectly, the local DM surface density), will select one of those lines (determining one degenerate direction in the plane (ρ s , r s )) and will translate into a determination of M χ (< R 0 ) and ρ 0 . The degeneracy is broken when we include the information on the velocity dispersion, which fixes the DM mass inside larger radii. On the other hand, in the (Σ d , R d ) plane, the lines of constant Σ d (R 0 , 1.1 kpc) are not parallel to those with constant M d (< R 0 ), so that the information on Σ * (R 0 ) and on Θ 0 (constraining the total amount of matter within R 0 ) act in a complementary way. The allowed region in the (Σ d , R d ) plane gets slightly larger when the terminal velocity data are included, since they prefer a less dense disk and a balance has to be made between the constraints on Σ * (R 0 ) and Θ 0 .
Microlensing and the proper motions of masers do not introduce significant addition information to the determination of the mass model. In the case of the microlensing data the baryonic component is already well determined by 4 http://www.ft.uam.es/personal/rruiz/superbayes/ the information on the local surface density and the rotation curve. While, for the masers, since we are marginalizing over the velocity of the peculiar motion of each maser, the scan practically has the effect of determining which values of (U s , V s , W s ) (for each object) provide a good fit to the data, given a rotation velocity fixed by the constraints on the local surface density and the rotation curve. This result was already discussed by Ref. [103].
The gravitational potential is dominated by the disk for R < R 0 , while the DM halo only becomes important around the Sun's radius. The bulge/bar is always subdominant: the probability distribution for ρ bb (0) is almost flat for values smaller than approximately 1 M ⊙ pc −3 and then goes rapidly to zero. This means that one should regard the upper end of the 68% and 95% credible contours for ρ bb as upper limits, while the lower end is practically determined by the range of values scanned. Fig. 1 shows the posterior pdf for the rotation curve Θ(r). The dark and light blue band indicate the 68% and 95% credible regions, respectively. The solid black line corresponds to the best-fit point and the fact that it falls within the 68% credible interval reflects the similar shape of the pdf and PL contours. The dashed (dotted) line shows the contribution of the DM halo (disk) to the MW rotation curve.
No tension is present between the results of the scans and the experimental data considered. At the best-fit point (independently of which profile is used for the DM halo) the χ 2 is dominated by the proper motion of masers, which is responsible for approximately half of the value of χ 2 .
The results of Tab. II for the case of a NFW DM halo are quite similar to the case of an Einasto halo, since the two parametrizations differ mainly in the inner region, where DM does not dominate the gravitational potential. On the other hand, a larger concentration and a smaller virial mass for the DM halo are obtained for the case of a Burkert profile. This is a consequence of the presence of the core: fixing the local DM density to a common value for the three halo models (as a result of the measurement of the local circular velocity and of the local surface density) implies a lower concentration for the Burkert DM halo. We stress that the goal of this paper is to infer local quantities (specifically the velocity distribution f (v)), not to determine whether a cuspy or cored profile provides a better fit to the data. Thus we leave a comparison of the goodness of fit of the different models for future work.
Our results are similar to those of Ref. [11]. The main difference is our value for the local circular velocity, which is ∼ 20 km/s smaller than theirs, independently of the profile chosen for the DM halo. This comes from the fact that Ref. [11] included a component of uncollapsed baryons, which they assumed follow the same density profile as the DM halo. It is thus expected that, for approximately the same values of other parameters, they find a larger circular velocity 5 . This is also related to the fact that the bulge/bar component in Ref. [  plays a significant role in the (very center of their) gravitational potential: in our case, the disk and DM components (fixed by the constraints on the local surface density) already saturates our smaller value of Θ 0 without leaving any room for a significant bulge/bar contribution. Ref. [15] also found a large value of Θ 0 , similar to that quoted in Ref. [11]. In the case of Ref. [15] we believe that this is due to their best-fit model having a larger contribution from the disk (in particular the thick one) which leads to a value of Σ * (R 0 ) larger than the measurement from Ref. [63], which we use. Note the mass of the bulge in Ref. [15] is fixed since this is one of their observational constraints.
Our local DM density is smaller than (but still compatible with) the value quoted in Ref. [18]: ρ 0 ≈ 0.5 GeV cm −3 . This is a direct consequence, again, of the value assumed for V φ,⊙ and the resulting Θ 0 .
The local DM density predicted for our mass models is also compatible with the values obtained by model-independent strategies, e.g., the solution of the equation of centrifugal equilibrium [12] or the Minimal Assumption method developed in Ref. [104] and applied to the study of approximately 2000 K dwarf stars in Ref. [105].
We now consider the quantities marked with * in Tabs. II,III and IV. We do not (for various reasons, as discussed in Sec. II B 8) include these quantities as data in our scans, instead we calculate the predicted values for our mass models. Our mass models predict that the mass within 50 kpc, M tot (r < 50 kpc), and the total mass, M tot , are lower than    found in Ref. [106] and Ref. [93] from the kinematics of satellite galaxies and halo stars. We believe this is a consequence of our smaller value of Θ 0 (or, equivalently, having assumed a large value of V RSR ⊙ ). Note also that our mass model ignores any interaction between the DM halo and the baryonic components. Hydrodynamic simulations of MW like galaxies, including DM and baryons, e.g. the Eris project [107], find that baryonic contraction, along with the formation of a DM disc, results in a larger local DM density than in DM-only simulations. However, the DM disc in the Eris simulations makes a relatively small contribution to the local DM density.

IV. SELF-CONSISTENT SOLUTION FOR THE PHASE-SPACE DENSITY OF DARK MATTER IN THE MILKY WAY HALO
In this section we summarize how we compute the DM phase-space density F(E, L) and, consequently, the DM velocity distribution f (v). We work under the following assumptions: • The system is in a steady state. In reality this assumption is unlikely to be satisfied completely. Simulated halos contain substructures and the high speed tail of the speed distribution has features [26] corresponding to incompletely phase-mixed DM, dubbed 'debris flow' [108]. However, at the Solar radius the dominant component of the DM distribution is expected to be smooth [31].
• The system is spherical and anisotropic, so that the phase-space distribution function is a function of only the energy and the angular momentum, which are integrals of motion. This also implies that the DM velocity tensor is diagonal in a spherical coordinate frame (i.e. the mixed terms σ i, j with i j are zero) and that the DM velocity moments satisfy the Jeans equation.
• The phase-space density is separable in the its two variables: F(E, L) = F E (E)F L (L). Ref. [35] tested this simplifying assumption qualitatively for simulated clustersized halos.
• F L (L) takes the form where the parameters β 0 , β ∞ and L 0 are defined in Sec. II C.
This ansatz was originally proposed in Ref. [35], which showed that the self-consistent solutions obtained from this assumption match the radial dependence of the anisotropy parameter, β(r), found in simulated clustersized halos. The clusters studied in Ref. [35] had an isotropic velocity distribution (i.e. β ≈ 0) close to their center, β ∼ 0.2 at the scale radius, increasing to a value as large as 0.4 at r ∼ 10 r s . Refs. [109,110] found similar behaviour in simulations of MW-like halos. Unlike the cluster-size halos studied by Ref. [35], in galaxy-sized halos β may decrease beyond r ∼ 5 r s [109,110]. However, this can be accommodated by the parametrization in Eq. (17). Thus we use Eq. (17) to parametrize the L-dependent part of the phase-space density for the MW DM halo.
Once the gravitational potential is fixed, and for a specific choice of the three parameters entering in Eq. 17, it is possible to determine F E (E) by inverting the following equation: This is what guarantees the self-consistency of our solutions, since, for a given F L (L), the phase-space density is completely determined by the gravitational potential of the system. Thus, reconstructing F(E, L) by inverting Eq. (18) can be considered as an extension of the Eddington formulism to the case of an anisotropic DM halo 6 . The solution has to be determined numerically since Eq. (18) is a Volterra integral equation. Details can be found in the appendix of Ref. [35]. F E (E) will depend on the gravitational potential of the system (this becomes evident as soon as one changes the integration variables in Eq. (18) to E and L) and on the choice made for β 0 , β ∞ and L 0 . The astronomical observations listed in Sec. II B constrain the gravitational potential without, however, providing any information on the parameters which appear in Eq. (17). Therefore we must marginalize over β 0 , β ∞ and L 0 , which can lead to large uncertainties in F(E, L). However, the relevant quantity for direct detection experiments is the local speed distribution f 1 (v), defined as 6 Note that the original Eddington formalism, for a spherical isotropic system, does not require any additional parameters and the speed distribution is completely determined by the gravitational potential [13,[32][33][34]. and f 1 (v) turns out to be much less dependent on the form of F L (L) than F E (E), and therefore it is possible to reconstruct f 1 (v) with reasonable uncertainties even when marginalizing over β 0 , β ∞ and L 0 . Fig. 2 shows how F E (E) depends on the gravitational potential of the system once the parameters entering in F L (L) are fixed. We only consider the effect of changing the characteristic density, ρ s , and the scale radius, r s , of the DM halo. For a given value of r s , decreasing ρ s corresponds to decreasing the amount of DM present in the halo and therefore the maximum value of the gravitational potential (i.e. the value at the center of the MW). This is why on decreasing the characteristic density, the range of E values over which F E (E) is defined becomes smaller. For ρ s below approximately 10 −1.5 M ⊙ pc −3 the baryonic component becomes important in particular for the largest values of energy, which correspond to the central region of the halo. Once the inner potential is dominated by baryons, the behaviour of F E (E) at high energies is independent of ρ s and the only remaining effect is at low energies, where F E (E) decreases as ρ s is decreased, simply because there is less DM.
The yellow bands show the behaviour of F E (E) as the scale radius is varied, with ρ s kept constant. Increasing r s leads to an increase in the amount of DM in the MW and therefore F E (E) moves to larger values.
Similarly Fig. 3 shows how f 1 (v) depends on the gravitational potential. For small values of ρ s the gravitational potential is completely determined by the baryonic components and the escape velocity is quite small, thus f 1 (v) goes to zero quite rapidly between 200 and 300 km s −1 . As the DM char- acteristic density is increased the range within which f 1 (v) is non zero gets larger and, therefore, the peak of the distribution is reduced since f 1 (v) is normalized to one. As in Fig. 2, increasing the scale radius (with ρ s fixed) increases the gravitational potential and the escape velocity, and therefore the peak in the speed distribution moves to higher speeds.
The yellow bands show the effect of changing β 0 and β ∞ . The effect is small and is localized at small/intermediate (intermediate/large) velocities for β 0 (β ∞ ). Note that the gravitational potential is fixed and therefore picking one specific velocity completely determines the energy of the DM particle. Small velocities correspond to large energies (remember that E = Φ(x ⊙ ) − v 2 /2) and to orbits localized close to the center of the halo. On the other hand, particles with large velocities will have small energies and will move on orbits that can reach large distances. Therefore the low (high)-velocity regime changes more as β 0 (β ∞ ) is varied.
For small negative values of L 0 increasing the transition scale has the same effect as decreasing β 0 , since it corresponds to reducing the portion of the halo with large β. For the same reason when the transition happens more or less at R 0 , the effect of increasing L 0 is the same as increasing β ∞ . Finally, for large L 0 , the velocity distribution becomes independent of L 0 since increasing the transition scale only affects large radii.
The ranges we take for the parameters defining F L (L) in Fig. 3 approximately match the ranges over which we carry out our scans (see Tab. I). The anisotropy at the center is constrained to be smaller than half the central slope of the DM halo profile [111]. This imposes an upper limit of 0.5 for the case of a NFW halo and only allows zero or negative β 0 for Einasto and Burkert profiles. However the halos formed in N-body simulations do not have negative β 0 , and therefore we only consider β 0 ≥ −0.5. We allow β ∞ to vary between -0.5 and 1.0, allowing for the possibility of β ∞ < β 0 .

V. RESULTS AND DISCUSSION
In this section we present our results for the probability distribution of the speed distribution f 1 (v), derived from the scans discussed in Sec. III. In Fig. 4, f 1 (v) is obtained through the Eddington formalism, assuming isotropy. The light (dark) blue bands show the 68% (95%) credible interval and the solid black line corresponds to the mean. The distribution is reconstructed with an uncertainty of a factor of 2 (at 68% credible interval) for speeds smaller than 500 km s −1 . This is consistent with the results of Refs. [13,33]. There seems to be a characteristic speed (around 300 km s −1 ) for which the speed distribution has a minimum uncertainty. This is particularly evident for the case of a NFW halo. We believe this is just a consequence of imposing the normalization of f 1 (v) to 1. Different mass models in the scan correspond to different gravitational potentials and, thus, to more peaked or extended speed distributions (see Fig. 3). However, increasing the amplitude of the peak has to lead to a less populated high-speed tail and, therefore, to a transition speed where smaller changes are observed 7 .
Our main result is displayed in Fig. 5 where we show the probability distribution of f 1 (v) obtained through the procedure summarized in the previous section (i.e. introducing a three-parameter form for F L (L), deriving F E (E) in a selfconsistent way and marginalizing over the three parameters). As in Fig. 4, the light (dark) green shows the uncertainty in the reconstruction by means of the 68% (95%) credible interval and the solid black line corresponds to the mean. As expected, the reconstruction is worse than the isotropic case in Fig. 4, due to the presence of the three parameters in F L (L), however the speed distribution can still be determined with an accuracy of a factor of 4-5 for velocities smaller than 500 km s −1 . The most uncertain part is, as before, the high-speed tail. The region of small uncertainty at around 300 km s −1 is more evident than in the isotropic case and an additional region with small uncertainty appears at around 100 km s −1 . These regions could already partially be seen in Fig. 3. The region at smaller (larger) speeds results from the effect of changing β 0 (β ∞ ) 8 .
The solid blue lines show the speed distribution obtained through the Eddington formalism as in Fig. 4. For a more physical comparison the two distributions should be normalized to the same value, and once that is done the difference is small. For all of the DM halo profiles considered the tail of the speed distribution, beyond ∼ 500 km s −1 , is larger in the anisotropic case, and consequently the peak of the distribution is lower. The effect is approximately a factor of 2 for the NFW halo, while it is less pronounced for the Einasto and Burkert halos (see Fig. 6). This is probably due to the marginalization over β ∞ and the uncertainty in the behaviour of the DM halo at large radii.
In the left panel of Fig. 5, for the case of the NFW halo, we also compare our results to the velocity distribution of the SHM (black dotted line) and the parametrization proposed by Mao et al. in Refs. [112,113] as a fit to the results of N-body simulations. This parameterization also provides a good fit to the DM speed distribution found in a hydrodynamical simulation of a MW-like galaxy containing baryons [107]. In each case the circular and escape speeds are set to the mean values from Tab. II and we set p, which parameterizes the shape of the high speed cut-off, to 1.5 as suggested by [113]. As already pointed out in Ref. [112], the self-consistent isotropic

VI. SUMMARY AND CONCLUSIONS
In this paper we developed a mass model for the MW, with the goal of studying the local DM velocity distribution in the case of a DM halo with an anisotropic velocity tensor. Our MW mass model, inspired by Refs. [11,15], assumes that the matter density of our own Galaxy can be written as a sum of four components: a disk, a combination of bulge and bar, a gaseous component and a DM halo. The free parameters entering in the modelling of these components are constrained by imposing a large set of astronomical observations (including the measurement of the proper motion of Sgr A * , the local surface density and terminal velocities, as well as information derived from the motion of masers and halo stars and from microlensing).
Such observations constrain the gravitational potential of the MW from which it is possible to reconstruct the DM phase-space density  Fig. 4 but for the speed distribution, f 1 (v), obtained using the procedure described in Sec. IV which allows the DM halo to be anisotropic. The light (dark) green band indicates the 68% (95%) credible region, while the solid black line corresponds to the mean. The solid blue line indicates the mean of the distribution of f 1 (v) obtained from the Eddington procedure (as in Fig. 4). In the left panel (for the NFW DM halo), the dotted line is the Standard Halo Model Maxwell-Boltzmann distribution, Eq. (1), while the dashed line shows the Mao et al. parametrization from Ref. [112]. In both cases the parameters are set to their mean values from Tab. II.
the assumption of an isotropic velocity tensor. The solution obtained will be self-consistent, since it depends entirely on the potential of the system. We extended this approach to the case of an anisotropic velocity tensor following the procedure introduced in Ref. [35], where the phase-space density F(E, L) is separable in the two variables F(E, L) = F E (E)F L (L). Ref. [35] parametrized the L-dependent component in terms of three quantities (β 0 , β ∞ and L 0 ) and showed that, once this is done, self-consistent solutions can be obtained for each given mass model simply by inverting a Volterra integral equation. The phase-space density obtained provides good fit to the radial dependence of the velocity anisotropy parameter measured in simulated DM halos surrounding galaxy clusters.
We apply the same procedure to the description of the MW DM halo, noting that this strategy extends the Eddington formalism to anisotropic scenarios (see also Ref. [20]). It fol-lows the same general approach employed when dealing with mass models of the MW: unknown quantities (e.g. the density profiles of the matter components of the MW) are parameterized and the parameters are constrained by a set of observations. There are no observations that directly constrain the three parameters introduced for F L (L), therefore they must be marginalized over. This could, in principle, spoil any hope to reconstruct the velocity distribution if f 1 (v) were to vary considerably as β 0 , β ∞ and L 0 are varied.
Our main conclusions are: • While F E (E) depends strongly on the values chosen for β 0 , β ∞ and L 0 , this is not the case for f 1 (v) and an acceptable reconstruction can still be achieved.
• The precision reached is, as expected, worse than when the Eddington formalism, which assumes isotropy, is used. However f 1 (v) is determined within a factor of a few (less than a factor of 5 below 500 km s −1 ), independently of which profile is chosen for the DM halo. The largest uncertainties are in the high-velocity tail, while there are two "sweet spots" (around 100 and 300 km s −1 ) where the precision is better than 10 − 15%. This is a consequence of the marginalization over β 0 , β ∞ and L 0 and the fact that f 1 (v) is, by definition, normalized to 1.
• The mean value of the distribution for the reconstructed f 1 (v) is very similar to that obtained assuming isotropy, apart from in the high speed tail. For high speeds the mean anisotropic speed distribution is larger (up to a factor of 2 for of a NFW halo) than the mean isotropic distribution derived with the Eddington formula (consequently the peak in the distribution is also lower).
• The parametrization of f 1 (v) introduced in Refs. [112,113], as a good fit to simulated DM halos, lies inside the band of our 68% credible interval.
Uncertainties in the local velocity distribution propagate in the analysis of direct detection data and may translate into large uncertainties in the reconstruction of the mass and interaction cross section of the DM particle, with the possibility of mis-reconstructions and biases (see, e.g., Refs. [23,25,34,36,[114][115][116], just to cite a few). Different experiments are sensitive to different regions of f 1 (v) so that they are affected differently by changes in the velocity distribution. Therefore realistic modelling of f 1 (v) is crucial when comparing results from different experiments.
The procedure presented here relaxes the assumption of isotropy in the analysis of mass models of the MW and allows, at the same time, a complete Bayesian assessment of the uncertainties involved and of the precision of the reconstruction of f 1 (v) (see also Ref. [20]). However, we have still assumed a spherical halo and a diagonal velocity tensor. N-body simulations and stellar dynamics have already demonstrated that these remaining assumptions are not true. The development of more general strategies, allowing a more flexible and realistic description of the MW, is therefore required. The present paper represents the first step towards this goal.
Also, realistic and theoretically-unbiased analysis techniques will be required in the near future, given the amount of information that will be delivered and the precision that will be reached by upcoming stellar surveys (e.g. the GAIA satellite [117,118]) or by the new data releases of already operational surveys, e.g. SEGUE [119], RAVE [120] and PanSTARRS [121].