On the estimation of the Local Dark Matter Density using the rotation curve of the Milky Way

TL;DR

This study assesses how different Galactic mass models influence the estimation of the local dark matter density using the Milky Way's rotation curve, highlighting the sensitivity to baryonic components and the limitations of the rotation curve method.

Contribution

It systematically analyzes the impact of baryonic model choices on dark matter density estimates, emphasizing the potential for unaccounted thin components to significantly alter results.

Findings

Estimated dark matter density is robust to halo model variations.

Uncertainty in baryonic components, especially the disk, affects density estimates.

Presence of a thin dark matter disk could increase density estimates by over a factor of 8.

Abstract

The rotation curve of the Milky Way is commonly used to estimate the local dark matter density $\rho_{{\rm DM},\odot}$. However, the estimates are subject to the choice of the distribution of baryons needed in this type of studies. In this work we explore several Galactic mass models that differ in the distribution of baryons and dark matter, in order to determine $\rho_{{\rm DM},\odot}$. For this purpose we analyze the precise circular velocity curve measurement of the Milky Way up to $\sim 25$ kpc from the Galactic centre obtained from Gaia DR2 [1]. We find that the estimated value of $\rho_{{\rm DM},\odot}$ stays robust to reasonable changes in the spherical dark matter halo. However, we show that $\rho_{{\rm DM},\odot}$ is affected by the choice of the model for the underlying baryonic components. In particular, we find that $\rho_{{\rm DM},\odot}$ is mostly sensitive to uncertainties in the disk components of the Galaxy. We also show that, when choosing one particular baryonic model, the estimate of $\rho_{{\rm DM},\odot}$ has an uncertainty of only about $10\%$ of its best-fit value, but this uncertainty gets much bigger when we also consider the variation of the baryonic model. In particular, the rotation curve method does not allow to exclude the presence of an additional very thin component, that can increase $\rho_{{\rm DM},\odot}$ by more than a factor of 8 (the thin disk could even be made of dark matter). Therefore, we conclude that exclusively using the rotation curve of the Galaxy is not enough to provide a robust estimate of $\rho_{{\rm DM},\odot}$. For all the models that we study without the presence of an additional thin component, our resulting estimates of the local dark matter density take values in the range $\rho_{{\rm DM},\odot} \simeq \text{0.3--0.4}\,\mathrm{GeV/cm^3}$, consistent with many of the estimates in the literature.

Tables

Table 1. Table 1: Fitted and derived quantities obtained from the analysis of the baryonic model B1 when their parameters are fixed, as in the analysis of ρ DM , ⊙ subscript 𝜌 DM direct-product \rho_{{\rm DM},\odot} presented by Eilers:1810.09466 , and an NFW dark matter halo (section  3.1.1 ). Top rows correspond to parameters of the fit, while the last rows are derived quantities. The values correspond to the maximum and 68% credible region of the marginal posteriors.

	without systematics	with systematics
$M_{200}\left[{10}^{11}{M}_{\odot }\right]$	$7.25\pm 0.26$	${6.5}_{-1.2}^{+1.6}$
$c_{200}$	$12.7\pm 0.3$	${12.9}_{-1.5}^{+1.7}$
${\rho }_{\mathrm{DM},\odot }\left[\mathrm{GeV}/{\mathrm{cm}}^3\right]$	${0.2981}_{-0.0023}^{+0.0022}$	${0.293}_{-0.010}^{+0.009}$
$r_{200}\left[\mathrm{kpc}\right]$	${189.4}_{-2.3}^{+2.1}$	${184}_{-12}^{+14}$
$r_s\left[\mathrm{kpc}\right]$	$14.9\pm 0.5$	${13.7}_{-2.3}^{+3.1}$

Table 2. Table 2: Dark matter related quantities obtained from the fit of the baryonic model B1, when the parameters of the baryonic components are allowed to vary, using an NFW, a gNFW gNFW {\rm gNFW} and an Einasto dark matter halo. For the last two columns γ 𝛾 \gamma and α 𝛼 \alpha are extra free parameters of the analysis. Top rows correspond to fitted variables and the last rows are derived quantities. The values correspond to the maximum and 68% credible region of the marginal posteriors.

	NFW	$\mathrm{gNFW}$	Einasto
$M_{200}\left[{10}^{11}{M}_{\odot }\right]$	${5.2}_{-1.1}^{+2.0}$	${5.5}_{-1.4}^{+3.1}$	${2.8}_{-1.2}^{+7.7}$
$c_{200}$	${15}_{-4}^{+5}$	$14\pm 5$	$12\pm 4$
Slope parameter	$\gamma =1$	$\gamma ={1.2}_{-0.8}^{+0.3}$	$\alpha ={0.11}_{-0.05}^{+0.20}$
${\rho }_{\mathrm{DM},\odot }\left[\mathrm{GeV}/{\mathrm{cm}}^3\right]$	${0.301}_{-0.025}^{+0.028}$	${0.300}_{-0.027}^{+0.028}$	$0.301\pm 0.027$
$r_{200}\left[\mathrm{kpc}\right]$	${173}_{-13}^{+19}$	${174}_{-15}^{+29}$	${182}_{-51}^{+43}$
$r_s\left[\mathrm{kpc}\right]$	${10}_{-3}^{+5}$	$9_{-8}^{+12}$	${11}_{-4}^{+10}$

Table 3. Table 3: Dark matter related quantities obtained from the fit of the baryonic model B2 with an NFW, a gNFW gNFW {\rm gNFW} and an Einasto dark matter halo. Top rows correspond to fitted variables and the last rows are derived quantities. As usual, for the NFW halo the inner slope is fixed to be γ = 1 𝛾 1 \gamma=1 . The values correspond to the maximum and 68% credible region of the marginal posteriors.

	NFW	$\mathrm{gNFW}$	Einasto
$M_{200}\left[{10}^{11}{M}_{\odot }\right]$	${7.4}_{-1.5}^{+1.8}$	${6.3}_{-1.3}^{+3.4}$	${3.0}_{-1.2}^{+5.7}$
$c_{200}$	${16}_{-3}^{+4}$	$17\pm 6$	${14}_{-4}^{+5}$
Slope parameter	$\gamma =1$	$\gamma ={1.3}_{-0.9}^{+0.3}$	$\alpha ={0.18}_{-0.09}^{+0.21}$
${\rho }_{\mathrm{DM},\odot }\left[\mathrm{GeV}/{\mathrm{cm}}^3\right]$	$0.376\pm 0.025$	${0.387}_{-0.036}^{+0.034}$	${0.384}_{-0.034}^{+0.038}$
$r_{200}\left[\mathrm{kpc}\right]$	${192}_{-13}^{+15}$	${184}_{-14}^{+29}$	${147}_{-19}^{+59}$
$r_s\left[\mathrm{kpc}\right]$	${11}_{-3}^{+4}$	${8.1}_{-7.8}^{+10.6}$	${9.2}_{-2.7}^{+5.3}$

