Which Milky Way masses are consistent with the slightly declining 5-25 kpc rotation curve?

Discoveries of extended rotation curves have suggested the presence of dark matter in spiral galaxy haloes. It has led to many studies that estimated the galaxy total mass, mostly by using the Navarro Frenk and White (NFW) density profile. We aim at verifying how the choice of the dark-matter profile may affect the predicted values of extrapolated total masses. We have considered the recent Milky Way (MW) rotation curve, firstly because of its unprecedented accuracy, and secondly because the Galactic disk is amongst the least affected by past major mergers having fully reshaped the initial disk. We find that the use of NFW profile (or its generalized form, gNFW) for calculating the dark-matter contribution to the MW rotation curve generates apparently inconsistent results, e.g., an increase of the baryonic mass leads to increase of the dark matter mass. Furthermore we find that NFW and gNFW profile narrow the total mass range, leading to a possible methodological bias particularly against small MW masses. By using the Einasto profile that is more appropriate to represent cold dark matter haloes, we finally find that the Milky Way slightly decreasing rotation curve favors total mass that can be as small as 2.6 $\times 10^{11}$ $M_{\odot}$, disregarding any other dynamical tracers further out in the MW. It is inconsistent with values larger than 18 $\times 10^{11}$ $M_{\odot}$ for any kind of CDM dark-matter halo profiles, under the assumption that stars and gas do not influence the predicted dark matter distribution in the MW. This methodological paper encourages the use of the Einasto profile for characterizing rotation curves with the aim of evaluating their total masses.


Introduction
Gaia DR2 provided accurate stellar proper motions to calculate the circular velocity curve of the Milky Way (MW) up to 25 kpc 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 ), 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 12 M 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(Navarro et al. , 2010Gao 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 Section 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 Section 3 we compare the χ 2 probability distribution for DM represented by the NFW or Einasto profiles. In Section 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. 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.

Rotation curve and error bars
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 Figure 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, 2020, private communication) 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).

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 promi-  (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, where r = √ 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) where Σ 0 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 Figure 1). Nevertheless, we consider it useful for testing the effect of an extremely high baryonic mass for the MW disk and bulge. 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), G2 : ρ bulge (x, y, z) = ρ 0 e −r 2 2 /2 , where ( Table 3): 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 Figure 1 or by comparing Table 1 with Tables 2 and 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), where r 0 is the scale radius, and ρ 0 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 0 and m NFW = 4 π ρ 0 r 3 0 , vary from 2 to 100 kpc and from 1 to 50 ×10 11 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 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 π ρ 0 h 3 n Γ(b E ), vary from 3 to 30, 0.05 to 3 and from 1 to 50 ×10 11 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 −7 M /pc 3 , which comes from Hinshaw et al. (2013). With this definition, the relation between virial radius and virial mass is Article number, page 3 of 10 A&A proofs: manuscript no. 41058corr_v2

Deriving the total MW mass and χ 2 probability
The total MW potential can be obtained through the Poisson equation, 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 We applied the χ 2 method to fit the RC and calculate its associated probability, for which we tested an extremely large parameter space. The χ 2 was calculated by the sum at each disk radius R i , 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 σ stat,i = (σ + v obs,i + σ − v 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 χ 2 probability can be expressed as 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 11 M , for instance. In  (right) profiles. The first profile shows that high χ 2 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 11 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 Figure 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).
The two panels in the middle row of Figure 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 χ 2 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 11 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 11 M , respectively. This suggests the following mechanism: for an increasing baryonic mass, the NFW DM scale radius (r 0 , 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.

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 Figure 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(χ 2 ), M tot ) plane as much as possible.
The solid lines in Figure 3 identify the envelop for each baryonic + DM model, which is defined as being the highest χ 2 probability calculated in mass slices with sizes of 0.3× 10 11 M . We assume that only χ 2 probabilities higher than 0.05 correspond to a good fit of the RC, which we verified after examining the latter. For comparison, Figure 3 also shows the averaged probabilities. 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(χ 2 ), 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 χ 2 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 χ 2 values (best fit, highest probabilities) or on averaging the χ 2 probabilities in each mass slice (average). As in Figure 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 11 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 11 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 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.
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 particu-lar the low-mass values that are favored when the Einasto profile is used.

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 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 1 . 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 Figure 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 . 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 χ 2 fitting because the error bars are very large in the outer disk because they account for the action of the vertical component (see Section 2.1 and Mackereth et al. 2019).
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 & 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.
Macciò (2014, see Figure 3 of Udrescu et al. 2019) values in the range of 10 11 − 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 χ 2 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 12 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 200 =M tot ) values lower than ∼ 5 ×10 11 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 Figure 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.
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 Notes. Models and associated baryonic mass (first and second columns), and estimated total mass using χ 2 probabilities for Einasto, NFW, and gNFW DM density profiles (third to eighth columns, all masses are given in units of 10 11 M ). The total mass and mass ranges are evaluated using the minimum χ 2 (best fit, columns 3 and 5) and by weighting the total masses by their χ 2 probabilities (average, columns 4 and 6), together with associated 1σ uncertainties. Uncertainties also account for systematics related to the Galactic distance and its motion, as well as to change in scale length (∼ 4% on masses, see Sect. 2.1), which have been added to the quoted error bars in this Table. 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 R 200 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. mass value that does not converge because it increases as the logarithm of the radius. Investigations by Nesti & Salucci (2013) of the internal r< 5kpc 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 11 M ?
The Einasto profile fit of the RC points toward low total mass values for the MW (see Table 4, Figure 2 and Figure 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 11 M (see Figures 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 11 M ) and the highest (total mass: 15 ×10 11 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 11 M and R 200 = 135 kpc (χ 2 ) probability =0.999) to 15 ×10 11 M and R 200 = 236 kpc (χ 2 probability =0.35), and even 18 ×10 11 M with a χ 2 probability =0.05. 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 12 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 11 M for the MW . 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 11 M 2 , with which they associated systematics up to 1.2 ×10 11 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. 2021, in preparation).
A&A proofs: manuscript no. 41058corr_v2  (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.

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 70s, 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 , 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 Figures  2, 3, and 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(Hammer et al. , 2009Hopkins 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 2020, 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 χ 2 probabilities for low-mass values (i.e., slightly below 3 ×10 11 M ) for the MW, although less probable higher values up to 18 ×10 11 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 11 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.
A&A proofs: manuscript no. 41058corr_v2 Appendix A: Table A.1 with 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 (σ − v c (km s −1 ) and σ + v c (km s −1 )), 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).