Bayesian reconstruction of the Milky Way dark matter distribution

TL;DR

This paper introduces a Bayesian method to accurately infer the Milky Way's dark matter distribution from rotation curve data, accounting for data consistency and baryonic morphology uncertainties.

Contribution

It presents a novel Bayesian approach that improves dark matter profile inference by selecting consistent data subsets and testing multiple baryonic models.

Findings

Local dark matter density measured as 0.43±0.02(stat)±0.01(sys) GeV/cm³

Method demonstrates high accuracy and robustness on simulated data

Results are stable against various systematic errors

Abstract

We develop a novel Bayesian methodology aimed at reliably and precisely inferring the distribution of dark matter within the Milky Way using rotation curve data. We identify a subset of the available rotation curve tracers that are mutually consistent with each other, thus eliminating data sets that might suffer from systematic bias. We investigate three different models for the mass distribution of the luminous (baryonic) component that bracket the range of likely morphologies. We demonstrate the statistical performance of our method on simulated data in terms of coverage, fractional distance, and mean squared error. Applying it to Milky Way data we measure the local dark matter density at the solar circle $\rho_0$ to be $\rho_0 = 0.43\pm 0.02(\rm{stat})\pm0.01(\rm{sys})$ GeV/cm$^3$, with an accuracy $\sim$ 6%. This result is robust to the assumed baryonic morphology. The scale radius and inner slope of the dark matter profile are degenerate and cannot be individually determined with high accuracy. We show that these results are robust to several possible residual systematic errors in the rotation curve data.

Figures

Tables

Table 1. Table 1 : Fiducial parameter values used in mock RC data generation. We use every possible combination of the parameters in the table, producing 25 fiducial points in our 5 dimensional parameter space (see Section  3.2 for details).

$r_s$ [kpc]	5	8.5	12.5	16.5	20
$\gamma$	0.1	0.4	0.8	1.2	1.5
${\rho }_0$, $\left[\mathrm{GeV}/{\mathrm{cm}}^3\right]$		0.4
${\Sigma }_{\ast }$, $\left[M_{\odot }/{\mathrm{kpc}}^2\right]$		$38\times {10}^6$
		$2.17\times {10}^{-6}$

Table 2. Table 2 : The maximum likelihood (ML) and the maximum a posteriori (MAP) estimates with uncertainties obtained from the 68% HPD region (except for γ 𝛾 \gamma in the BR morphology case, where we indicate the 95% upper limit as the 1D marginal posterior is compatible with 0). We also give the reduced χ 2 superscript 𝜒 2 \chi^{2} values for the best fit (i.e. ML) parameters in the last column.

Disk morph.	ML $r_s$	ML $\gamma$	ML ${\rho }_0$	HPD68 $r_s$	HPD68 $\gamma$	HPD68 ${\rho }_0$	${\chi }_{\mathrm{red}}^2$
$\mathrm{BR}$ [57]	6.4	0.01	0.43	${7.1}_{-1.1}^{+2.3}$	<0.93 (95%)	${0.42}_{-0.02}^{+0.02}$	0.91
HG [52]	8.7	0.80	0.44	${6.3}_{-1.4}^{+7.2}$	${1.20}_{-0.58}^{+0.25}$	${0.44}_{-0.03}^{+0.02}$	0.94
CM [61]	10.1	1.11	0.44	${7.9}_{-3.0}^{+6.9}$	${1.44}_{-0.48}^{+0.16}$	${0.44}_{-0.02}^{+0.02}$	$1.01$

Table 3. Table 3 : The maximum likelihood (ML) and the maximum a posteriori (MAP) estimates with uncertainties obtained from the 68% HPD region (except for γ 𝛾 \gamma in the BR disk morphology cases, where we indicate the 95% upper limit as the 1D marginal posterior is compatible with 0). We also give the reduced χ 2 superscript 𝜒 2 \chi^{2} values for the best fit (i.e. ML) parameters in the last column.