Table 4. Table 4: Comparison between the mass of the Milky Way obtained from the mass models of our analyses and estimates from other studies. Marginal 68% credible regions are shown. The first three columns correspond to our baryonic model B1 and columns 4–6 to the baryonic model B2. Different radial distances are considered.

	Baryonic model B1			Baryonic model B2
	NFW	$\mathrm{gNFW}$	Einasto	NFW	$\mathrm{gNFW}$	Einasto	Reference value	Ref.
	6.5–9.8	6.5–10.9	3.6–15.6	7.9–11.9	6.8–12.5	3.4–10.6	$14\pm 3$	Watkins:2010fe
	4.4–5.9	4.5–6.4	3.6–7.8	5.1–6.9	4.8–7.2	3.4–7.2	$4.1\pm 0.4$	Gibbons:2014ewa
	3.2–3.9	3.3–4.0	3.1–4.5	3.6–4.4	3.5–4.5	3.2–4.6	${3.7}_{-0.3}^{+0.4}$	Eadie:1810.10036
$M_{\mathrm{MW}}\left(R\lesssim 20\mathrm{kpc}\right)\left[{10}^{11}{M}_{\odot }\right]$	1.9–2.0	1.9–2.0	1.8–2.0	1.9–2.1	1.9–2.1	1.9–2.1	${2.1}_{-0.3}^{+0.4}$	Watkins:1804.11348

Equations

$$v_{\mathrm{c}}^2(R)=R{\frac{\partial \Phi }{\partial R}|}_{z=0}.$$

$${\nabla }^2{\Phi }_i=4\pi G{\rho }_i,$$

$${\Phi }_{\mathrm{Plummer}}(r)=-\frac{G{M}_{\mathrm{bulge}}}{\sqrt{r^2+r_b^2}},$$

$${\Phi }_{\mathrm{MN}}(R,z)=-\frac{G{M}_{\mathrm{disk}}}{\sqrt{R^2+{\left(R_d+\sqrt{z^2+z_d^2}\right)}^2}},$$

$$\rho (R,z)={\rho }_0\exp \left(-\frac{R}{R_d}-\frac{|z|}{z_d}\right),$$

$${\Phi }_{\mathrm{Hernquist}}(r)=-\frac{G{M}_{\mathrm{bulge}}}{r_b+r},$$

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

$${\rho }_{\mathrm{Ein}}(r)={\rho }_0\exp \left\{-\frac{2}{\alpha }\left({\left(\frac{r}{r_s}\right)}^{\alpha }-1\right)\right\},$$

$$M_{200}=\frac{4\pi }{3}{r}_{200}^3{\Delta }_{200}{\rho }_{\mathrm{crit}},$$

$$c_{200}=\frac{r_{200}}{r_{-2}},$$

$$M_{200}=4\pi {\int }_0^{r_{200}}{\rho }_{\mathrm{DM}}(r^{\prime },r_s,{\rho }_0,\eta )r^{\prime 2}\mathrm{d}{r}^{\prime }.$$

$$\ln \mathcal{L}=-\frac{1}{2}\sum _i{\left(\frac{v_{\mathrm{c}}^{\mathrm{d}}(R_i)-v_{\mathrm{c}}^{\mathrm{m}}(R_i)}{{\sigma }_i}\right)}^2,$$

$$v_{\mathrm{c}}^{\mathrm{m}}(R)={\left(v_{\text{c,DM}}^2(R)+v_{\mathrm{c},\mathrm{B}}^2(R)\right)}^{1/2},$$

$${\log }_{10}{c}_{200}=a+b{\log }_{10}\left(M_{200}/[10^{12}{h}^{-1}{\mathrm{M}}_{\odot }]\right),$$

Full text

1. 1 Introduction

2. 2 Galactic models

3. 2.1 Baryonic models

4. 2.1.1 Baryonic model B1

5. 2.1.2 Baryonic model B2

6. 2.2 Dark matter profiles

7. 3 Analyses

8. 3.1 Analysis of model B1

9. 3.1.1 Fixing baryonic parameters as in the case studied by Eilers:1810.09466

10. 3.1.2 Varying baryonic parameters

11. 3.1.3 Changing the dark matter profile

12. 3.2 Analysis of model B2

13. 3.2.1 NFW dark matter profile

14. 3.2.2 Changing the dark matter profile

15. 4 Discussion

16. 4.1 Estimated local dark matter density

17. 4.2 Sensitivity to disk uncertainties

18. 4.3 Studying the $c_{200}\text{-}M_{200}$ relation

19. 4.4 Estimated Galactic mass

20. 5 Conclusions

21. A Triangular plots of the dark matter parameters

1 Introduction

An accurate and robust estimate of the local dark matter density $(\rho_{{\rm DM},\odot})$ at the location of the Sun has been a goal in astrophysics for a long time. A correct value is imperative for the interpretations of dark matter signals in a variety of experimental efforts. For example, in direct detection experiments searching for Weakly Interacting Massive Particle (WIMPs), the parameter $\rho_{{\rm DM},\odot}$ is inversely proportional to the derived scattering cross section of the WIMP/nucleus interaction. Furthermore, since $\rho_{{\rm DM},\odot}$ is also linked to the normalization of the dark matter halo, a precise measurement of this quantity plays an important role for indirect detection searches, as well as in cosmology and astrophysics for those studies that rely on the dark matter distribution of the Milky Way.

One possible way of constraining $\rho_{{\rm DM},\odot}$ comes from the analysis of the rotation curve of the Milky Way. This method is based on the comparison of the observationally determined circular velocity (obtained from the measured rotation curve of our Galaxy111We note that the rotation curve and the circular velocity curve are often confused. The rotation curve is an observable quantity representing the mean rotational velocity at different radii. The measurement of the rotation curve of a galaxy is an essential part of determining the circular velocity curve, and this field is pioneered by the work of Vera Rubin Rubin:1980zd for external galaxies. with the radial Jeans equation) and its theoretical estimates at different radii, $v_{\rm c}(R)$, where $R$ is the radius in the galactocentric cylindrical coordinates. The theoretical estimates of the circular velocity are derived from the potential $\Phi$ of the Galaxy through the relation

[TABLE]

The potential $\Phi$ is the sum of the contribution of different components $\Phi_{i}$. In our study we consider a spherical dark matter halo, a spherical Galactic bulge and axisymmetric distributions for the Galactic disk. In those cases in which a distribution is given in terms of its energy density $\rho_{i}$, the potential is obtained through the Poisson equation

[TABLE]

where $G$ is the gravitational constant. For a comprehensive review of the rotation curve’s method and other methods to infer $\rho_{{\rm DM},\odot}$ we refer the reader to Read:2014qva .

Given our location in the Galaxy, determining the Milky Way’s rotation curve to a good precision has been a persistent challenge (Sofue:2011kw, ; Russeil:2017aa, ). This has been partly due to the complexities that come into play from the dense and obscuring inter-galactic gas and dust present in the Galaxy’s disk. However, the major drawback so far had been the lack of well measured 3D velocities and distances of the stars in the Milky Way. This handicap has now been overcome with the help of the second data release (DR2) of the ESA/Gaia mission (GaiaDR2_2018_Brown, ; GaiaCollab2018kinematics, ).

