Handling the Uncertainties in the Galactic Dark Matter Distribution for Particle Dark Matter Searches

TL;DR

This paper characterizes the Milky Way's dark matter distribution and its uncertainties, providing a comprehensive likelihood grid to improve indirect and direct dark matter detection analyses.

Contribution

It offers a detailed analysis of Galactic dark matter parameters with marginalized uncertainties and releases a likelihood grid for use in dark matter search studies.

Findings

Local DM density rho0 constrained to 0.3-0.8 GeV/cm3

Standard rho0=0.3 GeV/cm3 and R0=8.5 kpc are excluded at >4 sigma

Inner slope gamma weakly constrained, no preference for cored or cuspy profiles

Abstract

In this work we characterize the distribution of Dark Matter (DM) in the Milky Way (MW), and its uncertainties, adopting the well known "Rotation Curve" method. We perform a full marginalization over the uncertainties of the Galactic Parameters and over the lack of knowledge on the morphology of the baryonic components of the Galaxy. The local DM density rho0 is constrained to the range 0.3 - 0.8 GeV/cm3 at the 2 sigma level, and has a strong positive correlation to R0, the local distance from the Galactic Center (GC). The not well-known value of R0 is thus, at the moment, a major limitation in determining rho0. Similarly, we find that the inner slope of the DM profile, gamma, is very weakly constrained, showing no preference for a cored profile (gamma~0) or a cuspy one (gamma~[1.0,1.4]). Some combination of parameters can be, however, strongly constrained. For example the often used standard rho0=0.3 GeV/cm3, R0=8.5 kpc is excluded at more than 4 sigma. We release the full likelihood of our analysis in a tabular form over a multidimensional grid in the parameters characterizing the DM distribution, namely the scale radius Rs, the scale density rhos, the inner slope of the profile gamma, and R0. The likelihood can be used to include the effect of the DM distribution uncertainty on the results of searches for an indirect DM signal in gamma-rays or neutrinos, from the GC, or the Halo region surrounding it. As one example, we study the case of the GC excess in gamma rays. Further applications of our tabulated uncertainties in the DM distribution involve local DM searches, like direct detection and anti-matter observations, or global fits combining local and GC searches.

Equations

$${\Omega }_{g,\odot }=\frac{V_0+V_{\odot }}{R_0}.$$

$$\rho (r)={\rho }_s{\left(\frac{r}{R_s}\right)}^{-\gamma }{\left(1+\frac{r}{R_s}\right)}^{-3+\gamma }$$

$$\bar{\omega }=\frac{{\sum }_{i=1}^{N_{bin}}{\omega }_i/{\sigma }_{{\omega }_i}^2}{{\sum }_{i=1}^{N_{bin}}1/{\sigma }_{{\omega }_i}^2},$$

$${\sigma }_{\bar{\omega }}^2=\frac{{\sum }_{i=1}^{N_{bin}}(\bar{\omega }-{\omega }_i{)}^2/{\sigma }_{{\omega }_i}^2}{{\sum }_{i=1}^{N_{bin}}1/{\sigma }_{{\omega }_i}^2}+\frac{N_{bin}}{{\sum }_{i=1}^{N_{bin}}1/{\sigma }_{{\omega }_i}^2},$$

$${\chi }_{\mathrm{RC}}^2(R_s,{\rho }_s,\gamma ,R_0,{\mathcal{M}}_i)=\sum _j\frac{{\left({\bar{\omega }}_j-{\bar{\omega }}_j^{th}\right)}^2}{{\sigma }_{{\bar{\omega }}_j}^2}+{\chi }_{⟨\tau ⟩}^2(⟨\tau ⟩)+{\chi }_{{\Sigma }_{\ast }}^2({\Sigma }_{\ast })$$

$$\mathcal{J}={\int }_{\Delta \Omega }d\Omega {\int }_{l.o.s.}ds{\rho }^2(r(s,\psi ))$$

$$\frac{d\Phi }{dE}=\frac{⟨\sigma v⟩}{8\pi {m_{\chi }}^2}\frac{dN}{dE}\mathcal{J}$$

$${\chi }_{\mathrm{GCE}}^2=\sum _{i,j}(d_i-t_i){\left({\Sigma }_{ij}+{\delta }_{ij}({\sigma }_{\mathrm{rel}}{t}_i{)}^2\right)}^{-1}(d_j-t_j),$$

$${\chi }_{\mathrm{total}}^2={\chi }_{\mathrm{GCE}}^2(⟨\sigma v⟩,m_{\mathrm{DM}},\mathcal{J})+{\chi }_{\mathrm{RC}}^2(\mathcal{J},\gamma )+{\chi }_{\gamma }^2,$$

$${\chi }_{\mathrm{total}}^2={\chi }_{\mathrm{GCE}}^2(⟨\sigma v⟩,m_{\mathrm{DM}},\mathcal{J})+{\chi }_{\mathrm{RC}}^2(R_0,R_s,{\rho }_s,\gamma )+{\chi }_{\gamma }^2+{\chi }_{R_0,R_s,{\rho }_s}^2,$$

$$\rho (r)={\rho }_s{\left(\frac{r}{R_c}\right)}^{-1}{\left(1+\frac{r^2}{R_c^2}\right)}^{-1},$$

$$\rho (r)={\rho }_s\exp \left(-\frac{2}{\alpha }\left[{\left(\frac{r}{R_s}\right)}^{\alpha }-1\right]\right),$$

Full text

TTK-18-37

LAPTH-051/18

1 Introduction

The distribution of Dark Matter (DM) within the Galactic Halo is one crucial ingredient in direct and indirect particle DM searches. Typically, the DM Halo is assumed to be spherically symmetric and various analytical parameterizations for the radial profile are adopted, e.g. Navarro-Frenk-White (NFW [1]), Burkert [2], Einasto [3] or generalized NFW. The uncertainty arising from our lack of knowledge over the DM Halo shape is thus often taken into account by the bracketing of benchmark models: to account for the uncertainty in the DM distribution in the inner Galaxy a cored profile, like Burkert, and a cuspy one, like NFW, are often adopted as benchmark, and the results shown for the extreme cases. Still, further uncertainties are present, like the overall normalization of the profile, often expressed as the DM density $\rho(R)$ at $R_{0}$ –the Solar system distance from the Galactic Center– and referred to as the local DM density $\rho_{0}$, or the overall scale radius of the profile $R_{s}$. The quantity $R_{0}$ itself is also affected by observational uncertainties, which propagate into the global profile uncertainty.

Taking into account in a self-consistent way all the above uncertainties is not straightforward, also because strong correlations among the various parameters are present. Several analyses constraining the DM Halo distribution and the related parameters have been performed to date, e.g. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], but the correlation matrix is not typically provided, thus making the full results not promptly usable from the reader.

