Issue 
A&A
Volume 621, January 2019



Article Number  A56  
Number of page(s)  10  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201833355  
Published online  08 January 2019 
Mass and shape of the Milky Way’s dark matter halo with globular clusters from Gaia and Hubble
Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
email: posti@astro.rug.nl
Received:
3
May
2018
Accepted:
5
November
2018
Aims. We estimate the mass of the inner (< 20 kpc) Milky Way and the axis ratio of its inner dark matter halo using globular clusters as tracers. At the same time, we constrain the distribution in phasespace of the globular cluster system around the Galaxy.
Methods. We use the Gaia Data Release 2 catalogue of 75 globular clusters’ proper motions and recent measurements of the proper motions of another 20 distant clusters obtained with the Hubble Space Telescope. We describe the globular cluster system with a distribution function (DF) with two components: a flat, rotating disclike one and a rounder, more extended halolike one. While fixing the Milky Way’s disc and bulge, we let the mass and shape of the dark matter halo and we fit these two parameters, together with six others describing the DF, with a Bayesian method.
Results. We find the mass of the Galaxy within 20 kpc to be M(<20 kpc) = 1.91_{−0.17}^{+0.18}×10^{11} M_{⊙}, of which M_{DM}(<20 kpc) = 1.37_{−0.17}^{+0.18}×10^{11} M_{⊙} is in dark matter, and the density axis ratio of the dark matter halo to be q = 1.30 ± 0.25. Assuming a concentrationmass relation, this implies a virial mass M_{vir} = 1.3±0.3×10^{12} M_{⊙}. Our analysis rules out oblate (q < 0.8) and strongly prolate halos (q > 1.9) with 99% probability. Our preferred model reproduces well the observed phasespace distribution of globular clusters and has a disc component that closely resembles that of the Galactic thick disc. The halo component follows a powerlaw density profile ρ ∝ r^{−3.3}, has a mean rotational velocity of V_{rot} ≃ −14km s^{−1} at 20 kpc, and has a mildly radially biased velocity distribution (β ≃ 0.2 ± 0.07, which varies significantly with radius only within the inner 15 kpc). We also find that our distinction between disc and halo clusters resembles, although not fully, the observed distinction in metalrich ([Fe/H] > −0.8) and metalpoor ([Fe/H] ≤ −0.8) cluster populations.
Key words: Galaxy: kinematics and dynamics / Galaxy: structure / Galaxy: halo / globular clusters: general
© ESO 2019
1. Introduction
Giant leaps in the physical understanding of our Universe are often made when new superb datasets that peer into previously uncharted territory become available. The second data release of Gaia (Gaia Collaboration 2018a) has just arrived and is by all means a perfect instance of such a leap. The extent, accuracy, and quality of the data provided is unprecedented, and it is leading the field of Galactic Astronomy to a completely new era.
While most of the new discoveries and exciting results are probably still hidden in the vastness of the dataset, the increase in the number of objects observed and the greater accuracy, for example of the measured tangential motions on the sky, has the potential to immediately lead to groundbreaking results (Gaia Collaboration 2018b,c; Antoja et al. 2018). In this study we exploit this novel dataset by employing a wellestablished, but very sophisticated method to infer new tight constraints on fundamental parameters, such as the total mass and the shape of the gravitational potential of the inner Milky Way.
We do this by using the catalogue of absolute proper motions that Gaia Collaboration (2018b) estimated for 75 globular clusters (GCs) around the Milky Way. Thanks to the unprecedented performances of the Gaia astrometry, proper motions for these satellites could be measured with a typical accuracy of a few tens of μ as yr^{−1}, which roughly corresponds to an accuracy of 0.5−2kms^{−1} in tangential velocity for a cluster located at 10 kpc from us. Remarkable measurements of similar accuracy were also recently released by Sohn et al. (2018) using the Hubble Space Telescope (HST). These authors estimated proper motions of 20 clusters at large Galactocentric distances (> 10 kpc), using two HST epochs at least ∼6 yr apart.
With the dataset provided by these two new catalogues, we can use the GCs as tracers of the Galactic potential and hope to infer its total mass (e.g. Bahcall & Tremaine 1981; Watkins et al. 2010). Moreover, given the size and the unprecedented accuracy of the new catalogue, there may well be enough signal in the data to constrain simultaneously the mass and the axis ratio of the isodensity surfaces of the total mass distribution of the Galaxy, similarly to previous work using the kinematics of individual stars in the Milky Way halo (e.g. Olling & Merrifield 2000; Smith et al. 2009; Bowden et al. 2016; Posti et al. 2018). In order to have enough inference from the data, we can model the Milky Way GC system with equilibrium models in an arbitrary axisymmetric Galactic potential. One possibility is then to use distribution functions (DFs) to represent the phasespace distribution of the GC system and to measure the characteristic parameters of such DFs with the data from Gaia and HST. This would not only allow the determination of the fundamental parameters of the Galactic potential, but would also simultaneously constrain the full phasespace distribution of the GCs around the Milky Way.
In this paper we use this approach and follow closely Binney & Wong (2017, hereafter BW17), who described the GC system with two distinct DFs, one with halolike dynamics (typically representing the metalpoor clusters), and one with disclike dynamics (following the more metalrich clusters, see Zinn 1985). These DFs are constructed in the space of the action integrals of motion. We build upon the work of BW17 and improve on it in several ways: i) we use a much larger and more accurate dataset provided by very recent measurements presented above, ii) we allow for the dark matter halo mass and shape to vary in our fit, and iii) we simplify the description of the DFs by fixing some characteristic parameters while still reproducing remarkably well the observed distribution of GCs.
The paper is organised as follows: we introduce the data used in Sect. 2 and our modelling technique in Sect. 3; we describe our Bayesian approach to determine the distribution of the model parameters in Sect. 4 and we discuss our results in Sect. 5; finally, we summarise and conclude in Sect. 6. Throughout the paper we use a distance to the Galactic centre of R_{⊙} = 8.3 kpc and a peculiar motion of the Sun of (U_{⊙}, V_{⊙}, W_{⊙}) = (11.1, 12.24, 7.25) km s^{−1} (Schönrich et al. 2010). We note, however, that reasonable differences in R_{⊙} and in the solar peculiar motion (see Sect. 3.2 and Sect. 5.3.3 in BlandHawthorn & Gerhard 2016, respectively) are too small to significantly affect our conclusions.
2. Data
2.1. Globular cluster data from Gaia DR2 and the Hubble Space Telescope
Our dataset is primarily composed of the recently estimated proper motions and radial velocities for 75 GCs from the Gaia Collaboration (2018b). This catalogue provides data of the best quality which is ideal to study the motion in phasespace of GCs with unprecedented accuracy. The following observables are provided: (α, δ) are right ascension and declination, and (μ_{α*}, μ_{δ}) are the respective proper motions on the sky. Alongside these measurements, we also have standard uncertainties and the correlation between the two proper motions C_{μα*, μδ}, which we use for generating samples from the error distribution of each cluster. We supplement the Gaia DR2 catalogue with the recent measurements of 16 other GCs from Sohn et al. (2018). Hence, we work with an unprecedented total of 91 GCs with exquisite proper motion data quality.
For both sets of clusters we compute heliocentric distances s from the extinctioncorrected distance moduli in the Vband, taken from the latest version of the catalogue compiled by Harris (1996), and we assume the uncertainty to be 0.05 mag for each object. The heliocentric radial velocities v_{rad} and their uncertainties are also taken from this catalogue. Therefore, the set of observables that we use is
where the error distribution of each cluster is assumed to be a multivariate normal distribution with nonnull correlation C_{μα*, μδ} between the proper motions and negligible variances in the sky positions (α, δ) for the Gaia DR2 clusters. We also add an additional systematic 35 μ as yr^{−1} to the uncertainty budget of the cluster proper motion measured by Gaia, as advised in Gaia Collaboration (2018b).
Finally, to make our tracer sample as complete as possible, we also consider 52 GCs for which no proper motion data is available. For this sample, we again obtain distance moduli and radial velocities, together with their uncertainties from Harris (1996). This leaves us with a total sample of 91 + 52 = 143 GCs.
2.2. Sagittarius clusters
Several of the known GCs have been tentatively associated with the Sagittarius dwarf galaxy/stream by different authors in the past (see e.g. Law & Majewski 2010). The Sagittarius dwarf galaxy is currently merging with our Galaxy (Ibata et al. 1994) and it is possibly bringing in a few GCs. If a significant fraction of the GCs that we use as tracers are actually associated with the dwarf galaxy, this would imply a nonrandom sampling of Galactic phasespace which could lead to biased estimates of the mass and shape of the Galactic potential using our DFbased method.
For this reason we run our algorithm with the full sample of GCs, but also excluding the four clusters Arp 2, Palomar 12, Terzan 7, and Terzan 8 that have been associated with the Sagittarius dwarf galaxy according to the dynamical analysis of Law & Majewski (2010) and have been recently confirmed by Sohn et al. (2018). We retain M 54, the nuclear cluster in Sagittarius, as this now traces the orbit of the dwarf singly. In addition, we also remove from our analysis NGC 5053, which was recently found to be one of a pair of clusters (with NGC 5024, Gaia Collaboration 2018b).
3. Dynamical model
Here we introduce the technique that we use to model the distribution in phase space of the GC population given a Galactic gravitational potential. Our method borrows heavily from BW17, to which we refer the reader to for a more exhaustive description. The Galactic potentials, the mapping from positionvelocities (x, v) to actionangles (θ, J) and the DF models are constructed with the AGAMA code (Vasiliev 2019).
3.1. Galactic potential
We model the mass distribution of the Milky Way with five axisymmetric and analytic components: a gaseous disc, stellar thin and thick discs, an oblate bulge, and a dark halo. The discs are described by doubleexponential density distributions
where R_{d} and h_{d} are scalelength and height, Σ_{0} is a normalisation constant, and R_{hole} is the size of the central cavity (nonzero only for the gaseous disc), while the bulge and halo are described by
where , r_{0} is a scale radius, q is the axis ratio, and ρ_{0} is a normalisation constant. The parameters of the three discs and the bulge were fitted to a compilation of observational constraints by Piffl et al. (2014b), most notably to the kinematics of ∼200 000 stars in the solar neighbourhood as observed by the RAdial Velocity Experiment (RAVE, Steinmetz et al. 2006). These are kept fixed in our analysis and are summarised in Table 1.
Fixed parameters of the model (see also Piffl et al. 2014b).
For the dark matter halo instead, we make a cosmologically motivated ansatz fixing the slopes in Eq. (3) to those of a Navarro et al. (1996, hereafter NFW, a = 1 and b = 3) model, truncating the halo with an exponential cutoff at the virial radius, r_{cut} = r_{vir}, and imposing the concentrationmass relation found in ΛCDM Nbody simulations: log_{10} c = 1.025 − 0.097log_{10}(M/10^{12} h^{−1} M_{⊙}), where c = r_{vir}/r_{0} is the concentration and h = 0.68 is the Hubble parameter (see Dutton & Macciò 2014). Thus, our halo model is described by two free parameters: a mass and the axis ratio q.
3.2. Distribution function of the GC population
Assuming the GCs to be in dynamical equilibrium within the potential of the Galaxy, we can describe their distribution in phasespace with a function f(J) of the three action integrals J = (J_{R}, J_{ϕ}, J_{z}). This function, the DF, is normalised such that (2 π)^{3} f(J) dJ is the probability that a GC moves on the orbit specified by the action triplet J. In a general axisymmetric potential, the action integrals can be computed from positions and velocities, J = J(x, v), with the “Staeckel Fudge” algorithm (Binney 2012).
Following Zinn (1985), who suggested that the GC system is composed of two populations distinct in their metallicity and their phasespace distribution, we allow the DF to have two dynamically distinct components, one disclike and one halolike: f(J) = f_{disc} (J)+f_{halo}(J). For the disc component, we use the “quasiisothermal” DF, introduced by Binney (2010), as implemented in Vasiliev (2019)
where Ω, κ, and ν are the circular, radial, and vertical epicycle frequencies evaluated at the radius of the circular orbit with angular momentum ; the factor controlling the disc surface density is Σ = Σ_{0} exp[−R_{c}(J_{ϕ)}/R_{d}] and
makes the DF odd in the angular momentum J_{ϕ} and controls the net rotation of the disc component. The factors σ_{R} and σ_{z}, which determine the radial and vertical velocity dispersions, are chosen to mimic the Galactic thick disc model as in Piffl et al. (2014b), with an exponential vertical velocity dispersion with constant scaleheight, , and a radial velocity dispersion σ_{R} = σ_{R, 0}exp(−R_{c}/R_{σ}). In analogy with the thick disc, we fix R_{σ} = 13 kpc (as measured by Piffl et al. 2014b) and h_{d} = 0.2R_{d}. We are thus left with two free parameters to be fitted, the disc scalelength R_{d} and the central radial dispersion σ_{R, 0}, that completely characterise f_{disc}(J).
For the halo component, we use the double powerlaw DF, introduced by Posti et al. (2015, see also Williams & Evans 2015), as implemented in the AGAMA code. Here we fix the halo DF to have a constant density core in phasespace, thus making the 3D density distribution of the model a single power law with a central core; the halo DF then is
where J_{0} is a scale action defining the size of the constant density core, J_{cut} is a (large) cutoff action, g(J) = k_{R}J_{R} + (3−k_{R})(J_{ϕ} + J_{z})/2 is a homogeneous function of degree one and
is the factor controlling the net rotation of the system, with χ = 0 being the nonrotating case and χ = ±1 being the maximally rotating/counterrotating case. Here we fix two parameters: i) the cutoff action to J_{cut} = 10^{5} kpc km s^{−1}, so that the halo density distribution is cut off at distances much larger than that of the farthest GC (≳300 kpc) ensuring a finite halo mass, and ii) J_{ϕ, 0} = 100 kpc km s^{−1}, so that the possible rotation of the halo smoothly goes to zero to the centre (needed for continuity of the DF, see BW17). These two are the only nuisance parameters of the halo DF, and we are left with four free parameters to be fitted: the slope B, the extent of the constant density core J_{0}, the rotation parameter χ, and k_{R} which controls the system’s velocity anisotropy.
Our choice for the GC halo DF is relatively simple, describing the density profile of the halo GCs with a single power law. This is convenient because the small number of GCs at distances > 20 kpc limits our inference on radial trends. Therefore, it does not have the same level of complexity as seen in the stellar halo of the Galaxy (e.g. Carollo et al. 2007; Deason et al. 2011; Sesar et al. 2013; Xue et al. 2015; Das & Binney 2016; Iorio et al. 2018).
Finally, we note that the normalisation parameters Σ_{0} for f_{disc} and M_{h} for f_{halo} are unimportant in our case, since the GC system is treated as massless and simply traces the gravitational potential.
4. Parameter estimation
We estimate the posterior distribution of the model parameters Ξ that best represent the data w using standard Bayesian inference
where P(Ξ) is the prior, P(wΞ) is the likelihood, and P(w) is the evidence, which is unimportant for our purposes and hence is neglected.
Our model has eight free parameters. There are two for the disc DF: the disc scalelength R_{d} and the central velocity dispersion σ_{R, 0}); four for the halo DF: the slope B, scale action J_{0}, anisotropy coefficient k_{R}, and rotation parameter χ; and two for the dark matter halo potential: the axis ratio q of the density distribution and, instead of varying M_{h} in Eq. (6), we use the total mass of the Galaxy within 20 kpc, which is a characteristic radius probed by the tracers, M_{20} = M_{20,baryons} + M_{20,DM}, where M_{20,baryons} = 5.4×10^{10} M_{⊙} is a constant that is fixed by the baryonic model in Sect. 3.1. Following Bowden et al. (2016), we transform the axis ratio to the quantity
which varies in a finite range [0, 1]. A spherical model has u = 1/2, oblate halos have 0 < u < 0.5, while prolate ones have 0.5 < u < 1.
4.1. Prior
We adopt noninformative priors for all parameters. Three of these have finite domains, hence we use uniform priors for them:
The other five are intrinsically positive quantities, hence noninformative priors are uniform in the logarithm:
Thus, our prior P(Ξ) is the product of these eight terms.
4.2. Likelihood
The likelihood P(wΞ) that the N globular clusters in our dataset are moving in the gravitational potential Φ and are drawn from the phase space distribution f(JΦ) is
where the Jacobian factor is (see e.g. BW17)
In Eq. (12) we have convolved the DF of the model with the error distribution G_{j} of each cluster. G_{j} is a multivariate normal with mean given by the measurements and covariance matrix C, whose only nonnull offdiagonal element is given by the correlation coefficient C_{μα*, μδ}, which is measured only for the 75 clusters analysed with Gaia DR2. To compute the integral in Eq. (12) we use a similar strategy to that of BW17: we employ an importancesampling Monte Carlo method, with the same sampling density as in BW17 (their f_{S} as in their Sect. 3.3.3).
In this work we assume that our sample of GCs is complete (e.g. Harris et al. 2001), meaning that we neglect the possibility that a significant number of GCs is still hidden by dust close to the Galactic midplane. This approach is validated by BW17, who have shown that including a selection function depending on dust extinction (as an extra factor multiplying the integrand in Eq. (12)) does not alter the final inference on the model parameters.
4.3. Posterior
We sampled the posterior distribution P(Ξw) of the model parameters via a Markov chain Monte Carlo (MCMC) method. We used a MetropolisHastings sampler (Hastings 1970) with a multivariate proposal distribution. We ran 50 independent chains for 3000 steps, which were found to converge after a short burnin phase of about ∼600 steps, thus we effectively sampled the posterior distribution with ∼120 000 samples. All the chains turned out to be well mixed and we tuned the parameters of the proposal distribution to have ∼40 − 60% acceptance rate and small autocorrelations (after burnin).
5. Results
Figure 1 shows the 1D and 2D marginalised posterior distributions of the eight free parameters of our model. All the parameters are relatively well constrained by our analysis, with relative uncertainties limited to ∼10%. In what follows we always associate error bars with each measurement representing the 16th and 84th percentiles of the marginalised posterior distribution. Figure 1 depicts some nonzero correlations between the parameters which we discuss below.
Fig. 1. Posterior distributions for the eight model parameters. The black solid curves in the 2D marginalised subspaces enclose 68% and 95% of the total probability. The black vertical dashed lines in the 1D marginalised subspaces are the 16th and 84th percentiles. 
5.1. Distribution function
5.1.1. Halo
The halo DF of the best model has a shallow slope (B = 3.01 ± 0.05) and a smallscale action ( kpc km s^{−1}) and we find these two quantities to be correlated as models with a steeper slope have a larger scale action. This degeneracy is not surprising, and it is inevitable since J_{0} controls the physical scale at which the spatial density profile steepens. Our best model has a very small constant density core, of ∼0.1 kpc, and its density distribution is a power law close to ρ ∝ r^{−3.3}.
Contrary to previous claims (e.g. BW17), we find the halo component not to be significantly rotating (κ = −0.14 ± 0.14), if not mildly counterrotating. In fact, in our best model the halo component has a small net retrograde rotation of the order of V_{rot} ∼ −14 km s^{−1} at 20 kpc (roughly consistent with other independent estimates, e.g. Beers et al. 2012; Kafle et al. 2017; Koppelman et al. 2018), while the disc component is rotating at about V_{rot} ∼ 210 km s^{−1} at the solar radius. This difference compared to previous work is due to the dramatic increase in data quality.
We also constrain relatively well the coefficient k_{R} of the homogeneous function g(J) entering in the halo DF definition, and which controls the velocity anisotropy of the halo component. Our best model with has a mildly radially biased velocity distribution, with a rather constant , consistent with previous measurements (e.g. Eadie et al. 2017). We show in Fig. 2 the spherically averaged anisotropy profile of the halo component.
Fig. 2. Spherically averaged velocity anisotropy profile of the halo component of the GC system. The black solid line is for our best model and the grey band encompasses the 16th and 84th percentiles of the distribution. 
5.1.2. Disc
The disc DF is consistent with that of the Galactic thick disc fitted by Piffl et al. (2014b): the scalelength of the disc is kpc (see Table 1), while the central radial velocity dispersion is km s^{−1}. However, given that we have fixed two of the DF parameters, R_{σ} and h_{d}, to resemble the thick disc, it comes as no surprise that we find agreement with the estimates by Piffl et al. (2014b). This serves more as a sanity check that the MCMC chains are converging to a physical solution, instead of a proper constraint. The fact that we find the most likely models precisely in the region of the parameter space where the Galactic thick disc is, implies that our analysis is consistent with the possibility that some GCs were born in the thick disc or were associated with its formation (e.g. Zinn 1985).
5.1.3. Spatial and velocity distribution
We now compare the observed and predicted spatial and velocity distributions of the GC system. In Fig. 3 we show the heliocentric distance distribution of the full sample of 143 GCs compared to that of our best model. For visualisation purposes, we have cut both distributions at 50 kpc. The model represents very well the observed distances to the cluster: running a Kolmogorov–Smirnov test on the full sample (143 GCs) we find a probability^{1} that they are drawn from the same distribution. It turns out that most of the numerical noise in this comparison comes from the six clusters found at distances between 50 and 120 kpc: if we exclude this large and very poorly sampled volume and we restrict the comparison to the clusters within 50 kpc we find a probability of that the data and model come from the same underlying distribution. Within 50 kpc, the largest discrepancy is found between 10 and 30 kpc, but it is not significant (see inset in Fig. 3).
Fig. 3. Cumulative distribution of heliocentric distances of the 143 GCs (red solid) compared to that of the best model (black dashed). We cut these at 50 kpc for illustration purposes. The Kolmogorov–Smirnov test returns a 63% probability that the two samples are drawn from the same distribution. The inset shows the normalised distribution of distances within 30 kpc with Poissonian error bars. 
For the first time we have a sufficiently large dataset of accurate proper motions of 91 GCs, and we can also make a meaningful comparison with the velocity distribution predicted by the model. We show in Fig. 4 the cumulative distributions of three cartesian heliocentric velocities (V_{x}, V_{y}, V_{z}), where V_{x} points towards the Galactic centre and V_{Y} is parallel to the Galactic disc’s rotation. Figure 5 also shows the 2D projections of the clusters velocity distribution compared to the models.
Fig. 4. Cumulative distribution of heliocentric cartesian velocities of the GCs (red solid) compared to that of the best model (black dashed). The Kolmogorov–Smirnov test returns a 88%, 83%, and 42% probability that the two samples in V_{x}, V_{y}, and V_{z}, respectively, are drawn from the same distributions. 
Fig. 5. Velocity distribution in the V_{x} − V_{y} and V_{z} − V_{y} planes of the model (black) compared to that of the data (red circles). 
The agreement between data and model is generally very good; however, it should be noted that the model appears somewhat more symmetric with respect to the observed velocity distribution of GCs. On the one hand, it should be kept in mind that our equilibrium models will always produce symmetric distributions in V_{x} and V_{z}; on the other hand, observed deviations from symmetry may be due to the presence of substructures or may simply be stochastic, due to the small number of tracers (91).
5.1.4. Metallicity distribution
For each model we can use Eq. (12) to compute the probability that a cluster belongs to the disc component rather than to the halo by substituting f with f_{disc} or f_{halo}, respectively. The ratio of these two probabilities, which we call P_{disc} and P_{halo}, can be used to estimate the membership of each GC to these components. We now compare these membership probabilities with the metallicity [Fe/H] of each cluster, whose distribution is bimodal with a minimum at [Fe/H] ≃ −0.8 that has been classically used to distinguish metalpoor halo clusters from metalrich disc ones (e.g. Zinn 1985).
Figure 6 shows the metallicity of the 143 GCs as a function of the logarithm of the ratio P_{disc}/P_{halo} in our best model. Similarly to BW17, we find that metalrich clusters are typically much more likely to belong to the disc component rather than the halo, with only a couple of exceptions. In particular, all but two clusters with [Fe/H] > −0.8 have log_{10}(P_{disc}/P_{halo}) > −2, making it very unlikely that they are a part of the halo component (see histogram in the bottom panel of Fig. 6). At lower metallicities, instead, there are many more clusters with log_{10}(P_{disc}/P_{halo}) < −2 (67% below [Fe/H]≤ − 0.8), making them very likely halo clusters. These are typically objects that are either found at very large distances, that are counterrotating, or that have negligible angular momentum. Figure 6 shows that only a few clusters with heliocentric distances larger than 20 kpc have a higher probability of being part of the disc component.
Fig. 6. Upper panel: metallicity of the GCs as a function of the ratio of the probability of belonging to the disc or the halo component of the DF. The points are colourcoded according to the heliocentric distance of the cluster. The horizontal dashed line marks [Fe/H] = −0.8, the classical distinction between metalpoor and metalrich clusters. Lower panel: histogram of the ratio of probabilities in two bins of metalrich ([Fe/H] > −0.8) and metalpoor ([Fe/H] ≤ −0.8) clusters. 
In the lefthand panel of Fig. 7 we show the distribution of the 91 clusters with accurate proper motions in the plane of integrals of motion angular momentumenergy, colourcoded by log_{10}(P_{disc}/P_{halo}). Here we see that all clusters less bound than E ≳ −1.25×10^{5} km s^{−1} belong to the halo component. At lower (more negative) energies, halo clusters clearly prefer retrograde (J_{ϕ} < 0) or nonrotating (J_{ϕ} ∼ 0) orbits, while the probability of belonging to the disc component is maximal for prograde angular momenta. For the likely disc clusters (log_{10} P_{disc}/P_{halo} > 0) we find an average angular momentum of , while for the very likely halo clusters (log_{10} P_{disc}/P_{halo} < 0) we find , from which we conclude that the halo clusters exhibit a hint of retrograde rotation. If we divide the sample of 91 GCs according to the traditional separation into metalrich [Fe/H] > −0.8 (18 clusters) and metalpoor [Fe/H]≤ − 0.8 (73 clusters), we find a mean zangular momentum of respectively and for the two populations. This is because some metalrich clusters are on lowenergy, sometimes even retrograde orbits, while a few clusters on disclike orbits happen to have low metallicity.
Fig. 7. Left panel: distribution of the 91 GCs with accurate proper motions in the angular momentumenergy plane, colourcoded by the ratio of probability of belonging to the disc or to the halo. Right panel: same as the left panel, but in the metallicityenergy plane. 
The righthand panel of Fig. 7 shows energy as a function of metallicity, with the points being again colourcoded by their membership probability. It can be seen that the metalpoor halo clusters are mostly confined to less bound orbits, while disc clusters with high angular momentum span a wide range of metallicities.
5.2. Gravitational potential
With our modelling technique we are able to constrain the characteristic parameters of the Galaxy’s inner dark matter halo that fit the observed dynamics of GCs. The main novelty of our work is the fact that the new dataset released by the Gaia Collaboration (2018b) gives us strong inference on both the total mass and shape of the dark matter halo. In particular, in Fig. 8 we show the posterior distribution of the Galaxy’s mass within 20 kpc and the halo’s axis ratio q as sampled by our MCMC analysis. We find that log_{10} M_{20}/M_{⊙} = 11.33 ± 0.05 and and that they are strongly correlated with more massive halos being more prolate. While such a degeneracy is expected (e.g. Bowden et al. 2016), here we nonetheless find that spherical and light halos are disfavoured by our analysis, while halos more oblate than q ≤ 0.8 are ruled out with 99% probability by our experiment, as are strongly prolate halos with q ≥ 1.9.
Fig. 8. Posterior distributions of the two parameters describing the Galaxy’s dark halo: the total mass of the Galaxy within 20 kpc, M_{20}, and the axis ratio of the dark halo, q. The number of bins is increased with respect to Fig. 1. Black solid curves and vertical dashed lines are as in Fig. 1. 
The posterior distribution that we sample is not unimodal, but it has two distinct peaks; the most likely is at log_{10} M_{20} ≃ 11.33, q ≃ 1.37 and the other at log_{10} M_{20} ≃ 11.38, q ≃ 1.75. Although it is located in a region where the likelihood is considerably lower, the second peak cannot be neglected and deserves further investigation.
As discussed in the introduction, an important number of GCs in our sample are likely associated with the Sagittarius dwarf galaxy and have thus been accreted on similar orbits. In this case, our assumption that all the GCs are a random sample of our DF breaks down, and this could lead to biases in the results of the modelling. We therefore remove the clusters mentioned in Sect. 2.2, and repeat the analysis outlined in Sect. 4. As a result we obtain the posterior distribution for the halo parameters shown in Fig. 9. Now the second lower peak has completely disappeared and the distribution has a wellconstrained peak at log_{10} M_{20}/M_{⊙} = 11.28 ± 0.04 and q = 1.22 ± 0.23. From this we conclude that the second peak was most probably driven by the few clusters on very polar orbits dynamically associated with the Sagittarius dwarf galaxy.
Within 20 kpc, we find the total mass of the Galaxy to be , which implies a dark matter mass of . If we extrapolate this to the virial radius^{2}, assuming a concentrationmass relation as calibrated from numerical simulations (Dutton & Macciò 2014), we obtain the following final estimate of the Galaxy’s dark matter halo mass,
and thus a virial radius of kpc.
Many different studies have already estimated the mass of the dark matter halo of the Milky Way using very different techniques and tracers, which typically are sensitive to very different physical scales. Thus, a proper comparison taking into account all the possible biases given by the different assumptions or data used is needed, but it goes beyond the scope of the present study. By naively comparing our results with the numerical estimates compiled by BlandHawthorn & Gerhard (2016, and references therein) and Eadie & Harris (2016, and references therein) we conclude that our estimate of M_{20} and of the extrapolated M_{vir} lies somewhere between the very light (M_{vir} < 10^{12} M_{⊙}, e.g. Battaglia et al. 2005; Gibbons et al. 2014) and of the very heavy values (M_{vir} > 2×10^{12} M_{⊙} Li & White 2008; Watkins et al. 2010). In general, we find good agreement with other works that use satellite kinematics (provided that distant enough tracers are included, e.g. BoylanKolchin et al. 2013; Gonzalez et al. 2013; Eadie & Harris 2016; Patel et al. 2018; Sohn et al. 2018), estimates of the escape velocity (e.g. Smith et al. 2007; Piffl et al. 2014a), the velocity distribution of either fast moving stars (e.g. Gnedin et al. 2010), or the dynamics of nearby stars (e.g. McMillan 2011; Piffl et al. 2014b).
Similarly, there has also been some work on the shape of the total gravitational potential. While a proper comparison taking into account the systematics induced by tracers and methods is not straightforward, here we note that our estimate roughly agrees with other studies of the Sagittarius stream (Helmi 2004), the flaring of the HI disc (Banerjee & Jog 2011), or the kinematics of halo stars (Bowden et al. 2016), while being apparently inconsistent with studies modelling the GD1 stream (Koposov et al. 2010; Bowden et al. 2015, which, however, probes a region much closer to the Galactic centre) and studies of the tilt of the local velocity ellipsoid (Smith et al. 2009; Posti et al. 2018).
It may well be the case that the shape of the halo of our Galaxy is considerably more complicated than we have assumed here; for instance, i) it can be axisymmetric, but with a variable flattening as a function of radius (e.g. VeraCiro & Helmi 2013), being oblate close to the disc and more spherical/prolate at greater distances (e.g. as probed by the HI disc flare, Banerjee & Jog 2011); ii) it can be triaxial (e.g. as probed by the modelling of the Sagittarius stream Law & Majewski 2010); or iii) it can be triaxial, but with variable axis ratios/orientations as a function of radius (e.g. as expected from numerical simulations, Bailin & Steinmetz 2005; VeraCiro et al. 2011, and as probed by the Sagittarius stream, VeraCiro & Helmi 2013). Different results are thus found both because of different model assumptions and because of the various observables used: since different tracers are probably sensitive to distinct portions of the halo, the methodologies used to infer the halo shape may be subject to different kinds of limitations. This implies that to make further progress in our understanding of the shape of the dark halo of the Milky Way it is imperative to combine the signal coming from tracers of different nature.
6. Conclusions
We have used the kinematics of Milky Way GCs to simultaneously constrain the phasespace distribution of the GC system, the mass of the inner Galaxy, and the shape of its inner dark matter halo. This was possible thanks to the novel measurements by Gaia and HST of the proper motions of 91 GCs with unprecedented accuracy.
Our method is based on the assumption that the GC system can be described by axisymmetric equilibrium models with two distribution functions, which are analytic functions of the action integrals, representing the phasespace distribution of disc and halo clusters. Our models include an axisymmetric NFW dark matter halo, whose mass and density axis ratio, have been determined by fitting the kinematics of 143 GCs (91 with accurate 6D data and 52 with only radial velocities), using a Bayesian approach, imposing uninformative priors on the eight model characteristic parameters.
We summarise here our main results:

(i)
we find all the parameters of the halo component to be well constrained. The density distribution is well described by a single power law ρ ∝ r^{−3.3}; the velocity distribution is mildly radially biased, with a rather constant anisotropy β ≃ 0.2 ± 0.07, and consistent with either a small counterrotation (V_{rot} ≃ −15km s^{−1} at 20 kpc) or no net rotation;

(ii)
we find all the parameters of the disc component to be well constrained. Their phasespace distribution closely resembles that of the thick disc of the Galaxy (Piffl et al. 2014b);

(iii)
metalrich clusters ([Fe/H] > −0.8) very likely belong to the disc component, while roughly 67% of the metalpoor clusters ([Fe/H] ≲ −0.8) are more likely part of the halo component. All clusters with low binding energies belong to the halo component and have both prograde and retrograde orbits. The more bound halo clusters are typically found on counterrotating or nonrotating orbits. Most metalrich clusters ([Fe/H] ≳ −0.8) have high a probability of being associated with the disc component. Our findings are broadly consistent with the picture proposed by Zinn (1985);

(iv)
we measure the mass of the dark matter halo of the Galaxy within 20 kpc to be log_{10} M_{20,DM}/M_{⊙} = 11.14±0.05. This very accurate measurement implies a total virial mass for the Milky Way of M_{vir} = 1.3±0.3×10^{12} M_{⊙} after assuming a concentrationvirial mass relation;