Gaia DR2 has provided extremely well measured proper motions and 2D positions on the sky for a significantly large volume of our Galaxy, as well as precise parallaxes within $\sim\text{2--3}{\rm\,kpc}$ from the solar neighbourhood. Combining this exquisite data with photometric information from 2MASS and WISE, and spectroscopic information from APOGEE, Hogg:1810.09468 estimated precise spectrophotometric parallaxes for red-giant stars up to $25{\rm\,kpc}$ from the Galactic centre. Using this information, Eilers:1810.09466 determined the Milky Way’s rotation and circular velocity curves222Likewise based on Gaia DR2 data, Mroz:1810.02131 also measured the rotation curve of our Galaxy. Their result is consistent with that of Eilers:1810.09466 . for the distance range of $5{\rm\,kpc}\leq R\leq 25{\rm\,kpc}$. Analyzing their accurate measurement of the rotation curve of the Galaxy, they also provide an estimate of $\rho_{{\rm DM},\odot}$. In order to infer the value of $\rho_{{\rm DM},\odot}$ from the rotation curve, they use a Galactic mass model in which they fix the baryonic components and fit a spherical dark matter halo. With their assumptions, Eilers:1810.09466 quote a result of $\rho_{{\rm DM},\odot}=0.30\pm 0.03{\rm\,GeV/cm^{3}}$.

In this paper we perform a more general study. We explore a variety of different Galactic mass models, comprised of different baryonic and dark matter components, to study the corresponding effects on the estimate of $\rho_{{\rm DM},\odot}$ from rotation curve measurements. Specifically, we reanalyze the circular velocity curve data from Eilers:1810.09466 to show that uncertainties in the baryonic components can lead to a significant uncertainty in the best-fit value of $\rho_{{\rm DM},\odot}$.

The paper is arranged as follows: In section 2 we detail all the mass models that we use in our study; in section 3 we present the analyses of the circular velocity curve based on different mass models; the results are compared and discussed in section 4, and we draw our conclusions in section 5.

2 Galactic models

In order to study the influence that the choice of a specific mass model has in the determination of $\rho_{{\rm DM},\odot}$ from the rotation curve of our Galaxy, we analyze a selection of models with different mass distributions. All the mass models that we study contain two main components: the baryonic part (in the form of disks and a bulge) and a dark matter halo. Apart from axisymmetry, we also assume in our analysis that the Galaxy is in a state of dynamical equilibrium. Below we describe the baryonic and dark matter halo components that we investigate in our study and the motivation behind these choices. We first lay out all the considered baryonic models, followed by those of the dark matter halo.

2.1 Baryonic models

One of the fundamental tasks of Galactic Astronomy is to determine the luminosity distribution of the different baryonic components of the Milky Way—the stars, the gas and the dust. Their luminosity profiles effectively deliver their mass profiles. Various studies have been dedicated to this endeavour, but we yet lack a consensus in the baryonic model for the Milky Way (see e.g. Freeman:2002wq ; Bland-Hawthorn:1602.07702 ). For this reason, we iterate our analysis over two different baryonic models that are backed up by different observational studies. The choice of two different models also allows us to examine the effect that the baryonic matter distribution has on the estimates of $\rho_{{\rm DM},\odot}$.

2.1.1 Baryonic model B1

We refer to our first mass model of baryons as B1. It essentially corresponds to the baryonic part of the Galactic model used by Eilers:1810.09466 , that they fix in their analysis. However, we consider both the case of parameters fixed to agree with their work and also the more general case in which we allow parameters to vary. This model is based on Model I of Pouliasis:1611.07979 . It consists of three components, namely, a bulge and two disks (divided in terms of a thin and a thick distribution). The bulge is modeled as a Plummer potential (Plummer:1911zza, ) given as

[TABLE]

where $r$ is the distance in galactocentric Cartesian coordinates and $r_{b}$ is the cut-off radius. The thin and the thick disks are represented independently by Miyamoto-Nagai potentials (Miyamoto:1975zz, ) that are expressed as

[TABLE]

where $R$ is the radius in the galactocentric cylindrical coordinates, and $R_{d}$, $z_{d}$ are the characteristic scales of the profile.

Whenever we fix the baryonic components of this model in our analyses, we set all the parameters at the same values used by Eilers:1810.09466 , as described in Pouliasis:1611.07979 , with $M_{\rm bulge}=1.067\times 10^{10}{\rm\,M_{\odot}}$ and $r_{b}=0.3{\rm\,kpc}$ for the bulge, and $M_{\rm thin}=M_{\rm thick}=3.944\times 10^{10}{\rm\,M_{\odot}}$, $R^{\rm thin}_{d}=5.3{\rm\,kpc}$, $z^{\rm thin}_{d}=0.25{\rm\,kpc}$, $R^{\rm thick}_{d}=2.6{\rm\,kpc}$ and $z^{\rm thick}_{d}=0.8{\rm\,kpc}$ for the disks. The same parametric setting is used for the mean values of the Gaussian priors when we allow the baryonic components to vary in the fits. The chosen values of the baryonic components are motivated by several observations, among which Pouliasis:1611.07979 used the local baryonic surface and volume densities of the Milky Way. The used observations are based on different studies (Bensby:2011ae, ; Bovy:2011ux, ; Bland-Hawthorn:1602.07702, ; Kafle:2014xfa, ; McKee:2015hwa, ; Flynn:2006tm, ). As we show in section 3, we study both the cases where we fix these parameters and also those in which we allow some variation.

2.1.2 Baryonic model B2

Miyamoto-Nagai disks have density profiles that do not decline as rapidly with radius as observed disk galaxies, and are expected to overestimate the mass of the baryons towards the outer Galactic radii (Rix:2013bi, ). For this reason, we study another baryonic model that we refer to as B2.

The motivation for B2 comes mainly from the observational studies of Misiriotis:2006qq , where they analyzed COBE dust emission maps (Sodroski:1993mn, ; Fixsen:1996nj, ) to constrain the parameters of their Galactic model. Such model is comprised of axisymmetric distributions for the stars (composed of a bulge and a disk), the dust (cold and warm) and the gas (molecular $\mathrm{H}_{2}$ and atomic HI). The stellar disk, the two dust components and the $\mathrm{H}_{2}$ gas distribution are all modelled as double exponential profiles expressed as

[TABLE]

where $\rho_{0}=M/(4\pi z_{d}R_{d}^{2})$ is the normalization, $M$ is the corresponding mass, and $R_{d}$ and $z_{d}$ are the characteristic scale length and height, respectively. We set these parameters as $M_{\rm disk}=3.65\times 10^{10}\,\mathrm{M}_{\odot}$ (Portail:1502.00633, ; Bland-Hawthorn:1602.07702, ), $R_{d}^{\rm disk}=2.35{\rm\,kpc}$, $z_{d}^{\rm disk}=0.14{\rm\,kpc}$; $M_{\rm cold}=7.0\times 10^{7}{\rm\,M_{\odot}}$, $R_{d}^{\rm cold}=5.0{\rm\,kpc}$, $z_{d}^{\rm cold}=0.1{\rm\,kpc}$; $M_{\rm warm}=2.2\times 10^{5}{\rm\,M_{\odot}}$, $R_{d}^{\rm warm}=3.3{\rm\,kpc}$, $z_{d}^{\rm warm}=0.09{\rm\,kpc}$, and $M_{\rm H_{2}}=1.3\times 10^{9}{\rm\,M_{\odot}}$, $R_{d}^{\rm H_{2}}=2.57{\rm\,kpc}$, $z_{d}^{\rm H_{2}}=0.08{\rm\,kpc}$ (Misiriotis:2006qq, ). The atomic HI gas distribution is also modelled as a double exponential, as per Eq. (5), with its mass defined as $M=4\pi\rho_{0}z_{d}R_{d}(R_{t}+R_{d})\,\mathrm{e}^{-R_{t}/R_{d}}$, where $M_{\rm HI}=8.2\times 10^{9}{\rm\,M_{\odot}}$, $R_{d}^{\rm HI}=18.24{\rm\,kpc}$, $z_{d}^{\rm HI}=0.52{\rm\,kpc}$ and $R_{t}=2.75{\rm\,kpc}$, as obtained in Misiriotis:2006qq . This mass definition corresponds to that of a truncated disk, although we are using a double exponential shape without truncation for the sake of easier numerical analyses. This is feasible because the contribution that gas and dust terms adds to the global circular velocity curve is significantly smaller than that of the remaining components. Therefore, this minor approximation does not cause any marked effect. We however include these components in our total baryonic mass model for consistency and in order to have a more realistic distribution of baryons.

