© Y. Jiao et al. 2021

1. Introduction

Gaia DR2 provided accurate stellar proper motions to calculate the circular velocity curve of the Milky Way (MW) up to 25 kpc (Eilers et al. 2019; Mróz et al. 2019). The result was based on a thorough analysis of a very large sample of 26 000 RGB stars in the MW disk (Eilers et al. 2019), resulting in a slightly but robustly determined decrease in circular velocity from 5 to 25 kpc. While Eilers et al. (2019, see also Hogg et al. 2019) used spectrophotometric distances in their analysis, their finding was confirmed by Mróz et al. (2019) using 773 Classical Cepheids with precise distances. Subsequent analyses of these rotation curves (RCs) have led to a total MW mass near or well below 10¹² M_⊙ (Eilers et al. 2019; de Salas et al. 2019; Grand et al. 2019; Karukes et al. 2020). Karukes et al. (2020) have used a considerable number of baryonic matter distributions to derive the overall mass distribution, while the de Salas et al. (2019) have accounted for very large error bars after cumulating all the systematics described in details by Eilers et al. (2019).

The accuracy of the MW RC also allows testing different mass profiles for the dark matter (DM) distribution in the MW halo. Recent studies have shown that the three-parameter Einasto profile (Einasto 1965, see also Retana-Montenegro et al. 2012) provides a better description of the CDM halo density profile than the NFW profile (Navarro et al. 2004, 2010; Gao et al. 2008), and it is even than the three-parameter generalized gNFW (Klypin et al. 2016).

We propose to test the Einasto and NFW (Navarro et al. 1997) density profiles and their effect on the total mass estimates when spiral rotation curves are fit. We consider the MW RC because of its unprecedented accuracy, and also because the history of the MW is likely quiescent when compared to other spirals (Hammer et al. 2007) because the last MW major merger occurred ∼10 Gyr ago, as has recently been confirmed based on the resulting debris identified by Gaia DR1 (Belokurov et al. 2018) and as will soon be confirmed by Gaia DR2 (Haywood et al. 2018; Helmi et al. 2018).

In Sect. 2 we present our proposed treatment of the error bars for the Eilers et al. (2019) RC, and then describe the choice and mathematical descriptions of the baryon and DM models. In Sect. 3 we compare the χ² probability distribution for DM represented by the NFW or Einasto profiles. In Sect. 4 we discuss which mass range is consistent with the combined constraints provided by the fit of the MW RC and by adopting DM halo profiles from the cold dark matter (CDM) theory.

2. Methods

2.1. Rotation curve and error bars

Eilers et al. (2019) provided a thorough analysis of the possible systematic errors that may affect the MW RC and summarized (see their Fig. 4) four different types of systematics. The first type includes the neglected term in their Jeans equation (see their Eq. (3)), which is a cross-term made by the vertical density gradient of the product of the radial and vertical velocities. This term is found to be small but not negligible at large distances. For example, Mackereth et al. (2019) showed that vertical velocities are higher for young stars, which is expected because the gaseous disk is likely affected by (former) gas infall. This may affect the derived RC because Eilers et al. (2019) selected relatively young stars (< 4 Gyr) for the MW RC in order to avoid asymmetric drift effects.

However, the effect is expected to be smaller (< 5 km s⁻¹ at 12 kpc) than the RC amplitude. The second possible systematics is empirical, and it is an estimate of the error variations with radius after splitting the sample into two parts. We consider here only the first type of systematics because it likely includes the second.

Adding to this, Eilers et al. (2019) considered a third category of systematics with a quite different nature because it proportionally applies in the same way to all RC points. It is revealed by the three almost horizontal lines in Fig. 4 of Eilers et al. (2019). This last category of systematics includes the effect of changing the distance of the Sun to the Galactic center, the proper motion of the latter, and it can be extended to the change in scale length. These uncertainties have to be applied to the derived mass as a whole after the fitting analysis. Added together, they correspond to an additional systematic uncertainty of ∼2% on the velocity scale and ∼4% on the mass scale. We note in agreement with Christina Eilers (Eilers, priv. comm.) that summing all the errors of Fig. 4 of Eilers et al. (2019) (as it has been done by de Salas et al. 2019) would strongly overestimate the error bars (see the discussion above), which dilutes the significance of the MW RC.

In the following we adopt the same parameters for the position of the Sun and for the solar velocity as Eilers et al. (2019). Karukes et al. (2020) have shown that the choice of the solar velocity may significantly affect the determination of its mass, while it has been considered determined at a 2−3% level by Eilers et al. (2019).

2.2. Milky Way baryonic mass models

The contribution of the baryonic components to the MW mass or RC is still uncertain, and this may well affect the determination of the DM distribution. Following Karukes et al. (2020), we adopt here a large number of models from the literature to describe the MW baryonic component, as described below. The baryonic component and its distribution in the bulge, disk, thick disk, gas, and even halo gas is still debated (see the review by Bland-Hawthorn & Gerhard 2016), and some modeling also introduces an ionized gas component (Cautun et al. 2020). The basic idea is to cope with uncertainties on baryons by using a very large grid of possible models, although we are aware that some baryons models may not be fully consistent with other important constraints from vertical dynamics of the disk stars (Bovy & Rix 2013) or from microlensing (Wegg et al. 2016).

