Evidence for an Intermediate-mass Milky Way from Gaia DR2 Halo Globular Cluster Motions

, , , and

Published 2019 March 12 © 2019. The American Astronomical Society. All rights reserved.
, , Citation Laura L. Watkins et al 2019 ApJ 873 118 DOI 10.3847/1538-4357/ab089f

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/873/2/118

Abstract

We estimate the mass of the Milky Way (MW) within 21.1 kpc using the kinematics of halo globular clusters (GCs) determined by Gaia. The second Gaia data release (DR2) contained a catalog of absolute proper motions (PMs) for a set of Galactic GCs and satellite galaxies. We select from the catalog only halo GCs, identifying a total of 34 GCs spanning $2.0\leqslant r\leqslant 21.1\,\mathrm{kpc}$, and use their 3D kinematics to estimate the anisotropy over this range to be $\beta ={0.46}_{-0.19}^{+0.15}$, in good agreement with, though slightly lower than, a recent estimate for a sample of halo GCs using Hubble Space Telescope (HST) PM measurements further out in the halo. We then use the Gaia kinematics to estimate the mass of the MW inside the outermost GC to be $M(\lt 21.1\,\mathrm{kpc})={0.21}_{-0.03}^{+0.04}\times {10}^{12}\,{M}_{\odot }$, which corresponds to a circular velocity at ${r}_{\max }$ of ${v}_{\mathrm{circ}}(21.1\,\mathrm{kpc})={206}_{-16}^{+19}\,\mathrm{km}\,{{\rm{s}}}^{-1}.$ The implied virial mass is ${M}_{\mathrm{virial}}={1.28}_{-0.48}^{+0.97}\times {10}^{12}\,{M}_{\odot }$. The error bars encompass the uncertainties on the anisotropy and on the density profile of the MW dark halo, and the scatter inherent in the mass estimator we use. We get improved estimates when we combine the Gaia and HST samples to provide kinematics for 46 GCs out to 39.5 kpc: $\beta ={0.52}_{-0.14}^{+0.11}$, $M(\lt 39.5\,\mathrm{kpc})={0.42}_{-0.06}^{+0.07}\times {10}^{12}\,{M}_{\odot }$, and ${M}_{\mathrm{virial}}={1.54}_{-0.44}^{+0.75}\times {10}^{12}\,{M}_{\odot }$. We show that these results are robust to potential substructure in the halo GC distribution. While a wide range of MW virial masses have been advocated in the literature, from below 1012 ${M}_{\odot }$ to above 2 × 1012 ${M}_{\odot }$, these new data imply that an intermediate mass is most likely.

Export citation and abstract BibTeX RIS

1. Introduction

The mass of the Milky Way (MW) is one of its most fundamental parameters, and yet, despite decades of intense effort, our best estimates are significantly scattered, with some estimates agreeing very well, and others differing by more than their uncertainties (see Bland-Hawthorn & Gerhard 2016 for a thorough review). These estimates are very sensitive to assumptions made in the modeling, including, but not limited to, which of the MW's satellites are bound or unbound and for how long they have been bound, the velocity anisotropy of the MW halo and of its satellite system, the shape of the MW halo, and the particular method used for the analysis. Estimates typically range from as low as ∼0.5 × 1012 ${M}_{\odot }$ (e.g., Watkins et al. 2010; radial anisotropy with Leo I unbound) to as high as (2–3) × 1012 ${M}_{\odot }$ from abundance-matching studies (e.g., Boylan-Kolchin et al. 2010), the timing argument (Li & White 2008; van der Marel et al. 2012b), or studies of tracers (e.g., Watkins et al. 2010; tangential anisotropy with Leo I bound).

Accurate determination of the mass profile of the MW has implications for our understanding of the dynamical history of the Local Group (both past evolution and future interactions; e.g., van der Marel et al. 2012a, 2012b), the MW's satellite population, particularly the Sagittarius dwarf spheroidal (dSph) and its impressive tidal stream (Fardal et al. 2019), and the Magellanic Clouds (Kallivayalil et al. 2013).

Furthermore, the mass of a galaxy and its distribution (or shape) are intrinsically linked to the formation and growth of structure in the universe (Conselice 2014), so accurately determining these parameters for the MW will give us a clearer understanding of where our Galaxy sits in a cosmological context (for an excellent review see Freeman & Bland-Hawthorn 2002). In particular, we can know whether the MW is typical or atypical, and thus how much of what we learn about the MW can be safely assumed for other galaxies as well.

The MW is composed of a central nucleus which harbours a supermassive black hole at its heart, a bulge, a disk, and a halo (Ivezić et al. 2012; Bland-Hawthorn & Gerhard 2016). The first three components are all primarily baryonic in nature, and while many of their properties remain topics of some debate, their masses are reasonably well determined. The final component, the halo, is dominated by dark matter (DM)—only a few percent of the mass of the halo is baryonic (Helmi 2008), the exact percentage depends on the unknown total mass of DM in the halo—and it is our inability to see DM directly that gives rise to our present uncertainty in the mass.

As the majority of mass in the MW is "invisible," we cannot measure it directly; instead we can infer its presence by its influence on its surroundings. Typically, this is the purview of dynamical studies. Any mass distribution gives rise to a gravitational potential that causes objects to move: by studying measurements of the motions of the objects, we can work backwards to recover the underlying gravitational potential and, thus, the mass distribution.

There are some mass-estimation methods, notably the timing argument and abundance-matching studies, that estimate the total mass of a system. However, most dynamical methods work by using tracer objects to probe the properties of the whole system, and can only estimate the mass over the range for which tracer data are available. Thus different families of tracers provide crucial information at different points depending on the range they cover. This is particularly crucial in the MW where globular clusters (GCs) tend to probe the inner regions of the halo, while dSph satellite galaxies offer better coverage further out (e.g., Wilkinson & Evans 1999; Watkins et al. 2010; Patel et al. 2018). The modeling of multiple stellar streams may provide a promising alternative (e.g., Gibbons et al. 2014; Sanderson et al. 2017), though six-dimensional phase space data are only available for the GD-1 stream at present (Koposov et al. 2010; Bowden et al. 2015).

One key problem with mass estimation via kinematics is that we need to know the total velocity of each tracer, but we are seldom fortunate enough to have all three components of motion for a large sample of tracers. Typically, we only have line-of-sight (LOS) velocities. This is especially troublesome for studies of the MW as the Sun is very close to the Galactic center, and so for most objects LOS velocities predominantly probe only one component of the motion (the Galactocentric radial direction) and offer little information about the Galactocentric tangential motions of the tracers.4 With only LOS velocities, the masses we estimate depend very strongly on what assumptions we make for the tangential motions: the well-known mass–anisotropy degeneracy (e.g., Binney & Tremaine 2008).

Some methods attempt to overcome the lack of 3D velocity information; Eadie et al. (2015) introduced a Galactic mass estimator that includes unknown velocity components as free parameters in their models. However, the best constraints on mass will come from having complete phase space information and, with proper motions (PMs), we are able to break this degeneracy. First, we can make a direct estimate of the velocity anisotropy and, second, we can correctly include the total velocity of the tracers in our mass calculations instead of having to make assumptions.

Absolute PMs have been measured from the ground for a number of GCs (e.g., Casetti-Dinescu et al. 2013 and other papers in the series), although typically this is only possible with sufficient accuracy for objects within ∼10 kpc, and even then ground-based measurements often suffer from a number of systematic effects. Space offers a more stable environment for astrometry, so thanks to its excellent precision and long time baseline, the Hubble Space Telescope (HST) has proved extremely valuable for providing absolute PMs for dSphs (e.g., Piatek et al. 2016; Sohn et al. 2017) and GCs (e.g., Anderson & King 2003; Kalirai et al. 2007). The recent study by Sohn et al. (2018) that measured PMs using the HST for 20 outer halo GCs in the MW represents the largest sample of absolute PMs measured to date in a single study.