The main goal of this work is to fill this gap and to provide an account as comprehensive as possible of the uncertainties of both observational and modelization nature. We use the well-known method of Galactic “Rotation Curve” to derive the constraints on the parameters of the Galactic DM Halo taking the above–mentioned uncertainties and dependencies into account and we provide the full likelihood, so that the interested readers can use these constraints and apply them to their own analysis. Possible applications include analyses of a DM signal, or upper limits, from the Galactic Center, or its surrounding. Typical cases are searches for a gamma-ray signal, in the form of lines [19, 20, 21, 22], or smooth spectra either at GeVs (see e.g. [23, 24]), or at TeV with Cherenkov telescopes [25, 26, 27, 28]. Analogous searches for DM signal from the GC can be performed in neutrinos [29, 30, 31, 32]. In all these cases the amount of DM signal in the GC region is proportional to the so-called $\mathcal{J}$-factor, which is a function of the properties of the DM Halo. The uncertainty in the $\mathcal{J}$-factor can then be taken into account through the likelihood provided in this work111The likelihood profile from the full data-driven analysis presented in this work is available at https://github.com/mariabenitocst/UncertaintiesDMinTheMW.. Sometimes extra-constraints are also needed, for example one would like to restrict the slope of the DM Halo near the GC to a certain range of values, or apply a specific prior on the local DM density. The advantage of providing the full likelihood is that extra-priors can be easily included just introducing them as extra multiplication factors in the likelihood, so that to use a new likelihood which includes all the desired information. Similar considerations can be applied to local searches, for example DM searches with antiprotons, e.g., [33, 34, 35, 36, 37, 38] or antimatter in general (antinuclei, positrons) [39], whose result is mostly sensitive to the local DM density but also to the shape of the DM profile within few kpc of the Solar System.

Finally, it is worth noticing here that direct searches are sensitive to the local DM density, and to the local DM velocity distribution [40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51]. Our results can thus be also used to investigate the astrophysical uncertainties in direct DM searches, although in this study we will only touch upon the subject.

2 Setup and Data

To derive constraints on the distribution of DM in the Galaxy we use the method commonly referred to as the “Rotation Curve”. This method, largely adopted in the literature to determine the DM existence and distribution in rotation supported disc galaxies and in particular in our own e.g. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], relies on the use of two main elements: on one hand the actual, observed gravitational potential, which is inferred through the motion of appropriate tracers of the circular velocity of the disk (the Rotation Curve); on the other hand the potential expected from the baryonic component (gas and stars). The mismatch between the two (if any) is fit with the expectation (in Newtonian approximation) from a spherical DM profile, typically some form of broken power-law. This is opposed to the “local” method, which, instead relies on local observables only, typically the position and motions of stars within few kpc from the Solar system to infer the local gravitational potential and DM density, e.g. [52, 53, 54, 55, 56, 57, 58, 59]. For a review and discussion of both methods see [60].

2.1 The observed rotation curve

As tracer of the total gravitational potential, or Rotation Curve (RC), we adopt the data from the galkin compilation [6, 61]. The compilation contains data up to Galactocentric radii of $R\sim$ 20 kpc and includes the kinematics of gas, stars and masers for a total of 2780 measurements collected from the literature. More details are reported in the original publications [6, 61], with extensive studies about the possible source of systematics described in the Supplementary Material of [6]. galkin uses as input the local standard of rest $(U,V,W)_{\odot}$ (see below), $R_{0}$, the Solar system distance from the Galactic Center, and $V_{0}$, the local circular velocity to self-consistently transform determinations from different observations into a point on the RC $(R_{i}\pm\delta R_{i},\omega_{i}\pm\delta\omega_{i})$ with respective errors. The data-points on the RC can be expressed in term of circular velocity $V$, or, circular angular velocity $\omega$, which is typically more convenient since in the latter case $R$ and $\omega$ errors are uncorrelated, contrary to the former. Since $R_{0}$ is a parameter which we vary in our analysis, the RC curve is self-consistently updated when $R_{0}$ is changed.

To trace the RC from 20 kpc up to $\sim$100 kpc, stellar dispersion data are typically used. Nonetheless, since these stellar tracers are not in circular orbit, the use of these data requires further assumptions, for example about the virialization of the system and the velocity anisotropy. We prefer to conservatively limit ourselves to the use of tracers in circular motion such as those contained in the recent compilation galkin –though they imply an intrinsic limitation to the innermost Galaxy – and postpone an accurate study with stellar tracers to future work.

2.2 Baryonic morphology

The visible (baryonic) component of the Milky Way, is typically separated between a stellar bulge (highly asymmetric, dominating the potential from the center up to 3-4 kpc), a stellar disk (with possibly more than one component), extending up to $\sim$15 kpc, and a disk of gas, lying approximately in the same plane of the stellar disk, and mostly subleading in dynamical terms, which we include nonetheless for the sake of completeness.

In order to study the bulge, rather than relying on the spherical approximation common to many previous analysis, we adopt the approach described in [7], which takes into account a full three-dimensional density distribution for stars in the Galactic bulge, solves the potential, and then finds the component within the disk, also allowing a precise estimate of the lack of axisymmetry in that region.

The different morphologies for the Bulge and the Disc (separately) –collected and presented in [6] and then adopted by the same authors in [7]– are inferred from observations of different population of stars in different regions (see original references in [6, 7]), and therefore fully empirical, three-dimensional, alternative descriptions of the stellar component of Bulge and Disk(s). Following [6, 7], we adopt –separately– 6 models of Bulge (labeled a,b,c,d,e,f) and 5 models of Disc (labeled I,J,K,L,M), which are then individually combined (one disk and one bulge at the time) thus obtaining a total of 30 combinations of Bulge plus Disc. To each of these possible stellar morphologies, we add an observationally inferred morphology for the interstellar gas disk taken from [62] from the inner 3 kpc and [63] above 3 kpc, instead of bracketing two possible alternatives, as done in [6, 7] given the subdominant contribution of the gas component to the RC.

For our fit we will thus have one (discrete) parameter to describe the uncertainty related to the baryonic mass, namely the index of the baryonic morphology $\mathcal{M}_{i}$. The normalization of each morphology, corresponding the mass of the Disc and mass of the Bulge however, also has its own uncertainty. To take into account this uncertainty we normalize the morphology so that to agree with microlensing optical depth measurements towards $(\ell,b)=(1.50^{\circ},-2.68^{\circ})$, $\langle\tau\rangle=2.17^{+0.47}_{-0.38}\times 10^{-6}$ [64], and local total stellar surface density $\Sigma_{*}=38\pm 4\,{\rm M}_{\odot}/{\rm pc}^{2}$ [65]. See again [7] for more details. To take into account the uncertainty in $\langle\tau\rangle$ and $\Sigma_{*}$ we will add them to the total $\chi^{2}$ used to constrain the DM Halo and vary them in the range $\pm 2\sigma$. This is discussed in more details in section 3.2.