Pouliasis et al. (2017) generated a new axisymmetric model (Model I) including a spherical bulge and a thin and thick disk. This model satisfies a number of observational constraints: stellar densities at the solar vicinity, thin- and thick-disk scale lengths and heights, and the absolute value of the perpendicular force K_z as a function of distance to the Galactic center. Although the disk is made of a thin and a thick disk, the associated density profiles are both described by a Miyamoto–Nagai profile (Eq. (1)). Pouliasis et al. (2017) concluded that Model I supersedes the axisymmetric model (Model A&S) proposed by Allen & Santillan (1991) because there is growing evidence for a strong thick-disk component and because the bulge is less prominent and less classical than assumed in Model A&S. Model A&S consists of a stellar thin disk with a Miyamoto–Nagai profile (Miyamoto & Nagai 1975) and a central bulge with a Plummer profile (Binney & Tremaine 2011). The description of the bulge and disks for both Model I and Model A&S is expressed in the form (Pouliasis et al. 2017) for (R, z) cylindrical coordinates,

$$\begin{array}{cc}& {\rho }_{\mathrm{thin}}(R,z)=\frac{b_{\mathrm{thin}}^2{M}_{\mathrm{thin}}}{4\pi }\hfill \\ \hfill & \times \frac{(R^2{a}_{\mathrm{thin}}+3{(z^2+b_{\mathrm{thin}}^2)}^{1/2}){(a_{\mathrm{thin}}+{(z^2+b_{\mathrm{thin}}^2)}^{1/2})}^2}{{(R^2+{[a_{\mathrm{thin}}+{(z^2+b_{\mathrm{thin}}^2)}^{1/2}]}^2)}^{5/2}{(z^2+b_{\mathrm{thin}}^2)}^{3/2}}\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\rho }_{\mathrm{bulge}}(r)=\frac{3b_{\mathrm{bulge}}^2{M}_{\mathrm{bulge}}}{4\pi {(r^2+b_{\mathrm{bulge}}^2)}^{5/2}},\hfill \end{array}$$(2)

where $r=\sqrt{R^2+z^2}$, and M_thin, M_thick, M_bulge, a_thin, a_thick, b_thin, b_thick, b_bulge are the disks and bulge mass and scale constants, respectively (see Table 1).

Sofue (2015) presented a model (Model S) of the MW by attempting to fit a ‘grand rotation curve’, which defines the combination of the actual rotation curves (up to 20−25 kpc) with estimates based on orbital motions of objects beyond 25 kpc in the MW halo, e.g., distant globular clusters. The bulge was approximated by a de Vaucouleurs profile (de Vaucouleurs 1958). We chose to adopt a Plummer profile (Eq. (2)) for the bulge, and the disk was assumed to follow an exponentially thin density profile. The surface mass density of the disk is expressed as (Sofue 2015)

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{d}}(R)={\mathrm{\Sigma }}_{0}exp(-R/a_{\mathrm{thin}}),\end{array}$$(3)

where Σ₀ is the central value and a_thin is the scale radius (see Table 1). This model provides the highest baryonic mass when compared to other models in the literature (see Fig. 1). Nevertheless, we consider it useful for testing the effect of an extremely high baryonic mass for the MW disk and bulge.

Fig. 1. Top: contribution to the rotation curve of different baryonic models and model components. Red points indicate the rotation curve of the Milky Way from Eilers et al. (2019). The error bars are estimated by bootstrapping and include the systematic uncertainties from the neglected term (see text). Bottom: fit of the rotation curve by the best-fit model (solid blue curve, total mass of 2.8 × 1011 M⊙), and with the most massive MW model for which the χ2 probability reaches P = 0.05 (dash-dotted orange line, total mass of 18 × 1011 M⊙), both associated with the baryonic distribution from model I of Pouliasis et al. (2017).

A great addition to our choices of baryonic components was presented by Iocco et al. (2015), and they allowed several possible combinations of models for the bulge and the disk. For the bulge we chose the two triaxial mass density distributions E2 and G2 presented by Stanek et al. (1997),

$$\begin{array}{cc}& \mathrm{E}2:{\rho }_{\mathrm{bulge}}(x,y,z)={\rho }_0{e}^{-r_1}\hfill \end{array}$$(4)

$$\begin{array}{cc}& \mathrm{G}2:{\rho }_{\mathrm{bulge}}(x,y,z)={\rho }_0{e}^{-r_2^2/2},\hfill \end{array}$$(5)

with

$$\begin{array}{c}\hfill r_1^2=\frac{x^2}{x_{\mathrm{b}}^2}+\frac{y^2}{y_{\mathrm{b}}^2}+\frac{z^2}{z_{\mathrm{b}}^2},r_2^4={(\frac{x^2}{x_{\mathrm{b}}^2}+\frac{y^2}{y_{\mathrm{b}}^2})}^2+\frac{z^4}{z_{\mathrm{b}}^4},\end{array}$$(6)

