The Cumulative Mass Profile of the Milky Way as Determined by Globular Cluster Kinematics from Gaia DR2

Since the Gaia Data Release 2 (DR2) on 2018 April 25, a number of studies have presented new estimates or lower bounds for the mass of the Milky Way (MW) Galaxy (e.g., Gaia Collaboration et al. 2018; Malhan & Ibata 2018; Sohn et al. 2018b; Posti & Helmi 2019; Vasiliev 2019; Watkins et al. 2019). Approaches have included maximum likelihood and Bayesian analyses, and have used data for kinematic tracers such as globular clusters (GCs) and stellar streams provided by DR2.

These studies usually interpret results from the inferred gravitational potential profile or the circular velocity profile given an assumed model and the data. Results are often reported as a mass within a specific distance from the Galactic center, or as M₂₀₀ or the virial mass. While it is useful to have visualizations of the potential and circular velocity profiles in addition to individual estimates of the MW mass within specific distances, it would also be beneficial for studies to present a cumulative mass profile (CMP) of the MW with credible regions. A CMP makes it straightforward to compare mass results reported within different Galactocentric distances, and it can also be compared to the CMP of MW-type galaxies from cosmological, hydrodynamical simulations. Moreover, CMPs resulting from different model assumptions can be compared and used to assess differences in model fits at different distances.

In Eadie et al. (2015, Paper I), Eadie & Harris (2016, Paper II), and Eadie et al. (2017, Paper III), a hierarchical Bayesian method was developed for estimating the total mass and CMP of the MW. This Bayesian model was applied to the MW's GC population data, using the 2010 edition of the Harris (1996, 2010) catalog and supplemented with proper-motion measurements made by many other studies (e.g., Majewski & Cudworth 1993; Zoccali et al. 2001; Feltzing & Johnson 2002; Casetti-Dinescu et al. 2010, 2013; Fritz & Kallivayalil 2015; Rossi et al. 2015; see Paper II for a complete table).

In this paper, we provide the method to calculate a CMP given any model for the total gravitational potential and the posterior distribution of the model parameters acquired from a Bayesian analysis. As a motivating example, we employ a simple model for the total gravitational potential and CMP, and confront the model with Gaia DR2 data of GC kinematics to arrive at a new mass estimate and CMP for the MW. For reasons outlined in Section 2, we only use GCs that reside beyond 15 kpc. This limits the sample considerably, and also increases the percentage of GCs without proper-motion measurements.

The kinematic measurements of GCs released by the Gaia DR2 collaboration do not vastly improve the completeness of tracer data beyond 15 kpc, but the HSTPROMO project has helped contribute valuable measurements for GCs at large distances (Sohn et al. 2018b). Additionally, Vasiliev (2019) estimated the mean proper motions of 150 GCs using the Gaia DR2 data set, thereby increasing the number of proper-motion measurements at larger distances considerably.⁴ The catalog provided by Vasiliev is appealing not only because the data are more complete, but also because (1) all estimates of GC proper motions are measured in a consistent manner, and (2) the measurements agree very well with the available Gaia and HSTPROMO observations of 20 distant GCs by Sohn et al. (2018b).

Therefore, we estimate the mass and CMP of the MW using first the Gaia collaboration and the HSTPROMO project data, and then again using the catalog presented in Vasiliev (2019). To obtain the posterior distribution of model parameters that will be used to determine the CMP, we use the hierarchical Bayesian method presented in Papers I–III, with small adjustments in light of the results in Eadie et al. (2018; hereafter EKH).

The paper is organized as follows.

Section 2 outlines the hierarchical Bayesian model from Papers I to III and the methods used in this study. In Section 2.1, we describe how to calculate the CMP for the MW given a function M(R < r), and given samples from the posterior distribution of model parameters.

Section 3 discusses the details of and differences between the two data sets (Gaia DR2 + HSTPROMO, and Vasiliev).

Section 4 presents the Bayesian estimates of the MW's CMP, given each data set, and some discussion and interpretation of results.

Section 5 summarizes our findings and the impact on future studies.

Specific details about our hierarchical Bayesian model and sampling methods for the posterior distribution can be found in Papers I–III. Here, we provide a brief review for completeness.