Finally, the morphologies depend on $R_{0}$. Thus, when changing the used value of $R_{0}$ we self-consistently recalculate the morphology and its contribution to the RC.

2.3 Local Standard of Rest

Further uncertainty comes from the not well known peculiar motion of the Solar system with respect to the local standard of rest (LSR), the system comoving along a circular orbit around the GC with a velocity equal to the local RC velocity $V_{0}$. Recent measurements find values $(U,V,W)_{\odot}=(11.10,12.24,7.25)$ km/s [66], where $(U,V,W)_{\odot}$ are, respectively, the velocity orthogonal to the circle of the orbit and pointing outward the GC, the velocity tangential to the circle, and orthogonal pointing in the $z$ direction. The most relevant for the RC analysis is $V_{\odot}$ which in [66] is found to be $V_{\odot}=12.24$ km/s, but which has a quite larger scatter in the range $5-24$ km/s from different analyses in the literature [66, 18, 67]. $V_{\odot}$ can be used together with the precise determination of the total Solar system angular velocity $\Omega_{g,\odot}=30.24\pm 0.12\;{\rm km\,s^{-1}kpc^{-1}}$ [68] based on observations of the peculiar motion of the GC source Sagittarius A∗. They are linked by

[TABLE]

from which the local circular velocity $V_{0}$ can be derived once $R_{0}$ is also specified. For $R_{0}=8.0$ kpc, Eq. 2.1 gives $V_{0}=230$ km/s, which are commonly adopted values. In the following we will use as free parameter $R_{0}$ which we will vary in the range $[7.5,8.5]$ kpc. We will instead fix $V_{\odot}=12.24$ km/s , since the uncertainty in $V_{\odot}$ introduces a variation in $V_{0}$ similar or smaller than the one caused by $R_{0}$. In practice, the uncertainty in $V_{\odot}$ can be taken into account by considering a slightly more conservative range of variation for $R_{0}$. Recently, the GRAVITY collaboration [69] reported the very precise result $R_{0}=8.122\pm 0.031$kpc. If confirmed, this would essentially fix the value of $R_{0}$, so that in this case it would be convenient to consider explicitly $V_{\odot}$ as a parameter to vary in the analysis.

2.4 Dark Matter distribution

We parameterize the DM distribution as a spherically symmetric generalized NFW profile [1]

[TABLE]

where $r$ is the spherical distance from the Galactic center (GC), $R_{s}$ the scale radius of the profile and $\rho_{s}$ the scale density. The density behaves like $r^{-\gamma}$ toward the GC, and the case $\gamma=1$ denoted the standard NFW profile. The fit will thus have 3 parameters related to DM, $R_{s}$, $\rho_{s}$, and $\gamma$.

3 Methodology

For our analysis, we adopt the angular velocity rotation curve $\omega(R)$ instead of the linear velocity $V(R)$ in order to get rid of existing correlations between the uncertainty of the latter and that of the galactocentric distance $R$ [61].

3.1 Binning scheme

The binning of the data of the observed RC is a key point, as we adopt binned data instead of an unbinned analysis such as that performed in e.g.[6, 10]. An unbinned analysis would exploit the full constraining power of the data assuming that the underlying systematics uncertainties are under control. Systematics of diverse nature are possible. For example, single data points have sometimes very small errors, which, however, might not necessarily reflect the true uncertainty of the RC which can have systematic contributions, e.g. from peculiar motion of the gas tracers thus not following exactly a circular orbit, and more in general from deviation from axial symmetry of the motion of the tracers (see e.g. the Supplementary information of [6] and references therein). Whereas the systematic errors have been shown not to affect more general conclusions [6], they are important in the details of the determination of the DM profile [7]. Here we will thus employ a conservative view and use binned data.

We will hence consider binned data and as error the dispersion of the data in the bin rather than a formal weighted mean of the data. Furthermore, we bin the data in $x=R/R_{0}$ rather than $R$ itself. This because when varying $R_{0}$ the data ‘move’ along the $R$ axis. Binning in $x$ mitigates this problem, so that for different $R_{0}$ a given $x$ bin contains roughly the same unbinned data-points. We start from $x=2.5/8$. Data below this value of $x$ are not considered in the fit in order to avoid the inner Galaxy region for which there are significant deviations from an axysimmetric motion of the tracers of the gravitational potential. More in detail, we apply the following binning scheme:

• •

15 bins from $x=2.5/8$ to $x=10/8$ with a step of $\Delta x=0.5/8$,

• •

7 bins from $x=10/8$ to $x=18/8$ with a step of $\Delta x=1/8$,

• •

2 bins from $x=18/8$ to $x=22/8$ with a step of $\Delta x=2/8$,

for a total of 25 $x$ bins. A bin of $\Delta x=0.5/8$ or smaller for $x<10/8$ is necessary to properly follow the $\omega(x)$ curve which is quite steep in this $x$ range. At large $x$ the binning is larger also due to the scarcity of data. In order to assign a given data-point to a bin we consider only the central value $x_{i}$ and neglect the error $\delta x_{i}$. Within each bin the center of the binned data-point $\bar{\omega}$ and its uncertainty $\sigma_{\bar{\omega}}$ are constructed as follow

[TABLE]

[TABLE]

i.e., $\bar{\omega}$ is just the weighted mean of the data in the bin. $N_{bin}$ is the number of data points in the bin. The uncertainty $\sigma_{\bar{\omega}}$ is composed of two terms, the first is the weighted dispersion of the data, the second term gives the mean weighted error of the data, so that the final error of the binned data-point will be larger than the latter. The second term is however subdominant, i.e., the scatter in the data is much larger than the average error, typically by one order of magnitude of more, so the role of the latter is marginal, except for the last 2-3 bins, for which the two errors are comparable. Finally, before applying the above procedure we take unbinned data-points above $x=10/8$ with a relative error of less than 10% and we increase artificially their uncertainty to 10%. This is to avoid that when few data are present in a single bin the final result gets dominated by a single data-point with very small error. In practice, however, this affects only a single binned data-point, namely the 24th. An example of the binned and unbinned RC for $R_{0}=8$ kpc, $V_{0}=230$ km/s, and $(U,V,W)_{\odot}$ from [66] is given in Fig. 1.

In the Appendix we will show the effect of using a different binning scheme to study the impact on the final result.

3.2 Fitting procedure