where (x, y, z) are the coordinates along the major, intermediate, and minor axes. For the thin and thick disks, we adopted a double exponential of the three models (CM from Calchi Novati & Mancini 2011, dJ from de Jong et al. 2010 and J from Jurić et al. 2008) as described below (see Table 3):

$$\begin{array}{cc}& \mathrm{CM}:\rho (R,z)={\mathrm{\Sigma }}_{\mathrm{thin}}(\frac{1}{2H_1}exp(-\frac{R}{L_1}-\frac{z}{H_1})\hfill \\ \hfill & +f_{\mathrm{thick}}\frac{1}{2H_2}exp(-\frac{R}{L_2}-\frac{z}{H_2}))\hfill \end{array}$$(7)

$$\begin{array}{cc}& \mathrm{J},\mathrm{dJ}:\rho (R,z)={\rho }_{\mathrm{thin},\odot }(e^{R_{\odot }/L_1}exp(-\frac{R}{L_1}-\frac{z}{H_1})\hfill \\ \hfill & +f_{\mathrm{thick}}{e}^{R_{\odot }/L_2}exp(-\frac{R}{L_2}-\frac{z}{H_2}))·\hfill \end{array}$$(8)

In decreasing order of baryonic mass, Model S, Model A&S, and then Model I assume significantly higher mass baryonic components than the six combinations of bulge (G2, E2) and disk (CM, dJ, J), which is illustrated by Fig. 1 or by comparing Table 1 with Tables 2 and 3.

2.3. Milky Way dark matter models

We considered the NFW and Einasto profiles to describe the density profiles of DM halos in spherical coordinates (r). The generalized NFW profile (gNFW, see Zhao 1996) can be expressed as in de Salas et al. (2019),

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_0}{{(r/r_0)}^{\gamma }{(1+r/r_0)}^{3-\gamma }},\end{array}$$(9)

where r₀ is the scale radius, and ρ₀ is the characteristic dark matter density. For γ = 1, the profile becomes the NFW profile (Navarro et al. 1997) for which we investigate which parameters are able to fit the MW RC, after letting the two NFW parameters, r₀ and m_NFW = 4 πρ₀ $r_0^3$, vary from 2 to 100 kpc and from 1 to 50 × 10¹¹ M_⊙, respectively. For the gNFW profile we let the additional parameter, γ, vary from 0.1 to 3 (see also Karukes et al. 2020). For each tested mass configuration, we verified later that the investigated parameter space was sufficiently large to avoid having missed any solution.

Using the Retana-Montenegro et al. (2012) mathematical framework, the Einasto profile can be written as

$$\begin{array}{c}\hfill \rho (r)={\rho }_0\exp [-{\left(\frac{r}{h}\right)}^{1/n}],\end{array}$$(10)

where n can determine how fast the density decreases with r. To determine which models are able to fit the MW RC, we let the three Einasto parameters, b_E = 3 × n, h_red = h^1/n, and m_E = 4 πρ₀h³nΓ(b_E), vary from 3 to 30, 0.05 to 3 and from 1 to 50 × 10¹¹ M_⊙, respectively. For each tested mass configuration, we verified later that the investigated parameter space was sufficiently large to avoid any missing solution.

In order to determine a non-indefinite total MW mass, the DM halo mass has to be limited by the virial radius, R_vir, which enclosed M_vir, which is the virial mass. We define the virial radius as the radius of the sphere for which the average dark matter density equals 200 times the critical density of the Universe ρ_cr. We adopted a critical density of ρ_cr = 1.34 × 10⁻⁷ M_⊙ pc⁻³, which comes from Hinshaw et al. (2013). With this definition, the relation between virial radius and virial mass is

$$\begin{array}{c}\hfill M_{\mathrm{vir}}=200\times \frac{4\pi }{3}{\rho }_{\mathrm{cr}}{R}_{\mathrm{vir}}^3.\end{array}$$(11)

3. Results

3.1. Deriving the total MW mass and χ² probability

The total MW potential can be obtained through the Poisson equation,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2{\mathrm{\Phi }}_{\mathrm{tot}}=4\pi G{\displaystyle \sum _i{\rho }_i,}\end{array}$$(12)

after adding all the different MW mass components. The theoretical estimate of the circular velocity is derived at different disk radii (R) from the potential Φ_tot of the Galaxy through

$$\begin{array}{c}\hfill v_{\mathrm{c}}^2(R)=R\frac{\partial {\mathrm{\Phi }}_{\mathrm{tot}}}{\partial R}{|}_{z\approx 0}.\end{array}$$(13)

We applied the χ² method to fit the RC and calculate its associated probability, for which we tested an extremely large parameter space. The χ² was calculated by the sum at each disk radius R_i,

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i^N\frac{{(v_{\mathrm{mod},i}-v_{\mathrm{obs},i})}^2}{{\sigma }_i^2},}\end{array}$$(14)