For the bulge, we substitute the original profile used in Misiriotis:2006qq with a Hernquist potential that is given as

[TABLE]

where we set $M_{\rm bulge}=1.55\times 10^{10}\,\mathrm{M}_{\odot}$ (Portail:1502.00633, ; Bland-Hawthorn:1602.07702, )333The value of $M_{\rm bulge}$ estimated by Portail:1502.00633 was obtained from a triaxial modelling of the bulge region. and $r_{b}=0.7{\rm\,kpc}$. This change is also motivated by an easier numerical calculation, since the original profile used in Misiriotis:2006qq has a more complicated shape. We decided to use a Hernquist profile (Eq. (6)) instead of a Plummer potential (Eq. (3)) because a Hernquist profile with $r_{b}=0.7{\rm\,kpc}$ adjusts better, given the same mass, to the bulge’s shape used in Misiriotis:2006qq . In our analysis, the values of $M_{\rm bulge}$ and $M_{\rm disk}$ (or equivalently the central values of their Gaussian prior distributions, when these parameters are allowed to vary) were chosen such that $M_{\rm bulge}/(M_{\rm bulge}+M_{\rm disk})=0.3$ (Bland-Hawthorn:1602.07702, ).

2.2 Dark matter profiles

One of the long-standing problems of astrophysics and cosmology has been to determine the dark matter distribution around the Milky Way. Since a change of the dark matter halo can affect the estimate of $\rho_{{\rm DM},\odot}$ from the Milky Way’s rotation curve, we use two different parameterizations for the dark matter component, allowing for different shapes that deviate particularly towards the centre. Specifically, one of the profiles that we use presents a cusp in the centre (see Eq. (7) below) while the other one can exhibit a cored halo (see Eq. (8) below) for some values of its free parameter $\alpha$.

We describe our first dark matter halo as a generalized version of the Navarro-Frenk-White (NFW) profile, that we refer to as gNFW,

[TABLE]

where $\rho_{0}$ is the normalization constant, $r_{s}$ the scale radius and $\gamma$ the inner slope. For $\gamma=1$ the profile becomes the well-known NFW halo ($\rho_{\rm NFW}$, Navarro:1996gj ). This profile is a common approximation to dark matter densities found in cosmological simulations, and it represents a cuspy profile that diverges towards smaller $r$ values. We incorporate both cases in our study, one where $\gamma$ is a free parameter (varying in the range $\gamma=\text{0--2}$) and another where we set $\gamma=1$ (which then takes the form of the usual NFW). 2. 2.

The second dark halo that we consider is the Einasto profile (Einasto:1965czb, ), which is expressed as

[TABLE]

where $\alpha$ determines how fast the density distribution falls with $r$, making it cored towards the central regions of the Galaxy when its value is close to $\alpha=1$. We allow for values of this parameter in the range $\alpha=\text{0--1}$.

In the next section, the parameters of the dark matter halo that we will directly constrain in our analyses (comparing the Milky Way’s circular velocity curve data from Eilers:1810.09466 with the outcome of our Galactic mass models) are the virial mass $M_{200}$, the virial concentration $c_{200}$, and the $\gamma$ and $\alpha$ parameters. With this definition, $M_{200}$ is the mass contained within the radius $r_{200}$ such that the energy density is 200 times larger than the critical energy density $\rho_{\rm crit}$ of the Universe,

[TABLE]

where $\Delta_{200}=200$. Notice that this is just a matter of definition, so choosing this or another virial overdensity (like the one derived from the collapse of a spherical top-hat perturbation Bryan:1997dn ) does not affect our results.

One remark worth making is about the concentration parameter, defined as

[TABLE]

where $r_{-2}$ is the radius at which the slope $\mathrm{d}\,\mathrm{ln}\rho/\mathrm{d}\,\mathrm{ln}r=-2$. In the $\rho_{{\rm gNFW}}$ case, $r_{-2}=(2-\gamma)r_{s}$. This means that for very steep ${\rm gNFW}$ halos, with values of $\gamma$ close to 2, $r_{s}$ becomes extremely large. On the other hand, for $\rho_{\rm Ein}$ we have the relation $r_{-2}=r_{s}$.

Defining $\eta$ as the parameter $\gamma$ for the ${\rm gNFW}$ profile and $\alpha$ for the Einasto profile, the connection between the varied dark matter quantities of our analyses, $\left\{M_{200},c_{200},\eta\right\}$, and a given dark matter profile $\rho_{\rm DM}(r,r_{s},\rho_{0},\eta)$, is made through Eqs. (9), (10) and

[TABLE]

For each set of values $\left\{M_{200},c_{200},\eta\right\}$ we first determine $r_{200}$ using Eq. (9), then we obtain $r_{s}$ from Eq. (10), and finally we get the value of the normalization $\rho_{0}$ that determines $\rho_{\rm DM}(r,r_{s},\rho_{0},\eta)$ solving Eq. (11).

3 Analyses

In this section, we present our analyses of the Milky Way’s rotation curve based on the implementation of the aforementioned Galactic mass models. To this end, we perform Bayesian analyses, comparing the circular velocity measurements against the circular velocity values obtained from the mass models, in order to constrain $\rho_{{\rm DM},\odot}$. In particular, we run a Markov Chain Monte Carlo (MCMC) using the affine invariant sampler emcee (ForemanMackey:2012ig, ) implemented in Python, with the log-likelihood function defined by

[TABLE]

where the summation $i$ is done over all the data points. $v_{\rm c}^{\rm d}(R_{i})$ is the measurement of the circular velocity at a given radius $R_{i}$ in the galactocentric cylindrical coordinates, and the corresponding model value is given by $v_{\rm c}^{\rm m}(R_{i})$. The Gaussian dispersion $\sigma_{i}$ is quadratically summed over the statistical and systematic uncertainties associated to the $i$-th data point, where the statistical uncertainties are taken from Tab. 1 of Eilers:1810.09466 and the systematic uncertainties are extracted from their Fig. 4, assuming for the last bins, which are saturated in the figure, a 12% deviation that corresponds to the maximum deviation shown in that plot. When quadratically summing both statistical and systematic uncertainties we are simplifying the treatment of systematic uncertainties. In particular, we are adding extra freedom that is not necessarily allowed by the data, since in principle there are sources of systematic uncertainties that imply some correlation in the radial direction. Alternatively, a better approach could be to use a Gaussian process, or to repeat the analyses by systematically shifting the data points within the allowed range of systematic uncertainties, with $\sigma_{i}$ being only the statistical uncertainties in the second case. Although the estimated systematic error that we get might be larger than the real error caused by systematic uncertainties, our chosen treatment is simple from a computational point of view and it provides a reasonable result, since the size of systematic uncertainties is not large enough to produce a dangerous discrepancy between the different approaches.