We include in the fit 7 parameters, i.e., $R_{s}$, $\rho_{s}$, $\gamma$, $R_{0}$, $\mathcal{M}_{i}$, $\langle\tau\rangle$ and $\Sigma_{*}$. The number of parameters is still sufficiently small to use a grid scan rather than a Monte Carlo scan. In the following we will thus just use a discrete grid. More precisely, we use 50 values for $\rho_{s}$ linearly spaced in the range $[0.0,2.0]$ GeV/cm3, 50 values for $R_{s}$ logarithmically spaced in the range $[5.0,100.0]$ kpc, 15 values of $\gamma$ linearly spaced in the range $[0.0,1.5]$, 11 values of $R_{0}$ linearly spaced in the range $[7.5,8.5]$, and 30 morphologies $\mathcal{M}_{i}$. For $\langle\tau\rangle$ and $\Sigma_{*}$ and we use 10 values each, linearly spaced in the range $[-2\sigma,+2\sigma]$.

Having specified the above methodology to bin the data, we then compare them with the model using a simple $\chi^{2}$ statistics,

[TABLE]

and we evaluate $\chi^{2}_{\rm RC}$ over the grid defined above. $\bar{\omega}_{j}^{th}$ is the model prediction depending on $R_{s},\rho_{s},\gamma,R_{0},\mathcal{M}_{i}$, $\langle\tau\rangle$ and $\Sigma_{*}$ and is given by $\bar{\omega}_{j}^{th}=\sqrt{\left(\bar{\omega}_{j}^{b}(R_{0},\mathcal{M}_{i},\langle\tau\rangle,\Sigma_{*})\right)^{2}+\left(\bar{\omega}_{j}^{DM}(R_{s},\rho_{s},\gamma)\right)^{2}}$, i.e., by the sum of the DM and baryonic contribution. For $\chi^{2}_{\langle\tau\rangle}$ and $\chi^{2}_{\Sigma_{*}}$ we use the expressions $\chi^{2}_{\langle\tau\rangle}=(\langle\tau\rangle-2.17)^{2}/0.42^{2}$ and $\chi^{2}_{\Sigma_{*}}=(\Sigma_{*}-38)^{2}/4^{2}$ from section 2.2, where the error on $\tau$ has been made symmetric for simplicity. We verified that including $\langle\tau\rangle$ and $\Sigma_{*}$ in the $\chi^{2}$ does not crucially affect the analysis. Keeping $\langle\tau\rangle$ and $\Sigma_{*}$ fixed to their central values only reduces slightly the error in the determination of the other parameters of the analysis. This is likely due to the fact that the bulk of the uncertainty from the baryonic morphology is already taken into account considering the 30 different models $\mathcal{M}_{i}$. Nonetheless, for consistency of the analysis and for a more robust error determination, we include $\langle\tau\rangle$ and $\Sigma_{*}$ in the overall $\chi^{2}$. Furthermore, again for the above reason, just 10 grid values of $\langle\tau\rangle$ and $\Sigma_{*}$ are already enough to properly include the effect of their uncertainty on the analysis.

Another point to mention is that, formally the above definition is not fully self-consistent since the data change when $R_{0}$ is changed and we explore changes in $R_{0}$ of the order of 10%, which introduces changes in the data of the same order. This issue is unavoidable as soon as binned data are used. Nonetheless, as long as the induced variations in the data are smooth as function of $R_{0}$, as it is the case, the data variation can be thought as being reabsorbed into a redefinition of the model, so that the use of Eq. 3.3 should be approximately valid. Another minor inconsistency is given by the fact that once the set of 7 parameters is specified, the full rotation curve $V(R)$ is also specified and so is $V_{0}\equiv V(R_{0})$. The relation $V_{0}=V(R_{0})=R_{0}\sqrt{\omega_{b}^{2}(R_{0})+\omega_{DM}^{2}(R_{0})}$ should thus be enforced and used to remove one parameter. This, in practice, is not a big issue, since the fit will automatically prefer the region where this relation is satisfied. Furthermore, this extra freedom is, in practice, equivalent to not strictly assume Eq. 2.1 linking $V_{0}$ and $R_{0}$ but leaving some freedom in their relation to be constrained by the fit, which is a conservative choice.

We use a fully frequentist framework to derive constraints in sub-spaces of the 7 full dimensional space. Specifically we employ the commonly used method of profiling [70]. For example, when we build two-dimensional $\chi^{2}$ in two given parameters, for each 2-d grid point, we take the minimum $\chi^{2}$ over the remaining 5 parameters. We proceed similarly when building 1-d, 2-d, or 3-d profiled $\chi^{2}$.

We will show in the following, as an example, the constraints on the parameters from the above $\chi^{2}$. Nonetheless, these constraints are not necessarily the optimal ones. For example, the use of a flat prior on $R_{0}$ in the range $[7.5,8.5]$ kpc is perhaps too conservative and, instead, more stringent priors could be used, as for example the Gaussian prior $R_{0}=8.2\pm 0.1$ kpc based on Ref. [68]. Similar considerations apply to $\gamma$ or to the other parameters. This extra information can be easily included starting from the tables we provide. The main goal of this analysis is to provide results in a general form such that they can be used by the community together with complementary information, with the aim to simplify the use of a thorough data–driven approach on astrophysical uncertainties to analysis including direct and indirect DM searches, as well as collider probes.

4 Results

4.1 NFW $\gamma=1$ case

In Fig. 2, we show the result of the fit for $\gamma=1$, i.e. the canonical NFW case, and different values of $R_{0}$. As explained above we effectively fit the $\omega(x)$ data, but Fig. 2 shows the fit results and data-points in the $V(R)$ plane, which gives a more familiar representation. It can be also seen that since the data-points are at fixed $x$ values, when shown as function of $R$ they move along the $R$ axis for different $R_{0}$ values. The various curves show the best-fit contribution to the RC from the different baryonic components (Bulge, Disc and gas) and DM, as well as the total. The gas component, as explained previously, is the same since it is not varied in the fit. The morphology instead can be different for the various cases. Speficically, Bulge b is the preferred bulge morphology for all cases, except $R_{0}$= 7.5 kpc which is bulge e. All $R_{0}$ cases prefer Disc j, except for $R_{0}$ = 7.5 kpc which prefers Disc i. The plot also lists the best-fit $\chi^{2}$ which are in the range 7-10 for 25 data-points and 7 fitting parameters for a reduced $\chi^{2}$ $\simeq 8/(25-7)\simeq 0.5$. The value is slightly low and indicates that the binning procedure somehow overestimates the errors. Nonetheless, since this will give conservative results, we prefer not to modify our procedure.