(iv)
we measure a flattening q = 1.22 ± 0.23, and find oblate halos with q < 0.7 to be ruled out by our analysis (with 99% probability) and spherical models to be disfavoured in comparison to prolate halos.
During the completion of this work, Watkins et al. (2018) presented an estimate of the Milky Way’s dark halo mass by applying a simple “tracer mass estimator” on a similar dataset. Their values for the mass within 21.1 kpc (and within the virial radius) are closely compatible with ours, as is their measurement of the anisotropy (although they favour a slightly more radially biased ellipsoid).
Our work demonstrates that with a giant leap in data quality such as that provided by the Gaia DR2, constraints on fundamental physical quantities become significantly more precise, and several of these parameters can be determined at once if models with enough sophistication are employed. Of course, with this work we have simply touched the tip of the iceberg of the signal present in the Gaia data regarding the mass and shape of the Galactic potential. The next challenge will be to simultaneously fit the dynamics of different mass tracers such as satellites, streams, field stars in the halo, and hypervelocity stars. This work may be seen as a necessary first step in this direction, and merely confirms that we may have just entered the era of “Precision Galactic Astronomy”.
We adopt the Bryan & Norman (1998) definition of the virial radius, with a virial overdensity at redshift zero of Δ ≃ 102.5.
Acknowledgments
Thank the referee, Matthias Steinmetz, for a constructive report and Simon White for useful discussions. We acknowledge financial support from a VICI grant from the Netherlands Organisation for Scientific Research (NWO). This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
References
 Antoja, T., Helmi, A., RomeroGomez, M., et al. 2018, Nature, 561, 360 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., & Tremaine, S. 1981, ApJ, 244, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Bailin, J., & Steinmetz, M. 2005, ApJ, 627, 647 [Google Scholar]
 Banerjee, A., & Jog, C. J. 2011, ApJ, 732, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Battaglia, G., Helmi, A., Morrison, H., et al. 2005, MNRAS, 364, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Beers, T. C., Carollo, D., Ivezić, Z., et al. 2012, ApJ, 746, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J. 2010, MNRAS, 401, 2318 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J. 2012, MNRAS, 426, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton University Press) [Google Scholar]
 Binney, J., & Wong, L. K. 2017, MNRAS, 467, 2446 [NASA ADS] [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bowden, A., Belokurov, V., & Evans, N. W. 2015, MNRAS, 449, 1391 [NASA ADS] [CrossRef] [Google Scholar]
 Bowden, A., Evans, N. W., & Williams, A. A. 2016, MNRAS, 460, 329 [NASA ADS] [CrossRef] [Google Scholar]
 BoylanKolchin, M., Bullock, J. S., Sohn, S. T., Besla, G., & van der Marel, R. P. 2013, ApJ, 768, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Carollo, D., Beers, T. C., Lee, Y. S., et al. 2007, Nature, 450, 1020 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Das, P., & Binney, J. 2016, MNRAS, 460, 1725 [NASA ADS] [CrossRef] [Google Scholar]
 Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903 [NASA ADS] [CrossRef] [Google Scholar]
 Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [NASA ADS] [CrossRef] [Google Scholar]
 Eadie, G. M., & Harris, W. E. 2016, ApJ, 829, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Eadie, G. M., Springford, A., & Harris, W. E. 2017, ApJ, 835, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018a, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Helmi, A., et al.) 2018b, A&A, 616, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Katz, D., et al.) 2018c, A&A, 616, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2014, MNRAS, 445, 3788 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, O. Y., Brown, W. R., Geller, M. J., & Kenyon, S. J. 2010, ApJ, 720, L108 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2013, A&A, 552, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harris, W. E. 1996, AJ, 112, 1487 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, W. E. 2001, in SaasFee Advanced Course 28, eds. L. Labhardt, & B. Binggeli (Berlin: Star Clusters, SpringerVerlag), 223 [NASA ADS] [CrossRef] [Google Scholar]
 Hastings, W. K. 1970, Biometrika, 57, 97 [Google Scholar]
 Helmi, A. 2004, ApJ, 610, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194 [NASA ADS] [CrossRef] [Google Scholar]
 Iorio, G., Belokurov, V., Erkal, D., et al. 2018, MNRAS, 474, 2142 [NASA ADS] [CrossRef] [Google Scholar]
 Kafle, P. R., Sharma, S., Robotham, A. S. G., et al. 2017, MNRAS, 470, 2959 [NASA ADS] [CrossRef] [Google Scholar]
 Koposov, S. E., Rix, H.W., & Hogg, D. W. 2010, ApJ, 712, 260 [NASA ADS] [CrossRef] [Google Scholar]
 Koppelman, H. H., Helmi, A., & Veljanoski, J. 2018, ApJ, 860, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Law, D. R., & Majewski, S. R. 2010, ApJ, 718, 1128 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Y.S., & White, S. D. M. 2008, MNRAS, 384, 1459 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J. 2011, MNRAS, 414, 2446 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Olling, R. P., & Merrifield, M. R. 2000, MNRAS, 311, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Patel, E., Besla, G., Mandel, K., & Sohn, S. T. 2018, ApJ, 857, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Piffl, T., Scannapieco, C., Binney, J., et al. 2014a, A&A, 562, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piffl, T., Binney, J., McMillan, P. J., et al. 2014b, MNRAS, 445, 3133 [NASA ADS] [CrossRef] [Google Scholar]
 Posti, L., Binney, J., Nipoti, C., & Ciotti, L. 2015, MNRAS, 447, 3060 [Google Scholar]
 Posti, L., Helmi, A., Veljanoski, J., & Breddels, M. A. 2018, A&A, 615, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sesar, B., Ivezic, Ž., Stuart, J. S., et al. 2013, AJ, 146, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007, MNRAS, 379, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, M. C., Wyn, Evans N., & An, J. H. 2009, ApJ, 698, 1110 [NASA ADS] [CrossRef] [Google Scholar]
 Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, ApJ, 862, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E. 2019, MNRAS, 482, 1525 [NASA ADS] [CrossRef] [Google Scholar]
 VeraCiro, C., & Helmi, A. 2013, ApJ, 773, L4 [NASA ADS] [CrossRef] [Google Scholar]
 VeraCiro, C. A., Sales, L. V., Helmi, A., et al. 2011, MNRAS, 416, 1377 [NASA ADS] [CrossRef] [Google Scholar]
 Watkins, L. L., Evans, N. W., & An, J. H. 2010, MNRAS, 406, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2018, AAS, submitted [arXiv:1804.11348] [Google Scholar]
 Williams, A. A., & Evans, N. W. 2015, MNRAS, 448, 1360 [NASA ADS] [CrossRef] [Google Scholar]
 Xue, X.X., Rix, H.W., Ma, Z., et al. 2015, ApJ, 809, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Zinn, R. 1985, ApJ, 293, 424 [Google Scholar]