The method allows the user to supply (1) kinematic and position data of tracers (such as GCs or halo stars), (2) an analytic distribution function (DF) for the specific energies ${ \mathcal E }$ and angular momenta L of these tracers (based on a gravitational potential and tracer density profile), and (3) hyperprior distributions for the model parameters. The hierarchical Bayesian model includes incomplete and complete velocity data simultaneously, and also accounts for observational error.

The observations of the kinematic tracers are the Galactocentric distance r, line-of-sight velocity v_los, and proper motions (${\mu }_{\alpha }\cos \delta$, μ_δ), and their uncertainties. These measurements are assumed to be samples drawn from normal distributions centered on the tracers' true but unknown values with standard deviations equal to the measurement uncertainties. In other words, the true tracer velocity components and distances are treated as parameters. This measurement model defines the likelihood in the hierarchical Bayesian model.

The prior distribution on the parameters for the true values of r, ${v}_{\mathrm{los}},{\mu }_{\alpha }\cos \delta$, and μ_δ is a DF of the specific energies ${ \mathcal E }$ and angular momenta L, given a model for the total gravitational potential Φ(r) and density profile of the tracer population ρ(r). In Papers II and III, we used a simple power-law profile for both the gravitational potential,

$$\begin{eqnarray}&&{\rm{\Phi }}(r)={{\rm{\Phi }}}_{\circ }{r}^{-\gamma },\end{eqnarray} \tag{ 1 }$$

and the density profile of the tracer population,

$$\begin{eqnarray}&&\rho \propto {r}^{-\alpha },\end{eqnarray} \tag{ 2 }$$

where ${{\rm{\Phi }}}_{\circ },\gamma ,\mathrm{and}\,\alpha$ are parameters. At large distance, a spherically symmetric gravitational potential (Equation (1)) is a reasonable approximation; recent studies regarding the velocity ellipsoid of the stellar halo have suggested the potential is more spherical than flattened (Evans et al. 2016; Wegg et al. 2019).

Together, Equations (1) and (2), determine the analytic DF $f({ \mathcal E },L)$ first presented in Evans et al. (1997). The DF also includes the constant velocity anisotropy parameter β for the tracer population as a free parameter. The derivation of this DF is shown in Evans et al. (1997), with a condensed version shown in Paper II (although notations differ).

Hyperprior distributions are also defined for the four model parameters: Φ_◦, γ, α, and β. In Papers II–III, the parameters Φ_◦, γ, and β were given uniform prior distributions with upper and lower bounds of (1, 200), (0.3, 0.7), and (−0.5, 1), respectively, while the prior distribution on α was a Gamma distribution determined by data not used in the analysis.

The posterior distribution is the probability distribution of model parameters given the likelihood, prior, hyperpriors, and data. Instead of calculating the normalization constant in Bayes' theorem, samples from a distribution proportional to the posterior distribution are acquired. This is done by running multiple Markov chains in parallel, until they have reached a common, stationary distribution assumed to be proportional to the posterior distribution. Both visual inspection and statistical diagnostics are used to assess the mutual convergence of the chains. The code, called GME, was written in the ${\boldsymbol{R}}$ Statistical Software environment (see Papers I–III).

The hierarchical Bayesian model outlined above was used to estimate the total mass and CMP of the MW, given kinematic data for the MW GC population before Gaia DR2 became available. More recently, EKH tested this method on tracer data from simulated galaxies made in the McMaster Unbiased Galaxy Simulations 2 (MUGS2).

The study of MUGS2 simulated data was done as a strictly blind test. There were two main caveats to the study: (1) the tracer data from MUGS2 were GC analogs (i.e., star tracer particles instead of actual GCs), and (2) half of the galaxies created in the MUGS2 simulations did not follow the standard stellar-mass-to-halo-mass relation. Despite these caveats, the galaxies had the basic components of a bulge, disk, and dark matter halo, and it was therefore instructive to observe how well the single power-law gravitational potential model could describe the overall system.