Further results of the fit are shown in Fig. 3. The upper-left panel shows 2-$\sigma$ ($\Delta\chi^{2}=6.18$ from the minimum for two degrees of freedom) contours in the $\rho_{s}-R_{s}$ plane for fixed values of $R_{0}$ and profiled over $\mathcal{M}_{i}$, $\langle\tau\rangle$ and $\Sigma_{*}$, as well as for the case profiled over both $R_{0}$ and $\mathcal{M}_{i}$, $\langle\tau\rangle$ and $\Sigma_{*}$ (black contour), and shows a strong degeneracy between $\rho_{s}$ and $R_{s}$. Upper-right panel shows, however, that when we visualize the results in the $\rho_{0}-R_{s}$ plane, where $\rho_{0}$ is the local DM density (which is a derived parameter in our framework) the degeneracy basically disappears. To generate this plot we used as independent parameter to build the grid $\rho_{0}$, rather than $\rho_{s}$ (from Eq. 2.2 the two are related by $\rho_{0}=\rho_{s}\left(R_{0}/R_{s}\right)^{-\gamma}\left(1+R_{0}/R_{s}\right)^{-3+\gamma}$ ). Interestingly the plot also shows that the constraints on $\rho_{0}$ strongly depend on $R_{0}$. This is best seen in the lower-right panel where 1-$\sigma$ ($\Delta\chi^{2}=2.30$ from the minimum) and 2-$\sigma$ contours in the $R_{0}-\rho_{0}$ plane, profiled over the remaining parameters, are shown. The plot shows that the analysis is sensitive to $R_{0}$, although not strongly. This is reasonable, since $R_{0}$ is better constrained by different types of analysis than the RC ones (see [71] for a list of works on the determination of $R_{0}$). In the absence of strong priors on $R_{0}$, $\rho_{0}$ values between 0.3 and 0.8 GeV/cm3 at 2-$\sigma$ are allowed, which is in agreement with the conservative estimate provided in [56]. Interestingly, the combination $\rho_{0}=0.4$ GeV/cm3, $R_{0}=8$ kpc, often used in the literature, is in tension at 2$\sigma$ level with the fit result. The previous “standard” used until recently, $\rho_{0}=0.3$ GeV/cm3, $R_{0}=8.5$ kpc has a $\chi^{2}$ of 53.6 (in the profiled $\rho_{0}$, $R_{0}$ plane) and it’s excluded at more than 4 $\sigma$ confidence level.

Finally, the lower-left panel displays the 1-d $\chi^{2}$ plot of the morphology $\mathcal{M}_{i}$ profiled over the remaining parameters, and it shows that no single morphology is preferred by the analysis but all of them give similarly good best-fit $\chi^{2}$ at the level of 1-$\sigma$ or slightly more. This confirms the importance of considering different morphologies in order not to bias the final results. Comparing results obtained fixing the morphology to a single one we find that the systematic effect on $\rho_{0}$ is around $\pm 0.1$ GeV/cm3. In particular, the two morphologies which are found to give results which differ the most are the ones presenting a model of single disc vs the ones with a double disc (see [6, 7] for more details on the morphologies). An example of fit with a fixed morphology is discussed in the Appendix. The structure of the degeneracy in the $(\rho_{0},R_{0})$ plane is, instead, unchanged for each single morphology, i.e., the slope of the degeneracy remains the same and no significant sensitivity to $R_{0}$ is present.

4.2 Results as function of $\gamma$

In Fig. 4 we show the analogous of Fig. 3 for the case in which $\gamma$ is varied. It can be seen that the results are similar, except for the fact that when large values of $\gamma$ are used ($>1$), the largest values of $\rho_{0}$ (in the range $0.6-0.8$ GeV/cm3) are disfavored. Also, in general, no constraints can be inferred on $\gamma$. Including the uncertainty on $\gamma$, the couple $\rho_{0}=0.3$ GeV/cm3, $R_{0}=8.5$ kpc has now a $\chi^{2}$ of 42.3 (vs 53.6 when $\gamma$ is fixed to 1) which is still excluded at more than 4 $\sigma$ confidence level.

The best-fit RCs are shown in Fig. 5 for a fixed value of $R_{0}=8.0$ kpc and for different values of $\gamma$. It can be seen that despite strong differences in the RC contribution from DM in the cases $\gamma=0$ (i.e., cored profile) and $\gamma=1$ (i.e., cuspy profile), an equally good fit can be achieved in both cases to the measured RC. The main reason of this result is the degeneracy with the morphology. The uncertainties in the bulge and disc mass and morphology are large enough that can compensate in the two cases the large change from the DM contribution. The disc, in particular, seems to play a dominant role in this degeneracy, while the contribution from the bulge is slightly less prominent. This also means that in the future a more precise determination of the bulge and disc mass and morphology should be able to break this degeneracy and allow a reliable determination of the inner slope $\gamma$. A similar conclusion was reached by the work in [72], using a different analysis involving only observations within $\sim 2$ kpc from the Galactic Center.

4.3 Comparison with Other Results

As seen in the above sections, a general result of our analysis is that the single parameters are only weakly constrained by the fit (even for the case of fixed NFW profile, i.e., without varying $\gamma$). For example, at $2\sigma$ C.L. $\rho_{0}$ lies in the range 0.3-0.8 GeV/cm3, thus with an error of $\sim 0.25$ GeV/cm3, while $R_{0}$ is simply unconstrained by the analysis with respect to the prior range 7.5-8.5 kpc. What robustly constrained by the fit is, instead, the degeneracy and correlation among the parameters, like, noticeably, the one between $R_{0}$ and $\rho_{0}$. This is somewhat at odd with similar analyses performed in the past, which, typically, tend to find very small errors and strong constraints on the parameters. We attribute this difference to three main effects. First, the accurate statistical treatment, which explores and maps in details the degeneracies among the parameters. This is important, since if strong degeneracies are present, as in this case, and they are not well characterized, the error on the single parameters will be underestimated. Second, as detailed in Sec. 2.1, we don’t use stellar-tracer data up to $\sim$100 kpc, sometimes adopted in other analysis. The DM potential is the dominant component in the range 20-100 kpc, so these data could actually provide an important contribution in reducing the DM halo parameters errors, although this comes at the cost of adding further assumptions. Third, we use binned data, with an uncertainty estimated from the spread of the datapoints in the bin. This was already found to be an important point in [7] which shares the same dataset and similar methods as the present analysis. More precisely, when using the unbinned analysis, in [7] the authors report $\rho_{0}=0.420^{+0.021}_{-0.018}(2\sigma)\pm 0.025$ GeV/cm3 where the first error is statistical and the second comes from the baryonic morphology uncertainty. On the other hand, a test with binned data gives errors a factor of $\sim 5$ larger, thus more compatible with the analysis performed here. Other analyses giving small errors, like in [10] which reports $\rho_{0}=0.389\pm 0.025$ GeV/cm3 also use unbinned data. A noticeable exception is the work in [11], where a binned analysis is performed, with bin errors estimated in a similar way as in this work. We thus expect an uncertainty similar to the one of our analysis. The authors indeed find $\rho_{0}=0.471^{+0.048}_{-0.061}$ GeV/cm3, which has an uncertainty larger than [7, 10] but still smaller than our analysis. This is likely related to the larger dataset used in [11] which includes stellar velocity dispersion measurements at $R>20$ kpc, or to the simplified procedure used to estimate the errors, as explained in [11]. Finally, in [73], the authors perform a binned analysis using, for the inner Galaxy, the galkin dataset also employed here, in combination, for the outer Galaxy, with stellar tracers up to 100 kpc from [12]. The final uncertainties derived there are thus smaller. A further difference is that in [73] a fixed value of $R_{0}=8.34$ kpc is used, whereas one of the goals of this present analysis is indeed to estimate the very impact of the uncertainties on $R_{0}$ –which we therefore vary as discussed in the previous sections– on the determination of the DM distribution.