The Gaia mission's (Gaia Collaboration et al. 2016b) first data release (Gaia Collaboration et al. 2016a) contained PMs for around two million stars in the TychoGaia Astrometric Solution (TGAS) catalog (Michalik et al. 2015), which used Tycho2 (Høg et al. 2000) measurements to provide a first epoch of data and Gaia data for the second epoch, and has been used already for multiple PM studies of objects in the MW, including for a handful of Galactic GCs. Watkins & van der Marel (2017) identified member stars for five GCs in the TGAS data and used the stars to estimate the absolute PMs of their host clusters; comparing these Gaia PMs with previous estimates, they found excellent agreement with previous HST measurements, but some differences from previous ground-based values due to systematics inherent in such measurements. Massari et al. (2017) used archival HST data combined with Gaia DR1 data to estimate the PM of GC NGC 2419.

However, the second Gaia data release (Gaia Collaboration et al. 2018a) has greatly expanded our view of the local universe. This data release provides PMs for billions of stars and has made it possible to measure absolute PMs for 75 Galactic GCs out to a Galactocentric distance of ∼21 kpc, along with nine classical dSphs, a single ultrafaint dwarf, and both the Large and Small Magellanic Clouds (Gaia Collaboration et al. 2018b). This is by far the largest catalog of GC and dSph PMs to date. Combined with position and LOS velocity information from previous studies, these measurements have enabled analysis of the orbits of these objects.

In this paper, we use these motions to provide new mass estimates for the MW. In Section 2, we describe the Gaia catalogs, calculate the Galactocentric motions of the objects, and describe which objects we select for our analysis; in Section 3, we estimate the mass of the MW; in Section 4, we compare our results with previous estimates; in Section 5, we summarize our findings.

2. Overview of Data

2.1. Gaia Halo Cluster Sample

Gaia Collaboration et al. (2018b) used the the Gaia DR2 Catalog (Gaia Collaboration et al. 2018a) to identify member stars for a number of MW GCs and dSphs based on their positions, photometry, and PMs, from which their absolute PMs were calculated.5 Combined with distances and LOS velocities, it was then possible to calculate the orbit of each object within the Galaxy. The orbits derived do depend on certain assumptions made for the potential of the Galaxy which, as we have discussed, is still somewhat uncertain. To mitigate the effects of this uncertainty, orbits were calculated in three different potentials that span a range of possible MW shapes and masses. We will use both the absolute PMs and the orbital properties here.

The first step is to calculate Galactocentric positions and motions from the observed heliocentric values. We take right ascensions, declinations, and the PMs in these coordinates, along with the full covariance matrix for the PM uncertainties from Gaia Collaboration et al. (2018b).6 Distances and LOS velocities we take from Harris (1996, 2010 edition), which are primarily determined photometrically; these distances generally agree with kinematical distances determined from HST PM studies (Watkins et al. 2015). We assume distance uncertainties ΔD = 0.023 D, which are equivalent to uncertainties on the distance moduli of 0.05 mag, which is typical for GC distance uncertainties (Dotter et al. 2010). There are a few GCs in the Harris catalog for which there is an LOS velocity measurement but no uncertainty listed. For the GCs with uncertainties, the average is $\overline{{\rm{\Delta }}{v}_{\mathrm{LOS}}}=0.06{v}_{\mathrm{LOS}}$, which we thus adopt for the remaining GCs. We assume a distance from the Sun to the Galactic center of R = 8.29 ± 0.16 kpc and a circular velocity at the solar radius of V = 239 ± 5 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (McMillan 2011). For the solar peculiar velocity relative to the local standard of rest, we assume ${\boldsymbol{V}}$pec = (11.10 ± 1.23, 12.24 ± 2.05, 7.25 ± 0.62) $\mathrm{km}\,{{\rm{s}}}^{-1}$ (Schönrich et al. 2010). We will use these solar parameters throughout the paper.

We calculate the positions and velocities of the GCs in a spherical coordinate system (radius r, latitude θ, and longitude ϕ) centered on the Galactic center. We use Monte Carlo sampling, using 1000 samples and assuming Gaussian uncertainties, to propagate all of the observational uncertainties; this includes the full covariance matrix for the PMs, uncertainties of the distances and LOS velocities, and all the uncertainties on the position and velocity of the Sun. Although the initial distributions are assumed Gaussian, the resulting distributions of Galactocentric properties may not be, so we take medians and 15.9 and 84.1 percentiles of the distributions as the best estimate and uncertainties.7 These Galactocentric positions and velocities are provided in the Appendix.

One of the main contributors to the uncertainty in the mass of the MW is the paucity of tracer objects. The dSphs are limited in number and appear to have been accreted in groups (Gaia Collaboration et al. 2018b), which is an extremely interesting result but makes mass modeling tricky as certain key assumptions are thus invalidated. For the GCs, on the other hand, the improvements in PM accuracy offered by Gaia over HST (where such data exist) are modest, with only 22 months of Gaia data. Where Gaia can play a key role here is to provide Galactocentric motions for many more clusters than were previously available, greatly increasing our sample size. As such, we choose to concentrate on the GC sample henceforth.

To probe the anisotropy and mass of the Galactic halo, we require a sample of halo clusters, free of disk and bulge contaminants. Zinn (1993) showed that the disk and bulge clusters separate cleanly from the halo clusters in metallicity. We use the same cut and keep only clusters with [Fe/H] ≤ −0.8 dex (metallicities from Harris 1996, 2010 edition) so as to have a pure halo sample. Furthermore, the mass estimators we will use assume that the potential is scale-free over the region of interest, so we wish to limit ourselves to clusters for which this assumption reasonably holds, that is we do not wish to include clusters that spend most of their time in the innermost regions of the Galaxy where the disk is a significant contributor and the potential is non-spherical. We use the orbital parameters from Gaia Collaboration et al. (2018b) to extract only GCs with apocenters8 rapo ≥ 6 kpc;9 this leaves us with 34 GCs that span a radial range $2.0\leqslant r\leqslant 21.1\,\mathrm{kpc}$.

2.2. HST Halo Cluster Sample

Recently, Sohn et al. (2018) presented HST PMs for 20 halo GCs that extend further out into the halo than the Gaia cluster sample. Four of these are in common with the Gaia Collaboration et al. (2018b) sample: NGC 2298, NGC 5024, NGC 5053, and NGC 5466. The HST estimates are in good agreement with the Gaia estimates for three of the four clusters, which is reassuring news for both catalogs. The measurements for the fourth cluster differ by ∼48 $\mathrm{km}\,{{\rm{s}}}^{-1}$ but this is still well below the velocity dispersion of the halo (see Table 1). This also indicates that we can confidently combine the catalogs and increase the size of our sample and its range. This is an improvement on both analyses as the Gaia catalog probes further in and the HST catalog probes further out, but they also have a substantial region of overlap to serve as a solid anchor and consistency check.

Table 1.  Galactocentric Mean Velocities, Dispersions, and Anisotropies for the GC Subsamples

Sample $\overline{{v}_{{\rm{r}}}}$ $\overline{{v}_{\theta }}$ $\overline{{v}_{\phi }}$ σr σθ σϕ β
  ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$)