where v_mod is the modeled circular velocity for the cumulative baryons + DM profiles, v_obs is the observed circular velocity, and σ_stat is the statistical uncertainty of the measurement so that ${\sigma }_{\text{stat},i}=({\sigma }_{v_{\text{obs},i}}^{+}+{\sigma }_{v_{\text{obs},i}}^{-})/2$, to which we added the systematic uncertainty σ_sys, i to calculate σ_i (see Sect. 2.1 and the table in Appendix A). Hence the χ² probability can be expressed as

$$\begin{array}{c}\hfill \mathrm{Prob}(\frac{{\chi }^2}{2},\frac{N-\nu }{2})=\frac{\gamma (\frac{N-\nu }{2},\frac{{\chi }^2}{2})}{\mathrm{\Gamma }\left(\frac{N-\nu }{2}\right)},\end{array}$$(15)

where N is the number of independent observed velocity points in the Eilers et al. (2019) RC, and ν is the number of degrees of freedom.

To fit the MW RC, we investigated a very large parameter space, allowing for the total MW mass from 1 to 50 × 10¹¹ M_⊙, for instance. In Fig. 2 each point (P(χ²), M_tot) represents an investigated baryon + DM model. The top panels of Fig. 2 present the χ² probability for the Model I baryon profile (Pouliasis et al. 2017) when associated with either the Einasto (left), the NFW (middle), or the gNFW (right) profiles. The first profile shows that high χ² probabilities are reached for low MW masses. In contrast, there is no similar trend for the NFW profile, which selects a narrow range of MW masses to fit the RC. The situation is improved with the gNFW, although it does not recover the whole range of masses and especially misses total masses below 5 × 10¹¹ M_⊙. The bottom panels present the same for model S, for which the probabilities are very low when associated either with the Einasto or the NFW DM profiles. Examination of the RC fit shows that the baryonic mass is so high that its radial profile is setting up most of the expected RC (see also Fig. 1), leading to differences with the observed RC at almost every radius. This is expected because model S is clearly at odds for the MW; its disk plus bulge mass is higher than that of M31, while half this value is more likely (see, e.g., Hammer et al. 2007).

Fig. 2. χ2 probability associated with the combination of different baryonic models with Einasto, NFW, or gNFW mass profiles for the DM. From top to bottom, we show Model I, bulge G2+disk J, bulge E2+disk J, and model S. The red points and error bars are the average and 1σ uncertainties.

The two panels in the middle row of Fig. 2 compare the results when the bulge is changed from G2 (top) to E2 (bottom), both added to disk J. The first shows a similar behavior as Model I associated with either Einasto or NFW DM profiles. When we used the Einasto profile for the DM, we found that increasing the bulge mass by 15% (from G2 to E2) is sufficient to exclude high values of the total mass of the MW. This is expected because when the baryonic mass is increased, a smaller amount of DM mass is available to reproduce the MW RC. Moreover, a too large bulge may limit the number of possible solutions that can fit the RC at low radii. However, for the NFW profile we find that a bulge mass increase from G2 to E2 is sufficient to prevent an efficient reproduction of the MW RC, providing very low χ² values. We also find that the associated total (and DM) masses are higher than the mass for the G2 bulge, which disagrees with our expectations. We note that these two properties disappear when the three-parameter gNFW model is used, which might be because for γ < 1, this profile is less cuspy and is therefore less affected by changes in bulge mass.

The above motivates us to investigate further why adding an additional baryonic mass could lead to an increase of the DM mass when the later is modeled by the NFW density profile. We tested the effect of changing the amount of baryonic mass on the NFW DM mass. We considered a range of baryonic masses scaled on the mass of Model I, with scale factors f varying from 0.85 to 1.15. For f = 0.85, 1, and 1.15, this confirmed that by increasing the baryonic mass, the NFW DM model leads to a significant increase in DM mass from 5.9, 7.2, and 10.7 × 10¹¹ M_⊙, respectively. This is an unexpected behavior because the DM role is to compensate for the lack of baryonic mass when a given RC is fit. Our first explanation was to relate this to the two-parameter nature of the NFW profile. However, a similar (although less pronounced) behavior affects the gNFW profile. For f = 0.85, 1, and 1.15, the gNFW DM model also leads to an increase in DM mass from 4.8, 5.25, and 7.1 × 10¹¹ M_⊙, respectively. This suggests the following mechanism: for an increasing baryonic mass, the NFW DM scale radius (r₀, see Eq. (9)) has to increase to dilute the DM mass from 5 to 25 kpc (the latest point of the RC). Because outer density slope of the NFW and gNFW is almost constant and shallow (−3) at large radii, this automatically leads to increasing DM masses. This indicates a possible methodological problem of using the NFW (or gNFW) to fit the RC as and to estimate the mass of a galaxy from it.

3.2. Systematics due to the NFW and gNFW when the total mass is estimated

To evaluate the differences between Einasto and NFW DM density profiles in fitting the MW RC, we need to ensure that our method does not depend on the initial conditions. In particular, the parameter grid might affect our results because Fig. 2 shows that the three-parameter space (Einasto or gNFW) might be more difficult to be populated than the two-parameter space (NFW). We further performed for each model a combination of several Monte Carlo simulations that also accounted for the variance due to the RC error bars, which are assumed to follow a Gaussian distribution, in order to fill the high-probability space in the (P(χ²), M_tot) plane as much as possible.