Several other analysis exist in the literature, such as e.g. [12, 13, 14, 15, 16], to which we address the reader for a sample of the wide range of methodologies, datasets and assumptions employed.

5 Implications for Direct and Indirect Dark Matter Searches

In this section, we provide examples of how to use the results derived above. In particular we consider the example of the Galactic Center $\mathcal{J}$-factor uncertainty for indirect DM searches and the uncertainties in direct DM searches.

5.1 Galactic Center $\mathcal{J}$ factor

An immediate application of the above analysis is the derivation of the GC $\mathcal{J}$-factor and its uncertainty. The GC $\mathcal{J}$-factor is given by

[TABLE]

where $\psi$ is the angle from the GC, $s$ a coordinate along the line of sight, and $\Delta\Omega$ the region of interest over which the angular integration is performed. In particular we consider the case of the Galactic Center excess (GCE) as given in [23]222See also [74, 75, 76, 77], and [78, 79, 80, 81] for an astrophysical interpretation of the excess. where the authors consider a square of 40${}^{\circ}\times$40∘ around the GC, with a stripe of $\pm 2^{\circ}$ along the Galactic Plane excluded, for a total area of 0.43 sr. Fig. 6 shows the $\chi^{2}$ profile of the GCE $\mathcal{J}$-factor from our analysis for different cases described in the caption. On technical note, we mention that when a derived parameter as $\mathcal{J}$ is involved, our frequentist profiling methodology is slightly more involved. In practice, first, for each point in our 7-d grid we derive the corresponding $\mathcal{J}$ value. Then, to build, for example, the 1-d profile $\chi^{2}$ for $\mathcal{J}$ we bin all the derived $\mathcal{J}$ values in a new $\mathcal{J}$ grid. For each $\mathcal{J}$ bin we then take the minimum $\chi^{2}$ among the $\chi^{2}$ corresponding to the $\mathcal{J}$ values falling in that bin. The final $\chi^{2}$ profile is shown in Fig. 6. This procedure can be easily generalized to 2-d cases, and is the more accurate the denser the original 7-d grid from which we start. Incidentally, we can see from Fig. 6 that the profile tend to a flat plateau at low $\mathcal{J}$ with a $\Delta\chi^{2}$ with respect to the minimum of $\Delta\chi^{2}\sim 90-7=83$. This, in practice, corresponds to the overall significance of our analysis to the presence of DM in the Galaxy which is thus $\sqrt{\Delta\chi^{2}}\sim 9\,\sigma$.

The methodology can be easily extended to the calculation of $\mathcal{J}$-factors over other regions for different analyses like GC searches for gamma-ray lines  [19, 20, 21, 22] or DM searches at TeV with Cherenkov telescopes [25, 26, 27, 28] where the considered region is of only few degrees, and thus even more sensitive to the uncertainties in the DM distribution.

5.2 Galactic Center Excess

Given the $\mathcal{J}$-factor and its uncertainty (or $\chi^{2}$ profile), it is easy to include it in the GCE analysis. The gamma-ray flux from the GCE is given by

[TABLE]

where $\langle\sigma v\rangle$ is the thermally averaged DM annihilation cross section, $m_{\chi}$ the DM particle mass and $dN/dE$ is the spectrum of gamma-ray photons from a single DM annihilation. For this last quantity we will use as example the case of annihilation into $b\bar{b}$ quarks taking the spectrum from [82]. We calculate the $\chi^{2}$ relative to the GCE as

[TABLE]

where $d_{i}$ are the GCE fluxes in the 24 energy bins given in [23], $t_{i}$ is the model prediction from Eq. 5.2 (to be more precise, in [23] fluxes are normalized to the area of the region analyzed, so from Eq. 5.2 we further need to divide by 0.43 sr), and $\Sigma_{ij}$ is the covariance matrix among the energy bins, again given in [23]. Similarly to what explained in [83], we further add to the covariance matrix a diagonal error equal to $\sigma_{\rm rel}$ per cent of the model prediction $t_{i}$ to account for the model uncertainty in the annihilation spectrum $dN/dE$. In particular, as explained in [83], a choice of $\sigma_{\rm rel}=10\%$ is appropriate.

To include in the GCE analysis the $\mathcal{J}$-factor uncertainties, we consider a $\chi^{2}$ with 3 contributions:

[TABLE]

where $\chi^{2}_{\rm GCE}$ is given by Eq. 5.3, $\chi^{2}_{\rm RC}(\mathcal{J},\gamma)$ is build from Eq. 3.3 profiling over $R_{0}$, $R_{s}$, $\rho_{s}$, $\langle\tau\rangle$, $\Sigma_{*}$ and $\mathcal{M}_{i}$ (but not $\gamma$) and $\chi^{2}_{\gamma}=(\gamma-1.2)^{2}/0.08^{2}$ is a Gaussian prior on $\gamma$, with mean 1.2 and $\sigma=0.08$ again coming from the analysis of the morphology of the GCE in [23]. Constraints in the plane $(m_{\rm DM},\langle\sigma v\rangle)$ derived from the above total $\chi^{2}$, further profiled over $\gamma$ and $\mathcal{J}$, are shown in Fig. 7 and compared with the results of [23]. We also cross-checked, for consistency, that fixing the $\mathcal{J}$-factor to the value adopted in [23] we obtain their same $(m_{\rm DM},\langle\sigma v\rangle)$ contours, which, for comparison are also shown in the same plot. It can be seen that including the DM distribution uncertainties significantly enlarges the contours, both at large and small $\langle\sigma v\rangle$. In particular, while the preferred region of [23] is in (mild) tension with the null observations of a gamma-ray signal from local dwarf galaxies [84, 85], this tension disappear when considering the $\mathcal{J}$-factor uncertainties. It should be mentioned, nonetheless, that the tension with dwarfs constraints can also be further relieved if more conservative estimates of the DM content of the dwarfs is adopted [86], or, similarly, if a more conservative analysis of the gamma-ray background at the dwarfs positions is performed [87, 88].