Given a theoretical mass model, we evaluate $v_{\rm c}^{\rm m}(R)$ as

[TABLE]

where $v^{2}_{\text{c,DM}}(R)$ is the contribution to the circular velocity coming from dark matter and $v^{2}_{\rm c,B}(R)$ is the summed contribution of the different baryonic components, as described in section 2. The circular velocity for each component is obtained from its potential using Eq. (1).

Furthermore, the uncertainty in the estimate of $\rho_{{\rm DM},\odot}$ is related to the dependence of this quantity upon the solar distance from the Galactic centre, $R_{\odot}$. Since the value of $R_{\odot}$ also has some associated uncertainty, it is something that needs to be taken into account. Therefore, $R_{\odot}$ is included as an additional parameter over which we sample in our MCMC. In particular, we use $R_{\odot}=8.122\pm 0.031\,\mathrm{kpc}$ (GravityColl:2018, ), which is the same value adopted by Eilers:1810.09466 . Larger uncertainties are usually adopted for this quantity (see e.g. Bland-Hawthorn:1602.07702 ). However, throughout our study we observe that the estimates of $\rho_{{\rm DM},\odot}$ are mostly susceptible to changes in the underlying distribution of baryons.

In all our analyses we consider flat priors for the dark halo parameters: $M_{200}$, $c_{200}$, $\alpha$ and $\gamma$, that can vary in the ranges $M_{200}/[10^{11}{\rm\,M_{\odot}}]\in[0.01,30]$, $c_{200}\in[0,100]$, $\alpha\in[0,1]$ and $\gamma\in[0,2]$. We also tested that choosing a broader flat prior in $\log_{10}M_{200}$ does not affect our results.

3.1 Analysis of model B1

We first study our baryonic model B1, defined in section 2.1.1, in combination with different dark halo profiles, to find estimates for $\rho_{{\rm DM},\odot}$ from matching to the circular velocity curve data. We begin in section 3.1.1 by studying exactly the same case as done in Eilers:1810.09466 , using their values of baryonic parameters together with a standard NFW halo. Our resulting $\rho_{{\rm DM},\odot}$ agrees with theirs, but we find even smaller error bars for this particular case, as we will explain. Then in section 3.1.2 we generalize by allowing the values of the baryonic parameters to vary, and in section 3.1.3 we generalize further by considering a variety of dark matter halo profiles.

3.1.1 Fixing baryonic parameters as in the case studied by Eilers:1810.09466

In this section we first study the value of $\rho_{{\rm DM},\odot}$ obtained with exactly the same mass model as in the illustration case presented by Eilers:1810.09466 . We use the model B1 for the baryonic components with the parameters fixed to the values presented in section 2.1.1, and we use an NFW profile for the dark halo ($\rho_{{\rm gNFW}}$ with inner slope $\gamma=1$), thereby fitting for the NFW parameters $M_{200}$ and $c_{200}$. In addition, we distinguish two cases, one in which we do not include the systematic uncertainties of the data in the fit and another in which we include them. We make this distinction because we are able to reproduce the results of Eilers:1810.09466 only when systematics are not included.

The best-fit circular velocity curve when we include systematic uncertainties is shown in the left panel of Fig. 1, together with the individual contribution from each Galactic component. The corresponding fitted values for the parameters of the dark matter profile are provided in Tab. 1. Comparing these results with those of Eilers:1810.09466 , we find that we are successfully able to reproduce their analysis only when no systematics are included, with the single exception of the uncertainty associated with the derived $\rho_{{\rm DM},\odot}$, for which we find smaller error bars. We think that the reason is due to a correlation found in our analysis between the two fitted quantities $M_{200}$ and $c_{200}$ (see Fig. 2). Since $\rho_{{\rm DM},\odot}$ is a derived parameter that depends on $M_{200}$, $c_{200}$ and $R_{\odot}$, one needs to propagate the error of the virial quantities in order to infer the error of $\rho_{{\rm DM},\odot}$. If we neglect the correlation shown in Fig. 2 and compute the uncertainty of $\rho_{{\rm DM},\odot}$ from quadratic error propagation of $R_{\odot}$, $M_{200}$ and $c_{200}$, assuming that they are independent, we get an uncertainty that is consistent with Eilers:1810.09466 findings. It is important to highlight that we correct this aspect and account for the correlation in the rest of the study.

3.1.2 Varying baryonic parameters

In the previous section we fixed the baryonic components of the Galaxy, but our knowledge on the distribution of baryons in the Milky Way is not perfect. Therefore, as a first approach to test the dependence of $\rho_{{\rm DM},\odot}$ on the baryonic parameters, we free our model B1 to let it vary in our analysis. However, in order to retain the essence of the model, we impose a Gaussian prior on the parameters, taking as central values those discussed in section 2.1.1, with a standard deviation corresponding to 15% of the central values for the masses and 10% for the characteristic scales. These standard deviations enter in the range of uncertainties presented in Bland-Hawthorn:1602.07702 ; Pouliasis:1611.07979 . As in the previous section, we assume a standard NFW halo for the dark matter.

The circular velocity curves obtained from this analysis are shown in the right panel of Fig. 1 and the corresponding dark matter halo fitted parameters are listed in the first column of Tab. 2. Comparing these results with the values of Tab. 1 we can see that the error bars of all the dark matter parameters broaden when baryons are allowed to change.

In addition, we noticed that the resulting posterior distributions of the baryonic parameters are dominated by their priors. This is an expected behaviour for two reasons: a) the constraining power of $v_{\rm c}$ is not enough to fit all the baryonic parameters together with those of the dark matter, and b) the baryonic contribution to $v_{\rm c}$, for most of the range of $R$ values covered by the data, is smaller than the contribution of the dark matter component, as can be realized from the circular velocity curve figures. We also checked that doubling the uncertainties on the baryonic parameters does not make a significant difference in the $\rho_{{\rm DM},\odot}$ estimate, and neither to its uncertainty. However, when we imposed flat priors on our baryonic model, we found that the analysis preferred to fit the data with the dark matter distribution alone, ending up making the contribution from baryons too small to be realistic. This already points towards the dependence of the estimated $\rho_{{\rm DM},\odot}$ on the assumed baryonic distribution.

3.1.3 Changing the dark matter profile

We now study the effect that changing the dark matter profile has on the inferred value of $\rho_{{\rm DM},\odot}$. We use the same varying baryonic parameters for model B1 as in section 3.1.2; however, for the dark matter profile, we now generalize from standard NFW to ${\rm gNFW}$ with $\gamma$ being a free parameter. In other words we now let the inner slope $\gamma$ of the $\rho_{{\rm gNFW}}$ dark matter halo (Eq. (7)) become a free parameter of the mass model. The values of the dark matter fitted and derived quantities are shown in Tab. 2. One can see that the resulting value of $\rho_{{\rm DM},\odot}$ is essentially the same as before, despite the freedom in the dark matter profile. However, the additional inclusion of another free parameter in general broadens the obtained uncertainties of the velocity curves. The standard NFW profile with $\gamma=1$ is well contained within the 68% credible region when $\gamma$ is allowed to change.

On the other hand, the values of the baryonic parameters are not affected by this change in the dark matter profile, since they are still prior dominated (see the discussion in section 3.1.2).

