The cumulative mass profile of the Milky Way as determined by globular cluster kinematics from Gaia DR2

We present new mass estimates and cumulative mass profiles (CMPs) with Bayesian credible regions for the Milky Way (MW) Galaxy, given the kinematic data of globular clusters as provided by (1) the Gaia DR2 collaboration and the HSTPROMO team, and (2) the new catalog in Vasiliev (2018). We use globular clusters beyond 15kpc to estimate the CMP of the MW, assuming a total gravitational potential model $\Phi(r) = \Phi_{\circ}r^{-\gamma}$, which approximates an NFW-type potential at large distances when $\gamma=0.5$. We compare the resulting CMPs given data sets (1) and (2), and find the results to be nearly identical. The median estimate for the total mass is $M_{200}= 0.71 \times 10^{12} M_{\odot}$ and the 50\% Bayesian credible region bounds are $(0.63, 0.81) \times 10^{12} M_{\odot}$. However, because the Vasiliev catalog contains more complete data at large $r$, the MW total mass is better constrained by these data. In this work, we also supply instructions for how to create a CMP for the MW with Bayesian credible regions, given a model for $M(<r)$ and samples drawn from a posterior distribution. With the CMP, we can report median estimates and 50\% Bayesian credible regions for the MW mass within any distance (e.g. $M(r=25\text{kpc})= 0.26~(0.24, 0.29)\times 10^{12} M_{\odot}$, $M(r=50\text{\kpc})= 0.37~(0.34, 0.41) \times 10^{12} M_{\odot}$, $M(r=100\text{kpc}) = 0.53~(0.49, 0.58) \times10^{12} M_{\odot}$, etc), making it easy to compare our results directly to other studies.


INTRODUCTION
Since the Gaia Data Release 2 (DR2) in April 2018, a number of studies have presented new estimates or lower bounds for the mass of the Milky Way (MW) Galaxy (e.g. Helmi et al. 2018;Malhan & Ibata 2018;Posti & Helmi 2018;Sohn et al. 2018b;Watkins et al. 2018;Vasiliev 2018). 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 200 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 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 15kpc. 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 com-pleteness of tracer data beyond 15kpc, but the HST-PROMO project has helped contribute valuable measurements for GCs at large distances (Sohn et al. 2018b). Additionally, Vasiliev (2018) estimated the mean proper motion of 150 GCs using the Gaia DR2 dataset, thereby increasing the number of proper motion measurements at larger distances considerably. The catalog provided by Vasiliev is appealing not only because the data is 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 (2018). 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 1-3, 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 1-3 and the methods used in this study. In Section 2.1, we describe how to calculate the CMP for the Milky Way given a function M (R < r), and given samples from the posterior distribution of model parameters.
• Section 3 discusses the details and differences between the two data sets (Gaia DR2 + HST-PROMO, 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.

METHOD
Specific details about our hierarchical Bayesian model and sampling methods for the posterior distribution can be found in Papers 1-3. 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 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 simulataneously, 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 (µ α cos δ, µ δ ), 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 los , µ α cos δ, and µ δ is a distribution function (DF) of the specific energies E and angular momenta L, given a model for the total gravitational potential Φ(r) and density profile of the tracer population ρ(r). In Papers 2 and 3, we used a simple power law profile for both the gravitational potential, and the density profile of the tracer population, where Φ • , γ, and α are parameters. Together, equations 1 and 2, determine the analytic DF f (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 2 (although notations differ). Hyperprior distributions are also defined for the four model parameters: Φ • , γ, α, and β. In Papers 2-3, 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 the 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 R Statistical Software environment (see Papers 1-3). 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 globular cluster 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 stellarmass-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 mass (M 200 ) 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 towards 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 3 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 3 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 > 15 kpc.
EKH also found that using a uniform distribution as a prior on γ may have been too broad. The uniform distibution 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). 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μ measurements (Figure 2). Because we are excluding GCs within 15kpc, there are only three GCs rith r > 15 kpc that are lacking both line-of-sight and proper motion measurements: AM 4, Ko 1, and Ko 2.

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. 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 (Deason et al. 2012b). Thus, for a Markov chain of length m, each (Φ •,j , γ 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 kpc to r 100 = 200 kpc. Our chains are of length 90000, with effective sample sizes of > 2000 for all parameters. 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. 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, 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. 2010de la Fuente Marcos et al. 2015;Rossi et al. 2015).
The total number of GCs in our compiled data set is 157, but four 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.
Together, the Gaia collaboration and the HST-PROMO 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 (Papers 2-3), in which approximately 50% of the data were incomplete. Even more recently, Vasiliev (2018) 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 found by Vasiliev. The Vasiliev catalog contains some new line-of-sight velocity measurements that were previously missing in the Harris Catalog, and so we adopt these values.
The Vasiliev catalog also contains information for the GC Leavens 1 (Crater), but we exclude this data point because it is 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 2018). Although we do not include the cluster 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 + HST-PROMO and Vasiliev catalogs is completeness, with the latter containing complete measurements for 143 GCs, and incomplete measurements for 11 GCs (Figure 4).
As mentioned in Section 2, we use GCs with Galactocentric distances r > 15 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μ measurements, three 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 hyperparamters (see Section 2, Figure 2, and Papers 2-3). Figure 5 shows the proper motion components and line-of-sight velocities as a function of Galactocentric distance for GCs at r > 15kpc. The Gaia+HSTPROMO measurements are shown with open blue circles and the Vasiliev's measurements are grey 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.7kpc. 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 > 15 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. 2012a,c;Kafle et al. 2012;Gibbons et al. 2014;Eadie et al. 2015;Küpper et al. 2015), including five which also use Gaia DR2 data (Malhan & Ibata 2018;Posti & Helmi 2018;Sohn et al. 2018a;Watkins et al. 2018;Vasiliev 2018).