Simplified attempts to take into account the DM distribution uncertainties for the GCE excess have been performed in [89, 90, 91, 92]. In [89], similarly to here, variations with respect to the Galactic parameters are studied, but without performing a formal marginalization. In [90] the correlation among the Galactic parameters (in particular $\gamma$, $R_{s}$ and $\rho_{s}$) are extrapolated from [11], and the resulting GCE $\mathcal{J}$-factor is slightly overestimated with respect to our results (contours in Fig. 7 reach $\langle\sigma v\rangle$ values a factor of two lower). Finally, [91, 92] uses correlations among the Galactic parameters inferred from N-body simulations of MW-like DM haloes and priors on $\rho_{0}$ from local analyses.

Finally, we mention that further priors on $R_{0}$ or $R_{s},\rho_{s}$ or derived quantities like $\rho_{0}$ can be easily taken into account in our framework, if desired. In this case, one needs to use the more general form of the $\chi^{2}$

[TABLE]

where $\chi^{2}_{\rm RC}$ is now profiled only over $\langle\tau\rangle$, $\Sigma_{*}$ and $\mathcal{M}_{i}$ and $\chi^{2}_{R_{0},R_{s},\rho_{s}}$ contains the other priors to be implemented. The $\mathcal{J}$-factor, this time is to be intended more generally as a function of the DM parameters, $\mathcal{J}=\mathcal{J}(R_{0},R_{s},\rho_{s},\gamma)$. We derived the contours in the $(m_{\rm DM},\langle\sigma v\rangle)$ using this more general procedure and using, as in the previous case, only the Gaussian prior on $\gamma$ and, as expected, we obtained exactly the same results as using the simpler $\chi^{2}$ version in Eq. 5.4.

We provide the full $\chi^{2}_{\rm RC}(R_{0},R_{s},\rho_{s},\gamma)$ table on a grid of the four parameters. This represent the full information needed to reproduce the results of this section and to specialize to specific cases of need. For completeness, we also provide a table containing $\chi^{2}_{\rm RC}(\mathcal{J},\gamma)$ although this can be derived from the first table using the procedure described above. This second table can only be used for GCE analyses, though, since the $\mathcal{J}$-factor refers to the GCE ROI, while the first table is completely general and can be used in any analysis involving uncertainties in the Galactic DM distribution.

5.3 Direct Detection

The dependence of the results of direct searches for DM on the uncertainties of the properties of the Galactic DM Halo are an active topic of research, see e.g. [40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51]. Results of direct searches depends on the DM Halo in two ways. The main dependence is typically from $\rho_{0}$, the local DM density, which enters linearly in the expected DM detection rate (see e.g. [93]). Thus, the typical exclusion limits of the DM-proton scattering cross-section as function of DM mass (see, e.g., the recent XENON1T experiment results [94]), can just be linearly rescaled to a new $\rho_{0}$ value and weighted according to the $\rho_{0}$ likelihood. The second dependence is from the velocity distribution of DM particles $f(v)$. Assuming this distribution is a Maxwelliann one can simply express the velocity dispersion parameter, $\sigma_{v}$, entering in the Maxwellian as $\sigma_{v}=\sqrt{3/2}V_{0}$. Thus, under the Maxwellian approximation, our results can be used to take into account the uncertainties in $f(v)$ through the uncertainty in $V_{0}$. We do not attempt here to include the effect of the variation of $V_{0}$, although this can be implemented starting from our provided $\chi^{2}_{\rm RC}$ table.

Dropping the assumption of Maxwellian distribution can have major effects on the direct DM constraints, the larger the deviation of $f(v)$ from a Maxwellian (see the recent study [95]). However, under the reasonable assumption of isotropic DM velocity and of system at equilibrium $f(v)$ cannot be arbitrary but has to satisfy constraints given by the DM spatial distribution and the Boltzmann equation. This strategy to constrain $f(v)$ has been indeed pursued in various studies [40, 43, 44]. Again, these kind of studies can be in principle gereralized including the DM distribution uncertainties tabulated in the present study to derive the related uncertainties on the reconstructed $f(v)$, i.e., in practice, propagating the DM distribution uncertainties into the velocity distribution.

5.4 Combined Fits

Finally, a more subtle effect can appear when performing combined fits of GC constraints or hints of signals like the GCE and local observations like the above direct detection constraints or, for example, antiproton constraints (see e.g., [90]). In this case correlations might appear between the two observables. This is illustrated in Fig. 8 which shows that indeed there is a degeneracy between the GCE $\mathcal{J}$-factor and $\rho_{0}$, mainly produced by the variation of $R_{0}$. Again this can be taken into account using our tabulated likelihood.

6 Summary and Conclusion

We have used the observed Rotation Curve of the Milky Way up to $\approx 25$ kpc in Galactocentric radius to constrain the parameters of a generalized NFW Dark Matter density profile. We have improved with respect to previous analyses in several ways. First, we have adopted a systematic statistical approach, scanning the relevant parameter space and accurately exploring the various degeneracies present. This last point is particularly important since several degeneracies exist, and precisely mapping them is necessary to have reliable final error estimates. Second, we use an accurate treatment of the systematic uncertainties arising from the modeling of the visible components of the MW, by both considering different baryonic morphologies for the Disc and the Bulge, and allowing for each morphology mass variations within the uncertainties given by microlensing and stellar surface-density measurements. These baryonic uncertainties are fully marginalized (profiled) away within our statistical framework. We find that the local DM density $\rho_{0}$ is constrained to the range $0.3-0.8\;{\rm GeV/cm^{3}}$ at the $2\sigma$ level, showing a strong positive correlation with the Sun’s Galactocentric distance $R_{0}$. The inner slope of the DM profile, $\gamma$, is very weakly constrained and both core ($\gamma\sim 0$) and cusp ($\gamma\sim 1$) DM density profiles are allowed. Some combination of parameters can be, however, strongly constrained. For example the often used standard $\rho_{0}=0.3$ GeV/cm3, $R_{0}=8.5$ kpc is disfavored at more than 4 $\sigma$. We release the likelihood of our analysis, namely a 4-dimensional table listing $\chi^{2}$ values over a grid in $\gamma$, $R_{0}$, $R_{s}$, $\rho_{s}$, the latter two parameters being the scale radius and density scale of the generalized NFW profile. In the above likelihood the baryonic physics and related uncertainties have been already profiled away. We have provided some example for the use of the likelihood, in particular we have employed it in the analysis of the Galactic Center gamma-ray excess. We have found that the uncertainties in the DM profile significantly enlarge the allowed cross-section range, by a factor 3 to 4. Other contexts in which our tabulated likelihood can be employed involve Galactic Center or Galactic Halo DM searches in gamma rays at GeV energies or TeV with Cherenkov telescopes, DM neutrinos searches, direct DM searches, local DM searches with antimatter, and combined local and GC searches.

Acknowledgements