The results of the blind tests were mixed; the total masses (M₂₀₀) of some galaxies were estimated well within the 95% Bayesian credible regions, while others were not. In particular, there was a tendency for the mass to be underestimated. EKH suggested a number of compounding reasons for the bias (e.g., the physical location of incomplete data and the spatial distribution of tracers in the simulation), but a main factor seemed to be the model for the total gravitational potential. The model had difficulty predicting both the inner and outer regions of the true CMPs for the simulated galaxies when tracers at all Galactocentric radii were included.

Of the galaxies whose masses were poorly estimated, the posterior distribution was pushed toward and reached the edge of allowable parameter space. Posterior distributions that are at an extreme end of parameter space allowed by the prior distribution usually indicate a poor fit to the data, and can be a red flag that the results should be interpreted with caution. Moving forward, we keep this in mind as we apply the method of Paper III to the new Gaia data.

EKH also showed that using only the outermost tracers from the simulated galaxies improved the mass estimates and increased the chance of containing the true mass profile within the Bayesian credible regions. These results echoed the sensitivity test in Paper III, which also used a power-law potential. In the latter case, the removal of inner GC data from the analysis caused a slight increase in the mass estimate for the MW. Thus, EKH recommended that in future studies, only outer kinematic tracers be used in the analysis when a single power-law model for the gravitational potential is employed. For this study, we use GCs at $r\gt 15\,\mathrm{kpc}$.

EKH also found that using a uniform distribution as a prior on γ may have been too broad. The uniform distribution was centered on 0.5, but also allowed values as low and as high as 0.3 and 0.7. When γ = 0.5, Equation (1) approaches a Navarro–Frenk–White-type halo at large distance, and this is the prior information we intended to include. Thus, for this study, we replace the uniform distribution for p(γ) with a normal distribution with mean 0.5 and standard deviation 0.06 (Figure 1).

Zoom In Zoom Out Reset image size

In this study, we continue to use a Gamma distribution for the p(α) with hyperparameters defined by GC data not used in the analysis because of their lack of v_los and $\bar{\mu }$ measurements (Figure 2). Because we are excluding GCs within 15 kpc, there are only three GCs with $r\gt 15\,\mathrm{kpc}$ that are lacking both line-of-sight and proper-motion measurements: AM 4, Ko 1, and Ko 2.

Zoom In Zoom Out Reset image size

2.1. Calculating the CMP with Bayesian Credible Regions

To calculate Bayesian credible regions for a CMP, given samples from the posterior distribution, perform the following steps:

• 1.  
  Create a sequence of n closely spaced r values for the desired range from the Galactic center.
• 2.  
  At each r_i for i in (1, ..., n), calculate M(<r_i) for every sample from the posterior distribution. For each r_i value, you should have m values of M(<r_i), where m is the number of samples from the posterior distribution (e.g., the length of the Markov chain).
• 3.  
  At each r_i, sort the m values of M(<r_i), and calculate the credible intervals of interest (e.g., inner 50%, 75%, and 95%).

In our study, the total gravitational potential assumed in Equation (1) leads to a CMP of the form

$$\begin{eqnarray}&&M(\lt r)=\displaystyle \frac{\gamma {{\rm{\Phi }}}_{o}}{G}{\left(\displaystyle \frac{r}{\mathrm{kpc}}\right)}^{1-\gamma }\end{eqnarray} \tag{ 3 }$$

(Deason et al. 2012c). Thus, for a Markov chain of length m, each $({{\rm{\Phi }}}_{\circ ,j},{\gamma }_{j})$ pair in the chain is used to calculate the mass within a given radius r_i. This results in a vector of masses within r_i, which are sorted and used to calculate the marginal distribution of the mass within r_i. We use n = 100 values logarithmically spaced from ${r}_{1}=0.1\,\mathrm{kpc}$ to ${r}_{100}=200\,\mathrm{kpc}$. Our chains are of length 90,000, with effective sample sizes of >2000 for all parameters.

As part of Gaia DR2, Gaia Collaboration et al. (2018) published kinematic data for 75 GCs in the MW, including proper motions, positions, distances, and line-of-sight velocities.⁵ We use these data to replace the relevant GC measurements in the Harris (1996, 2010 edition) catalog, with some exceptions.