Aa $-{15.9}_{-24.0}^{+23.4}$ ${22.1}_{-14.9}^{+15.6}$ $-{29.8}_{-19.0}^{+18.6}$ ${138.9}_{-15.8}^{+20.1}$ ${91.3}_{-10.7}^{+13.0}$ ${109.5}_{-12.4}^{+15.9}$ ${0.46}_{-0.19}^{+0.15}$
Bb $-{26.7}_{-23.5}^{+23.2}$ ${20.8}_{-14.5}^{+14.8}$ $-{18.8}_{-16.3}^{+16.0}$ ${153.3}_{-14.9}^{+18.0}$ ${100.7}_{-10.3}^{+12.0}$ ${110.4}_{-11.0}^{+13.8}$ ${0.52}_{-0.14}^{+0.11}$
Cc $-{21.5}_{-26.2}^{+26.9}$ ${22.2}_{-21.0}^{+21.0}$ $-{41.0}_{-24.2}^{+24.0}$ ${132.9}_{-16.7}^{+22.6}$ ${103.7}_{-13.4}^{+17.2}$ ${124.0}_{-15.5}^{+21.2}$ ${0.24}_{-0.31}^{+0.23}$
Dd $-{27.1}_{-25.0}^{+23.9}$ ${7.9}_{-11.7}^{+12.1}$ $-{27.5}_{-12.5}^{+11.9}$ ${159.1}_{-16.4}^{+19.8}$ ${76.4}_{-8.3}^{+9.8}$ ${78.4}_{-8.4}^{+10.0}$ ${0.76}_{-0.08}^{+0.06}$

Notes.

aGaia GCs. bGaia and HST GCs. cGaia GCs with possible merger group removed. dGaia and HST GCs with high ${v}_{\tan }$ clusters removed.

Download table as:  ASCIITypeset image

We follow our approach from Sohn et al. (2018) and exclude NGC 2419 as its distance is much greater than the rest of the sample,10 and three of the four clusters associated with the Sagittarius dSph as they represent a group, not a well mixed population (and including all four can again lead to biases). We also choose to use the Gaia values for the clusters in common. Overall, this gives us an extra 12 clusters in our sample, bringing the total to 46, and increasing the radial range of the sample to 2.0 ≤ r ≤ 39.5 kpc.

In what follows, we will provide analysis using only the Gaia cluster sample (Sample A), and the combined Gaia and HST cluster samples (Sample B).

Figure 1 shows the distribution of the Galactocentric velocities of the halo GC sample as a function of Galactocentric distance. The upper three panels show the radial vr, latitudinal vθ, and longitudinal vϕ velocity components. The next panel shows the tangential velocity where ${v}_{{\rm{t}}}=\sqrt{{v}_{\theta }^{2}+{v}_{\phi }^{2}}$, and the bottom panel shows the total velocity v. The red circles, orange diamonds, and green stars show the halo GCs in the Gaia DR2 sample. The blue triangles and cyan squares show the halo GCs from HST measurements. The different groupings identify various subsamples that we will use later to verify the robustness of results against substructure.

Figure 1.

Figure 1. Galactocentric velocity distributions as a function of Galactocentric distance of the halo GCs, from top to bottom: Galactocentric radial velocities, Galactocentric latitudinal velocities, Galactocentric longitudinal velocities, Galactocentric tangential velocities, Galactocentric total velocities. In all panels, red circles, orange diamonds, and green stars show the Gaia measurements and the cyan squares and blue triangles show the HST measurements. Orange diamonds (Gaia) and cyan squares (HST) highlight the GCs with ${v}_{\tan }$ ≥ 250 $\mathrm{km}\,{{\rm{s}}}^{-1}$, and the green stars highlight the clusters tentatively identified as part of a recent merger—these subsamples are discussed in later sections. Uncertainties are shown in the figure but are too small to be visible in most cases.

Standard image High-resolution image

3. MW Mass

We can use the Galactocentric positions and motions of the clusters to estimate the mass of the MW within the radius of the outermost cluster in our sample, which lies at ${r}_{\max }$ = 21.1 kpc. Watkins et al. (2010) introduced a family of simple, yet effective, tracer mass estimators (TMEs), which we will use here. Subsequent extensions and applications of this method were presented in Annibali et al. (2018) in a study of NGC 4449, and in Sohn et al. (2018) in a similar study of MW halo GCs.

The estimators work with different types of distance and velocity data, depending on what is available. As we have full 6D phase-space information for our cluster sample, we are able to use the estimator that uses distances and total velocities (Equation (24) of Watkins et al. 2010). That is,

Equation (1)

The estimators assume that the underlying potential is a power law with index α over the region of interest, that the tracer objects have a number density distribution that is a power law with index γ over the region of interest, and that the velocity anisotropy of the tracer sample is a constant β over the region of interest. Before we can proceed, we need to estimate α, β, and γ.

3.1. Anisotropy

In Section 2, we calculated Galactocentric motions vj of all the clusters in our sample in a spherical coordinate system (j, k = {r, θ, ϕ }), along with a full covariance matrix for their uncertainties δj and the correlations between the uncertainties ρjk. That is, for each cluster i, we have velocities

Equation (2)

with uncertainties

Equation (3)

We assume that the cluster population has a mean velocity

Equation (4)

and covariance

Equation (5)

where $(\overline{{v}_{{\rm{r}}}},\overline{{v}_{\theta }},\overline{{v}_{\phi }})$ are the mean velocities and (σr, σθ, σϕ) the velocity dispersions for each coordinate. In setting the cross-terms of the covariance to zero, we have assumed that the axes of the velocity ellipsoid are aligned with the spherical coordinate system. We further assume that the velocity distributions are Gaussian. Then the likelihood of the observed measurements for mean $\overline{{\boldsymbol{v}}}$ and covariance ${\boldsymbol{C}}$ is

Equation (6)

We use flat priors for the mean velocities in each coordinate, that is

Equation (7)

We insist that the dispersions must be positive, but otherwise use a flat prior for positive dispersion values, that is

Equation (8)

Finally, the posterior is the product of the likelihood and the priors.

To estimate the means and dispersions that best describe the data, we use the affine-invariant Markov Chain Monte Carlo package emcee (Foreman-Mackey et al. 2013) to find the region of parameter space where the posterior is maximized and to sample that region. We draw 10,000 points from the final posterior distribution for our final sample. For each parameter, we adopt the median as the best estimate and use the 15.9 and 84.1 percentiles for the uncertainties.

From these, we are then able to estimate the anisotropy

Equation (9)

of the system. For each of the 10,000 points in our final sample, we calculate the anisotropy β. The distribution of these anisotropy values is shown in Figure 2, and has a median and 15.9- and 84.1-percentile uncertainties of $\beta ={0.46}_{-0.19}^{+0.15}$. We discuss this value in the context of previous work in Section 4.

Figure 2.

Figure 2. Posterior distribution of anisotropy β for the Gaia DR2 halo cluster sample.

Standard image High-resolution image

Table 1 summarizes our fits to the velocity ellipsoid of the halo. We provide estimates of the mean velocities and velocity dispersions for each velocity component and the inferred anisotropy. The Gaia sample is Sample A in the table. The other samples are described later in the text. We see that the Gaia sample has a mean radial velocity consistent with 0 within its uncertainty, but shows hints of net tangential motion, which we will address later.

3.2. Density

To estimate the power-law index γ of the halo cluster number density profile, we start with Galactocentric distances from the Harris (1996, 2010 edition) catalog of 157 Galactic GCs. This catalog contains both clusters that move on bulge-like and disk-like orbits and are found in the inner regions of the Galaxy, and clusters that move on halo-like orbits and are found further out. All the clusters in our PM sample were deliberately chosen to be part of the halo cluster population, so that is the number density profile of interest for this analysis. We use a least-squares fitting algorithm to fit (with uncertainties) a broken power law to the data that has an index ${\gamma }_{\mathrm{in}}$ in the inner regions and an index ${\gamma }_{\mathrm{out}}$ in the outer regions, with a break radius of ${r}_{\mathrm{break}}$. We assume that ${\gamma }_{\mathrm{in}}$ describes the density profile of the bulge and disk clusters and that ${\gamma }_{\mathrm{out}}$ describes the density profile of the halo clusters.