RESULTS AND DISCUSSION
The two CMPs in Figure 6 are nearly identical, but the Vasiliev catalog provides a better constraint on the mass because the data are more complete. The median estimates of M 200 given each catalog are also extremely similar at 0.70 (0.47, 1.14)×10 12 M (Gaia + HSTPROMO) (4) and 0.71 (0.52, 1.07) × 10 12 M (Vasiliev), 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 20kpc of the Galactic centre ( 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 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 kpc through r = 150 kpc, by steps of 25kpc, to show the marginal posterior distributions for (M < r) at those distances ( Figure 7). Marginal distributions for the mass within any r are also easily calculated. 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 100kpc, Eadie et al. (2015) at 125kpc, and Vasiliev (2018) at 80kpc. 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 100kpc.
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 200 or the virial mass M vir . Instead, we compare three of these studies' results to our marginal distribution for M 200 given the Vasiliev catalog (Figure 8). Patel et al. (2017) performed a Bayesian analysis in conjunction with cosmological simulations of MWtype galaxies to infer the Galaxy's mass. They used the kinematics of one massive satellite, and inferred M 200 = 1.223 × 10 12 M and a 95% Bayesian credible of (1.207, 1.248) × 10 12 M . 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 ≈ 1.4 × 10 12 M . As a final example, Monari et al. (2018) used counterrotating halo stars in Gaia to estimate the escape velocity of the MW between 5-10kpc. Their results also implied a heavy MW with a dark matter virial mass of M 200 = 1.55 +0.64 −0.51 × 10 12 M . All three of these mass estimates are almost completely ruled out by our estimate of M 200 (Figure 8), which did not include high/extreme velocity stars or dwarf galaxies.
It is less clear why the estimate at 80kpc by Vasiliev (2018) is significantly higher than our estimate. 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 15kpc. 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.5kpc 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 find M (r < 39.5 kpc) = 0.61 +0.18 −0.12 × 10 12 M . Our estimate for the mass within 39.5kpc is M (r < 39.5kpc) = 0.33 (0.26, 0.44) × 10 12 M , (6) which does not agree with their results within the uncertainties. However, more recently Watkins et al. (2018) combined Gaia DR2 data with HSTPROMO data, and used the same estimator as Sohn et al. (2018b). They found M (r < 39.5 kpc) = 0.44 0.07 −0.06 × 10 12 M . 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 (2018) also use 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 use a two-component DF for the GC system of the MW. The assumption is that the GC population dynamics can be decomposed into disklike and halo-like components, with the latter containing parameters for the mean rotational velocity and velocity anisotropy of the GC system. They also include a shape parameter for the dark matter halo. Because their study uses a Bayesian paradigm, like ours, their method includes measurements uncertainties and incomplete data directly in the analysis. Despite using a more complicated DF than this study, their mass estimate within 20kpc (M (r < 20kpc) = 0.191 +0.17 −0.15 × 10 12 M ) is within the bounds of our 95% Bayesian credible regions: M (r < 20 kpc) = 0.24 (0.18, 0.32) × 10 12 M .
Vasiliev (2018) also estimates the mass of the MW using their GC catalog in a maximum likelihood analysis. Their analysis assumes a single-component DF for the GC population, but a three-component gravitational potential to account for the bulge, disc, and halo components. They find M (r < 50 kpc) = 0.6 +0.14 −0.09 × 10 12 M The uncertainties almost overlap with the upper bound of the 95% Bayesian credible for our mass estimate within 50kpc: given by our method.
Overall, it appears that studies are reaching a consensus on the mass of the MW within 40kpc 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.
Mass estimates of the Galaxy out to larger distances, however, are still rather uncertain. For example, Sohn et al. (2018b) and Watkins et al. (2018)

CONCLUSIONS & FUTURE WORK
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 40kpc, and in general they are in good agreement with one another. Our CMP of the MW within r < 50 kpc is also in agreement with many recent results.
Our CMP overlaps more than half of the mass esti-mates 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 HST-PROMO data have greatly improved our understanding of the mass within 40kpc 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 author would like to extend graitude to Dr. Vasiliev for sharing their catalog of new proper motions measurements.