All Tables
All Figures
Fig. 1. Posterior distributions for the eight model parameters. The black solid curves in the 2D marginalised subspaces enclose 68% and 95% of the total probability. The black vertical dashed lines in the 1D marginalised subspaces are the 16th and 84th percentiles. 

In the text 
Fig. 2. Spherically averaged velocity anisotropy profile of the halo component of the GC system. The black solid line is for our best model and the grey band encompasses the 16th and 84th percentiles of the distribution. 

In the text 
Fig. 3. Cumulative distribution of heliocentric distances of the 143 GCs (red solid) compared to that of the best model (black dashed). We cut these at 50 kpc for illustration purposes. The Kolmogorov–Smirnov test returns a 63% probability that the two samples are drawn from the same distribution. The inset shows the normalised distribution of distances within 30 kpc with Poissonian error bars. 

In the text 
Fig. 4. Cumulative distribution of heliocentric cartesian velocities of the GCs (red solid) compared to that of the best model (black dashed). The Kolmogorov–Smirnov test returns a 88%, 83%, and 42% probability that the two samples in V_{x}, V_{y}, and V_{z}, respectively, are drawn from the same distributions. 

In the text 
Fig. 5. Velocity distribution in the V_{x} − V_{y} and V_{z} − V_{y} planes of the model (black) compared to that of the data (red circles). 

In the text 
Fig. 6. Upper panel: metallicity of the GCs as a function of the ratio of the probability of belonging to the disc or the halo component of the DF. The points are colourcoded according to the heliocentric distance of the cluster. The horizontal dashed line marks [Fe/H] = −0.8, the classical distinction between metalpoor and metalrich clusters. Lower panel: histogram of the ratio of probabilities in two bins of metalrich ([Fe/H] > −0.8) and metalpoor ([Fe/H] ≤ −0.8) clusters. 

In the text 
Fig. 7. Left panel: distribution of the 91 GCs with accurate proper motions in the angular momentumenergy plane, colourcoded by the ratio of probability of belonging to the disc or to the halo. Right panel: same as the left panel, but in the metallicityenergy plane. 

In the text 
Fig. 8. Posterior distributions of the two parameters describing the Galaxy’s dark halo: the total mass of the Galaxy within 20 kpc, M_{20}, and the axis ratio of the dark halo, q. The number of bins is increased with respect to Fig. 1. Black solid curves and vertical dashed lines are as in Fig. 1. 

In the text 
Fig. 9. Same as Fig. 8, but obtained excluding the GCs associated with the Sagittarius dwarf galaxy. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.