Figure 3 shows the cumulative number density profile of all MW GCs in blue. The solid black line is the best-fitting broken power law to the data, and the dashed black line marks the break radius at which the power-law index changes. The parameters of the fit are shown in the top-left corner, with uncertainties that are purely statistical and do not encompass systematic effects. As we assume that the outer profile describes the halo clusters, we adopt this value γ = 3.53 ± 0.01 for our mass analysis.11 This is the same value we adopted in a recent study of halo GCs with HST PMs by Sohn et al. (2018) and agrees well with previous studies (e.g., Harris 2001).

Figure 3.

Figure 3. Power-law fit to the cumulative number profile of Milky Way (MW) globular clusters. To this we fit a broken power-law model that has an index γin in the inner regions, which we assume to be dominated by disk clusters, and an index γout in the outer regions, which we assume to be dominated by halo clusters, and is the value of interest for our analysis. The cumulative histogram in blue shows the data, the solid black line shows the best-fitting broken power law, and the dashed black line marks the break radius. The parameters for the best-fitting broken power law are shown in the top-left corner of the figure. The index for the outer power law is the value of interest and is highlighted in blue. The uncertainty in this value is the uncertainty in the fit, and does not account for uncertainty in the estimated positions of the MW clusters or the functional form of the fit. For guidance, we also mark as blue dotted lines the radius of the innermost cluster in the Gaia sample at 2.0 kpc, the outermost clusters in the Gaia sample at 21.1 kpc, and the outermost cluster in the HST sample at 39.5 kpc.

Standard image High-resolution image

3.3. Potential

To estimate the power-law index α of the gravitational potential, we assume that the MW is composed of a nucleus, bulge, disk, and halo. We adopt the nucleus, disk, and bulge prescriptions from Price-Whelan (2017), that is, we assume a Hernquist nucleus with mass ${M}_{\mathrm{nucleus}}=1.71\times {10}^{9}$ ${M}_{\odot }$ and scale length lnucleus = 0. 07 kpc, a Hernquist bulge with mass ${M}_{\mathrm{bulge}}=5\times {10}^{9}$ ${M}_{\odot }$ and scale length lbulge = 1 kpc, and a Miyamoto–Nagai disk ${M}_{\mathrm{disk}}=6.8\times {10}^{10}$ ${M}_{\odot }$, scale length ldisk = 3 kpc and scale height hdisk = 0. 28 kpc.

The shape and mass of the halo is uncertain—indeed, this uncertainty is the primary motivation for the current analysis—and so we choose to sample a range of possible halos to see what α values they imply. To do this, we assume that the halo is spherical and NFW (Navarro et al. 1996) in form, but with unknown virial radius ${r}_{\mathrm{virial}}$ and a concentration c, which we will sample.12

The oft-cited study by Klypin et al. (2002) favors an NFW halo with virial mass ${M}_{\mathrm{virial}}$ = 1 × 1012 ${M}_{\odot }$ and scale radius ${r}_{\mathrm{scale}}$ = 21.5 kpc (virial radius ${r}_{\mathrm{virial}}$ = 258 kpc and a concentration c = 12). More recent halo estimates have been less massive: Bovy (2015) favors an NFW halo with virial mass ${M}_{\mathrm{virial}}$ = 0.8 × 1012 ${M}_{\odot }$ and scale radius ${r}_{\mathrm{scale}}$ = 16.0 kpc (virial radius ${r}_{\mathrm{virial}}$ = 245 kpc and a concentration c = 15.3), and Price-Whelan (2017) favors an NFW halo with virial mass ${M}_{\mathrm{virial}}$ = 0.54 × 1012 ${M}_{\odot }$ and scale radius ${r}_{\mathrm{scale}}$ = 15.62 kpc (virial radius ${r}_{\mathrm{virial}}$ = 214 kpc and a concentration c = 13.7). However, other studies have found somewhat larger total masses for the MW (e.g., Watkins et al. 2010) or concentrations (e.g., Deason et al. 2012) than these halo parameters would imply. As such, we choose to sample a range of virial radii 200 ≤ ${r}_{\mathrm{virial}}$ ≤ 400 kpc, sampled at 1 kpc intervals, and a range of concentrations 8 ≤ c ≤ 20, sampled at 0.1 intervals.

To further narrow down the list of allowed halos, we insist that the circular velocity at the solar radius must be consistent with observed values. Estimates for the circular velocity typically span 220–250 $\mathrm{km}\,{{\rm{s}}}^{-1}$; for each halo, we estimate the circular velocity for the best-fitting power-law model at R and reject halos with velocities outside of this range.

For each halo model, we calculate the total potential profile from the nucleus, bulge, disk, and halo components, and use a least-squares fitting algorithm to fit a power law across the range spanned by our cluster sample $2.0\leqslant r\leqslant 21.1\,\mathrm{kpc}$. The index of the power-law fit α is the quantity we require for our models. Figure 4 shows the range of potentials sampled and the range of the power-law fits to those halos.

Figure 4.

Figure 4. Range of power-law fits to the MW potential. We assume that the MW consists of a Hernquist nucleus (orange), a Hernquist bulge (yellow), a Miyamoto–Nagai disk (light green), and an NFW halo (dark green). Their sum (the total potential) is shown in gray. We assume that the nucleus, bulge, and disk components are fixed. The disk potential shows some slight broadening around 2 kpc as we have plotted potential as a function of spherical radius but the disk is not spherical. We sample a range of NFW halo parameters so, for the halo and total potentials, the solid lines show the median profiles, the dark shaded regions show the range between the 25th and 75th percentiles of the profiles, and the light shaded regions show the range of profiles for all halos sampled. The dotted lines show the region $2.0\leqslant r\leqslant 21.1\,\mathrm{kpc}$ spanned by our halo cluster sample. The solid black lines show the extent of the best-fitting power laws in this region, and the dashed lines show the extent of the best-fitting power laws extended outside of the region of interest.

Standard image High-resolution image

Figure 5 shows the variation in α (upper panel) and ${v}_{\mathrm{circ}}$ (lower panel) across our halo sample, as indicated by the respective color bars. The white regions are halos with circular velocities inconsistent with observations. The corresponding distribution of α values is shown in Figure 6.

Figure 5.

Figure 5. Variation in potential slope (upper panel) and circular velocity at the solar radius (lower panel) for a grid of halo models assumed to be NFW in shape and defined by a concentration c and a virial radius rvir. We only show halos with circular velocities at the solar radius R = 8.29 kpc of 220 ≤ ${v}_{\mathrm{circ}}$ ≤ 250 $\mathrm{km}\,{{\rm{s}}}^{-1}$ to force consistency with observations. The white regions are halos with circular velocities outside of this range.

Standard image High-resolution image
Figure 6.

Figure 6. Distribution of α values fitted to the grid of sample halos; halos for which the circular velocity at the solar radius is inconsistent with observations have been removed.

Standard image High-resolution image

3.4. Monte Carlo Simulations

The variation in the mass estimates for different halos assesses our uncertainty in the particular density profile of the halo, but does not assess how well the mass estimators themselves are able to recover the true mass of the MW using 34 tracers drawn from the underlying distribution. We do this using the set of Monte Carlo simulations described in Watkins et al. (2010). Briefly, the simulations create a set of tracer objects drawn from a power-law density with index γ, with velocities consistent with a power-law potential with index α and with anisotropy β. We then use the TME to estimate the mass ${M}_{\mathrm{TME}}$ within ${r}_{\max }$, and compare this estimated mass with the known true mass of the simulation ${M}_{\mathrm{true}}$ within ${r}_{\max }$.