The distance measurements to GCs as measured by Gaia show a systematic difference that should improve with the next data release (Gaia Collaboration et al. 2018). In the meantime, we follow the Gaia Collaboration et al. (2018) guidelines and continue to use the distances provided in Harris (2010). The Gaia data are also not entirely complete—only 57 of the 75 GCs have line-of-sight velocity measurements. For those that are missing v_los measurements, we use the ones listed in Harris (2010) when available.

The HSTPROMO team recently provided proper motions for 20 GCs at larger distances (Sohn et al. 2018b). Four GCs measured by HSTPROMO are also measured by Gaia (NGC 362, NGC 2298, NGC 2808, NGC 3201), and in these cases we use the DR2 data. In all other cases, we let the HSTPROMO measurements of GCs supersede any previous measurements in the literature. We then supplement other missing proper motions measurements with those reported in other studies (Majewski & Cudworth 1993; Zoccali et al. 2001; Casetti-Dinescu et al. 2010, 2013; de la Fuente Marcos et al. 2015; Rossi et al. 2015).

The total number of GCs in our compiled data set is 157, but 4 GCs—Arp 2, Pal 12, Terzan 7, and Terzan 8—were recently confirmed by Sohn et al. (2018b) to be associated with the Sagittarius (Sgr) dwarf galaxy. Thus, Sohn et al. (2018b) used only one of these GCs (Arp 2) in their study of the MW's mass. They also report no significant change in the mass estimate if Pal 12, Terzan 7, or Terzan 8 is used instead. We follow their lead and exclude all but Arp 2 from our analysis. Law & Majewski (2010) listed a few other GCs that may be associated with Sgr, but these have yet to be confirmed and so we leave these GCs in our analysis. Thus, our sample consists of 154 GCs, with a breakdown of complete and incomplete data shown in Figure 3.

Zoom In Zoom Out Reset image size

Together, the Gaia collaboration and the HSTPROMO team have greatly increased the number of proper-motion measurements of GCs in the MW. Only 52 out of 154 GCs in our updated Gaia+HSTPROMO catalog have incomplete measurements. This is a significant improvement from the data set used previously (Paper II–Paper III), in which approximately 50% of the data were incomplete.

Even more recently, Vasiliev (2019) used the Gaia DR2 data to calculate the mean proper motions of 150 GCs in the MW and substantially increase the percentage of complete data. The new Vasiliev catalog is an appealing data set not only because of its completeness, but also because the proper motions are calculated in a consistent manner. Thus, we use this new data set to form a second catalog of GCs and estimate the MW mass and CMP for comparison.

We generate the second catalog from the Vasiliev data such that it contains the same GCs as our Gaia+HSTPROMO catalog. We again start with the Harris catalog and then replace all GC measurements of proper-motion and line-of-sight velocity with the new measurements reported by Vasiliev (the line-of-sight velocities originate from Baumgardt & Hilker 2018).

The Vasiliev catalog also contains information for the GC Leavens 1 (Crater) and recently discovered GC FSR 1716 (Minniti et al. 2017; Koch et al. 2017), but we exclude these data because they are not in the Gaia+HSTPROMO catalog. Crater is also the most distant GC known to date (Laevens et al. 2014), and its proper-motion measurement is quite uncertain (Vasiliev 2019). Although we do not include Crater in this study, its distance could make it an important GC for future MW mass estimates once its proper motion is known with more certainty.

Overall, the main difference between the Gaia+HSTPROMO and Vasiliev catalogs is completeness, with the latter containing complete measurements for 145 GCs, and incomplete measurements for 9 GCs (Figure 4).

Zoom In Zoom Out Reset image size

As mentioned in Section 2, we use GCs with Galactocentric distances $r\gt 15\,\mathrm{kpc}$, thus limiting our sample size to 35 in both our Gaia+HSTPROMO and Vasiliev catalogs. In this reduced Gaia+HSTPROMO catalog, 16 GCs are missing $\bar{\mu }$ measurements, 3 of which are also missing v_los measurements. In contrast, the latter three GCs (AM 4, Ko 1, and Ko 2) are the only incomplete data in the Vasiliev catalog. Because only the distances and positions of AM 4, Ko 1, and Ko 2 are known, these GCs are removed from both catalogs and used to define the α prior distribution's hyperparameters (see Section 2, Figure 2, and Papers II–III).