The solid lines in Fig. 3 identify the envelop for each baryonic + DM model, which is defined as being the highest χ² probability calculated in mass slices with sizes of 0.3 × 10¹¹ M_⊙. We assume that only χ² probabilities higher than 0.05 correspond to a good fit of the RC, which we verified after examining the latter. For comparison, Fig. 3 also shows the averaged probabilities.

Fig. 3. Maximum (solid lines) and averaged (dotted lines) χ2 probabilities for the different baryonic models. Model names are labeled in each panel, with Einasto, NFW, and gNFW mass predictions in green, magenta, and orange, respectively. The two panels associated with baryonic model E2+dJ and E2+J show no histogram for the NFW because this density profile fails to reproduce the MW RC. The horizontal dotted lines indicates the χ2 probability limit of 0.05 below which a model is found to be unable to fit the MW RC.

Figure 3 shows that for all baryonic models, a narrower range of total masses is found for the MW when an NFW or gNFW instead of a Einasto profile is adopted for the DM. Conversely, using the Einasto profile suffices to sample most of the points generated by the NFW profile in the (P(χ²), M_tot) plane. We find that the total mass solutions based on the NFW and gNFW profiles are often included in those from the Einasto profile, while using an NFW does not match the highest χ² probabilities found by the Einasto model (compare the peaks of the solid magenta and green lines). However, in the case of a massive bulge (E2, especially when associated with disk J), the three-parameter gNFW may sample total MW mass values that cannot be reached by the Einasto model.

Table 4 gives the estimated total masses based on the minimum χ² values (best fit, highest probabilities) or on averaging the χ² probabilities in each mass slice (average). As in Fig. 3, the rows are sorted from high- to low-mass baryonic models. This indicates that the best fit of the MW RC for all baryonic models, except for A&S, are unavoidably related to low total masses (from 2.3 to 3.3 × 10¹¹ M_⊙) if a Einasto profile is chosen for the DM. Conversely, adopting an NFW (or gNFW) profile for the DM leads to much higher total mass values by a factor of 2 to 4.

Of the models we studied, Model A&S possesses the second highest baryonic mass, close to 10¹¹ M_⊙, and we investigated why the behavior it shows is so different from that of other baryonic models, especially Model I. In addition to baryonic masses that differ by 11%, the main difference between the two models is the presence of a thick disk incorporating half the disk mass in Model I, with a scale length that is half that of the thin disk of Models I and A&S. By modifying the thick-disk scale length of Model I to a higher value, we find that this suffices to provide a similar behavior to Model A&S for the normalized cumulative probabilities of both NFW and Einasto DM profiles. As previously noted for model S, this suggests that an extended and relatively massive baryonic disk determines a significant part of the RC shape.

Considering the averaged total masses slightly improves the similarities between predictions based on NFW and Einasto DM mass profiles. This is true for Models A&S and I, which lead to almost consistent NFW and Einasto values of the total masses. However, for lighter baryonic models, the NFW profile for DM still leads to a mass that is higher by factors from 1.5 to 3 when compared to that resulting from the Einasto profile. The NFW (and to a lesser extent, the gNFW) profile appears to preferentially select a narrow range of total masses, excluding in particular the low-mass values that are favored when the Einasto profile is used.

4. Discussion

4.1. Limitations of this study and comparison with other works

The goal of this paper is mostly methodological, that is, we search for the range of total MW masses that reproduces the MW RC, and then evaluate which mass density profile is the most suitable for estimating the DM mass. We focus on the rotation curve provided by Gaia DR2 alone (Eilers et al. 2019; Mróz et al. 2019) because its accuracy is several times better than those of any former studies (see Fig. 3 of Eilers et al. 2019). This is also because disk stars correspond to dynamical points that are well anchored in the stellar disk, which is assumed to be well in equilibrium with the MW potential. In this context, our study broadens the recent work of de Salas et al. (2019) and Karukes et al. (2020) because here we consider a wider range of baryonic matter models of the MW to fit the Gaia DR2 RC¹. Our resulting total masses for the baryonic Model I are indeed quite similar to the values in Table 2 of de Salas et al. (2019), thus confirming that using Einasto profile will predict significantly lower total MW masses than when NFW or gNFW profiles are used. Small differences between the two works are probably due to the different schemes in interpreting the systematics of the Eilers et al. (2019) RC. We also retrieved similar results by Karukes et al. (2020), who also studied the effect of changing the DM density profile. While it goes in the same direction (the Einasto profile predicts lower total masses than the gNFW), their results have not been applied on the accurate Gaia DR2 MW RC, which prevents a detailed comparison.