For these simulations, we use density slope γ = 3.53 (as calculated in Section 3.2), and anisotropy β = 0.46 (the median value found in Section 3.1). For the potential slope α, we select a "typical" halo using the halo grid described in Section 3.3: we choose to use ${r}_{\mathrm{virial}}$ = 300 kpc (the middle value of the range sampled), and then use the grid to find the value of concentration c for which the circular velocity at the solar radius is closest to the observed value of V, which results in c = 16.8. We use the value of α = 0.26 calculated for this halo. We further use the observed solar radius and the circular velocity at the solar radius for the fiducial radius and fiducial velocity in the potential power law. Although our analysis does use a range of values for both α and β, using fixed values for our simulations is sufficient for our purposes here, which is to assess how far our mass estimate based on a single cluster sample may be from the true value, as we do not expect the performance of the TME with 34 tracers to change with α or β.

We generate 1000 simulations of 34 clusters, from which we find that the ratio of the estimated mass to the true mass $f={M}_{\mathrm{TME}}/{M}_{\mathrm{true}}={1.01}_{-0.16}^{+0.17}$, suggesting that the estimators are able to recover the true mass on average, assuming our models are a good description of nature. The full distribution of f values is shown in Figure 7 with the median and 15.9 and 84.1 percentiles marked as solid and dashed lines. We will later use the results of these simulations to ensure that the scatter in the recovery of the true mass from a single sample is accurately accounted for in our uncertainties.

Figure 7.

Figure 7. Performance of the TME using 34 tracers. We ran 1000 Monte Carlo simulations (see text for details). Here we show a histogram of the ratio ${M}_{\mathrm{TME}}$/${M}_{\mathrm{true}}$ of the estimated mass to the true mass for the simulations, with the median value marked with a solid line and the 15.9 and 84.1 percentiles shown as dashed lines. These values are all given in the upper-right corner. We see that on average the estimator does indeed recover the true mass. We incorporate the scatter in the recovered values into the uncertainties for our final mass estimate.

Standard image High-resolution image

3.5. Tracer Mass Estimates

Now that we have estimated the anisotropy β, density power-law index γ, and potential power-law indices α, we can combine this information with the cluster position and velocity data, and finally estimate the mass of the MW within ${r}_{\max }$ = 21.1 kpc using Equation (1). From this mass, we can also estimate the circular velocity of the MW at ${r}_{\max }$ as

Equation (10)

We begin with the set of halos generated in Section 3.3 with circular velocities at the solar radius consistent with observations. We previously estimated α values for every halo. Additionally, we draw an anisotropy β at random from the posterior distribution of anisotropy values calculated in Section 3.1 for each halo. The value of γ is fixed from Section 3.2 and is the same for every halo. Finally, we draw positions and velocities from the distributions for each cluster and use Equation (1) to estimate the mass ${M}_{\mathrm{TME}}$ for each halo in the grid.

We know from the Monte Carlo simulations in Section 3.4 that there is some scatter in the performance of the estimator for a single sample of 34 clusters and that the true mass is related to the estimated mass via ${M}_{\mathrm{true}}$ = ${M}_{\mathrm{TME}}$/f. For each sample halo, we draw a value of f at random from the posterior distribution calculated in Section 3.4 and use that to infer the true mass from the TME mass. We adopt the median of the resulting distribution as the best mass estimate, and use the 15.9 and 84.1 percentiles to estimate uncertainties. Thus, we estimate the mass of the MW within 21.1 kpc to be

Equation (11)

and the circular velocity of the MW at this distance to be

Equation (12)

Note that the scatter in the performance of the estimator for a single sample and our uncertainties on both α, β, and the cluster properties have been naturally folded into our results using this method.

The TME only allows us to estimate the mass inside the radius of the outermost cluster in our sample. However, we can use our halo grid to predict what the virial mass13 of the MW might be, given the value of $M\left(\lt {r}_{\max }\right)$ we have estimated.

For each of the sample halos in the grid, we can estimate the mass ${M}_{\mathrm{grid}}\left(\lt {r}_{\max }\right)$ of the MW inside ${r}_{\max }$ for the halo model. We then define the probability density Phalo of the model value, given the distribution of $M\left(\lt {r}_{\max }\right)$ we estimate from the data. The variation of Phalo over the halo grid is shown in Figure 8. Halos on the left side of the swath in the diagram are most consistent with the observations. For each halo in the grid, we can also estimate the mass inside the virial radius ${M}_{\mathrm{grid},\mathrm{virial}}$. Now we ask: what are the virial masses of the halos for which ${M}_{\mathrm{grid}}\left(\lt {r}_{\max }\right)$ agrees with our measured value $M\left(\lt {r}_{\max }\right)$? That is, we look at the distribution of virial masses over the grid, weighted by the Phalo values that assess the consistency of that grid point with our measured value. From this, we can calculate the expected value and 15.9 and 84.1 percentiles of the virial mass implied by our tracer mass estimates, taking these probabilities into account. Thus, we find the virial mass of the MW to be

Equation (13)

Again, the scatter inherent in the estimator and the uncertainties on the α and β have been propagated into this result. Our estimate for ${M}_{\mathrm{virial}}$ has a larger fractional uncertainty than $M\left(\lt {r}_{\max }\right)$ owing to the uncertainty in the extrapolation of the dark halo mass profile to large radii.

Figure 8.

Figure 8. Probability Phalo of a set of NFW models defined by concentration c and virial radius rvir, based on the ${M}_{\mathrm{grid}}\left(\lt {r}_{\max }\right)$ predicted for the model and the actual $M\left(\lt {r}_{\max }\right)$ measured from the data. White regions are halos previously rejected as they have circular velocities inconsistent with observed values. The vertical lines mark, from left to right, halos with virial masses of (0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5) × 1012 ${M}_{\odot }$.

Standard image High-resolution image

Figure 9 summarizes our results. The left panel shows the distribution of mass estimates M (<21.1 kpc), the middle panel shows the distribution of ${v}_{\mathrm{circ}}$ (21.1 kpc) estimates, and the right panel shows the resulting distribution of predicted virial masses ${M}_{\mathrm{virial}}$.

Figure 9.

Figure 9. Summary of results using the Gaia GC sample for which ${r}_{\max }$ = 21.1 kpc. Left: distribution of tracer mass estimates $M\left(\lt {r}_{\max }\right)$. Middle: distribution of estimates of the circular velocity ${v}_{\mathrm{circ}}\left({r}_{\max }\right)$ made using the mass estimates in the left panel. Right: virial mass ${M}_{\mathrm{virial}}$ estimates inferred from the grid of halos sampled, weighted by the match to distribution of masses in the left panel. In each case, we adopt the median of the distribution and the 15.9 and 84.1 percentiles as the estimate and its uncertainties; these are given in the top-right corner of each panel.

Standard image High-resolution image

3.6. Expanded GC Sample

Our analysis so far has been limited to 21.1 kpc, and only probes the very inner parts of the halo. Furthermore, the uncertainties on our virial mass are large due to extrapolation of our results out to significantly larger radii. More clusters further out in the halo would thus be extremely beneficial. Increased sample size would also help to decrease uncertainties in all estimates, assuming our model assumptions are reasonable, or highlight problems if they are incorrect.

As described, in Section 2.2, we augment our sample with 12 extra clusters with HST PMs from Sohn et al. (2018). We repeat the analysis as laid out above. The density slope γ remains the same, but all other parts of the analysis necessarily change with the new sample. We find an anisotropy,

Equation (14)

The Monte Carlo simulations are run with the same α and γ as before, but with the revised β and with an increased radial range and give $f={1.00}_{-0.13}^{+0.14}$. For the new ${r}_{\max }$ = 39.5 kpc, we then estimate