On a similar basis, running the analysis by implementing an Einasto profile to model the dark matter distribution (Eq. (8)) also has a negligible effect on the estimate of $\rho_{{\rm DM},\odot}$ (see Tab. 2). Since this profile is naturally less peaked at the centre, the virial mass is reduced by about a factor of two with respect to the standard and generalized NFW studies. However, their values are compatible within their uncertainties.

In Fig. 7 of the appendix A we show the plot including the 2D marginal credible intervals for the parameters of Tab. 2.

3.2 Analysis of model B2

We now repeat the study presented in the previous section, although this time modelling the baryonic components with those of the model B2 (see section 2.1.2 for the details). The fact that the resulting $\rho_{{\rm DM},\odot}$ is driven by our choice of the baryonic model will become clearer in this section.

We again allow baryons to change with Gaussian priors as described before (see section 3.1.2), with the exception of the parameters for the gas distribution, for which we consider standard deviations corresponding to 25% and 20% from the central values of their masses and characteristic scales, respectively. The reason for this broadening is that these components are distributed in a clumpier way in the Galaxy, so the axisymmetric profiles provide a worse fit to the COBE maps in the analysis of Misiriotis:2006qq .

As we did in our study of the baryonic model B1, we start by assuming an NFW dark halo in section 3.2.1 and change it later in section 3.2.2.

3.2.1 NFW dark matter profile

The first dark halo that we study is the standard NFW halo (${\rm gNFW}$ with $\gamma=1$). The resulting values of the dark matter parameters are presented in Tab. 3 (first column of results) and the associated circular velocity curves are shown in Fig. 3.

Partly because the total baryonic mass of the B1 model is larger than the baryonic mass of B2, the contribution to the circular velocity curve of baryons from the B2 model is smaller than the contribution of baryons from the B1 model. Therefore, a slightly more massive and concentrated dark halo is needed to fit the data well (compare the values of Tabs. 2 and 3). As a result, the estimated value of $\rho_{{\rm DM},\odot}$ in the analysis of the baryonic model B2 is about 30% larger than in the analysis of the B1 model.

In principle, even though the mass of the stellar disk of model B2 is smaller than the mass of each of the two disks of model B1, the contribution to the circular velocity curve coming from the disk of model B2 is larger than the individual contributions of the baryonic disks of model B1, as can be appreciated comparing Figs. 1 and 3. The difference is due to the larger density in the Galactic plane (at a height $z=0$) of the double exponential disk from the model B2, compared to the same density from the Miyamoto-Nagai disks of the B1 model.

The different choice of the bulge profile in the two baryonic models affects mostly the circular velocity at distances smaller than $5{\rm\,kpc}$ from the Galactic centre, having little effect on the estimated value of $\rho_{{\rm DM},\odot}$. Furthermore, observationally we can not associate properly a value of the circular velocity to measurements in the region $R\lesssim 5{\rm\,kpc}$ (Chemin:1504.01507, ). In addition, it is known that the bulge is not spherical at such distances (see e.g. Portail:1502.00633 ), so a spherical approximation can be used as long as we do not include data at $R\lesssim 5{\rm\,kpc}$.

3.2.2 Changing the dark matter profile

As we did when analyzing the model B1, we test how robust the estimated value of $\rho_{{\rm DM},\odot}$ is to the change in shape of the dark halo. We analyze both $\rho_{{\rm gNFW}}$ and $\rho_{\rm Ein}$. The estimated quantities are presented in Tab. 3 and the corresponding circular velocity curves are similar to those obtained for an NFW (see Fig. 3) and are presented in Fig. 4.

Using a ${\rm gNFW}$ distribution, with free inner slope $\gamma$, we find a consistent measurement for all the parameters with respect to the case with $\gamma=1$ (the usual NFW profile). However, there is a noticeable increase in the error bars because of the presence of an additional free parameter (the inner slope $\gamma$). On the other hand, an Einasto profile produces a smaller virial mass $M_{200}$ and virial radius $r_{200}$. The maximum of their posteriors are nonetheless fully compatible with the studies using an NFW or a ${\rm gNFW}$ halo.

The plot including the 2D marginal credible intervals for the parameters of Tab. 3 is shown in Fig. 8 of the appendix A.

Similarly to what we obtained when analyzing the baryonic model B1, the estimate of $\rho_{{\rm DM},\odot}$ from studying the model B2 is robust under the change of the dark halo, although it presents a larger variation in the maximum and credible regions of its posterior.

The main source of change in the estimate of $\rho_{{\rm DM},\odot}$ is then the different baryonic distributions considered in the model B1 (analyses of section 3.1) and B2 (analyses of this section 3.2).

4 Discussion

In order to understand better the consistency of the results obtained in our analyses, we compare them with other estimates in the literature, both regarding $\rho_{{\rm DM},\odot}$ and other derived results. First, we comment on the local dark matter density in section 4.1. Then we test in section 4.2 the constraining power of the rotation curve method to a change in a disk structure. Since $N$-body simulations provide values for the $c_{200}\text{-}M_{200}$ relation, in section 4.3 we compare the 2D posterior distribution from our analyses with the $c_{200}\text{-}M_{200}$ trend found in simulations. Finally, we compare the estimated mass of the Milky Way from our mass models with measurements at different distances from the Galactic centre in section 4.4.

4.1 Estimated local dark matter density

Our study shows that the estimated dark matter parameters are affected by the choice of the distribution of baryons. In particular, as can be seen in Fig. 5 (where we present the marginalized posterior distributions of $\rho_{{\rm DM},\odot}$ for the different mass models analyzed), and in Tabs. 2 and 3, we find a larger value of $\rho_{{\rm DM},\odot}$ in the case of model B2 (that comprises double exponential disk profiles) than in the case of B1 (that is based on Miyamoto-Nagai disk profiles). This happens not only because of the larger mass of the B1 model, which contains around 50% more baryonic mass than the model B2, but also because of the different shape assumed in the components, especially the disks. Assuming the same mass, the double exponential stellar disk of model B2 contributes more to the circular velocity curve at intermediate-to-large Galactocentric distances $(\text{3--10}{\rm\,kpc})$ than the Miyamoto-Nagai disks of the model B1. However, the presence of two disks in the B1 model (with different scales but same mass) enhances the contribution of baryons to the circular velocity curve with respect to model B2. Therefore, less dark matter is needed in order to fit the measured circular velocities in the case of the baryonic model B1, and the resulting $\rho_{{\rm DM},\odot}$ value is $\sim 30\%$ smaller than its estimate from the analysis of the model B2.

The difference in the two baryonic models can also be seen comparing the distance at which the dark matter contribution to the circular velocity curve becomes larger than the contribution of baryons. In particular, we obtain that for the analysis using the model B1 this distance corresponds to $R_{v_{\rm c},{\rm eq}}\sim\text{10--16}{\rm\,kpc}$, where the values are the minimum and maximum of the marginal 68% posteriors for the three dark matter profiles (NFW, ${\rm gNFW}$ and Einasto). Similarly, the dark matter mass enclosed in a sphere of radius $r$ overcomes the baryonic mass inside the sphere at a radius $r_{M,{\rm eq}}\sim\text{8--14}{\rm\,kpc}$. On the other hand, when a baryonic B2 model is used, the distances at which dark matter contributes more than baryons are $R_{v_{\rm c},{\rm eq}}\sim\text{4--11}{\rm\,kpc}$ and $r_{M,{\rm eq}}\sim\text{3--9}{\rm\,kpc}$. These distances are closer to the Galactic centre when the model B2 is used, which is consistent with a larger $\rho_{{\rm DM},\odot}$ with respect to the model B1.