We are aware that using the RC up to 25 kpc to constrain the mass density profile of the MW is a limited exercise because it needs to be extrapolated to larger radii (see Fig. 4). Extrapolations of the total mass from a rotation curve is incorrect, although it has been used very often in the literature either for giant spirals such as the MW (see, e.g., Eilers et al. 2019 and references therein) of for dwarfs (see, e.g., Read et al. 2016). Other works used different mass tracers such as globular clusters (Vasiliev 2019), massive and very bright stars (Deason et al. 2021), or dSph galaxies assumed to be satellites of the MW (Callingham et al. 2019). These methods have the advantage of sampling objects much farther out in the MW halo, although their virial equilibrium with the MW potential is less guaranteed than for rotating disk stars (Eilers et al. 2019). A warp and flare that occur at radii larger than 12 kpc may also limit our study, especially in the outer disk. However, the effect is possibly limited for our χ² fitting because the error bars are very large in the outer disk because they account for the action of the vertical component (see Sect. 2.1 and Mackereth et al. 2019).

Fig. 4. Mass model derived from the MW RC (left panel) and extrapolated to larger radii (right panel), using Model I for the MW baryonic mass. The best-fit low- and high-mass models are shown as solid, dotted, and dash-dotted lines from Einasto (green), NFW (magenta), and gNFW (orange) models, respectively. Areas showing the possible mass ranges are shaded using the same color code. This shows how the NFW and gNFW bias the mass determination from RCs.

There are two other limitations of our study. The first is linked to the adoption of a spherical halo, although constraints on the dark matter halo shape in the Milky Way are still weak (see Read 2014). The second limitation is linked to our choice of initial (flat) priors for DM halo profiles, and this might alter the validity of our results. We compared our initial halos with the Dutton & Macciò (2014) CDM simulations, in particular, through the relation between concentration and total mass (M200). Our very broad range of parameters encompasses all the Dutton & Macciò (2014, see Fig. 3 of Udrescu et al. 2019) values in the range of 10¹¹ − 10^12.5 M_⊙ and in the (c, M200) plane. The solutions that fit the MW rotation curve are also well within the range of halos simulated by Dutton & Macciò (2014).

Interestingly, our mass boundaries for the χ² fitting of the MW RC encompass all these values using other mass tracers. The question remains which mass density profile is the most suitable to properly evaluate the DM contribution to the MW RC. During the submission of this paper, a study by Cautun et al. (2020) was published. They provide a detailed analysis of the effect of baryons on the DM profile. It results in a contracted halo in the spatial region in which the RC is determined. While the total mass is assumed to be consistent within the error bars with the Callingham et al. (2019) value (M_tot = 1.17 ± 0.18 × 10¹² M_⊙), Cautun et al. (2020) succeeded to fit the MW RC that provided most of the constraints, given its accuracy. Together with our study, this leads us to three important remarks:

• When the MW RC alone is used as a constraint, we find that the Einasto mass density profile leads to the largest range of MW total masses that can reproduce its RC, while both NFW and gNFW profiles lead to a narrow mass range, in particular, by excluding total mass (M₂₀₀ = M_tot) values lower than ∼5 × 10¹¹ M_⊙.

• The contracted halo density profile might be difficult to reproduce by NFW or by gNFW profiles (Cautun et al. 2020), while it is part of the solutions of this paper using an Einasto profile combined with the baryonic model of Cautun et al. (2020) (see Fig. 5).

• We find that both NFW and gNFW profiles provide total masses that increase with baryonic masses (see rows four to nine in Table 4), which contradicts expectations that the DM compensates for a lack of mass from baryons in a galaxy. This contradicts the Einasto predictions, according to which the DM mass is higher when the assumed baryonic mass is lower.

Fig. 5. Dark-matter enclosed mass vs. radius for the mass profile of the contracted halo of Cautun et al. (2020), from which points and error bars are given in red. When the same baryon content is assumed, the black curves show the result from 24 Einasto models that fit both the RC and the contracted halo. The total mass are very similar within a few percent, and the only small difference is that R200 ranges from 200 to 213 kpc instead of 218 kpc for the contracted halo of Cautun et al. (2020). The inset shows a zoom of the mass distribution below 25 kpc to show the similarity of the Einasto DM and the contracted halo near the range of radii of the RC.

It might have been envisioned that these limitations of the NFW profile are related to its two-parameter nature, but this seems to be ruled out by the (almost) similar behavior of the three-parameter gNFW profile. Alternatively, this might be attributed to the density profile of the models in the outskirts. Both have an analytical form that imposes a constant slope of the density profile reaching −3 at large distances, leading to an enclosed mass value that does not converge because it increases as the logarithm of the radius. Investigations by Nesti & Salucci (2013) of the internal r < 5 kpc MW kinematics showed that a cusp-like NFW (or gNFW for γ > 0) profile may also experience some difficulties when combined to baryonic mass.

4.2. Can the MW has a total mass as low as 2.6 × 10¹¹ M_⊙?