Equation (15)

which corresponds to a circular velocity at ${r}_{\max }$ of

Equation (16)

We again extrapolate out to estimate the virial mass and find

Equation (17)

We summarize our results in Table 2. Sample A is the set of Gaia measurements and Sample B is the set of combined Gaia and HST measurements. Samples C and D are discussed in Section 4.

Table 2.  Summary of Mass Results for the Halo GC Subsamples

Sample N ${r}_{\min }$ ${r}_{\max }$ $M\left(\lt {r}_{\max }\right)$ ${v}_{\mathrm{circ}}\left({r}_{\max }\right)$ ${M}_{\mathrm{virial}}$
    (kpc) (kpc) (1012 ${M}_{\odot }$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) (1012 ${M}_{\odot }$)
Aa 34 2.0 21.1 ${0.21}_{-0.03}^{+0.04}$ ${206}_{-16}^{+19}$ ${1.28}_{-0.48}^{+0.97}$
Bb 46 2.0 39.5 ${0.42}_{-0.06}^{+0.07}$ ${214}_{-15}^{+17}$ ${1.54}_{-0.44}^{+0.75}$
Cc 26 2.0 21.1 ${0.21}_{-0.03}^{+0.04}$ ${209}_{-17}^{+21}$ ${1.34}_{-0.50}^{+1.02}$
Dd 41 4.1 39.5 ${0.36}_{-0.06}^{+0.07}$ ${199}_{-16}^{+19}$ ${1.22}_{-0.36}^{+0.63}$

Notes.

aGaia GCs. bGaia and HST GCs. cGaia GCs with possible merger group removed. dGaia and HST GCs with high ${v}_{\tan }$ clusters removed.

Download table as:  ASCIITypeset image

4. Discussion

4.1. Anisotropy

Using only the Gaia halo GCs we estimated the anisotropy $\beta ={0.46}_{-0.19}^{+0.15}$ over $2.0\leqslant r\leqslant 21.1\,\mathrm{kpc}$, and using the expanded sample we found $\beta ={0.52}_{-0.14}^{+0.11}$ over $2.0\leqslant r\leqslant 39.5\,\mathrm{kpc}$. Both values indicate that the halo over the range of the GC sample is radially anisotropic. Sohn et al. (2018) reported $\beta ={0.609}_{-0.229}^{+0.130}$ over the range 10.6 < r < 39.5 kpc. All of these values are consistent within their uncertainties, but the trends suggest that the halo becomes more radially anisotropic in its outer regions, in line with predictions from cosmological simulations (e.g., Diemand et al. 2007).

There are a number of other estimates for the anisotropy of the halo over radial ranges that overlap that of our sample. Our radial β is inconsistent with the estimates for individual halo star samples of Sirko et al. (2004) and Cunningham et al. (2016), both of which favor an isotropic or even tangentially anisotropic halo (although the uncertainties in the latter were large enough that that results cannot be called truly discrepant), but in good agreement with the radial estimates from Bond et al. (2010) and Deason et al. (2012).

4.2. Masses

As we discussed in Section 1, MW mass estimates can vary markedly based on the types of data used, the techniques used, and the assumptions that go into the mass estimate, with estimates for the total mass of the MW varying between ∼0.5 and 3 × 1012 ${M}_{\odot }$. Furthermore, only the timing argument and abundance-matching studies actually estimate the total mass of the MW; most other estimates can only measure the mass within the extent of the data set being used, as we have done here, so most estimates are given at different radii, which makes them hard to compare.

Indeed, our mass estimates at 21.1 kpc and 39.5 kpc cannot easily be compared directly. We do note that the latter mass is about twice as large for about twice the size of the enclosed radius, consistent with naive expectations for a nearly isothermal sphere. Our values are instead most usefully compared to other estimates near these radii.

Kafle et al. (2012) estimated M (<25 kpc) ≈ 0. 21 × 1012 ${M}_{\odot }$, unfortunately without uncertainties, by analyzing blue horizontal branch stars in the halo, and Küpper et al. (2015) estimated M (<19 kpc) = (0. 21 ± 0.04) × 1012 ${M}_{\odot }$ by analyzing the orbit of the Palomar 5 tidal stream, both in good agreement with our estimate of $M(\lt 21.1\,\mathrm{kpc})={0.21}_{-0.03}^{+0.04}\times {10}^{12}\,{M}_{\odot }$. Since completing this work, there have been further mass estimates from Gaia measurements: Posti & Helmi (2019) applied a Bayesian estimator to the Gaia sample as we used here to estimate $M(\lt 20\,\mathrm{kpc})={0.191}_{-0.015}^{+0.017}\times {10}^{12}$ ${M}_{\odot }$; and Eadie & Jurić (2018) applied an alternative Bayesian estimator to the expanded sample from Vasiliev (2019) to estimate $M(\lt 25\,\mathrm{kpc})\,={0.26}_{-0.02}^{+0.03}\times {10}^{12}$ ${M}_{\odot }$, also in good agreement with our results.

Sohn et al. (2018) estimated $M(\lt 39.5\,\mathrm{kpc})={0.60}_{-0.11}^{+0.17}\,\times {10}^{12}$ ${M}_{\odot }$ using the same a set of HST clusters we used to augment our sample, slightly higher than, but still consistent with, our estimate here of $M(\lt 39.5\,\mathrm{kpc})={0.42}_{-0.06}^{+0.07}\,\times {10}^{12}\,{M}_{\odot }$.

There are a number of previous estimates of the mass within 50 kpc including ${0.54}_{-0.36}^{+0.02}\times {10}^{12}$ ${M}_{\odot }$ (Wilkinson & Evans 1999), ${0.55}_{-0.02}^{+0.00}\times {10}^{12}$ ${M}_{\odot }$ (Sakamoto et al. 2003), and 0.42 ± 0.04 × 1012 ${M}_{\odot }$ (Deason et al. 2012). Remembering that these estimates were made ∼10 kpc further out than our sample, the first two estimates are in extremely good agreement with our estimate given expectations for a nearly isothermal sphere, and the last is slightly lower, but still consistent with our estimate.

More recently, Vasiliev (2019) estimated $M(\lt 50\,\mathrm{kpc})\,={0.6}_{-0.09}^{+0.14}\times {10}^{12}$ ${M}_{\odot }$ and Eadie & Jurić (2018) estimated $M(\lt 50\,\mathrm{kpc})={0.37}_{-0.03}^{+0.04}\times {10}^{12}$ ${M}_{\odot }$, both using an expanded sample of PMs measured with Gaia from the former work. Our estimate is consistent with the Vasiliev result and marginally inconsistent with the Eadie & Jurić result; our estimate further in at ∼21 kpc was consistent with the Eadie & Jurić estimate at 25 kpc, suggesting that their method predicts a steeper mass density slope than favored by our models.

We also note that Gibbons et al. (2014) estimated a mass (0.41 ± 0.04) × 1012 ${M}_{\odot }$ inside 100 kpc by fitting to the Sagittarius tidal stream, clearly at odds with our results. These authors only fitted the locations of the apocenters and pericenters of the trailing and leading arms of Sagittarius, and not the intervening locations. The data set is now much richer, with Gaia DR2 likely to add to our knowledge of the PMs along the stream. It would be interesting to re-visit the work of Gibbons et al. in light of this.

Now let us consider our virial mass estimates. Using only the Gaia GCs, we find ${1.28}_{-0.48}^{+0.97}$ × 1012 M and, using the expanded sample, we find ${1.54}_{-0.44}^{+0.75}\times {10}^{12}\,{M}_{\odot }$. Both of these values are in good agreement, but note that the error bars on the latter are smaller than those on the former. The increased sample size may have a small effect here, but the main contribution to this decrease comes from the fact that there is less extrapolation involved in estimating a virial mass from data at ∼40 kpc than there is from data at ∼20 kpc. This means that there is less variation in the allowed halos, which in turn allows us to constrain the virial mass much better.