F. I. acknowledges support from the Simons Foundation and FAPESP process 2014/11070-2. This research was supported by resources supplied by the Center for Scientific Computing (NCC/GridUNESP) of the São Paulo State University (UNESP). This work has been possible following the visit of A.C. to ICTP-SAIFR in São Paulo, supported by a Theodore von Kármán fellowship from RWTH University of Aachen, Germany. We wish to thank Gabrijela Zaharias, Francesca Calore, Roberto Trotta and Jan Heisig for useful comments on the manuscript.

Appendix A: Additional Tests

In this appendix we study the effect of different binning and fitting methods on the results presented in the main text.

($x$,$v$) vs ($x$,$\omega$) fit

As explained in section 3.1, we perform the main fit in the ($x$,$\omega$) plane where the unbinned data-points have uncorrelated $x$ and $\omega$ uncertainties. Here, we test the effect of fitting, instead, in the ($x$,$v$) plane. Unbinned data in the ($x$,$v$) plane are binned according to the procedure described in section 3.1. Results are shown in Fig.9 for the example case of the $(\rho_{0},R_{0})$ plane . It can be seen that the results of the ($x$,$v$) fit are compatible with the ($x$,$\omega$) fit, although some differences can be seen, as a slight shift toward lower $\rho_{0}$ values of about 0.05 GeV/cm3, which is anyway small with respect to the overall width of the contours, and slightly different slope of the $(\rho_{0},R_{0})$ degeneracy. Despite the overall agreement, the small differences among the two fits suggest nonetheless that performing the analysis in the ($x$,$\omega$) plane is a more robust procedure since the properties of the errors (i.e., uncorrelated) are more straightforward.

Fit with a larger number of bins

We have, furthermore, tested the effect of increasing the number of bins used in the fit in the default $(x,\omega)$ procedure. In the present case the bins are chosen again starting from $x=2.5/8$ and ending at $x=22/8$, but we take i) 40 bins (up to 10.5/8) with a width of $\Delta x=0.2/8$ (up to 10.5/8), ii) 9 bins with $\Delta x$=0.5/8 (up to 15/8), iii) 3 bins with $\Delta x=1/8.0$ (up to 18/8) and iv) 2 last bins with $\Delta x=2/8$, for a total of 54 bins. This is roughly double with respect to the default setup, which has 25 bins. In each bin we use the procedure outlined in section 3.1 to derive the central value and the error. Results of the fit for this case are shown in Fig. 10 for the example case of the $(\rho_{0},R_{0})$ plane. As expected, the main effect is a reduction of the errors, as also discussed in section 4.3 and in ref. [7]. In particular, the allowed range for $\rho_{0}$ at 2$\sigma$ is now in the range $0.3.-0.6$ GeV/cm3. The analysis also becomes more sensitive to $R_{0}$, and values above 8.2 kpc are disfavored at 2$\sigma$.

Frequentist vs Bayesian

Another possible uncertainty is given by the use of the statistical methodology. To test this effect we compare our default methodology, which makes use of a grid in the parameter space and frequentist formalism, which a fully Bayesian analysis. To this purpose we use a simplified framework where we fix $R_{0}=8.34$ kpc and a single baryonic morphology, in particular the one labeled bJ (which assumes the E2 bulge given in [96] and the stellar disc from [97]). In this case, we thus have only five parameters, instead of the seven ones explored in the main analysis. To perform the Bayesian analysis we use a Monte Carlo scan of the parameter space with the emcee tool [98], and use flat priors on the parameters. The results are shown in the triangle plot of Fig. 11. The purple lines show the 1 and 2 $\sigma$ frequentist contours build with the method described in the main text, while the black lines show the analogous Bayesian result. The triangle diagonal shows the 1d Bayesian posterior for the single parameters. The triangle plot focus on the three DM Halo parameters $\gamma$, $R_{s}$ and $\rho_{0}$, and it does not show $\langle\tau\rangle$, $\Sigma_{*}$, although the fit is five-dimensional. As can be seen the Bayesian and frequentist contours are in excellent agreement. The only clear difference is that the frequentist contours are slightly larger, and thus more conservative. This is a typical result, especially when some of the parameters is not well constrained, as in this case. In the case when all the parameters are well constrained typically the agreement between the two methods is even closer (for example, see [99].)

Appendix B: Burkert Profile

In this appendix we discuss the results of the fit when the Burkert (BUR) profile [2] is adopted, i.e.,

[TABLE]

instead of the gNFW used in the main text. A simplification in this case arises from the fact that only two parameters define the model, i.e., $R_{c}$ the core radius, and, $\rho_{s}$ the scale density, instead of the three of the gNFW case. The left panel of Fig. 12 is the BUR analogue of the lower-right panel of Fig. 3 in the main text. The results for BUR and gNFW are fully compatible, with the BUR case giving a slightly tighter degeneracy between $R_{0}$ and $\rho_{0}$. The right panel shows the $(R_{c},\rho_{0})$ plane. An emerging interesting feature is that a minimum core size of about $\sim 5$ kpc is present. This appears to be a peculiarity of the BUR profile, while smaller core sizes should be possible if different profile parameterization are employed. The best-fit $\chi^{2}$ for the BUR case is $\sim 6$, similar to the gNFW case, indicating the two profiles can provide equally good fits to the Galactic rotation curve. As for the gNFW case, the tabulated likelihood in $R_{0},R_{c},\rho_{s}$, profiled over $\mathcal{M}_{i}$, $\langle\tau\rangle$ and $\Sigma_{*}$ is provided at https://github.com/mariabenitocst/UncertaintiesDMinTheMW.

Appendix C: Einasto Profile

Finally, we also derive results for another commonly employed profiled, i.e., the Einasto profile  [3] ,

[TABLE]

which is defined in terms of $R_{s}$, $\rho_{s}$ and $\alpha$, which is a shape parameters which plays a role similar to $\gamma$ for gNFW case, although with ‘opposite’ values, i.e., when $\alpha=1$ the profile is cored, while small values of $\alpha$ give a cuspy profile. Fig. 13 is the analogue of Fig. 4 in the main text, and it shows that the same degeneracies of the gNFW profile are present for the Einasto one. In particular, both cuspy ($\alpha\ll 1$) and cored ($\alpha=1$) profiles are compatible with the data. The best-fit $\chi^{2}$ for the Einasto case is $\sim 6$, similar to the gNFW and BUR case, indicating the all the profiles considered can provide equally good fits to the Galactic rotation curve. As for the gNFW and BUR case, the tabulated likelihood in $R_{0},R_{s},\rho_{s}$, $\alpha$ profiled over $\mathcal{M}_{i}$, $\langle\tau\rangle$ and $\Sigma_{*}$ is provided at https://github.com/mariabenitocst/UncertaintiesDMinTheMW.