© 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 ground-breaking results (Gaia Collaboration 2018b,c; Antoja et al. 2018). In this study we exploit this novel dataset by employing a well-established, 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⁻¹, which roughly corresponds to an accuracy of 0.5−2kms⁻¹ 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 phase-space 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 phase-space 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 halo-like dynamics (typically representing the metal-poor clusters), and one with disc-like dynamics (following the more metal-rich 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⁻¹ (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 Bland-Hawthorn & 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 phase-space 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 extinction-corrected distance moduli in the V-band, 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

(1)

where the error distribution of each cluster is assumed to be a multi-variate normal distribution with non-null 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⁻¹ 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 non-random sampling of Galactic phase-space which could lead to biased estimates of the mass and shape of the Galactic potential using our DF-based 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 position-velocities (x, v) to action-angles (θ, 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 double-exponential density distributions

(2)

where R_d and h_d are scale-length and -height, Σ₀ is a normalisation constant, and R_hole is the size of the central cavity (non-zero only for the gaseous disc), while the bulge and halo are described by

(3)

where , r₀ is a scale radius, q is the axis ratio, and ρ₀ 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.

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 cut-off at the virial radius, r_cut = r_vir, and imposing the concentration-mass relation found in ΛCDM N-body simulations: log₁₀ c = 1.025 − 0.097log₁₀(M/10¹² h⁻¹ M_⊙), where c = r_vir/r₀ 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 phase-space 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 π)³ 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 phase-space distribution, we allow the DF to have two dynamically distinct components, one disc-like and one halo-like: f(J) = f_disc (J)+f_halo(J). For the disc component, we use the “quasi-isothermal” DF, introduced by Binney (2010), as implemented in Vasiliev (2019)

(4)

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 Σ = Σ₀ exp[−R_c(J_ϕ)/R_d] and

(5)

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 scale-height, , and a radial velocity dispersion σ_R = σ_R, 0exp(−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 scale-length R_d and the central radial dispersion σ_R, 0, that completely characterise f_disc(J).

For the halo component, we use the double power-law 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 phase-space, thus making the 3D density distribution of the model a single power law with a central core; the halo DF then is

(6)

where J₀ is a scale action defining the size of the constant density core, J_cut is a (large) cut-off action, g(J) = k_RJ_R + (3−k_R)(|J_ϕ| + J_z)/2 is a homogeneous function of degree one and

(7)

is the factor controlling the net rotation of the system, with χ = 0 being the non-rotating case and χ = ±1 being the maximally rotating/counter-rotating case. Here we fix two parameters: i) the cut-off action to J_cut = 10⁵ kpc km s⁻¹, 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⁻¹, 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₀, 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 Σ₀ for f_disc and M_h for f_halo are unimportant in our case, since the GC system is treated as mass-less 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

(8)

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 scale-length R_d and the central velocity dispersion σ_R, 0); four for the halo DF: the slope B, scale action J₀, 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₂₀ = M_20,baryons + M_20,DM, where M_20,baryons = 5.4×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

(9)

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 non-informative priors for all parameters. Three of these have finite domains, hence we use uniform priors for them:

(10)

The other five are intrinsically positive quantities, hence non-informative priors are uniform in the logarithm:

(11)

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

(12)

where the Jacobian factor is (see e.g. BW17)

(13)

In Eq. (12) we have convolved the DF of the model with the error distribution G_j of each cluster. G_j is a multi-variate normal with mean given by the measurements and covariance matrix C, whose only non-null off-diagonal 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 importance-sampling 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 mid-plane. 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 Metropolis-Hastings sampler (Hastings 1970) with a multi-variate proposal distribution. We ran 50 independent chains for 3000 steps, which were found to converge after a short burn-in 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 auto-correlations (after burn-in).

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 non-zero 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 small-scale action ( kpc km s⁻¹) 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₀ 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 counter-rotating. In fact, in our best model the halo component has a small net retrograde rotation of the order of V_rot ∼ −14 km s⁻¹ 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⁻¹ 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 scale-length of the disc is kpc (see Table 1), while the central radial velocity dispersion is km s⁻¹. 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¹ 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 Vx, Vy, and Vz, respectively, are drawn from the same distributions.

Fig. 5. Velocity distribution in the Vx − Vy and Vz − Vy 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 metal-poor halo clusters from metal-rich 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 metal-rich 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₁₀(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₁₀(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 counter-rotating, 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 colour-coded according to the heliocentric distance of the cluster. The horizontal dashed line marks [Fe/H] = −0.8, the classical distinction between metal-poor and metal-rich clusters. Lower panel: histogram of the ratio of probabilities in two bins of metal-rich ([Fe/H] > −0.8) and metal-poor ([Fe/H] ≤ −0.8) clusters.

In the left-hand panel of Fig. 7 we show the distribution of the 91 clusters with accurate proper motions in the plane of integrals of motion angular momentum-energy, colour-coded by log₁₀(P_disc/P_halo). Here we see that all clusters less bound than E ≳ −1.25×10⁵ km s⁻¹ belong to the halo component. At lower (more negative) energies, halo clusters clearly prefer retrograde (J_ϕ <  0) or non-rotating (J_ϕ ∼ 0) orbits, while the probability of belonging to the disc component is maximal for prograde angular momenta. For the likely disc clusters (log₁₀ P_disc/P_halo > 0) we find an average angular momentum of , while for the very likely halo clusters (log₁₀ 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 metal-rich [Fe/H] > −0.8 (18 clusters) and metal-poor [Fe/H]≤ − 0.8 (73 clusters), we find a mean z-angular momentum of respectively and for the two populations. This is because some metal-rich clusters are on low-energy, sometimes even retrograde orbits, while a few clusters on disc-like orbits happen to have low metallicity.

Fig. 7. Left panel: distribution of the 91 GCs with accurate proper motions in the angular momentum-energy plane, colour-coded by the ratio of probability of belonging to the disc or to the halo. Right panel: same as the left panel, but in the metallicity-energy plane.

The right-hand panel of Fig. 7 shows energy as a function of metallicity, with the points being again colour-coded by their membership probability. It can be seen that the metal-poor 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₁₀ M₂₀/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, M20, 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₁₀ M₂₀ ≃ 11.33, q ≃ 1.37 and the other at log₁₀ M₂₀ ≃ 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 well-constrained peak at log₁₀ M₂₀/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.

Fig. 9. Same as Fig. 8, but obtained excluding the GCs 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², assuming a concentration-mass relation as calibrated from numerical simulations (Dutton & Macciò 2014), we obtain the following final estimate of the Galaxy’s dark matter halo mass,

(14)

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 Bland-Hawthorn & Gerhard (2016, and references therein) and Eadie & Harris (2016, and references therein) we conclude that our estimate of M₂₀ and of the extrapolated M_vir lies somewhere between the very light (M_vir < 10¹² M_⊙, e.g. Battaglia et al. 2005; Gibbons et al. 2014) and of the very heavy values (M_vir > 2×10¹² 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. Boylan-Kolchin 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 GD-1 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. Vera-Ciro & 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; Vera-Ciro et al. 2011, and as probed by the Sagittarius stream, Vera-Ciro & 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 phase-space 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 phase-space 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 counter-rotation (V_rot ≃ −15km s⁻¹ at 20 kpc) or no net rotation;

• (ii)

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

• (iii)

  metal-rich clusters ([Fe/H] > −0.8) very likely belong to the disc component, while roughly 67% of the metal-poor 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 counter-rotating or non-rotating orbits. Most metal-rich 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₁₀ 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¹² M_⊙ after assuming a concentration-virial 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”.

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

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 Vx, Vy, and Vz, respectively, are drawn from the same distributions.
In the text

	Fig. 5. Velocity distribution in the Vx − Vy and Vz − Vy 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 colour-coded according to the heliocentric distance of the cluster. The horizontal dashed line marks [Fe/H] = −0.8, the classical distinction between metal-poor and metal-rich clusters. Lower panel: histogram of the ratio of probabilities in two bins of metal-rich ([Fe/H] > −0.8) and metal-poor ([Fe/H] ≤ −0.8) clusters.

In the text

	Fig. 7. Left panel: distribution of the 91 GCs with accurate proper motions in the angular momentum-energy plane, colour-coded by the ratio of probability of belonging to the disc or to the halo. Right panel: same as the left panel, but in the metallicity-energy 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, M20, 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