Figure 5 shows the proper-motion components and line-of-sight velocities as a function of Galactocentric distance for GCs at $r\gt 15\,\mathrm{kpc}$. The Gaia+HSTPROMO measurements are shown with open blue circles and the Vasiliev's measurements are gray diamonds. There is excellent agreement between the catalogs' measurements. The only outlier is the μ_δ measurements for Pal 3, which can be seen in the second panel at 95.7 kpc.

Zoom In Zoom Out Reset image size

Figure 6 shows the CMP of the MW with Bayesian credible regions (50%, 75%, and 95%) given the Gaia+HSTPROMO catalog (left) and the Vasiliev catalog (center). These CMPs were calculated using the posterior distribution of model parameters (Section 2.1), and GCs at $r\gt 15\,\mathrm{kpc}$. The points with error bars in Figure 6 are mass estimates reported in other studies, which are identified in the legend (Kochanek 1996; Wilkinson & Evans 1999; Sakamoto et al. 2003; Battaglia et al. 2005; Xue et al. 2008; Gnedin et al. 2010; Watkins et al. 2010; McMillan 2011; Deason et al. 2012b, 2012a; Kafle et al. 2012; Gibbons et al. 2014; Eadie et al. 2015; Küpper et al. 2015), including five that also use Gaia DR2 data (Malhan & Ibata 2018; Sohn et al. 2018a; Posti & Helmi 2019; Vasiliev 2019; Watkins et al. 2019).

Zoom In Zoom Out Reset image size

The two CMPs in Figure 6 are nearly identical, but the Vasiliev catalog provides a slightly better constraint on the mass because the data are more complete. The median estimates of M₂₀₀ given each catalog are also extremely similar at

$$\begin{eqnarray}&&0.70\,(0.47,1.14)\times {10}^{12}\,{M}_{\odot }\,\,(\mathrm{Gaia}+\mathrm{HSTPROMO})\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&0.70\,(0.51,1.10)\times {10}^{12}\,{M}_{\odot }\,\,\,(\mathrm{Vasiliev}),\end{eqnarray} \tag{ 5 }$$

where numbers in brackets are the lower and upper bounds of the 95% Bayesian credible region. The quantiles for the model parameters calculated from the posterior distributions are also extremely similar between the two data sets. Henceforth, we refer to the Vasiliev results in the discussion.

A large number of GCs in our catalogs lie between 15 and 20 kpc of the Galactic center (Figure 5), and yet we extrapolate the results out to larger distances with the CMP. To test how sensitive our method is to these inner GCs, we perform the analysis again on the Vasiliev catalog, using a new cut of $r\gt 20\,\mathrm{kpc}$.

Figure 7 compares the marginal posterior distributions and medians for M₂₀₀, given r_cut values of 15 kpc and 20 kpc. The marginal distributions are very similar, but the median estimate for M₂₀₀ and the upper bound of the 95% credible interval are slightly higher:

$$\begin{eqnarray}&&0.77\,(0.51,1.40)\times {10}^{12}\,{M}_{\odot }\,\,({r}_{\mathrm{cut}}=20\,\mathrm{kpc}).\end{eqnarray} \tag{ 6 }$$

The new CMP is virtually unchanged barring the larger upper values for the Bayesian credible regions, indicating that the method is relatively insensitive to the tracers between 15 and 20 kpc.

Zoom In Zoom Out Reset image size

Recently, Myeong et al. (2018) found evidence that up to ten GCs (NGCs 362, 1261, 1851, 1904, 2298, 2808, 5286, 6779, 6864, and 7089) were brought into the Galaxy during a massive merger event (Belokurov et al. 2018; Kruijssen et al. 2018). All of these GCs reside within 20 kpc of the Galaxy, and thus the sensitivity test just mentioned has excluded them (Figure 7). Four of Myeong's GCs, however, do reside between 15 and 20 kpc, so testing their possible influence on the mass estimate when ${r}_{\mathrm{cut}}=15\,\mathrm{kpc}$ is a worthwhile exercise.

Interestingly, when the four Myeong GCs between $15\,\mathrm{kpc}\lt r\lt 20\,\mathrm{kpc}$ are excluded from the analysis, the median estimate of M₂₀₀ increases slightly from Equation (5) to