The Einasto profile fit of the RC points toward low total mass values for the MW (see Table 4, Figs. 2 and 3), disregarding any other dynamical tracers farther out in the Milky Way. However, the main result of this paper is provided by the combination of RC fitting with either an Einasto or NFW profile for the MW DM halo, leading to a range of the total MW mass of between 2.5 to 18 × 10¹¹ M_⊙ (see Figs. 1–4). This range is consistent with many studies, including that based on other mass indicators, although they generally disagree with our lowest mass range. Figure 6 compares the orbital energy of globular clusters (GCs) from Vasiliev (2019) with that expected from the most likely (total mass: 2.6 × 10¹¹ M_⊙) and the highest (total mass: 15 × 10¹¹ M_⊙) MW mass model that could reproduce the MW RC when combining Model I for baryons and the Einasto profile for DM. Both are consistent with the scenario that GCs are gravitationally bounded to the MW except for one, Pyxis, which appears to disagree significantly for the lighter model. However, the Pyxis eccentric orbit, metallicity, and age indicate an extragalactic origin of Pyxis (Fritz et al. 2017). This indicates that in absence of other precise mass indicators from 25 to 70 kpc, it may be premature to conclude on the total MW mass value from 2.6 × 10¹¹ M_⊙ and R₂₀₀ = 135 kpc (χ² probability = 0.999) to 15 × 10¹¹ M_⊙ and R₂₀₀ = 236 kpc (χ² probability = 0.35), and even 18 × 10¹¹ M_⊙ with a χ² probability = 0.05.

Fig. 6. Kinetic energy (KE) of the GCs from Vasiliev (2019) (crosses with red error bars) compared to the blue and green thick lines that indicate the potential energy (PE, absolute values) expected from the most likely and the heaviest Einasto model when associated to Model I for baryons, respectively. Error bars have been estimated with Monte Carlo randomly sampling by considering the errors in distance and radial velocity, as well as errors in proper motion and their covariance. The (small) thickness of the potential lines is due to the presence of the axisymmetric disk component.

We remark that a low value for the MW mass would have considerable consequences on the orbits of many dSph galaxies, for instance. For example, Boylan-Kolchin et al. (2013) convincingly showed that an MW mass significantly higher than 10¹² M_⊙ is necessary to bound Leo I. Using the Boylan-Kolchin et al. (2013) phase space plot, Hammer et al. (2020) showed that Gaia DR2 orbits might indicate a passage more recent than 4 Gyr ago for many dSphs, assuming a total mass of 8.66 × 10¹¹ M_⊙ for the MW (Eilers et al. 2019). Because MW dSphs also have a peculiar planar alignment (Pawlowski et al. 2014), Deason et al. (2021) opted to use halo stars. After a thorough analysis of the possible recent accretions based on phase-space diagrams, they derived a total mass within 100 kpc of 6.07 × 10¹¹ M_⊙², with which they associated systematics up to 1.2 × 10¹¹ M_⊙. This is only marginally consistent with a very low MW mass and makes a future study of Gaia EDR3 results promising that combines the MW RC and GC motions (Wang et al., in prep.).

5. Summary

Rotation curves are major tools for determining the dynamical mass distribution in the Milky Way and spiral galaxies (Sofue 2013). They are also historically at the root for the requirement of DM in galactic halos (Bosma 1978; Rubin et al. 1980), especially when they have been derived from the HI gas, which often extends far beyond the optical disk. Since the end of the 70 s, many estimates of the DM content in many spiral galaxies were derived, generally through extrapolations of the observed rotation curves of spiral galaxies.

We have tested the most frequently used density profile to perform numerous analyses of galaxy RCs, namely the NFW density profile (Navarro et al. 1997), and its generalization to three parameters, the gNFW profile. We considered the MW RC because it is one of the most accurately determined RCs (Eilers et al. 2019), and also because the MW has not had a major merger since ∼10 Gyr (Hammer et al. 2007; Helmi 2020). This supports the idea that its disk is dynamically virialized to at least 30 kpc because Gnedin & Ostriker (1999) showed that it takes more than three dynamical times for a system to virialize after a perturbation.

In contrast to the NFW (or gNFW) profile for DM, the three-parameter Einasto profile (Einasto 1965, see also Retana-Montenegro et al. 2012) may account for many types of outer slopes, and it provides a much better fit of the simulated DM properties (Dutton & Macciò 2014, and references therein), including for the physically motivated contracted halo (Cautun et al. 2020). It also shows consistent results that can fit the MW RC with most combinations of baryonic mass models, generating a plausible wide range of possible total masses (see Figs. 2–4).

Methodological problems due to the use of a too analytically constrained density model may affect the current estimates of the MW mass such as were reported by Eilers et al. (2019). Perhaps this also applies to the numerous galaxies for which the RC has been analyzed. Other galaxy RCs have yet to be analyzed using a three-parameter density model for the DM as we did here, although see Chemin et al. (2011) for their promising results. These future investigations should focus first on galaxies that did not experience a recent major merger during which most of the disk was resettled or rebuilt (Hammer et al. 2005, 2009; Hopkins et al. 2009). For example, an event like this might complicate the interpretation of the M31 RC, whose recent major merger 2−3 Gyr ago has had a more serious impact (see Hammer et al. 2018) than that of the Sagittarius passage near the MW. The Sagittarius passage is thought to have created vertical waves within the MW disk, although this is still disputed (see Bennett & Bovy 2021, and references therein), while the recent merger in M31 has completely destroyed the thin disk of M31 for stars with ages older than 2 Gyr (see the modeling by Hammer et al. 2018, which reproduced the anomalous age-velocity dispersion discovered by Dorman et al. 2015). In addition, it is also possible that other two-parameter models are affected in a similar manner, for instance, the isothermal model, which renders comparisons of the validity of these profiles for fitting RCs somewhat obsolete.