Sohn et al. (2018) estimated a virial mass of ${1.87}_{-0.47}^{+0.67}\,\times {10}^{12}$ ${M}_{\odot }$, slightly higher than but in good agreement with both of our values.

Comparison with other studies is again tricky as the definition of "virial" can change from study to study. The recent review by Bland-Hawthorn & Gerhard (2016) put a number of different estimates at a different radii onto the same scale for comparison. Our estimates are best compared against their M100 values, that is the mass within the radius for which the mean overdensity is 100, which is very close to the mean overdensity value we use in this work. Our estimates are inconsistent with the most massive and least massive of these, and agree best with the intermediate values. (See also Figure 1 of Wang et al. (2015) and Figure 1 of Eadie et al. (2017) for further comparisons.)

More recently, Vasiliev (2019) estimated a virial mass ${0.8}_{-0.2}^{+0.5}\times {10}^{12}$ ${M}_{\odot }$ using a large sample of GCs with Gaia PMs. This is considerably smaller than our value, and the virial radius they estimate (∼160 ± 20 kpc) is much smaller than those our method favors, implying that the outer halo density is steeper than the NFW value of −3. Many of the additional clusters in the Vasiliev data set are further out in the halo than the sample we used, which in principle should give stronger constraints on the mass at larger distances. However, these outer clusters are more likely to have been recently accreted and may not be fully phase-mixed, which could invalidate any mass estimates. To fully explore these differences and their implications is beyond the scope of this paper.

4.3. Halo Density Profile

We can also consider the implications of our results for the density profile of the halo. Figure 8 shows the agreement between the mass enclosed at ${r}_{\max }$ for the model and from the data for the set of allowed halos with consistent circular velocities at the solar radius. On this, we mark the virial radii for which the halo virial masses are, from left to right, (0.5–3.5) × 1012 ${M}_{\odot }$ at 0.5 × 1012 ${M}_{\odot }$ intervals for reference. It is clear that the data can be explained by either low concentration, high virial radius (and thus high virial mass) halos—similar to those favored by cosmological simulations (Klypin et al. 2011)14 —or by high concentration, low virial radius/mass halos. The latter seem to be favored over the former, as there are more such halos that can adequately fit the data, but we cannot rule out either.

To derive our results, we have made a number of assumptions, one of which is that the halo and the GC distributions are spherical, neither of which is necessarily true. Indeed, there are some hints of non-sphericity in that σθ < σϕ (see Table 1). To lowest order, the predictions of a spherical model, when applied to non-spherical distributions, can be interpreted as an estimate of the spherically averaged quantities of the actual distribution. In reality, some bias may be introduced. The exact size of this needs to be quantified through analysis of either, e.g., triaxial equilibrium models, or cosmological simulations with non-spherical distributions, both of which are outside the scope of the present paper.

4.4. Substructure

One of the assumptions in our work is that the halo GCs comprise a statistically independent, well-mixed tracer population in the MW's gravitational potential. It is for this reason that we excluded all but one of the GCs associated with the Sagittarius dSph. However, there have been recent claims that up to two thirds of the local stellar halo may have been deposited by the single encounter of a massive dwarf galaxy (e.g., Deason et al. 2013; Belokurov et al. 2018; Myeong et al. 2018a). The evidence for this rests mainly on the highly radially anisotropic, non-Gaussian distribution of velocities of metal-rich halo stars in both the SDSS–Gaia and Gaia DR2 catalogs. If true, then such a massive satellite will also have been accompanied by its own retinue of GCs. Very recently, Myeong et al. (2018b) have tentatively used Gaia DR2 to identify eight halo GCs associated with this merger event (NGCs 1851, 1904, 2298, 2808, 5286, 6779, 6864, and 7089).

It is prudent to check the robustness of our results to these claims. Re-running the calculations for our DR2-only sample with these eight GCs removed, we obtain

Equation (18)

very close to our earlier result in Equation (11). We note that while the mass enclosed hardly changes, the inferred anisotropy does:

Equation (19)

This change is understandable, as the most eccentric GCs have been excised from the sample. Indeed, Myeong et al. (2018b) identified their eight GCs from clustering in radial action, and so the removed GCs do make a prominent contribution to the anisotropy parameter. These results are added to Tables 1 and 2 as Sample C.

Alternatively, recently accreted GCs or young halo GCs that have not yet phase-mixed with the rest of the GC population could be on highly tangential orbits. As we noted in Section 3.1, we see that the Gaia (and HST) GCs show net tangential motion; this is also apparent in Figure 1 where the high ${v}_{\tan }$ GCs have been colored in orange and cyan. To assess the dependence of our results on these outliers, we re-ran the analysis of the combined Gaia+HST sample with these clusters removed. We now obtain an anisotropy estimate of

Equation (20)

The value of β is now obviously more radially anisotropic than before as we have removed high ${v}_{\tan }$ clusters. The enclosed mass within ${r}_{\max }$ becomes

Equation (21)

and the inferred virial mass is

Equation (22)

These results are added to Tables 1 and 2 as Sample D. The inferred masses are somewhat decreased, but still within the uncertainty ranges of the previously quoted values. These tests with modified samples add confidence that our inferred MW mass estimates are pleasingly robust to potential substructure in the halo GC distribution.

5. Conclusions

We have used Galactocentric motions for a set of 34 halo GCs estimated using PM data from the second Gaia data release to estimate the anisotropy of the halo GC population and then to estimate the mass of the MW inside 21.1 kpc, the position of the most distant cluster in our sample. Combined with a catalog of clusters from a recent HST study, we were able to estimate the anisotropy and mass of the MW inside 39.5 kpc.

Using only the Gaia sample, we find an anisotropy $\beta ={0.46}_{-0.19}^{+0.15}$ in the range of the clusters $2.0\leqslant r\leqslant 21.1\,\mathrm{kpc}$. With the addition of the HST clusters, we find $\beta ={0.52}_{-0.14}^{+0.11}$ over $2.0\leqslant r\leqslant 39.5\,\mathrm{kpc}$. This suggests that the halo is radially anisotropic, consistent with a number of previous estimates and predictions from cosmological simulations and in good agreement with a number of similar studies, including the study of 16 halo GCs with HST by Sohn et al. (2018).

We estimate the masses $M(\lt 21.1\,\mathrm{kpc})={0.21}_{-0.03}^{+0.04}\,\times {10}^{12}\,{M}_{\odot }$ and $M(\lt 39.5\,\mathrm{kpc})={0.42}_{-0.06}^{+0.07}\times {10}^{12}\,{M}_{\odot }$. These masses correspond to circular velocities of ${v}_{\mathrm{circ}}(21.1\,\mathrm{kpc})\,={206}_{-16}^{+19}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${v}_{\mathrm{circ}}(39.5\,\mathrm{kpc})={214}_{-15}^{+17}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ respectively which, compared with estimates for the circular velocity at the solar radius and to each other, is consistent with a rotation curve that is not falling rapidly over the radial range of our data.

From these, we are also able to place constraints on the virial mass of the MW. We favor the results for the combined sample here, as it is more reasonable to perform this extrapolation for a larger sample size with a broader radial range and, more importantly, greater reach. We find ${M}_{\mathrm{virial}}={1.54}_{-0.44}^{+0.75}\,\times {10}^{12}\,{M}_{\odot }$, again of intermediate size compared with previous estimates. All of our mass estimates are intermediate in value when compared to the range of values found in the literature, with both low-mass (<1012 ${M}_{\odot }$) and very high-mass (≳2.5 × 1012 ${M}_{\odot }$) MWs generally disfavored.