$$\begin{eqnarray}&&0.77\,(0.53,1.23)\times {10}^{12}\,{M}_{\odot }.\end{eqnarray} \tag{ 7 }$$

Yet, when we put these GCs back into the analysis and remove four other GCs between $15\,\mathrm{kpc}\lt r\lt 20\,\mathrm{kpc}$ at random, the median estimate is

$$\begin{eqnarray}&&0.70\,(0.49,1.11)\times {10}^{12}\,{M}_{\odot }.\end{eqnarray} \tag{ 8 }$$

This result seems to suggest that the Myeong's GCs push our mass estimate to smaller values. Although this is intriguing, the effect is minor when we consider the width of the Bayesian credible regions. We summarize the results of this sensitivity test in Table 1.

Table 1.  Mass Estimates and 95% Credible Regions for M₂₀₀ (${10}^{12}\,{M}_{\odot }$) Given Different Subsets of the Vasiliev Data

Data used		M200 (${10}^{12}\,{M}_{\odot }$)
	Lower 95% c.r.	Median	Upper 95% c.r.
r > 15 kpc	0.51	0.70	1.10
r > 15 kpc, Myeong GCs removed	0.53	0.77	1.23
r > 15 kpc, 4 random GCs removed	0.49	0.70	1.11
r > 20 kpc	0.51	0.77	1.40

Download table as:  ASCIITypeset image

The main advantage of our CMP for the MW is the ability to easily compare our mass estimate to many other studies both visually (i.e., Figure 6) and quantitatively, within any distance from the Galactic center. To illustrate the quantitative comparisons, we present slices of our CMP M(<r) taken at $r=25\,\mathrm{kpc}$ through $r=150\,\mathrm{kpc}$, by steps of 25 kpc, to show the marginal posterior distributions for (M < r) at those distances (Figure 8). Marginal distributions for the mass within any r are also easily calculated.

Zoom In Zoom Out Reset image size

Our CMP agrees with many of the results reported in the literature within their uncertainties, except for three higher-mass estimates in Figure 6: Watkins et al. (2010) (with Draco) at 100 kpc, Eadie et al. (2015) at 125 kpc, and Vasiliev (2019) at 100 kpc. Notably, the first two of these studies included dwarf galaxies in their samples, whereas we do not. There is evidence suggesting that some dwarf galaxies may inflate the mass estimate of the MW. Case in point, the result of Watkins et al. (2010) without Draco is in good agreement with our estimate at 100 kpc.

Other studies that use dwarf galaxies, or even high/extreme velocity stars, also tend to infer "heavier" mass estimates of the MW. Many of these studies' results cannot be included on the CMP (Figure 6) because they are reported as M₂₀₀ or the virial mass M_vir. Instead, we compare three of these studies' results to our marginal distribution for M₂₀₀ given the Vasiliev catalog (Figure 7).

Patel et al. (2017) performed a Bayesian analysis in conjunction with cosmological simulations of MW-type galaxies to infer the Galaxy's mass. They used the kinematics of one massive satellite, and inferred ${M}_{200}=1.223\times {10}^{12}\,{M}_{\odot }$ and a 95% Bayesian credible of $(1.207,1.248)\times {10}^{12}\,{M}_{\odot }$. Hattori et al. (2018) also deduced a high mass using newly discovered extreme velocity stars in the Gara DR2 data. If these stars are bound to our Galaxy, then they imply ${M}_{200}\approx 1.4\,\times {10}^{12}\,{M}_{\odot }$. As a final example, Monari et al. (2018) used counter-rotating halo stars in Gaia to estimate the escape velocity of the MW between 5 and 10 kpc. Their results also implied a heavy MW with a dark matter virial mass of ${M}_{200}={1.55}_{-0.51}^{+0.64}\,\times {10}^{12}\,{M}_{\odot }$. All three of these mass estimates are almost completely ruled out by our estimate of M₂₀₀ (Figure 7), which did not include high/extreme velocity stars or dwarf galaxies.

It is less clear why the estimate at 100 kpc by Vasiliev (2019) is slightly higher than our estimate, although the uncertainties do overlap with our 95% c.r.. Two major differences exist between our studies. The first difference is that they used all GCs in their catalog, whereas we only used GCs beyond 15 kpc. The second difference is the choice of gravitational potential model. While we chose a simple power-law model and excluded inner GCs in our analysis, they used a three-component model and included GCs at all distances.