Baryonic morph.	ML $r_s$	ML $\gamma$	ML ${\rho }_0$	HPD68 $r_s$	HPD68 $\gamma$	HPD68 ${\rho }_0$	${\chi }_{\mathrm{red}}^2$
G2 [50]BR [57]	6.5	0.01	0.43	${7.4}_{-1.4}^{+1.9}$	<0.91 (95%)	${0.43}_{-0.02}^{+0.02}$	0.94
E2 [50]BR [57]	6.4	0.01	0.43	${7.1}_{-1.1}^{+2.3}$	<0.93 (95%)	${0.42}_{-0.02}^{+0.02}$	0.91
V [66]BR [57]	6.4	0.0	0.43	${7.3}_{-1.3}^{+2.2}$	<0.99 (95%)	${0.43}_{-0.02}^{+0.02}$	0.94
BG [67]BR [57]	6.5	0.0	0.43	${7.1}_{-1.0}^{+2.0}$	<0.86 (95%)	${0.42}_{-0.02}^{+0.02}$	0.94
Z [68]BR [57]	6.4	0.01	0.43	${7.4}_{-1.5}^{+1.9}$	<0.93 (95%)	${0.43}_{-0.02}^{+0.02}$	0.94
R [69]BR [57]	6.3	0.01	0.43	${7.1}_{-1.0}^{+2.2}$	<0.91 (95%)	${0.43}_{-0.02}^{+0.02}$	0.94
G2 [50]HG [52]	10.4	1.0	0.44	${9.3}_{-4.0}^{+6.8}$	${1.44}_{-0.54}^{+0.11}$	${0.43}_{-0.02}^{+0.02}$	0.96
E2 [50]HG [52]	8.7	0.80	0.44	${6.3}_{-1.4}^{+7.2}$	${1.20}_{-0.58}^{+0.25}$	${0.44}_{-0.03}^{+0.02}$	0.94
V [66]HG [52]	10.6	1.05	0.43	${9.6}_{-3.8}^{+6.8}$	${1.39}_{-0.46}^{+0.17}$	${0.43}_{-0.02}^{+0.02}$	0.95
BG [67]HG [52]	10.0	0.96	0.43	${9.2}_{-4.3}^{+6.1}$	${1.28}_{-0.47}^{+0.23}$	${0.42}_{-0.02}^{+0.03}$	0.95
Z [68]HG [52]	10.5	1.03	0.44	${9.6}_{-4.6}^{+5.8}$	${1.25}_{-0.42}^{+0.28}$	${0.43}_{-0.02}^{+0.02}$	0.94
R [69]HG [52]	9.9	0.96	0.44	${8.7}_{-3.7}^{+5.6}$	${1.25}_{-0.51}^{+0.23}$	${0.44}_{-0.03}^{+0.02}$	0.94
G2 [50]CM [61]	14.3	1.39	0.43	${8.2}_{-2.1}^{+8.9}$	${1.46}_{-0.33}^{+0.16}$	${0.43}_{-0.02}^{+0.02}$	1.00
E2 [50]CM [61]	10.1	1.11	0.44	${7.9}_{-3.0}^{+6.9}$	${1.44}_{-0.48}^{+0.16}$	${0.44}_{-0.02}^{+0.02}$	1.01
V [66]CM [61]	14.2	1.38	0.43	${11.5}_{-4.2}^{+8.7}$	${1.50}_{-0.25}^{+0.16}$	${0.43}_{-0.02}^{+0.02}$	1.00
BG [67]CM [61]	11.0	1.19	0.44	${9.3}_{-4.5}^{+6.0}$	${1.29}_{-0.28}^{+0.30}$	${0.44}_{-0.02}^{+0.02}$	1.00
Z [68]CM [61]	11.6	1.25	0.44	${10.7}_{-5.6}^{+6.7}$	${1.51}_{-0.41}^{+0.12}$	${0.44}_{-0.03}^{+0.01}$	1.00
R [69]CM [61]	11.3	1.20	0.44	${10.7}_{-5.7}^{+4.6}$	${1.35}_{-0.35}^{+0.24}$	${0.43}_{-0.02}^{+0.02}$	1.00
G2 [50]dJ [70]	13.1	1.29	0.43	${10.2}_{-3.8}^{+9.6}$	${1.47}_{-0.32}^{+0.19}$	${0.43}_{-0.02}^{+0.03}$	0.99
E2 [50]dJ [70]	9.9	1.07	0.44	${7.8}_{-2.5}^{+7.3}$	${1.32}_{-0.41}^{+0.23}$	${0.43}_{-0.02}^{+0.02}$	0.98
V [66]J [70]	14.7	1.38	0.43	${11.8}_{-5.3}^{+6.6}$	${1.45}_{-0.33}^{+0.17}$	${0.43}_{-0.03}^{+0.02}$	0.99
BG [67]dJ [70]	11.5	1.20	0.43	${10.3}_{-5.0}^{+6.4}$	${1.41}_{-0.36}^{+0.21}$	${0.43}_{-0.02}^{+0.02}$	0.99
Z [68]dJ [70]	12.3	1.24	0.43	${11.1}_{-4.8}^{+6.9}$	${1.48}_{-0.38}^{+0.13}$	${0.43}_{-0.02}^{+0.02}$	0.99
R [69]dJ [70]	11.4	1.17	0.44	${8.6}_{-3.8}^{+7.2}$	${1.31}_{-0.33}^{+0.29}$	${0.44}_{-0.02}^{+0.02}$	0.99
G2 [50]J [71]	15.2	1.31	0.42	${11.1}_{-4.7}^{+8.0}$	${1.41}_{-0.33}^{+0.18}$	${0.42}_{-0.02}^{+0.02}$	0.96
E2 [50]J [71]	11.0	1.05	0.43	${7.8}_{-2.6}^{+8.0}$	${1.36}_{-0.55}^{+0.20}$	${0.43}_{-0.03}^{+0.02}$	0.95
V [66]J [71]	14.8	1.31	0.42	${12.5}_{-6.1}^{+6.9}$	${1.39}_{-0.32}^{+0.20}$	${0.42}_{-0.02}^{+0.02}$	0.96
BG [67]J [71]	12.3	1.14	0.42	${10.2}_{-4.8}^{+6.7}$	${1.39}_{-0.42}^{+0.18}$	${0.42}_{-0.02}^{+0.02}$	0.96
Z [68]J [71]	13.3	1.23	0.42	${9.4}_{-3.7}^{+7.2}$	${1.48}_{-0.51}^{+0.11}$	${0.42}_{-0.02}^{+0.03}$	0.96
R [69]J [71]	11.9	1.14	0.43	${8.1}_{-1.8}^{+8.9}$	${1.32}_{-0.43}^{+0.21}$	${0.43}_{-0.02}^{+0.02}$	0.96

Equations

$$\begin{gathered} \bar{V}=\frac{{\sum }_{i=1}^{N_{bin}}{V}_i/{\sigma }_i^2}{{\sum }_{i=1}^{N_{bin}}1/{\sigma }_i^2}\mathrm{and} \\ {\sigma }_{\bar{V}}^2=\frac{{\sum }_{i=1}^{N_{bin}}(\bar{V}-V_i{)}^2/{\sigma }_i^2}{{\sum }_{i=1}^{N_{bin}}1/{\sigma }_i^2} \end{gathered}$$

$${\Lambda }_{\mathcal{M}}=\frac{p({\mathrm{d}}_1,{\mathrm{d}}_2|\mathcal{M})}{p({\mathrm{d}}_1|\mathcal{M})\mathrm{p}({\mathrm{d}}_2|\mathcal{M})},$$

$${\rho }_{gNFW}(r)={\rho }_0{\left(\frac{R_0}{r}\right)}^{\gamma }{\left(\frac{r_s+R_0}{r_s+r}\right)}^{(3-\gamma )},$$

$$M_{DM}(<r)=4\pi {\int }_0^r{\rho }_{gNFW}(r)R^{2}dR.$$

$${\omega }_{DM}^2=\frac{G{M}_{DM}(<r)}{r^3}.$$

$${\omega }_c^2(\Theta ,{\Sigma }_{\ast },⟨\tau ⟩)={\omega }_{disk}^2({\Sigma }_{\ast })+{\omega }_{gas}^2+{\omega }_{bulge}^2({\Sigma }_{\ast },⟨\tau ⟩)+{\omega }_{DM}^2(\Theta ),$$

$$0<r_s/[\mathrm{kpc}]<40,$$

$$0<\gamma <3,$$

$$0<{\rho }_0/[\mathrm{GeV}\mathrm{c}{\mathrm{m}}^{-3}]<1,$$

$$19\times 10^6<{\Sigma }_{\ast }[{\mathrm{M}}_{\odot }\mathrm{kp}{\mathrm{c}}^2]<57\times 10^6,$$

$$0.1\times 10^{-6}<⟨\tau ⟩<4.5\times 10^{-6}.$$

$$\begin{array}{rl}P(d|\Theta ,{\Sigma }_{\ast },⟨\tau ⟩)& =\prod _{i=1}^m\left\{\frac{1}{\sqrt{2\pi }{\sigma }_{\bar{\omega },i}}\exp \left[-\frac{1}{2}\frac{{\left({\omega }_c(r_i,\Theta ,{\Sigma }_{\ast },⟨\tau ⟩)-{\bar{\omega }}_i\right)}^2}{{\sigma }_{\bar{\omega },i}^2}\right]\right\}\\ & \times \frac{1}{\sqrt{2\pi }{\sigma }_{⟨\tau ⟩}}\exp \left[-\frac{1}{2}\frac{\left(⟨\tau ⟩-⟨\tau {⟩}^{obs}\right)}{{\sigma }_{⟨\tau ⟩}^2}\right]\\ & \times \frac{1}{\sqrt{2\pi }{\sigma }_{{\Sigma }_{\ast }}}\exp \left[-\frac{1}{2}\frac{\left({\Sigma }_{\ast }-{\Sigma }_{\ast }^{obs}\right)}{{\sigma }_{{\Sigma }_{\ast }}^2}\right],\end{array}$$