Using the Einasto profile, we find that the MW mass is mostly constrained by its slightly declining RC, which leads to higher χ² probabilities for low-mass values (i.e., slightly below 3 × 10¹¹ M_⊙) for the MW, although less probable higher values up to 18 × 10¹¹ M_⊙ cannot be excluded. This causes a revision of the available total mass range of the MW down to values that can be as low as 2.6 × 10¹¹ M_⊙, which are also consistent with the kinetic energy distribution of globular clusters. Further improvements of the accuracy of the MW RC will be invaluable to support or reject these low total masses. They would be invaluable in particular for determining precise orbits for the MW dSphs, for which, given the Gaia EDR3 precision, most uncertainties now come from our insufficient knowledge of the total MW mass.

Acknowledgments

We are very grateful for the useful and insightful discussions with Christina Eilers about the MW RC and the treatment of the systematics. We warmly thank Piercarlo Bonifacio for his careful advices and comments on the manuscript, and also Frederic Arenou and Carine Babusiaux for their contributions in the meetings during which the methodology of the paper has been adopted. We are acknowledging the support of the International Research Program Tianguan led by the French CNRS and the Chinese NAOC and Yunnan University.

References

Appendix A: Data of the MW RC and adopted error bars

Table A.1 provides the data for the MW RC given by Eilers et al. (2019), for which they defined the statistical errors (${\sigma }_{{\mathit{v}}_{\mathrm{c}}}^{-}$ (km s⁻¹) and ${\sigma }_{{\mathit{v}}_{\mathrm{c}}}^{+}$ (km s⁻¹)), and to which we added the systematic error (see Sect. 2.1) in the last column as a fraction of the observed velocity, following the definition made in Figure 4 of Eilers et al. (2019).

All Tables

All Figures

Fig. 1. Top: contribution to the rotation curve of different baryonic models and model components. Red points indicate the rotation curve of the Milky Way from Eilers et al. (2019). The error bars are estimated by bootstrapping and include the systematic uncertainties from the neglected term (see text). Bottom: fit of the rotation curve by the best-fit model (solid blue curve, total mass of 2.8 × 1011 M⊙), and with the most massive MW model for which the χ2 probability reaches P = 0.05 (dash-dotted orange line, total mass of 18 × 1011 M⊙), both associated with the baryonic distribution from model I of Pouliasis et al. (2017).

In the text

	Fig. 2. χ2 probability associated with the combination of different baryonic models with Einasto, NFW, or gNFW mass profiles for the DM. From top to bottom, we show Model I, bulge G2+disk J, bulge E2+disk J, and model S. The red points and error bars are the average and 1σ uncertainties.
In the text

Fig. 3. Maximum (solid lines) and averaged (dotted lines) χ2 probabilities for the different baryonic models. Model names are labeled in each panel, with Einasto, NFW, and gNFW mass predictions in green, magenta, and orange, respectively. The two panels associated with baryonic model E2+dJ and E2+J show no histogram for the NFW because this density profile fails to reproduce the MW RC. The horizontal dotted lines indicates the χ2 probability limit of 0.05 below which a model is found to be unable to fit the MW RC.

In the text

Fig. 4. Mass model derived from the MW RC (left panel) and extrapolated to larger radii (right panel), using Model I for the MW baryonic mass. The best-fit low- and high-mass models are shown as solid, dotted, and dash-dotted lines from Einasto (green), NFW (magenta), and gNFW (orange) models, respectively. Areas showing the possible mass ranges are shaded using the same color code. This shows how the NFW and gNFW bias the mass determination from RCs.

In the text

Fig. 5. Dark-matter enclosed mass vs. radius for the mass profile of the contracted halo of Cautun et al. (2020), from which points and error bars are given in red. When the same baryon content is assumed, the black curves show the result from 24 Einasto models that fit both the RC and the contracted halo. The total mass are very similar within a few percent, and the only small difference is that R200 ranges from 200 to 213 kpc instead of 218 kpc for the contracted halo of Cautun et al. (2020). The inset shows a zoom of the mass distribution below 25 kpc to show the similarity of the Einasto DM and the contracted halo near the range of radii of the RC.

In the text

Fig. 6. Kinetic energy (KE) of the GCs from Vasiliev (2019) (crosses with red error bars) compared to the blue and green thick lines that indicate the potential energy (PE, absolute values) expected from the most likely and the heaviest Einasto model when associated to Model I for baryons, respectively. Error bars have been estimated with Monte Carlo randomly sampling by considering the errors in distance and radial velocity, as well as errors in proper motion and their covariance. The (small) thickness of the potential lines is due to the presence of the axisymmetric disk component.

In the text