Previous mass estimates have often been limited by either the mass–anisotropy degeneracy for LOS velocity samples, or the small sample sizes for distant objects with 3D motions. Given the new results from Gaia and HST, these are now both resolved. Various kinds of systematics may now become the dominant source of uncertainty. Nevertheless, further progress will come from having yet larger sample sizes. Gaia Collaboration et al. (2018b) measured PMs for only 75 GCs out of the 157 known in the MW (Harris 1996, 2010 edition) by making extremely conservative cuts on the number of member stars identified. It is likely that PMs and, hence, Galactocentric motions can be measured for many more Galactic GCs, both using DR2 and future data releases. Such measurements will further refine our understanding of the MW mass.

L.L.W. wishes to thank Alis Deason, Mark Fardal, Elena Pancino, and Mark Gieles for very useful conversations related to this work. We also thank Amina Helmi for making supporting data for Gaia Collaboration et al. (2018b) available electronically.15 We thank the referee for the suggestions that improved the presentation of our results. Support for this work was provided by grants for HST program GO-14235 provided by the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555.

This work has made use of data from the European Space Agency (ESA) mission Gaia,16 processed by the Gaia Data Processing and Analysis Consortium (DPAC).17 Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This research made use of Astropy,18 a community-developed core Python package for Astronomy. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

This project is part of the HSTPROMO (High-resolution Space Telescope PROper MOtion) Collaboration,19 a set of projects aimed at improving our dynamical understanding of stars, clusters and galaxies in the nearby universe through measurement and interpretation of PMs from HST, Gaia, and other space observatories. We thank the collaboration members for the sharing of their ideas and software.

Facilities: Gaia - , HST. -

Software: astropy (Astropy Collaboration et al. 2013, 2018), emcee (Foreman-Mackey et al. 2013).

Appendix: Galactocentric Positions and Motions

In Section 2, we calculated Galactocentric positions and velocities for all 75 GCs in Gaia Collaboration et al. (2018b). We provide these positions and motions in both spherical and Cartesian coordinates in Table 3 along with their uncertainties and the correlations between velocity components. Note that the astrometric measurements and heliocentric Cartesian positions and velocities are provided in Gaia Collaboration et al. (2018b).20

Table 3.  Galactocentric Positions and Velocities in Spherical and Cartesian Coordinates for the Full Gaia GC Sample

Name r θ ϕ vr vθ vϕ Crθ Crϕ Cθϕ
  (kpc) (deg) (deg) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$)      
NGC 0104 7.6 ± 0.1 −24.7 ± 0.7 202.0 ± 0.8 −12.0 ± 2.8 44.8 ± 1.5 −187.9 ± 5.2 −0.006 0.004 −0.001
NGC 0288 12.2 ± 0.2 −46.8 ± 0.8 179.7 ± 0.0 −31.9 ± 1.1 41.0 ± 1.0 46.0 ± 8.6 −0.001 0.001 −0.001
NGC 0362 9.5 ± 0.1 −40.6 ± 0.7 224.5 ± 1.4 143.9 ± 4.4 30.2 ± 7.3 4.3 ± 4.8 0.007 0.008 −0.044
NGC 1851 16.9 ± 0.3 −24.3 ± 0.3 215.5 ± 0.5 131.5 ± 3.1 −30.3 ± 3.3 7.0 ± 4.7 −0.001 0.013 −0.039
NGC 1904 19.0 ± 0.3 −19.4 ± 0.2 207.4 ± 0.4 43.9 ± 2.7 20.4 ± 2.6 −8.3 ± 5.9 0.000 −0.023 −0.068
x y z vx vy vz Cxy Cxz Cyz v
(kpc) (kpc) (kpc) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$) ($\mathrm{km}\,{{\rm{s}}}^{-1}$)       ($\mathrm{km}\,{{\rm{s}}}^{-1}$)
−6.4 ± 0.2 −2.6 ± 0.1 −3.2 ± 0.1 −77.6 ± 2.3 171.3 ± 5.9 45.7 ± 0.9 −0.000 0.000 −0.000 193.6 ± 5.0
−8.4 ± 0.2 0.0 ± 0.0 −8.9 ± 0.2 −8.3 ± 1.3 −46.0 ± 8.6 51.3 ± 0.7 0.007 −0.000 −0.000 69.7 ± 5.7
−5.2 ± 0.2 −5.1 ± 0.1 −6.2 ± 0.1 −89.0 ± 4.2 −93.5 ± 7.1 −70.8 ± 2.0 0.002 −0.001 −0.001 147.3 ± 5.7
−12.5 ± 0.2 −8.9 ± 0.2 −6.9 ± 0.2 −83.4 ± 1.3 −68.1 ± 5.8 −81.6 ± 2.3 −0.000 0.000 −0.001 135.2 ± 3.0
−15.9 ± 0.2 −8.3 ± 0.2 −6.3 ± 0.1 −46.6 ± 2.0 −14.8 ± 6.3 4.6 ± 2.4 −0.008 −0.016 0.107 49.6 ± 2.1

Note. This table is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Footnotes

  • The degree to which the LOS and Galactocentric radial velocities are similar depends on the geometry of the system, specifically the position of the object relative to both the Sun and the Galactic center, and there are a few objects with more favorable geometry to give some insight on the tangential motions. However, anisotropy measurements rely on averaging over many objects, so favorable geometry for a few objects provides limited benefit overall.

  • We have not included the Gaia systematic errors of ∼0.035 mas yr−1 as, at the distances we are considering, these errors translate to a few km s−1 and will be very much smaller than the velocity dispersion of the halo.

  • For Gaussian distributions, the 15.9 and 84.1 percentiles enclose the 1σ confidence interval. The posterior distributions we derive thoroughout are generally not Gaussian, so that these uncertainties should be interpreted as actual percentiles ranges, with any analogy to Gaussian errors only being approximate at best.

  • Orbits were calculated in three different potentials; we insisted that at least two of the three apocenter estimates must pass the cut.

  • As we will see later, this value is equivalent to two disk scale lengths for our adopted disk model. We experimented with different halo samples, and found that our results were not sensitive to the particular choice of apocenter cut that we used.

  • 10 

    As mass estimates depend strongly on r, a single cluster far out in the halo can unduly bias any mass estimates.

  • 11 

    The fitted double power law provides a better fit to the cumulative number profile inside ∼20 kpc than outside, as there are far fewer clusters at large radii to constrain the outer fits. Mass estimates correlate mildly with γ such that a change in γ of 0.2 will change the mass estimate inside 21.1 kpc by ∼5%, but we do not include this in our calculations.

  • 12 

    Our goal here is to approximate the slope of the potential over the region of interest. The precise functional form of the potential is unknown. We are assuming here that picking a single functional form and varying its parameters will give a similar distribution of α values as sampling a variety of functional forms, especially as we insist that that circular velocity at the solar radius must be consistent with observations.

  • 13 

    We calculate the virial mass at each grid point directly given the NFW parameters of the point. The virial mass is defined as the mass inside the virial radius, and the virial radius is defined as the radius at which the mean overdensity of the halo relative to the critical density is Δvir. We adopt the prescription for Δvir from Bryan & Norman (1998). Note that we do not estimate the virial mass from the power-law fits, as the power laws are assumed to hold over the range of the GC data, not out into the outskirts of the halo.

  • 14 

    Note though that these simulations do not fully consider the effect of baryons in shaping the halo.

  • 15 
  • 16 
  • 17 
  • 18 
  • 19 
  • 20 
Please wait… references are loading.
10.3847/1538-4357/ab089f