Regarding a change in the dark matter profile, we find that choosing a spherical NFW, ${\rm gNFW}$ or Einasto halo does not greatly affect the $\rho_{{\rm DM},\odot}$ estimates, given a baryonic model, as can be inferred from Fig. 5.

The obtained values of $\rho_{{\rm DM},\odot}$ from the analyses of the B1 and the B2 baryonic models are around $\rho_{{\rm DM},\odot}\sim 0.30\,\mathrm{GeV/cm}^{3}$ and $\rho_{{\rm DM},\odot}\sim 0.38\,\mathrm{GeV/cm}^{3}$, respectively. These values are consistent with most of the previous studies (some recent examples are McKee:2015hwa ; Xia:2015agz ; Pato:2015dua ; Sivertsson:2017rkp ; Karukes:2019jxv ; Benito:2019ngh ), but smaller than others (such as the recent estimates of Buch:2018qdr ; Widmark:2018ylf ). One explanation could be the different method used to estimate $\rho_{{\rm DM},\odot}$, since Buch:2018qdr ; Widmark:2018ylf exploit the local $z$-Jeans equation method, which is based on the vertical movement of stars in a region close to the Solar System (see e.g. Binney:2008-book ; Read:2014qva ). This method is likely to overestimate $\rho_{{\rm DM},\odot}$ when the assumption of equilibrium breaks down. Observations of asymmetries in the densities and velocities of stars in the disk show that the disk is likely experiencing vertical oscillations, probably because of the passage of a massive satellite (see e.g. Antoja:1804.10196 ; Heines:1903.00607 ). In addition, the Solar System could be placed at an overdense region, compared to the surroundings, and it could also be possible that the baryonic models for gas and dust are incorrect. We further comment on this possibility in the next section.

4.2 Sensitivity to disk uncertainties

In the previous section we discussed that the main difference between the best-fit values of $\rho_{{\rm DM},\odot}$ that we obtain from the analyses of the two baryonic models, B1 and B2, arise from the difference in their disk structure. Thus, the estimate of $\rho_{{\rm DM},\odot}$ using the rotation curve method is particularly sensitive to differences in the shape of the Galactic disk. In this section we study by how much the determined value of $\rho_{{\rm DM},\odot}$ could be affected by disk uncertainties.

In principle, disk uncertainties could move $\rho_{{\rm DM},\odot}$ either up or down. However, it is rather difficult to obtain a smaller value of $\rho_{{\rm DM},\odot}$, because it would imply that the distribution of baryons decline even less rapidly, both vertically and in radius, than in the Miyamoto-Nagai disks of the model B1. Such a spread distribution seems to be incompatible with other observations (Rix:2013bi, ; Bland-Hawthorn:1602.07702, ). On the contrary, an additional thin disk distribution would increase $\rho_{{\rm DM},\odot}$ and could explain the rather large values measured by Buch:2018qdr ; Widmark:2018ylf . Remember, however, that these authors used the $z$-Jeans equation method in order to infer $\rho_{{\rm DM},\odot}$, and it is known that local methods can be affected by disequilibria Banik:2016yqm .

If the overdensity found in local studies is not a feature of disequilibrium, it could be argued that the excess is produced by a dark matter disk (motivated by beyond the Standard Model physics, Read:2008fh ; Purcell:2009yp ; Fan:2013yva ), while another possible explanation considers that it is made of baryonic gas Widmark:2018ylf . Independently of the reason for such an excess, we can assume for a moment that there really is an additional, extremely thin component in the Galaxy, that will increase $\rho_{{\rm DM},\odot}$, and study whether our analysis using the rotation curve method can determine the presence of the new component.

Therefore, we performed a simple test. We used our baryonic B1 model configuration with varying baryons plus an NFW spherical halo. To these components we added an extremely thin double exponential disk, with fixed parameters $R_{d}=3.93{\rm\,kpc}$, $z_{d}=0.038{\rm\,kpc}$ and $M=7.3\times 10^{9}{\rm\,M_{\odot}}$, that increases $\rho_{{\rm DM},\odot}$ by an amount compatible with the local observations of Buch:2018qdr ; Widmark:2018ylf . This additional component contributes to the circular velocity curve of the Galaxy in a similar way as gas from the B2 model (see e.g. Fig. 3). We find that the varied parameters of the fit adapt very well, acquiring values consistent with the analysis without the extra thin disk. Thus, the rotation curve method is unable to notice the presence of such a thin component. However, the new component adds some extra energy density, in particular close to the Galactic plane. If we interpret the energy density of the new thin disk as a dark, non visible source, we should add an extra $\Delta\rho_{{\rm DM},\odot}=\{4.7,2.4,0.3\}{\rm\,GeV/cm^{3}}$ for a distance from the Galactic plane of $z=\{0,25,100\}\,\mathrm{pc}$. Since the Earth is at around $z_{\odot}\approx 25\,\mathrm{pc}$ from the plane (Bland-Hawthorn:1602.07702, ), this would imply a local dark matter density 8 times larger than the one obtained in the same analysis without the extra thin component.

We insist that we are not arguing in favor of the existence of the additional thin disk. But it is important to show that the rotation curve method is insensitive to the existence of such a component, and the obtained values of $\rho_{{\rm DM},\odot}$ from our analyses might be smaller than the real $\rho_{{\rm DM},\odot}$ value if there was indeed a new component whose distribution is too spread out to be noticed from the analysis of the Milky Way’s rotation curve.

4.3 Studying the $c_{200}\text{-}M_{200}$ relation

In order to investigate further the consistency with other studies, here we explore the relation between the two main virial quantities, $c_{200}$ and $M_{200}$. Note, however, that the relation that we want to study in this section is independent of the correlation found in our analyses between the two fitted dark matter quantities; such correlation between $c_{200}$ and $M_{200}$, depicted in Fig. 2, is a result of the way that both parameters affect the fit of the dark matter distribution to the data. Given the values of $v_{\rm c}$ determined in Eilers:1810.09466 and for a fixed baryonic mass distribution, a more concentrated dark halo (with larger $c_{200}$) needs to have a smaller virial mass $M_{200}$ in order to fit the data, as shown in Fig. 2. However, what we are going to study in this section is the comparison of the relation between the two virial parameters, $c_{200}$ and $M_{200}$, found in our analyses with respect to the same relation obtained in $N$-body simulations (see for instance Dutton:2014xda ; Klypin:2014kpa ; Bose:2019hwt ).

Note that in our analyses we have used flat priors on the values of the two main virial quantities, $c_{200}$ and $M_{200}$. However, $N$-body simulations reveal a characteristic relation between these two quantities (which we reiterate is different than the relation shown in Fig 2), that can be expressed as

[TABLE]

where $a=0.520$, $b=-0.101$ for an NFW halo, and $a=0.459$, $b=-0.130$ for an Einasto profile (Dutton:2014xda, ).

In Fig. 6 we show the comparison between our analyses and $N$-body simulations. In the left panel the fitted dark matter halo is an NFW (or a ${\rm gNFW}$) and in the right panel we show the same when an Einasto profile is considered.

All marginal credible regions shown in the left panel are perfectly centered at the $N$-body trend, which assumes an NFW shape of the halo. The studies of the baryonic model B1 prefer a slightly less massive halo and the 2D posteriors of our analyses are wider for the ${\rm gNFW}$ profile since it has an extra free parameter.