$$P(\Theta ,{\Sigma }_{\ast },⟨\tau ⟩|\mathrm{d},\mathcal{M})=\frac{P(\mathrm{d}|\Theta ,{\mathrm{\Sigma }}_{\ast },⟨\tau ⟩,\mathcal{M})\mathrm{P}(\mathrm{\Theta },{\mathrm{\Sigma }}_{\ast },⟨\tau ⟩|\mathcal{M})}{P(\mathrm{d}|\mathcal{M})},$$

$$\mathrm{P}\left({\bar{\omega }}^{ij}|\Theta ,{\Sigma }_{\ast },⟨\tau ⟩,ϵ\right)=\frac{\sqrt{{ϵ}_j}}{\sqrt{2\pi }{\sigma }_{\bar{\omega },ij}}\exp \left[-\frac{1}{2}\frac{{\left({\omega }_c(r_i,\Theta ,{\Sigma }_{\ast },⟨\tau ⟩)-{\bar{\omega }}^{ij}\right)}^2}{{\sigma }_{\bar{\omega },ij}^2/{ϵ}_j}\right],$$

$${\rm P}(\log{\epsilon_{j}})=\left\{\begin{array}[]{ll}2/3\hskip 17.07182pt{\rm for-3/2\leq\log\epsilon_{j}}\leq 0\\ 0\hskip 28.45274pt{\rm otherwise}\end{array}\right.$$

$$\mathrm{P}({ϵ}_j)\propto 1/{ϵ}_j.$$

$$\mathrm{FSE}=\frac{\sqrt{\frac{1}{n_i}{\sum }_{i=1}^{n_i}({\hat{\theta }}_i-{\theta }_{\text{true}}{)}^2}}{{\theta }_{\text{true}}}.$$

Full text

1 Introduction

Dark matter is a crucial ingredient for the formation of structures in our universe. Without dark matter it would be impossible for small density perturbations to grow into the huge potential wells that host the galaxies we observe today. Dark matter thus constitutes the invisible gravitational backbone of our universe. It cannot, however, be accommodated within the Standard Model of particle physics and large experimental and theoretical efforts are underway to identify its nature. At the same time astronomical observations of increasing refinement are characterizing the distribution of dark matter within the visible galaxies. On one hand, such observations may be useful to characterize ancillary quantities required to interpret dark matter detection experiments. On the other, the distribution of dark matter within galaxies of different sizes is a prediction of the $\Lambda$CDM model, and measuring these distributions through astrophysical analyses provide an important test of consistency of the cosmological framework.

Different methods are adopted in the literature to probe the dark matter distributions in galaxies and clusters of diverse sizes, ranging from gravitational lensing for clusters [1, 2, 3], Jeans analysis of the motions of stars for dwarf and elliptical galaxies [4, 5, 6], and the motion of tracer objects in circular orbits in disk galaxies [7, 8, 9]. The latter is often referred to as the Rotation Curve (RC) method, and relies on the assumption that the stellar disk of a galaxy is rotation supported, and that the tracers are in circular motion around the gravitational center of the galaxy. These assumptions typically hold in disk galaxies over a large range of masses including our very own, the Milky Way.

The dark matter distribution (often referred to as the “density profile”) is inferred by fitting the observed velocity distribution (as a function of the galactocentric distance) with the one expected from the observed, luminous matter plus the “missing component” [10].

The values (and uncertainties) resulting from this procedure are usually adopted, and a dark matter density profile for the Galaxy at hand assumed from this. Studies in the literature [11, 12, 13, 14, 15] usually assume spherical distribution for dark matter, and one (or at most few) possibilities for the underlying visible component.

In this paper we revisit this well known method. First we propose a new Bayesian procedure to characterize the Milky Way’s dark matter density profile based on the current state-of-the-art for both the observed rotation curve and the luminous (hereafter, baryonic) component of the Milky Way. In this context, we also study the possible incompatibility of different data sets using a Bayesian procedure in order to rule out the possibility of systematic discrepancies between data sets.

Second, we study the accuracy of such a method: we propose for the first time the validation of the RC method in the Milky Way by generating mock rotation curves constructed on the basis of a known, controlled dark matter profile, and holding the same statistical properties as the observed ones. We test their reconstruction through the RC method against the known truth. The two separate procedures that we have developed and present here can be easily applied to new incoming data sets for the Milky Way, as well as to other galaxies.

The paper is structured as follows: in Section 2 we present our astrophysical setup, namely the set of observations adopted for the rotation curve and for the luminous (baryonic) component of the Milky Way; in Section 3 we present our methodology for reconstruction and describe our statistical and numerical tools. Section 4 contains an extended discussion on our main findings, and we summarize our conclusions in Section 5. In Section 4.2 we directly compare the results of our procedure for three extreme baryonic morphologies, while in Appendix B we show the results for the whole set of different baryonic morphologies collected in the literature and presented in [16].

2 Astrophysical setup

The rotation curve (RC) of a galaxy provides important constraints on its mass distribution, including its dark matter content (see e.g. [17, 18, 11]). Paradoxically, the RC of the Milky Way is much harder to measure than those of external galaxies. This is due to our interior position that complicates some measurements, such as the extended RC of the gas in its disk. Hence, in order to measure the RC of both the inner and outer regions one can use two different classes of tracers that we summarize below.

In what follows we adopt the following fixed Galactic parameters $R_{0}=8.34$ kpc, $V_{0}=239.89\;\rm{km/s}$, ($U_{\odot}$, $V_{\odot}$, $W_{\odot}$)=(7.01, 12.20, 4.95) $\rm{km/s}$. These are, respectively, the galactocentric distance to the Sun, the circular velocity at the Sun’s position, and the local standard of rest. Given these fixed parameters, we explore several important sources of uncertainty, including systematic errors in the rotation curve data, the dependency on the morphological models we assume for the disk component, and the observational uncertainties on the normalizations of the baryonic components.

It is to be noticed that the uncertainty existing in the literature about the Galactic parameters does not impact the general conclusions on the presence of dark matter even within the solar circle (see Supplementary Information of [16] for details). However, such uncertainties can affect the actual determination of the dark matter distribution (see e.g. [19]), and therefore they must be taken into account by dedicated analysis, also in light of its effects on searches for the very nature of the dark matter [20, 21]. A quantitative investigation of the effect of Galactic parameter uncertainties will be addressed in a future work  [22].

2.1 The observed rotation curve

For the innermost region of the Galaxy (galactocentric radii in the range $2.5-22$ kpc) we adopt the galkin data compilation, first presented as a whole in [16]. This compilation includes kinematics of gas, masers, and stars for a total of 2780 tracers in the disk, separated into 25 individual datasets, collected from almost four decades of literature. We address the reader to the public release of this data (and of the online tool to manage them) [24] for further details on the individual datasets, their nature, and original literature sources.

A close inspection of the data (plotted as gray data points in Fig. 1) shows that, due to their heterogeneous nature, some individual data sets show marked discrepancies with others. For instance, in the region around $R=10\leavevmode\nobreak\ \rm kpc$ we observe data points that are clearly discrepant (by many standard deviations) from the bulk of the observations. Such data points might suffer from distinctive systematic effects, or other data-collection artifacts, and it would be unwise to simply bin all the data together lest the outliers exert a disproportionate (and incorrect) influence on the bin average. In particular, most data sets that we have excluded are based on observations of such tracers as HI regions and CO emission associated with HII regions. The offsets between different data sets based on similar type of tracers might arise due to several reasons, as for instance: differences in the method used to estimate the maximum line-of-sight velocities (see discussion in [25]), differences in the Galactic longitudes at which data were taken (the long-known asymmetry between the southern and the northern rotation curves, see e.g. [26, 27]). Therefore it is necessary to invoke a statistical procedure to decide which subset of the 25 data sets are mutually compatible. The aim of the procedure is to weed out highly discrepant data sets, whose inclusion would systematically bias the bin value. To this end, we use a method that uses the Bayesian evidence (or model likelihood) to determine the degree of overlap in parameter space between parameter constraints (in a given parametric model) obtained using two different data sets. The idea is that if the data sets are compatible, constraints obtained by using them jointly will overlap significantly with those obtained when using each data set individually (see e.g. [28, 29]). The details on the procedure are given in Section 2.1.1 below.

After determining the subset of data within galkin that is mutually compatible via the above procedure, we carry out error-weighted binning of the data in radius. This is necessary in order to use a binned likelihood approach across the entire range of radial distances. We require an approximately equal number of data points per bin in order to obtain well-balanced bins errors. The number of bins is chosen such that the size of it is approximately equal to the bin size of [23]. This yields to 30 bins in total. The mean value in each bin is calculated using the weighted mean, and the corresponding standard deviation by using the weighted standard deviation:

[TABLE]

where $\bar{V}$ is the weighted mean of the observed circular velocity, $\sigma_{i}$ is the uncertainty in the circular velocity measurement $V_{i}$, and $N_{bin}$ is the number of data points of the $i^{\rm th}$ bin. We have verified that this procedure leads to a scatter within each bin that closely follows a Gaussian distribution. This is a necessary condition for our likelihood function (defined later on) to be a good representation of the data in each bin.

For the intermediate region of the Milky Way we use a data set obtained from the data compilation of Huang et al. [23]. The authors use $\sim$ 16,000 primary red clump giants (PRCGs) to derive a binned determination of the RC of the outer stellar disk (galactocentric distances in the range $8.34-14.95$ kpc) (see Fig. 1).

For the outer Milky Way (with radius in the range $\sim 15$ to $\sim$ 100 kpc) we use the binned RC measurements obtained by Huang et al. [23] from $\sim$5,700 so-called halo K giant (HKG) stars. These stars, differently from the tracers in the disk, do not follow circular orbits. Therefore, in order to derive the RC one needs to apply the Jeans analysis (see [23] for more details). We plot the resulting RC measurement on Fig. 1.

2.1.1 Compatibility of data sets via the Bayesian evidence

Our method to select mutually compatible data sets uses a Bayesian model comparison framework. For further details about the methodology, see [28, 30, 29, 31, 32, 33, 34, 35]. Given a data set $\rm d_{1}$, we consider whether a second data set $\rm d_{2}$ is compatible, assuming that the two data sets are statistically independent. We perform a model comparison and evaluate the model likelihood ratio $\Lambda_{\mathcal{M}}$, defined as

[TABLE]

where $p(\rm{d}|\mathcal{M})$ is the Bayesian evidence (or model likelihood) for data $\rm{d}$ given model $\mathcal{M}$. In the above, $\rm{d}_{1}$, $\rm{d}_{2}$ are the two data sets being tested for compatibility; model $\mathcal{M}$ has parameters $\Theta$ (as described below in Section 3), and in the numerator $p(\rm{d_{1}},\rm{d_{2}}|\mathcal{M})$ is the probability of obtaining both data sets if they both come from the same underlying parameters. The denominator is the probability of obtaining both data set when each of them is separately described by its own set of parameters, $\Theta$. The model likelihood ratio $\Lambda_{\mathcal{M}}$ thus represents the change in odds ratio induced by the data when comparing a model where both data sets are described by a single set of parameters and are thus compatible, to a disjoint model where each of the data sets is described by its own set of parameters. If $\Lambda_{\mathcal{M}}>1$, then the two data sets are deemed to be compatible (as their induced parameter space constraints overlap), and data set $\rm{d}_{2}$ is added to the list of compatible data sets. If $\Lambda_{\mathcal{M}}<1$ data set $\rm{d}_{2}$ is deemed incompatible with the data sets previously accepted and it is discarded. The details on the model we use and the form of the likelihood are given in Section 3.

The procedure begins by selecting an initial data set from galkin. We choose the measurements of Malhotra et al. [36] since by inspection it appears smooth with no outliers and covers the inner $\sim$ 8 kpc of the Milky Way RC. Then, we consider additional data sets one by one and use Eq. 3.4 to decide whether to include them. We end up with 12 compatible data sets out of 25, which we denote galkin12. The latter contains the following data sets: gas kinematics [36, 37, 38], star kinematics [39, 40, 41, 42] and masers [43, 44, 45, 46, 47]. However, it is worth noticing that the order in which the data sets are tested can be important. Indeed, if one starts from (or includes at the first steps) a data set that is incompatible with the rest (or with most of the others), then the final number of the compatible data sets is smaller111Notice that we did not perform this test for all the possible combinations of the 25 data sets since we deal with time-consuming computations. than 12. We have checked that among the 12 selected data sets the order of inclusion is not important, demonstrating the robustness of the procedure. We then bin the selected subsample in the manner described in Subsection 2.1. The resulting binned RC is shown in Fig. 1. We check to see whether the derived RC is compatible with the RC measurements we adopt for the outer Milky Way [23]. However, notice that the outermost data of [23] separately are not compatible with its inner part as well as with our binned subsample galkin12. Therefore, we perform an error bar re-scaling analysis that has the ability to correct for any systematic underestimation of the measurement uncertainties that might be present in our adopted data sets (see Section 4.2).

The model comparison procedure weeds out a sizable fraction of data from the galkin compilation, and the results of the analysis with the whole galkin compilation versus the “trimmed” one differ by a sizable factor of approximately 50 % in $\rho_{0}$. We therefore use the “trimmed” galkin12 compilation for our analysis even if the power of the Huang et al. [23] data -once included- drives the final results thus making the “full galkin + Huang” virtually indistinguishable from the “galkin12 + Huang”.

We also note that the Huang et al. [23] data have a dip in circular velocity between 10-20 kpc that is also noticeable in the galkin compilation. The origin of the dip is currently unclear. The authors of [23] state that this dip can be either an artifact due to high measurement uncertainties or it might reflect some actual feature of the kinematic tracers. Therefore, in order to check whether this dip in the data has any impact on the model fitting results, we perform a test where we exclude from our analysis the data between 10-20 kpc (see Section 4 for the results).

We have also carried out a consistency check by adopting different — and extreme — baryonic morphologies: for the central, heaviest, and lightest disks (namely those that provide the central, highest, and lowest RC) we have performed the procedure to select mutually compatible galkin data sets. The compatible data sets resulting from the procedure are always the same 12, and the inferred value of the local density $\rho_{0}$ remains constant, its precision staying virtually unchanged (the same does not apply to the anyway poorly constrained scale radius $r_{s}$ and density profile inner slope $\gamma$, see Section 4 and Appendix B for details), thus validating the procedure as insensitive to the various baryon morphologies.

2.1.2 The visible component (baryonic morphology)

In order to constrain the Milky Way’s dark matter content and its distribution, we need to include a model of the visible (baryonic) component of the Milky Way. Despite the improvements in the amount and quality of observations, controversies about the actual distribution of visible matter in the Milky Way remain. The uncertainties on the distribution of baryonic components do in turn impact the determination of the dark matter distribution, as addressed in detail by Iocco et. al. [16].

Following the approach of [16, 24, 48], we account for the stellar bulge, stellar disk(s), and gaseous disk using a fully observation-driven approach. The mass density of each one of the components is described by a three-dimensional function inferred from the observation of the stellar and gas target populations. Different functional shapes are available in the literature for both the disk(s) and the bulge component, with the latter presenting the highest degree of variation (it is to be noted that the differences are most important at small galactocentric distances). The diversity of observationally inferred morphologies of the baryonic component, available to date in the literature, must necessarily be accounted for as it impacts final conclusions on the matter distribution, as shown in previous work, e.g. [19]. We therefore collect a vast array of observationally inferred morphologies, separately for the bulge and the disk, and describe their treatment to our purposes in the following.

It is noteworthy that whilst all disks (thin, and thick) have necessarily cylindrical symmetry, the bulge morphologies show no degree of symmetry, presenting bar(s) at different angles (we devote to the bulges a fully consistent 3–dimensional treatment). The array of disk(/bulge) morphologies that we adopt here is taken from [19], where a summary of the different morphologies and of their normalization has been originally presented, and the interested reader may find appropriate citation to the original sources of observations for each of them. As of the use of such large array of possibilities for the visible component, rather than opting for choosing a single choice of bulge and disk, which would necessarily force us into a strong assumption on the baryonic components’ distributions, we follow the wake of previous analysis [19, 48, 20, 22], and adopt each possible combination of bulge and disk morphologies independently, thus constructing a catalogue of alternative “baryonic morphologies” out of all the possible individual combinations of the two (plus the gas disk). By treating each baryonic combination as individual and independent we are able to bracket the systematic uncertainties that might alter the inferred dark matter density distribution [16, 19]. Along these lines, we first adopt one single morphology to test our methodology and then test the generality of our conclusions on the methodology for the most extreme configurations allowed in the literature, namely those producing the “heaviest” and “lightest” rotation curves, see Section 4.2.

In what follows we introduce the reference baryonic morphology that we adopt throughout this study. Additionally, we present posterior results when adopting two other baryonic morphologies that are expected to bracket the possibilities for the disk and bulge. Finally, in Appendix B we introduce the details and we show the model fitting results for all baryonic morphologies presented in [16], thus bracketing the systematic uncertainty due to the observational ignorance on the morphology of the visible component of the Milky Way.

Stellar bulge. The bulge dominates the inner 2-3 kpc of the Milky Way and presents a triaxial shape with a bar extending at positive Galactic longitudes. This general picture is consistently painted by different observations, but the morphological details are rather uncertain. Following [49] we implement their E2 triaxial mass density model of the bulge that is based on COBE photometric maps. This model and the further calculations of the corresponding circular velocity have been presented in detail elsewhere [13]. We follow the same type of calculation as in [13] and all the parameters of the E2 fitting formula except the total normalization are fixed according to [50]. The normalization of the mass density of the stars in the bulge, or more concretely the microlensing optical depth $\left<\tau\right>$ is kept as a free parameter (but constrained with data from MACHO [51]). The microlensing optical depth towards the galactic bulge defines the total mass density of stars in the bulge (see for more details [13] on the relation of the microlensing optical depth and the bulge density normalization) and it is a measured quantity with an uncertainty of about 20%.

Stellar disk. It is well-known that the Milky Way stellar disk consists of two components: a thin and a thick disk that correspond to distinct stellar populations. In the present analysis we adopt a double exponential disk mass profile where, according to [52], the thickness of the thin (thick) disk and the exponential disk length scale are $\rm{z_{h}^{thin}}=0.27$ kpc $(\rm{z_{h}^{thick}}=0.44$ kpc) and $R_{d}=2.9$ kpc, respectively222For consistency, we provide the values that have been recalculated adopting $R_{0}$ = 8.34 kpc. Meanwhile, in our analysis, we keep free the local stellar surface mass density $\Sigma_{*}$.

Gas. The gas component is a sub-dominant component of the Milky Way. It takes the form of molecular, atomic and ionized hydrogen and its distribution along with its total mass density are relatively well-known. We follow the observationally inferred morphology for the density distribution of the gaseous disk taken from [53] for the inner 3 kpc and from [54] for the remaining range. We do not take into account the uncertainties on the density distribution profile of the gaseous disk. This is due to both its dynamically sub-dominant nature and its relatively well measured density distribution.

Finally, in Fig. 1 we plot circular velocities that correspond to the above-described baryonic components. The shaded regions in Fig. 1 illustrate the observational uncertainties on the normalization factors $\left<\tau\right>$ and $\Sigma_{*}$ (see Section 3 for details) in the inferred bulge and disk RCs.

2.2 The dark matter component

As a model of the dark matter distribution within the Galaxy we adopt the spherical generalized Navarro, Frenk, and White (gNFW) profile [55, 56]:

[TABLE]

where $r_{s}$ is the characteristic radius of the halo, $\rho_{0}$ is the dark matter density at Sun’s location $R_{0}=8.34$ kpc, and $\gamma$ is the inner slope of the density profile333The value $\gamma=1$ corresponds to the standard NFW profile.. The enclosed mass within a sphere of radius $R$ is

[TABLE]

In the assumed spherical dark matter distribution the angular velocity for a circular orbit is

[TABLE]

3 Statistical framework

The RC curve obtained in Sections 2.1–2.2 reflects the total (baryonic + dark matter) kinematics of the Milky Way. We now use this RC, in combination with information on the distribution of gas, stars and dark matter as described in Section 2, to perform the global mass modeling and to gather insights on the underlying dark matter distribution. To do so we fit a global model of the Galaxy that consists of four components: stellar disk, gaseous disk, stellar bulge and dark matter halo,

[TABLE]

where the first three baryonic components are presented in Section 2.1.2 and the dark matter component (see Section 2.2) depends on parameters $\Theta=(\rho_{0},r_{s},\gamma)$. Note that each term in Eq. 3.1 is implicitly a function of galactocentric radius $r$ and that angular velocity $\omega_{c}$ is used instead of linear velocity $V_{c}\equiv r\omega_{c}$.

3.1 Priors, likelihood, and posterior sampling

For a given morphology, our model has five free parameters: the scale radius of the dark matter halo $r_{s}$, the inner dark matter density slope $\gamma$, the local dark matter density $\rho_{0}$, the microlensing optical depth $\left<\tau\right>$, and the stellar surface density $\Sigma_{*}$. We explore the model parameter space using uniform priors in the following variables and ranges:

[TABLE]

We use fairly wide priors, which encompass the support of the likelihood. The last two parameters ($\Sigma_{*}$ and $\left<\tau\right>$), are nuisance parameters which are each independently constrained by a Gaussian likelihood. For the mean and standard deviation of this likelihood we adopt the values of the stellar surface density at the Sun’s position provided by [57], i.e. $\Sigma_{*}^{obs}=(38\,\pm\,4)\times 10^{6}\,\rm{M_{\odot}/kpc^{2}}$, as well as the measurement of the microlensing optical depth provided by the MACHO collaboration in Popowski et al. [51], namely $\langle\tau\rangle^{obs}=2.17^{+0.47}_{-0.38}\times 10^{-6}$. For simplicity, we symmetrize the error in the microlensing optical depth and adopt a standard deviation of $\sigma_{\langle\tau\rangle}=\pm 0.47$ (which is conservative, as it uses the larger of the two errors).

The likelihood function takes the form

[TABLE]

where $\Theta=(r_{s},\gamma,\rho_{0})$, $\bar{\omega}_{i}$ is the measured angular velocity and $\sigma_{\bar{\omega},i}$ is the corresponding uncertainty, here $i$ runs over the radial RC bins. The posterior from Bayes theorem is given by

[TABLE]

where $\mathcal{M}$ represents the assumed model and morphology (see Eq. 3.1), the likelihood $P(\rm{d}|\Theta,\Sigma_{*},\langle\tau\rangle,\mathcal{M})$ is given by Eq. 3.3, the prior $P(\Theta,\Sigma_{*},\langle\tau\rangle|\mathcal{M})$ is uniform in the variables and ranges indicated above, and the Bayesian evidence $P(\rm{d}|\mathcal{M})$ is an irrelevant normalization constant in this context.

The posterior is sampled using the open source affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler emcee [58]. In emcee we use 50 walkers and $2.5\times 10^{5}$ posterior samples (in the case of mock data we use $2\times 10^{5}$ posterior samples). We find it sufficient to discard the first 300 samples from each walker as burn-in. For all the runs, including different baryonic morphologies, the acceptance fraction varies between $0.28$ and $0.4$.

We further allow for the presence of residual undetected systematic effects between the different data sets (this is, between galkin12 and the two Huang et. al. [23] data sets) that we are using in the observed RC compilation by using the procedure described in, e.g., [59, 29]. The procedure has the ability to detect systematic offsets in the data as well as underestimated measurement uncertainties. Specifically, it allows to check whether the reported error bars capture the intrinsic noise in the data. For each data set $j$ and bin $i$, we rescale the variances of the binned data by a factor $1/\epsilon_{j}$, and replace the variance in the likelihood by $\sigma^{2}_{ij}/\epsilon_{j}$. The nuisance parameters $\epsilon_{j}$ describe possible error rescaling due to undetected systematics in each data set. If there are no residual systematics, the observed variance in the RC is fully explained by the statistical error in our binned data and the data sets are mutually compatible (as they should be, for in compiling them we used an appropriate compatibility test), then $\epsilon_{j}\approx 1$. We expect residual systematics to yield $\epsilon_{j}\ll 1$, which would inflate the binned data variance to account for sources of variance not captured by the statistical errors.

The first term of the likelihood function in equation (3.3) needs to be modified to include the $\epsilon$ reascaling parameters. With this, the likelihood takes the form

[TABLE]

where $i$ runs for the number of bins in data set $j$.

For the rescaling parameter $\epsilon_{j}$, following [29] we use Jeffreys’ priors which are uniform in $\log{\epsilon_{j}}$ and take the following form,

[TABLE]

which correspond to a prior on $\epsilon_{j}$ of the form,

[TABLE]

We further investigate whether our analysis is robust with respect to the $\sim$ 20 kpc dip apparent in our data compilation but whose physical origin is unclear. In order to do so, we perform an alternative analysis excluding from the data sets bins between 10 and 20 kpc. The results are presented in Section 4.

3.2 Performance of parameter reconstruction

We wish to investigate the sampling errors of our parameter reconstruction procedure (i.e. over repeated data realizations) in order to establish, for example, coverage properties of our Bayesian posterior regions under different procedures for creating them. As these properties in general depend on the true value of the parameters, we investigate 25 different choices for such fiducial parameter values, as given by the possible combinations in Table 1. For each set of fiducial parameters, we generate 100 mock data sets, drawn from the likelihood given by Eq. (3.3), with measurement error in each bin equal to the estimated standard deviation within that bin. The underlying angular velocity are obtained using Eq. 3.1, where $r_{s},\gamma,\rho_{0},\Sigma_{*}$, and $\left<\tau\right>$ are set to the fiducial values. For the baryonic components we adopt the HG disk [52], E2 bulge [49], and gas [53, 54] as our fiducial reference morphologies (see discussion in Section 2.1.2).

4 Results

4.1 Mock data reconstruction

Using the priors described in Section 3.1, we have fit 100 mock RCs for each fiducial point in the parameter space in order to test our reconstruction procedure. Figure 2 shows, as a first example, the full posterior distribution given mock data generated by two different choices of the fiducial parameters. In both cases, the local dark matter density $\rho_{0}$ is very well recovered while the inner slope of the density profile $\gamma$ and the scale radius $r_{s}$ are highly degenerate. These features are present in the posteriors for all 25 fiducial points.

We further estimate the accuracy of the reconstruction for the dark matter parameters of interest: $r_{s}$, $\rho_{0}$, and $\gamma$. In order to do so we define the ‘fractional standard error’ where in the definitions below $n_{i}$ is the number of mock data realizations for each fiducial point in the parameter space, $\theta_{\mathrm{true}}$ is the value of the fiducial parameters used and $\hat{\theta}$ is a point estimate. We consider four such point estimates: the mean, median, maximum a posteriori (MAP), and the maximum likelihood (ML). These estimates (except for the ML, which is obtained directly from the likelihood) are derived from the marginal 1D posterior distribution for the parameter in question.

The fractional standard error (FSE) is defined as the square root of the mean squared error normalized by the true value:

[TABLE]

The FSE is a measure of the (fractional) precision and accuracy of the reconstructed point estimate.

The detailed plot for the FSE with respect to different estimators for each of the 3 parameters is shown on Fig. A.9 of Appendix A where it is presented in form of color maps. We can summarize these results as follows:

• •

The values for the FSE of $\rho_{0}$ varies between 0.03 and 0.06 (see Fig. A.9) independently of the chosen point estimator. This indicates that the recovered value of the local dark matter density can be expected to have better than $6\%$ precision for any reasonable values of the dark matter halo parameters.

• •

The accuracy and precision of the $\rho_{0}$ estimators slightly improve for larger values of $\gamma$ (except for the ML).

• •

The accuracy and precision of the estimators for $\gamma$ and $r_{s}$ (except for MAP and ML estimators of $r_{s}$) is better for higher values of $\gamma$ and $r_{s}$.

• •

In general the estimates of $\gamma$ and $r_{s}$ have poor accuracy. This is due to the curved degeneracy between these two parameters as illustrated in Fig. 2. While a particular combination of $\gamma$ and $r_{s}$ may be constrained, marginalizing the posterior to 1D yields poor estimates of these quantities individually.

We now turn our attention to the coverage properties of the reconstructed intervals for each parameter. Coverage is the fraction of the time, over repeated data realizations, that the inferred interval contains the true value of the parameter. We consider different definitions for our 1D intervals: the highest posterior density (HPD) interval from the 1D marginal posterior (which by definition is the shortest credible interval containing a given fraction of the total probability); median credible regions (or equal-tailed interval) (CR), that are defined as the 16th and 84th percentile of the posterior; and the profile likelihood444The profile likelihood function is defined for each value of $\Theta$ by maximization of the likelihood function over the nuisance parameters. interval (PL), the region where the log-likelihood is within 1 of its maximum value (this is a classical way of constructing confidence intervals, based on Wilks’ theorem [60]).

Coverage plots from 100 mock data realization as a function of the assumed fiducial parameter values are presented in Fig. A.10 of Appendix A. These results can be summarized as follows:

• •

In the case of the local density $\rho_{0}$, all intervals overcover (i.e., they are conservative), with the most severe overcoverage occurring for the PL interval (bottom left panel of Fig. A.10).

• •

For the scale radius $r_{s}$ the situation is similar to the previous point, but the HPD and CR Bayesian intervals suffer from undercoverage in some parts of the parameter space, notably for small values of $\gamma$ and more severely so for the CR intervals.

• •

For $\gamma$ we observe both over- and under-coverage for the Bayesian intervals (both HPD and CR), while the PL interval always overcover.

The reasons for the behaviours described in the last two points can be traced back to the 2D degeneracy between $r_{s}$ and $\gamma$: once marginalized to 1D the resulting 1D Bayesian intervals can have poor coverage depending on where in the 2D degenerate region the true value is located. The over-coverage observed for all the fiducial parameter values for the PL is a consequence of the PL intervals being typically larger than both the Bayesian ones. For this reason, for $\gamma$ and $r_{s}$ the 1D coverage of Bayesian intervals can stray from being exact, while being much better behaved in the 2D plane spanning both. We leave further investigation of this effect for future work. Notice that the two nuisance parameters are always well recovered within the observed $1\sigma$ uncertainties, independently of where we are in the parameter space, as they are constrained by an independent piece of the likelihood.

4.2 Reconstruction of parameters using real data

In this section we present the inference results when using the data sets for the Milky Way presented in Section 2.1. We have argued above that our Bayesian procedure for selecting compatible data sets robustly eliminates data sets that are likely to be systematically offset. An inspection of the galkin12 compilation selected in this way shows that the data within each bin are nicely Gaussian distributed around their mean value. This is another indication of compatibility, and it validates our usage of a Gaussian likelihood function. Here, we further investigate the possibility of remaining undetected systematic effects within our compilation of RC data by applying the Bayesian formalism explained in Section 3, where we introduced the variance scaling parameters $\epsilon_{j}$ (see Eq. 3.5).

We look for systematic differences between the galkin12 and the two data sets of Huang et al. [23] by introducing three epsilon scaling parameters, one for each data set. Notice that if there were any systematic errors within the galkin12 data, the corresponding $\epsilon_{j}$ parameter would prefer values smaller than unity to compensate for them. The 1D marginal posterior distributions for the $\epsilon_{j}$ parameters are shown in Fig. 3, where we see that the posteriors peak at $\log(\epsilon_{j})=0$ for all three parameters, thus indicating that there is no need to further increase their uncertainties. In the following, we thus remove the $\epsilon_{j}$ parameters from the analysis (having checked that the ensuing posterior distributions for the parameters of interest do not exhibit any significant shift when including or excluding the $\epsilon_{j}$ parameters). In Fig. 4 we also show the comparison of the MCMC analysis with and without including error rescaling parameters. This plot indicates practically no difference in the resulting posteriors.

We now check whether the $\sim 20$ kpc dip in the data (see Fig. 1) affects the parameter reconstruction. In order to do so, we re-run the MCMC analysis (without including error rescaling parameters) for both the whole binned RC data and omitting the bins between 10 and 20 kpc. The results of these two runs are compared on Fig. 5, where we see essentially no difference in the resulting posteriors (the ML point does however suffer a shift, negligible in the case of $\rho_{0}$). This exercise, together with the fact mentioned in Section 2.1.1 that binned whole galkin data compilation once combined with the outer RC data show almost no change in the resulting Bayesian inference with that of galkin12 + extended Huang et al. RC, indicates the importance of the RC data beyond the $\sim$ 20 kpc.

Finally in Fig. 6 we superimpose RC based on the maximum likelihood parameter values over the observed RC data, and indicate the 68% and 95% HPD regions from our analysis. The latter is derived from the corresponding pdf of the circular velocity which, in turn, is calculated based on the pdfs of the emcee samples.

The quality of the best-fit RC when adopting different morphologies is presented in Appendix B, in terms of its reduced chi-squared (Table 3), showing a good quality of fit for all morphologies. Instead, here we summarize details of how the Bayesian posterior inferences change when adopting three different and critical disk morphologies. In particular, we make use of two additional disk morphologies (demonstrated in Fig. 7) that are obtained with the lightest disk (also referred to as CM [61]) and heaviest disk (also referred to as BR [57]) so as to present a “bracketing” of the actual systematic uncertainty due to ignorance of the actual morphology of the visible component. Notice that we keep the bulge morphology fixed and the same as in our reference morphology555If we vary bulge morpholgies while keeping a disk morphology fixed, the resulting dark matter parameters remain almost unchanged with mild variations well within the corresponding 68% HPD regions, see Appendix B..

The results of the MCMC analysis for the three morphologies are presented on Fig. 8 and summarized in Table 2. In this table we quote the reduced $\chi^{2}$ that is calculated for 49 degrees of freedom (54 bins - 5 free parameters). From this figure and the table it is evident, that whilst posteriors on $\gamma$ and $r_{s}$ do change with different morphologies, the local density is unchanged with a maximum likelihood value of $\sim 0.44\,\rm{GeV/cm^{3}}$. We also find that it is impossible, at the moment, to rule out any of them using the goodness of the fit. Finally, we emphasize that in our analysis we fix the Sun’s galactocentric distance, whose exact value might impact the determination of the dark matter distribution. For example, reference [22] shows that there is a correlation between the local dark matter density and the value of the Sun’s galactocentric distance. However, [22] only makes use of the inner Milky Way RC data (differently from us), and adopts a different statistical approach than ours. Furthermore, [22] does not assess the compatibility of the galkin sub-sets of data, as we do here, and use instead the original galkin compilation. Therefore, the impact of varying the Sun’s galactocentric distance on the determination of the local dark matter density in our study cannot be directly compared with what obtained in [22]. Despite the differences in the methodology, we notice that our value of the local dark matter density for the Sun’s galactocentric distance of 8.34 kpc is compatible within 95 % HPD region with the result of [22].

Although our derived value of the local dark matter density is in good agreement with recent measurements (see e.g. [62, 63] and references therein)666Notice, that when we apply our analysis with our reference baryonic morphology to the recently published Gaia DR2 data [64], the local dark matter estimate is compatible with our estimate (when using galkin12 data assuming the local standard of rest and the values of $\rm V_{0},\,R_{0}$ as in [64]) and with the estimate from [65] when the authors use their B2 baryonic model. Ref. [65] appeared in arXiv while our paper was submitted to this Journal., it is in tension with the value derived in [23]. The reason for that is the difference in the modelling of the total Milky Way gravitational potential between our work and [23]. In particular, in the latter the authors add two additional components (with respect to ours) in order to model a hypothetical caustic ring of dark matter in the Galactic plane. Doing so the derived value of the local dark matter density (not accounting for that in the ring) decreases by $\sim\leavevmode\nobreak\ 25$ %. In order to derive this number we have included the dark matter ring in our model in the same manner as in [23], we then fixed all the related parameters to their best-fit values and we obtained the value of the local dark matter density consistent within 68 % HPD region with that of [23].

5 Conclusions

We have developed a novel Bayesian methodology to reliably and precisely infer the distribution of dark matter within the Milky Way using the rotation curves method. Based on the most recent data sets for the rotation curve, we first selected a subset of tracers that are mutually compatible with each other, thus excluding the ones that may be suffering from systematic bias. We then demonstrated the statistical performance of our rotation curve method (over many data realizations) by applying it to simulated data, which we generated respecting the statistical properties of the observed rotation curve. We generalized our procedure against (systematic) uncertainties on the visible component of the Milky Way, by applying it to different possible baryonic morphologies, in order to bracket the current uncertainties on the shape and normalization of the Galactic stellar disk(s) and bulge, as well as of the (subleading) interstellar gas component (see Section 4.2 and Appendix B).

Our most relevant findings can be summarized as follows:

• •

Based on the kinematic data and the Galactic parameters used in this work, we infer a local dark matter density is $\rho_{0}=0.43\pm{0.02}\pm{0.01}\,\rm{GeV/cm^{3}}$, where the first errors are statistical and the second systematic (adopting different baryonic morphologies). This number is in fair agreement with previous determinations in the literature (see e.g. [12, 13, 19]).

• •

The determination of $\rho_{0}$ is very precise (68 % credible region always between 4% and 7%), and extremely accurate: when tested against mock data, the value of $\rho_{0}$ inferred through the procedure does not fall more than $\sim$6% away from the true value. This is of course predicated on our model for the RC data to be a complete and accurate representation of how the true data are generated, in that the reconstruction accuracy is only valid for mock data generated from the model itself.

• •

In the gNFW formulation, the scale radius $r_{s}$ and the slope of the inner density profile $\gamma$ are degenerate, thus making the separate reconstruction of the two parameters more challenging. Whereas the combination of the two is relatively well constrained in a 2-dimensional space, the individual posteriors, once marginalized to 1D, are not very informative.

• •

The last two statements are robust against uncertainties in the baryonic morphologies (i.e. they are unaffected by the systematic ignorance of the baryonic morphology), and the shape of the underlying dark matter profile (precision and accuracy are little sensitive of the parameter space region where the underlying true value lies).

From the above, we conclude that the inference of the local dark matter density $\rho_{0}$ is extremely robust: it is important to stress that the central value obtained with our method is both very precise (with the 68% credible interval being $\sim$ 4%) and highly accurate (never falling farther than $\sim$ 6% from the true value in simulated data). We also found that the inferred value and error bars are stable with respect to the morphology of the visible, baryonic component (with $<3\%$ shift in the inferred $\rho_{0}$). Our results are also stable with respect to various tests for the presence of systematic errors in the RC data used.

The methodology presented here is very flexible, as it allows the effortless incorporation of further data sets and accompanying nuisance parameters. We therefore expect that it will prove very useful in analyzing current and upcoming data about the structure and properties of the Milky Way mass distribution.

Acknowledgements

We thank the referee for his/her comments that helped improve the paper. EK work was supported by the São Paulo Research Foundation (FAPESP) under Grant No. 2016/26288-9 and by the grant AstroCeNT: Particle Astrophysics Science and Technology Centre is carried out within the International Research Agendas programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund. F. I. acknowledges support from the Simons Foundation and FAPESP process 2014/11070-2. AG-S and RT are supported by Grant ST/N000838/1 from the Science and Technology Facilities Council (UK). RT was partially supported by a Marie-Sklodowska-Curie RISE (H2020-MSCA-RISE-2015-691164) Grant provided by the European Commission. This work has been possible through the FAPESP/Imperial College London exchange grant “SPRINT”, process number 2016/50006-3. This research was supported by resources supplied by the Center for Scientific Computing (NCC/GridUNESP) of the São Paulo State University (UNESP).

Appendix A Mock data reconstruction: accuracy and precision

On Fig. A.9 we demonstrate the fractional standard error from 100 mock data realization as a function of the assumed fiducial parameter values, where the four columns correspond to the four different definitions of point estimate, and rows correspond to three parameters of interest. The color bars next to the plots decode the corresponding values of FSE. The results of this figure are summarized in Section 4.1.

On Fig. A.10 we demonstrate the coverage plots from 100 mock data realization as function of the assumed fiducial parameter values, where the three columns correspond to the three interval definitions, and rows correspond to three different parameters of interest. In these plots, exact coverage would appear as a yellow square, with overcoverage shown in shades of green and undercoverage in orange/red. The color bars next to the plots indicate the fraction of the time, over 100 data realizations, that the inferred interval of interest contains the true value of the parameter. The results of this figure are also summarized in Section 4.1.

Appendix B Dark matter distribution adopting different baryonic morphologies

Here we present the results of the MCMC analysis for all combinations of 5 disk morphologies and 6 bulge morphologies.

In the first column of Table 3 we list the bulge/disk morphology combinations and references on it in the literature. The details on these morphologies can be found in [16, 19]. The main conclusions here are very same as summarized in Section 4 for three different disk morppholgoes, but generalized to all the possible disk and bulge morphologies and its combinations:

• •

the value of the local density $\rho_{0}$ is independent of the assumed baryonic morphology

• •

the values of the reduced $\chi^{2}$ for all morphologies do not rule out any of them.