We also compare our results to studies that use GC data from Gaia and/or HSTPROMO. For example, Sohn et al. (2018b) estimated the MW mass within 39.5 kpc using 16 GCs with proper motions measurements by HSTPROMO. Utilizing the estimator from Watkins et al. (2010) and Monte Carlo simulations to estimate the uncertainty, they found $M(r\lt 39.5\,\mathrm{kpc})\,={0.61}_{-0.12}^{+0.18}\times {10}^{12}\,{M}_{\odot }$. Our estimate for the mass within 39.5 kpc is

$$\begin{eqnarray}&&M(r\lt 39.5\,\mathrm{kpc})=0.33\,(0.26,0.45)\times {10}^{12}\,{M}_{\odot },\end{eqnarray} \tag{ 9 }$$

which does not agree with their results within the uncertainties. However, more recently Watkins et al. (2019) combined Gaia DR2 data with HSTPROMO data, and used the same estimator as Sohn et al. (2018b). They found $M(r\lt 39.5\,\mathrm{kpc})\,={0.44}_{-0.06}^{0.07}\times {10}^{12}\,{M}_{\odot }$. Reassuringly, the latter agrees with our estimate at this distance within the uncertainties, despite their use of a different model for M(<r).

Posti & Helmi (2019) also used the kinematic data of 75 GCs provided by Gaia DR2 and HSTPROMO to estimate the mass of the MW in a Bayesian analysis. In contrast to this study, they used a two-component DF for the GC system of the MW. The assumption is that the GC population dynamics can be decomposed into disk-like and halo-like components, with the latter containing parameters for the mean rotational velocity and velocity anisotropy of the GC system. They also included a shape parameter for the dark matter halo. Because their study used a Bayesian paradigm, like ours, their method included measurements uncertainties and incomplete data directly in the analysis. Despite using a more complicated DF than this study, their mass estimate within 20 kpc ($M(r\lt 20\,\mathrm{kpc})={0.191}_{-0.15}^{+0.17}\,\times {10}^{12}\,{M}_{\odot }$) is within the bounds of our 95% Bayesian credible regions:

$$\begin{eqnarray*}&&M(r\lt 20\,\mathrm{kpc})=0.24\,(0.18,0.32)\times {10}^{12}\,{M}_{\odot }.\end{eqnarray*}$$

Vasiliev (2019) estimates the mass of the MW using his GC catalog in a maximum likelihood analysis. Vasiliev's analysis assumes a single-component DF for the GC population, but a three-component gravitational potential to account for the bulge, disk, and halo components. He finds $M(r\lt 50\,\mathrm{kpc})\,={0.54}_{-0.08}^{+0.11}\times {10}^{12}\,{M}_{\odot }$. The uncertainties overlap with the upper bound of the 95% Bayesian credible for our mass estimate within 50 kpc:

$$\begin{eqnarray*}&&M(r\lt 50\,\mathrm{kpc})=0.37\,(0.29,0.51)\times {10}^{12}\,{M}_{\odot }.\end{eqnarray*}$$

Finally, Malhan & Ibata (2018) used Gaia data for the GD-1 stellar stream and found $M(r\lt 14.5\,\mathrm{kpc})={0.175}_{-0.05}^{+0.06}\,\times {10}^{12}\,{M}_{\odot }$, which is in agreement with

$$\begin{eqnarray}&&M(r\lt 14.5\,\mathrm{kpc})=0.20\,(0.15,0.28)\times {10}^{12}\,{M}_{\odot }\end{eqnarray} \tag{ 10 }$$

given by our method.

Overall, it appears that studies are reaching a consensus on the mass of the MW within 40 kpc since the release of Gaia DR2 and HSTPROMO. Even when different analysis methods and models are confronted with the same data, similar conclusions are reached. Albeit, beyond 15 kpc the broken power-law models used in Vasiliev (2019) and Posti & Helmi (2019) are similar to our single power-law model, so the agreement here is perhaps less surprising.