For those cases that include an Einasto profile the comparison is less attractive, but the result of our analyses and the $N$-body trend are still fully compatible. As shown in Tabs. 2 and 3, the preferred $M_{200}$ value is smaller for an Einasto halo than for a ${\rm gNFW}$. However, an Einasto profile fits more loosely to data and the 68% and 95% posteriors extend beyond those found with a ${\rm gNFW}$, except for the largest allowed values of $c_{200}$ (see Fig. 6).

Overall, our results respect what is found in $N$-body simulations, indicating that our Galactic model behaves as expected. However, it is worth mentioning that the $c_{200}\text{-}M_{200}$ relations of Eq. (14) were obtained from dark matter only simulations, and the lines in Fig. 6 from Dutton:2014xda could change if baryons are included (see e.g. Sawala:2015cdf ; Kelso:2016qqj ).

4.4 Estimated Galactic mass

Another quantity that we can compare with observations is the total dynamical mass of the Galaxy at different distances from the Galactic centre. In Tab. 4 we present the values obtained in our main studies, with varying baryonic parameters and flat priors on the virial quantities, and compare them with a selection of observational estimates. As can be seen in the table, our results are compatible with other studies in a wide range of values; however, since the reach of the data that we have used only covers up to $\sim 25{\rm\,kpc}$, the error of our estimated masses grow at larger distances, making our estimate at the largest distance of $300{\rm\,kpc}$ comparatively smaller than the value found by Watkins:2010fe .

This difference is another indication that the constraining power of circular velocities is limited. Although current well measured $v_{\rm c}$ data have increased their reach up to $25{\rm\,kpc}$, this distance is not enough to put robust constraints on the virial mass and concentration parameters of the halo, since its real shape, in particular beyond $25{\rm\,kpc}$, could be very different from the assumed profiles of this work. However, we decided to present the estimated Galactic mass up to such large distances to facilitate the comparison of the properties of our fitted mass models with those that can be found in the literature.

5 Conclusions

Our goal in this work was both to estimate the value of $\rho_{{\rm DM},\odot}$ from the rotation curve of the Galaxy and to study the robustness of the determination. This last point was addressed using different Galactic mass models, with different dark matter and baryonic density distributions. In order to estimate $\rho_{{\rm DM},\odot}$, we fitted the parameters of the mass models to the precise measurements (using Gaia DR2 data) of the Milky Way’s circular velocity presented in Eilers:1810.09466 .

We examined two baryonic models, referred to in the text as B1 and B2 (see section 2 for the description of the models), and tested three different dark matter spherical halos (NFW, ${\rm gNFW}$ and Einasto). The main difference between the baryonic models stands in the fact that B1 is denser than the model B2 in the intermediate Galactic radii regions and out of the Galactic plane.

Given a baryonic mass model, we found the value of $\rho_{{\rm DM},\odot}$ to be robust under the change of the spherical dark matter halo. However, $\rho_{{\rm DM},\odot}$ was more sensitive to a change in the baryonic model. In particular, in the analyses of the model B1 we found a value of $\rho_{{\rm DM},\odot}\simeq 0.30{\rm\,GeV/cm^{3}}$ (see Tab. 2), while $\rho_{{\rm DM},\odot}\simeq 0.38{\rm\,GeV/cm^{3}}$ was found instead for the baryonic model B2 (see Tab. 3). Notice that the $68\%$ region for the baryonic model B1 goes up to $\rho_{{\rm DM},\odot}\approx 0.33{\rm\,GeV/cm^{3}}$, while the $68\%$ region for the B2 model goes down to $\rho_{{\rm DM},\odot}\approx 0.35{\rm\,GeV/cm^{3}}$. One can see that the regions do not overlap at the $68\%$ level, but they do at $95\%$. Thus the baryonic models are not totally incompatible.

Although the difference between models B1 and B2 is not extremely large, it is important, in particular given that usually people only cite the error bars for one baryonic model. Within one single model the largest $1\sigma$ error bars on $\rho_{{\rm DM},\odot}$ correspond to $\pm 0.036{\rm\,GeV/cm^{3}}$ (for the case of baryonic model B2 and Einasto dark matter profile), see Tab. 3. On the other hand, when we take into account the variety of mass models considered in this paper, in particular the uncertainty in the baryonic profile by considering both models B1 and B2, the $1\sigma$ uncertainty range reaches $0.149{\rm\,GeV/cm^{3}}$. This error is at least twice that from any of the individual mass models.

We also tested the robustness of the rotation curve method to estimate $\rho_{{\rm DM},\odot}$ assuming the existence of an hypothetical thin disk, together with an NFW dark halo and the baryonic components of the B1 model. In this scenario the fit to the $v_{\rm c}$ data is not able to perfectly distinguish the presence of such thin disk, since it barely contributes to the Galactic circular velocity curve. Thus, uncertainties in the Galactic disk can significantly increase the uncertainty in the estimate of $\rho_{{\rm DM},\odot}$, that could be more than a factor 8 larger if an additional thin disk is present in the Milky Way.

Furthermore, we compared the resulting virial quantities of our analyses with the $c_{200}\text{-}M_{200}$ trend from $N$-body simulations, finding them to be consistent. We additionally compared the total Galactic mass resulting from our analyses with respect to observational estimates at different distances from the Galactic centre. This comparison shows that the constraining power of circular velocities is limited. The reach of current $v_{\rm c}$ measurements have extended its precision up to $25{\rm\,kpc}$, but this distance is still not enough to put robust constraints on the parameters (virial mass and concentration) of the halo. Indeed the halo shape could be different from the assumed profiles of this work, in particular at larger distances from the Galactic centre.

Overall, we proved that the estimate of $\rho_{{\rm DM},\odot}$ from Milky Way’s rotation curve measurements depends on the assumed shape for the mass distribution, regardless of the precision of the $v_{\rm c}$ measurements.

In order to constrain $\rho_{{\rm DM},\odot}$ in a more robust way, more data apart from circular velocities need to be taken into account, like phase-space information of different type of tracers such as halo stars that travel far from the disk plane. We plan to do that in future works, also combining different methods to that of the rotation curve (for instance we can explore a combination with local methods, such as the vertical Jeans equation method) in a search for an estimate of $\rho_{{\rm DM},\odot}$ as precise and consistent as possible.

Acknowledgements.

PFdS thanks Justin Alsing for useful discussions about our Bayesian analyses, and Sofia Sivertsson, Eric F. Bell and Oleg Gnedin for their comments and discussions. PFdS, KM and KF acknowledge support by the Vetenskapsrådet (Swedish Research Council) through contract No. 638-2013-8993 and the Oskar Klein Centre for Cosmoparticle Physics. KF acknowledges support from DoE grant DE-SC007859 and the LCTP at the University of Michigan. MV and KH are supported by NASA-ATP award NNX15AK79G to the University of Michigan. PFdS thanks the LCTP at the University of Michigan for the hospitality received while this work was been finalized.

Appendix A Triangular plots of the dark matter parameters

In this appendix we present the triangular plots corresponding to the 2D marginal 68% and 95% credible regions of the dark matter fitted and derived parameters of our analyses. Figure 7 shows the triangular plot for the analyses including the baryonic model B1 (section 3.1). Figure 8 shows the triangular plot for the analyses including the baryonic model B2 (section 3.2).