Mass estimates of the Galaxy out to larger distances, however, are still rather uncertain. For example, Sohn et al. (2018b) and Watkins et al. (2019) reported virial mass estimates of ${M}_{\mathrm{vir}}\,={2.05}_{-0.79}^{+0.97}\times {10}^{12}\,{M}_{\odot }$ and ${M}_{\mathrm{vir}}$ = ${1.41}_{-0.52}^{+1.99}\times {10}^{12}\,{M}_{\odot }$, respectively. Posti & Helmi (2019) provided a tighter constraint on the virial mass at ${M}_{\mathrm{vir}}=1.3\pm 0.3\times {10}^{12}\,{M}_{\odot }$, as does our Bayesian estimate of M₂₀₀ (Equation (5)). Although there is some tension between results, our marginal credible regions for M₂₀₀ (Figure 7) do overlap with those of Posti & Helmi (2019).

The state of knowledge about the MW's total mass is still uncertain, but the state of data for tracers has vastly improved. Multiple studies have now estimated the mass within 40 kpc, and in general they are in good agreement with one another. Our CMP of the MW within $r\lt 50\,\mathrm{kpc}$ is also in agreement with many recent results.

Our CMP overlaps more than half of the mass estimates at different radii, suggesting a less massive MW. However, mass estimates that disagree with our results tend to be at larger distances. These discrepancies may be attributed to the inclusion/exclusion of tracers such as dwarf galaxies and/or high-velocity stars in the analyses, but it could also be due to differing model assumptions. Thus, the Galaxy's mass within larger distances is still up for debate. Until more complete measurements for tracers at larger distances are acquired, the total mass of the MW will remain uncertain.

Comparing and verifying mass estimates at large distances will remain challenging without more data, because of the variety of models and data analysis techniques. As a community, we compare results by looking at the uncertainties presented in each others' studies, with the caveat that everyone's methods and data are different. We have seen how the Gaia DR2 and HSTPROMO data have greatly improved our understanding of the mass within 40 kpc of the MW within this context. Thus, it is imperative to obtain complete velocity measurements of tracers at large distances to reach a similar consensus about the Galaxy's total mass.

At the same time, studies continue to report mass estimates within different distances from the Galactic center, making comparisons difficult. The CMP presented in this work instead provides a holistic view of our understanding and estimate of the MW's mass at all distances. Thus, we hope CMPs with Bayesian credible regions become a standard way to present and compare estimates of the MW's mass in future studies, and openly share instructions to create these profiles.

In summary, future work on the mass of the MW will benefit greatly from complete velocity measurements of tracers at greater distances, and by comparing CMPs from different studies.

The authors thank the anonymous referee for their insightful and helpful report. G. Eadie would like to extend gratitude to E. Vasiliev for sharing his catalog of new proper motions measurements and for useful discussions.

This research was funded by an eScience Institute Postdoctoral Fellowship to G. Eadie, made possible from generous donations from the three Moore, Sloan, and Washington Research Foundations at the eScience Institute, University of Washington (UW). G. Eadie and M. Juric acknowledge support from the University of Washington Institute for Data Intensive Research in Astrophysics and Cosmology (DIRAC). The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. M. Juric acknowledges the support of the Washington Research Foundation Data Science Team Chair fund, and the UW Provost's Initiative in Data-Intensive Discovery. The authors would like to thank members of the DIRAC Institute for insightful conversations regarding this work.

This work uses the data from the European Space Agency mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (https://www.cosmos.esa.int/web/gaia/dpac/consortium).

Software: R Statistical Software Environment (R Development Core Team 2012), including the following packages: stats (part of base R), CODA: Convergence Diagnosis and Output Analysis for MCMC (Plummer et al. 2006a; Plummer et al. 2006b), emdbook: Ecological Models and Data in R (Bolker 2008, 2018 version), ggplot2 (Wickham 2016), limma (vennDiagram function) (Ritchie et al. 2015), MASS: Modern Applied Statistics with S Venables & Ripley (2002), moments: Moments, cumulants, skewness, kurtosis and related tests (Komsta & Novomestky 2015), pracma: Practical Numerical Math Functions (Borchers 2017), RColorBrewer: ColorBrewer Palettes, and SNOW: Simple Network of Workstations (Tierney et al. 2013) (Neuwirth 2014).