Modelling the Milky Way's globular cluster system

TL;DR

This paper develops a detailed dynamical model of the Milky Way's globular cluster system using a Bayesian approach, revealing insights into its structure, rotation, and relation to stellar components.

Contribution

It introduces a novel action-based distribution function model for the globular cluster system, integrating disc and halo components with Bayesian parameter estimation.

Findings

Model reproduces density profile and rotation velocity of clusters

Disc and halo parameters align with stellar component estimates

Identifies potential issues with cluster distance estimates

Abstract

We construct a model for the Galactic globular cluster system based on a realistic gravitational potential and a distribution function (DF) analytic in the action integrals. The DF comprises disc and halo components whose functional forms resemble those recently used to describe the stellar discs and stellar halo. We determine the posterior distribution of our model parameters using a Bayesian approach. This gives us an understanding of how well the globular cluster data constrain our model. The favoured parameter values of the disc and halo DFs are similar to values previously obtained from fits to the stellar disc and halo, although the cluster halo system shows clearer rotation than does the stellar halo. Our model reproduces the generic features of the globular cluster system, namely the density profile, the mean rotation velocity. The fraction of disc clusters coincides with the observed fraction of metal-rich clusters. However, the data indicate either incompatibility between catalogued cluster distances and current estimates of distance to the Galactic Centre, or failure to identify clusters behind the bulge. As the data for our Galaxy's components increase in volume and precision over the next few years, it will be rewarding to revisit the present analysis.

Tables

Table 1. Table 1: Parameters of our gravitational potential, which are fixed in our dynamical model.

	Thick disc	Thin disc	Gas disc
${\Sigma }_0$	274.5	532.4	87.3	$M_{\odot }{\text{pc}}^{-2}$
$R_{\text{d}}$	2.58	2.58	5.16	kpc
$z_{\text{d}}$	0.67	0.20	0.04	kpc
$R_{\text{hole}}$	0	0	4	kpc
	Stellar bulge	Dark halo
${\rho }_0$	94.9	0.0196	$M_{\odot }{\text{pc}}^{-3}$
$r_0$	0.075	15.5	kpc
$r_{\text{cut}}$	2.1	0	kpc
$a$	0	1
$b$	1.8	3
$q$	0.5	0.8

Table 2. Table 2: Expectation values and standard deviations of the parameters of the models in the MCMC chain.

Parameter	expectation	s.d.
$\alpha$	0.77	0.36
$\beta$	8.83	2.01
$J_0/\mathrm{kpc}\mathrm{km}{s}^{-1}$	5650	3140
$k_h$	0.33	0.20
$F_{\text{disc}}$	0.32	0.07
${\sigma }_r/\mathrm{km}{s}^{-1}$	94.3	22.5
${\sigma }_z/\mathrm{km}{s}^{-1}$	130.3	26.7
${ℛ}_d/\mathrm{kpc}$	1.90	0.39

Equations

$${\rho }_{\text{disc}}(R,z)=\frac{{\Sigma }_0}{2z_{\mathrm{d}}}\exp \left(-\frac{R}{R_{\mathrm{d}}}-\frac{R_{\text{hole}}}{R}-\frac{|z|}{z_{\text{d}}}\right),$$

$${\rho }_{\text{halo}}(R,z)=\frac{{\rho }_0}{m^a(1+m{)}^{b-a}}\exp \left[-(m{r}_0/r_{\text{cut}}{)}^2\right],$$

$$m(R,z)=\sqrt{(R/r_0{)}^2+(z/q{r}_0{)}^2}.$$

$$H=\frac{1}{2}\left(v_R^2+J_{\varphi }^2/R^2+v_z^2\right)+\Phi (R,z).$$

$$\mathbf{u}\equiv (l,b,s,v_{\mathrm{los}},{\mu }_{\alpha }^{\ast },{\mu }_{\delta }),$$

$$(\mathbf{x},\mathbf{v})\equiv (R,\varphi ,z,v_R,v_{\varphi },v_z),$$

$$f(\mathbf{J}){\mathrm{d}}^3\mathbf{\theta }{\mathrm{d}}^3\mathbf{J}=f(J(\mathbf{x},\mathbf{v})){\mathrm{d}}^3\mathbf{x}{\mathrm{d}}^3\mathbf{v}.$$

$$(2\pi {)}^3\int {\mathrm{d}}^3\mathbf{J}f(\mathbf{J})=1,$$

$$f(\mathbf{J}|\Pi )=F_{\text{disc}}{f}_{\text{disc}}(\mathbf{J})+(1-F_{\text{disc}})f_{\text{halo}}(\mathbf{J}),$$

$$f_{-}(\mathbf{J})=k\tanh (J_{\varphi }/L)f_{+}(\mathbf{J}).$$

$$f_{\text{disc+}}(\mathbf{J})=\frac{\Omega \nu \Sigma }{2{\pi }^2\kappa {\sigma }_r^2{\sigma }_z^2}\exp \left(-\frac{\kappa J_r}{{\sigma }_r^2}-\frac{\nu J_z}{{\sigma }_z^2}\right)$$

$$\Sigma (J_{\varphi })={\Sigma }_0\exp \left[-\frac{R_{\text{c}}(J_{\varphi })}{{\mathcal{R}}_{\text{d}}}\right].$$

$${\sigma }_i={\sigma }_{i0}\exp \left[(R_0-R_{\text{c}})\gamma \right],$$

$$h(\mathbf{J})=J_{\text{core}}+J_r+\frac{\Omega }{\kappa }(|J_{\varphi }|+J_z),$$

$$\begin{array}{rl}{f}_{\text{halo+}}(\mathbf{J})=& N_0\frac{[1+J_0/h(\mathbf{J}){]}^{\alpha }}{[1+h(\mathbf{J})/J_0{]}^{\beta }}\exp \left[-{\left(\frac{h(\mathbf{J})}{J_{\text{max}}}\right)}^4\right]\end{array}$$

$$\text{Pr}(\Pi |\text{Data})\propto \text{Pr}(\text{Data}|\Pi )\times \text{Pr}(\Pi ),$$

$$\text{Pr}(\beta )=\frac{B}{\Gamma (A)}(B\beta {)}^{A-1}{\mathrm{e}}^{-B\beta },$$

$$\textup{Pr}(\Pi)\propto\begin{cases}{\displaystyle\frac{\textup{Pr}(\beta)\,\textup{Pr}(\gamma)}{L_{\textup{d}}L_{\textup{h}}J_{0}\sigma_{r}\sigma_{z}\mathcal{R}_{\textup{d}}}}&{\raggedright\vbox{\noindent$k_{\textup{d}},k_{\textup{h}}\in[-1,1]$, $F_{\textup{disc}}\in[0,1]$, $L_{i},\alpha,\beta,J_{0},\sigma_{i0},\mathcal{R}_{\textup{d}}>0$,}\@add@raggedright}\\ 0&\textup{otherwise,}\end{cases}$$

$$G(u,\bar{u},{\sigma }_u)=\frac{1}{\sqrt{2\pi {\sigma }_u^2}}\exp \left[-\frac{(u-\bar{u}{)}^2}{2{\sigma }_u^2}\right],$$

$$G_L(s,\bar{s},{\sigma }_s)=\frac{1}{s}\frac{1}{\sqrt{2\pi {\sigma }_s^2}}\exp \left[-\frac{(\ln s-\mathcal{S}{)}^2}{2{\sigma }_s^2}\right],$$

$$\begin{array}{rl}G(\mathbf{u},{\bar{\mathbf{u}}}^{\eta },{\mathbf{\sigma }}^{\eta })=& \delta (l-{\bar{l}}^{\eta })\times \delta (b-{\bar{b}}^{\eta })\\ & \times G_L(s,{\bar{s}}^{\eta },{\sigma }_s^{\eta })\times G(v_{\mathrm{los}},{\bar{v}}_{\mathrm{los}}^{\eta },{\sigma }_{v_{\mathrm{los}}}^{\eta })\\ & \times G({\mu }_{\alpha }^{\ast },{\overline{{\mu }_{\alpha }^{\ast }}}^{\eta },{\sigma }_{{\mu }_{\alpha }^{\ast }}^{\eta })\times G({\mu }_{\delta },{\bar{\mu }}_{\delta }^{\eta },{\sigma }_{{\mu }_{\delta }}^{\eta })\end{array}$$

$$\mathcal{L}=\prod _{\eta }{\mathcal{L}}_{\eta }=\prod _{\eta }\int {\mathrm{d}}^6\mathbf{u}G(\mathbf{u},{\bar{\mathbf{u}}}^{\eta },{\mathbf{\sigma }}^{\eta })\mathrm{Pr}(\mathbf{u}|\text{Model}),$$

$$\mathrm{Pr}(\mathbf{u}|\text{Model})=AS(\mathbf{u})f(\mathbf{J})\left|\frac{\partial (\mathbf{\theta },\mathbf{J})}{\partial (\mathbf{u})}\right|,$$

$$\left|\frac{\partial (\mathbf{\theta },\mathbf{J})}{\partial (\mathbf{u})}\right|\equiv s^6\cos b.$$

$$1/A=\int {\mathrm{d}}^3\mathbf{x}{\mathrm{d}}^3\mathbf{v}f[\mathbf{J}(\mathbf{x},\mathbf{v})]S[\mathbf{u}(\mathbf{x},\mathbf{v})].$$

$$S(\mathbf{u})=\{\begin{array}{ll}1& \text{if }{\mathrm{E}(B-V)|}_{(l,b)}<\mathrm{E}(B-V{)}_{\text{max}}\\ 0& \text{otherwise.}\end{array}$$

$$1/A\simeq \frac{1}{N_{\text{S}}}\sum _{k=1}^{N_{\text{S}}}\frac{f(\mathbf{J}({\mathbf{x}}_k,{\mathbf{v}}_k))}{f_{\text{S}}({\mathbf{x}}_k,{\mathbf{v}}_k)}S[\mathbf{u}({\mathbf{x}}_k,{\mathbf{v}}_k)],$$

$${\mathcal{L}}_{\eta }\simeq \frac{A}{N_{\eta }{C}_{\eta }}\sum _{k=1}^{N_{\eta }}\frac{f(\mathbf{J}({\mathbf{x}}_k,{\mathbf{v}}_k))}{f_{\text{S}}({\mathbf{x}}_k,{\mathbf{v}}_k)}S(\mathbf{u}),$$

$$\xi (\mathbf{u}|{\bar{\mathbf{u}}}^{\eta })=C_{\eta }G(\mathbf{u},{\bar{\mathbf{u}}}^{\eta },{\mathbf{\sigma }}^{\eta })f_{\text{S}}(\mathbf{x},\mathbf{v})\left|\frac{\partial (\mathbf{\theta },\mathbf{J})}{\partial (\mathbf{u})}\right|,$$

$$L_{\text{d}}=L_{\text{h}}=100\mathrm{kpc}\mathrm{km}{\mathrm{s}}^{-1},k_{\text{d}}=1,\gamma =0,$$

Full text

††pubyear: 2017††pagerange: Modelling the Milky Way’s globular cluster system–Modelling the Milky Way’s globular cluster system

1 Introduction

In recent years, it has become possible to construct sophisticated models for the Milky Way. These models assume a gravitational potential that is axisymmetric, and components of the Galaxy are approximated by distribution functions (DFs) chosen to be analytic in three isolating integrals of motion. This allows for the construction of equilibrium DFs via Jeans’ (1915) theorem. Candidate DFs for the stellar discs were proposed by Binney (2010) and refined and extended by Binney & McMillan (2011) and Sanders & Binney (2015). A DF for the dark halo was implemented by Piffl, Penoyre & Binney (2015), while Posti et al. (2015) and Williams & Evans (2015) describe a wide range of DFs for spheroidal systems, some of which have been applied and extended by Das & Binney (2016) and Das et al. (2016). In this paper, we explore how this technique of constructing action-based DFs can be used to understand the Galactic globular cluster (GC) system.

Harris & Canterna (1979) noted that the distribution in metallicity of the Galactic globular clusters is bimodal. Zinn (1985) showed that the metal-rich and metal-poor sub-populations have distinct phase-space distributions, the metal-rich clusters being more strongly concentrated to the Galactic centre and forming a more rapidly rotating body. Now one divides the 157 GCs in the Harris 1996 catalogue (2010 edition)111http://www.physics.mcmaster.ca/Globular.html into 44 metal-rich GCs with $\hbox{[Fe/H]}>-0.8$ and 113 metal-poor GCs with $\hbox{[Fe/H]}<-0.8$. Zinn (1993) later suggested that the metal-poor subpopulation could be further divided into ‘young’ and ‘old’ based on the relation between their metallicities and the morphologies of their horizontal branches (see e.g. Fig. 6 of Mackey & Gilmore, 2004). Moreover, he argued that the old metal-poor GCs are more concentrated closer to the Galactic centre than the young GCs. However, very precise ages of 55 GCs extracted from HST data do not confirm a relationship between age and Galactocentric radius (VandenBerg et al., 2013). Mackey & Gilmore (2004) found that the metal-rich GCs in their sample were all within a Galactocentric radius $r\sim 6\,\mathrm{kpc}$, with 60 per cent of them situated at $r<3.0\,\mathrm{kpc}$. Unlike the metal-poor GCs, this system is very flattened, with all metal-rich GCs at $|z|<3.6\,\mathrm{kpc}$, and all but four located at $|z|<2.0\,\mathrm{kpc}$ ($z$ is the vertical distance from the plane of the Galaxy).

Studying the velocity distribution of the GC system proves to be more difficult, as 53 per cent of GCs do not have any proper motion data. Frenk & White (1980) provide a method to estimate the mean rotation velocity $v_{\textup{rot}}$ of the GC system based only on position and line-of-sight velocities. Using this method, and assuming a circular speed of $220\,\mathrm{km}\,\mathrm{s}^{-1}$ at the solar neighbourhood, Thomas (1989) found that the GC system has a mean rotation of $v_{\textup{rot}}\sim 65\,\mathrm{km}\,\mathrm{s}^{-1}$. The metal-rich GCs in his sample systemically rotate with $v_{\textup{rot}}\sim 113\,\mathrm{km}\,\mathrm{s}^{-1}$, while the metal-poor GCs have $v_{\textup{rot}}\sim 43\,\mathrm{km}\,\mathrm{s}^{-1}$. Zinn (1985) expands on this result, finding that the old metal-poor GCs rotate with $v_{\textup{rot}}\sim 70\,\mathrm{km}\,\mathrm{s}^{-1}$, while the younger metal-poor GCs have a mean rotation consistent with zero, although the uncertainty in the latter result is large.

These considerations are all consistent with Zinn’s (1985) identification that the metal-rich GCs exhibit disc-like kinematics, while the metal-poor GCs exhibit halo-like kinematics. As the GC system is thought to consist of these two distinct subpopulations, our DF will also be constructed using distinct disc and halo components. However, in fitting the DF to the data we make no assumption regarding the metallicities of the components, and in this way investigate whether the GC system can be effectively divided using phase-space data alone.

Our DF is described in Section 2. Section 3 explains how the posterior distribution for our model parameters was determined. Section 4 describes models favoured by the data, including the spatial and kinematic properties of the two components. Section 5 sums up and suggests some directions for further work.

2 Model

2.1 The Galactic potential

Piffl et al. (2014) sought an axisymmetric Galactic potential that is consistent with a wide range of observational data. Specifically, they required the potential to reproduce gas terminal velocities at various longitudes, the kinematics of stellar masers with very precise astrometry, the proper motion of Sgr A*, the run of stellar density with distance from the Galactic plane near the Sun, and the kinematics of $\sim 200\,000$ stars in the RAdial Velocity Experiment (RAVE; Steinmetz et al., 2006; Kordopatis et al., 2013). The potential assumes that the Galaxy’s mass is dominated by a gas disc, thin and thick stellar discs, a stellar bulge and a dark halo. The disc components are described by density distributions of the form

[TABLE]

where $(R,z)$ are cylindrical coordinates in the Galactocentric frame, $R_{\mathrm{d}}$ and $z_{\mathrm{d}}$ are characteristic scales, and $\Sigma_{0}$ is the surface density of the disc. A non-zero value of $R_{\textup{hole}}$ allows for a central cavity.

The stellar bulge and dark halo have

[TABLE]

where

[TABLE]

Here, $\rho_{0}$ is a normalisation constant, $r_{0}$ is a scale radius, and $q$ is the axis ratio for surfaces of constant density. The exponential term permits the spheroids to extend only to a finite distance set by $r_{\textup{cut}}$, and also ensures that the density profile has a finite mass for all sensible values of $a$ and $b$. The Galactic potential $\Phi(R,z)$ is then given by solving the Newton–Poisson equation.

As the constraints adopted by Piffl et al. (2014) are probes of the vertical profile of gravitating matter, both baryonic and dark, their data are consistent with a model with a spherical dark halo and a heavier baryonic disc, or a more flattened dark halo and a lighter baryonic disc. Consequently, the axis ratio $q\leq 1$ of the dark halo is not tightly constrained, although Piffl et al. (2014) argue that a comparison with the results of Bienaymé et al. (2014) favours an axis ratio $q\simeq 0.8$. We have chosen to use $q=0.8$ for the dark halo, although a larger value has no significant effect on our results. The full set of parameter values we used for the gravitational potential are shown in Table 1. Section 2.3 of Binney & McMillan (2016) gives details of downloadable code that evaluates $\Phi(R,z)$ and its derivatives given the numbers in Table 1.

Note that the GC system does not enter into this model because its mass is negligible. We treat the GCs as a system of 157 identical, non-interacting point particles moving in the static potential $\Phi(R,z)$, each with Hamiltonian

[TABLE]

The velocity components $v_{R}$ and $v_{z}$ are in the radial and vertical directions respectively, and $J_{\phi}\equiv Rv_{\phi}$ is the conserved angular momentum about the axis of Galactic rotation. Positive velocity $v_{\phi}$ is in the direction of Galactic rotation.

2.2 Angle–action variables

We define the set of observables for a GC as

[TABLE]

where $(l,b)$ are the Galactic longitude and latitude, $s$ is the heliocentric distance, $v_{\mathrm{los}}$ is the line-of-sight velocity, and $\bm{\mu}\equiv(\mu_{\alpha}^{*}=\dot{l}\cos b,\mu_{\delta}=\dot{b})$ is the proper motion vector.

In order to transform these into the phase space coordinates

[TABLE]

we have assumed that the Sun is located at $(R,\phi,z)=(R_{0},0,0)$, where we take $R_{0}=8.3\,\mathrm{kpc}$ (Schönrich, 2012). We further assume that the local circular speed is $240\,\mathrm{km}\,\mathrm{s}^{-1}$ (Schönrich, 2012), and that the Sun has velocity $\bm{(}U,V,W)_{\odot}=(11.1,12.24,7.25)\,\mathrm{km}\,\mathrm{s}^{-1}$ (Schönrich, Binney & Dehnen, 2010) relative to the local standard of rest (LSR). Positive velocities $U$, $V$ and $W$ point in the direction of the Galactic centre, Galactic rotation and Galactic north pole respectively.

It proves useful to transform these phase space coordinates into angle–action variables $(\bm{\theta},\bm{J})$. These are a set of canonical coordinates where the momenta $J_{i}$ are integrals of motion. It then follows from Hamilton’s equations that $H\equiv H(\bm{J})$ and the angles $\theta_{i}$ increase linearly with time (Binney & Tremaine, 2008). The orbit of a GC is specified by $\bm{J}$, and if at time $t=0$ it is at position $\bm{x}(\bm{\theta}(0),\bm{J})$, then we simply increase the angles linearly to evolve the position forwards in time. Despite its advantages, this formalism has until recently been little used on account of the difficulty in evaluating $(\bm{\theta},\bm{J})$. Recent technical progress now makes this possible. We make extensive use of the Stäckel Fudge (Binney, 2012; Sanders & Binney, 2016), which yields $(\bm{\theta},\bm{J})$ given $(\bm{x},\bm{v})$ in an axisymmetric potential like that of Section 2.1. We identify $\bm{J}\equiv(J_{r},J_{z},J_{\phi})$, where $J_{r},J_{z}\geq 0$ can be thought of as quantifying oscillations in the radial and vertical directions.

2.3 Distribution function

Given that the GC system is in dynamical equilibrium, Jeans’ (1915) theorem allows us to assume that the DF is a function $f(\bm{J})$, so the probability that a randomly chosen GC has phase-space coordinates in $\mathrm{d}^{3}\bm{\theta}\,\mathrm{d}^{3}\bm{J}$ is

[TABLE]

We normalize $f(\bm{J})$ such that

[TABLE]

so the quantity $(2\pi)^{3}f(\bm{J})\;\mathrm{d}^{3}\bm{J}$ gives the probability that a randomly selected GC moves on the orbit specified by $\bm{J}$.

Although we do not assume that the metal-rich clusters form a disc-like component and the metal-poor clusters form a spheroid, our DF is a linear combination of a DF for a disc-like population and a DF for a halo population. We write

[TABLE]

where $F_{\textup{disc}}\in[0,1]$ is the fraction of disc GCs and $\Pi$ is a set of parameters for the model. We will often write $f(\bm{J}|\Pi)$ as simply $f(\bm{J})$ in what follows, but they will mean the same thing. It turns out that favoured models assign roughly as much probability to the disc-like component as the fraction of observed clusters that are metal rich, but this is an empirical result rather than an assumption.

The DF of an axisymmetric system is usually best considered to be the sum of a part $f_{+}$ even in $J_{\phi}$ and a part $f_{-}$ odd in $J_{\phi}$. The latter does not contribute to the density of the system, but is largely responsible for the system’s rotation. If we wish to avoid discontinuities in the DF, $f_{-}$ must vanish with $J_{\phi}$. A convenient way to satisfy this condition and obtain a non-negative and physically reasonable DF is to posit

[TABLE]

Here $-1\leq k\leq 1$ is a constant that controls the sign and intensity of any rotation, and $L$ is a constant that controls the steepness of the system’s central rotation curve. We have used this ansatz for the odd parts of the DFs of both disc and halo.

2.3.1 Disc

The disc component is described by the ‘quasi-isothermal’ DF introduced by Binney & McMillan (2011). We take the part even in $J_{\phi}$ to be

[TABLE]

where

[TABLE]

$R_{\textup{c}}(J_{\phi})$ is the radius of the circular orbit with angular momentum $J_{\phi}$, $\Sigma(J_{\phi})$ is approximately the surface density of the disc at that radius, $\mathcal{R}_{\textup{d}}$ is a characteristic scale length, and $\Sigma_{0}$ is a normalisation constant that ensures that $(2\pi)^{3}\int{\rm d}^{3}{\bm{J}}\,f_{\rm disc}=1$. The functions $\Omega(J_{\phi})$, $\kappa(J_{\phi})$ and $\nu(J_{\phi})$ are the circular, radial and vertical epicycle frequencies of the potential evaluated at the radius $R_{\textup{c}}(J_{\phi})$. Following Binney & McMillan (2011) we let the velocity-dispersions parameters $\sigma_{r}$ and $\sigma_{z}$ vary with $J_{\phi}$ as

[TABLE]

where $\sigma_{i0}$ and $\gamma$ are free parameters.

2.3.2 Halo

Our halo DF is inspired by Posti et al. (2015). Das & Binney (2016) and Das et al. (2016) recently extended DFs of this type to depend on [Fe/H] and age in addition to $\bm{J}$. However, given the small number of halo clusters, we doubt the ability of the data to constrain dependence on [Fe/H]. So we model the halo clusters with a metallicity-blind and age-blind DF.

The DFs of Posti et al. (2015) are constructed from a function $h({\bm{J}})$ that is a homogeneous function of degree one in the $J_{i}$. The resulting stellar system has a plausible structure near the $z$ axis only when the coefficients in $h$ of $J_{\phi}$ and $J_{z}$ become equal as $J_{\phi}\to 0$. We ensure satisfaction of this condition in the simplest possible way, namely by making $h$ a function of $|J_{\phi}|+J_{z}$. In fact we define $h$ to be

[TABLE]

where the small constant $J_{\rm core}=10^{-3}\,\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$ ensures that the DF remains finite at $\bm{J}=0$ at the cost of making $h$ not quite a homogeneous function of the $J_{i}$. The epicycle frequencies are evaluated at $R_{\textup{c}}(J_{\textup{tot}})$, where $J_{\textup{tot}}=J_{r}+J_{z}+|J_{\phi}|$. This choice is made to prevent the epicycle frequencies becoming large for eccentric or highly inclined orbits (when $|J_{\phi}|$ is small) (Piffl et al., 2015).

For the halo clusters, the part of the DF of that is even in $J_{\phi}$ is

[TABLE]

with $N_{0}$ a normalisation constant chosen so $(2\pi)^{3}\int{\rm d}^{3}{\bm{J}}\,f_{\textup{halo}}=1$. The exponents $\alpha$ and $\beta$ in equation (15) control the inner and outer power-law slopes of the DF, and the two regimes are separated by a break scale $J_{0}$. Posti et al. (2015) show that DFs with $(\alpha,\beta)=(2,5)$ self-consistently produce models with radial density profiles that closely approximate that of Jaffe (1983), while DFs with $(\alpha,\beta)=(5/3,5)$ and $(5/3,3)$ self-consistently produce models with radial density profiles very similar to those of Hernquist (1990) and Navarro, Frenk & White (1996) models, respectively. Unfortunately, since the GC system does not generate the potential that confines it, the relations given by Posti et al. (2015) between $\alpha,\beta$ and the slopes of the density profile in real space do not apply. Lastly, the exponential term with $J_{\textup{max}}=10^{6}\,\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$ ensures that the DF has a finite mass for all choices of $\alpha$ and $\beta$; this has little effect on the DF provided we choose $J_{\textup{max}}$ to be sufficiently large.

3 Bayesian inference

To determine the values and associated uncertainties of the model parameters $\Pi$ that best fit the GC data, we compute the posterior distribution

[TABLE]

where $\textup{Pr}(\Pi)$ is an appropriately chosen prior and $\textup{Pr}(\textup{Data}|\Pi)\equiv\mathcal{L}$ is the data’s likelihood. We follow the approach to the application of Bayesian inference to Galactic structure described in McMillan & Binney (2013) (hereafter MB13).

3.1 Prior distribution

As we are interested in how the GC data constrain the parameters of this model, we opt for the least informative prior possible. The variables $k_{\textup{d}},k_{\textup{h}}\in[-1,1]$ and $F_{\textup{disc}}\in[0,1]$ have finite ranges, so we take the prior to be uniform in these ranges. The scale parameters $L_{\textup{d}}$, $L_{\textup{h}}$, $J_{0}$, $\sigma_{i0}$, and $\mathcal{R}_{\textup{d}}$ are intrinsically positive, so the unbiased prior is uniform in the logarithms of these quantities (Jeffreys, 1961). The indices $\alpha,\beta$ can, in principle take arbitrarily large positive values, but if we require only $\beta>0$, the data yield a degeneracy in which both $\beta$ and $J_{0}$ increase. From the work of Posti et al. (2015) and Das & Binney (2016) we expect $\beta\sim 5$. To incorporate this knowledge into our prior, we take the prior on $\beta$ to be the gamma distribution

[TABLE]

where $A$ and $B$ are hyperparameters and $\Gamma(A)$ is the gamma function. We set $A=4$ and $B=1$ such that the distribution has a mean $\beta=4$, a mode $\beta=3$, and a variance of $4$. With this prior on $\beta$, implausibly large values of $\beta$ and $J_{0}$ acquire low values of the posteriori probability, but it might be argued that a weaker prior would be preferable. We have investigated two priors on the inverse distance $\gamma$: uniform and uniform in its logarithm.

In summary, our prior is

[TABLE]

where, as discussed in Section 4 $\textup{Pr}(\gamma)$ is successively (i) uniform in $\gamma$, (ii) uniform in $\log\gamma$, and (iii) $\delta(\gamma)$.

3.2 Data

Data on the position variables $(l,b,s)$ of the 157 known GCs were obtained from the 2010 edition of the Harris 1996 catalogue. This catalogue also contains line-of-sight velocities for 143 of these GCs. Proper motion data exist for only 64 GCs; these have been obtained from Massari et al. (2013) and a catalogue by Casetti–Dinescu et al. (2013)222The proper motions catalogue can be accessed at http://www.astro.yale.edu/dana/gc.html..

The Galactic coordinates $(l,b)$ are measured to high precision, so the effects of their errors will be neglected. The errors in the remaining observable quantities are assumed to be independent. If $u\in\{v_{\mathrm{los}},\mu^{*}_{\alpha},\mu_{\delta}\}$ denotes a component of velocity, then we assume that its measured value $\bar{u}$ is related to its true value $u$ by the Gaussian distribution

[TABLE]

where $\sigma_{u}$ is the uncertainty in $u$. Where possible, the values of $\bar{u}$ and $\sigma_{u}$ are obtained from the data sources listed above. When no such data are available, following MB13 the value of $\sigma_{u}$ is taken to be sufficiently large that the Gaussian distribution is effectively constant for all sensible values of $u$.

Since heliocentric distance $s$ is an intrinsically positive quantity, we take its probability distribution to be a lognormal distribution

[TABLE]

where following (Casetti–Dinescu et al., 2013) the values of $\sigma_{s}$ and $\mathcal{S}(\bar{s},\sigma_{s})$ are chosen such that the distribution has a mean equal to the measured distance $\bar{s}$ and a variance of $(0.1\bar{s})^{2}$. Actually the uncertainties in cluster distances are surely sometimes smaller and sometimes larger than this estimate, so this work could be upgraded by obtaining errors on individual cluster distances. We shall see also that it would be very desirable to take pains to ensure that all cluster distances are assigned using the same value of $R_{0}$ as that adopted for the model.

For brevity, we use the notation

[TABLE]

to describe the six-dimensional error distribution for the observables of the $\eta^{\textup{th}}$ GC.

3.3 Likelihood

While MB13 include the apparent magnitude $m$ as an observable, we do not because we assume the completeness of our sample is independent of magnitude. The likelihood $\mathcal{L}$ that the GC data are drawn from the model $f(\bm{J})$ is the product

[TABLE]

where $\mathrm{Pr}(\bm{u}|\textup{Model})$ is the probability that a randomly chosen GC has true observables $\bm{u}$:

[TABLE]

with

[TABLE]

Here the coordinates $(\bm{\theta},\bm{J})$ are evaluated at the given point $\bm{u}$, and the selection function $S(\bm{u})$ gives the probability that if a GC with true observables $\bm{u}$ exists, it has been included in the Harris and Casetti–Dinescu catalogues. The normalisation constant $A$ is given by the condition $\int{\rm d}^{6}{\bm{u}}\,\mathrm{Pr}(\bm{u}|\textup{Model})=1$, so

[TABLE]

3.3.1 Selection function

Given that globular clusters have been discovered over many decades and in a range of wavebands, it is impossible to characterise the incompleteness of our sample with any precision. It is, however, believed that nearly all our Galaxy’s GCs have been observed – it is estimated that the total number of Galactic GCs lies between 160 to 180 (Harris, 2001; Kurtev et al., 2008). If this is accepted, we do not need a sophisticated selection function. It is likely that that any GCs that have still eluded astronomers are likely to lie close to the Galactic plane, so they are hidden by dust. To take this effect into account, we adopt the selection function

[TABLE]

This simply states that a GC located at Galactic coordinates $(l,b)$ will be unobservable if the extinction in that direction is above a threshold value $\mathrm{E(B-V)}_{\textup{max}}$.

We can obtain estimates for $\mathrm{E(B-V)}$ at any Galactic coordinate using the dust map produced by Schlegel, Finkbeiner & Davis (1998)333Dust map values were obtained from http://irsa.ipac.caltech.edu/applications/DUST/.. We find that we must choose $\mathrm{E(B-V)}_{\textup{max}}\geq 4$ to ensure that all Harris catalogue GCs that were observed in the visible have $S(\bm{u})=1$. An important subtlety is that Schlegel et al. (1998) produce estimates for extinctions from the observer to infinity, rather than from the observer to a given heliocentric distance $s$. This is only an issue for GCs that might be located in the Galactic plane but have sufficiently small values of $s$ that they can nevertheless be detected. Such instances are probably rare we ignore them.

It turns out that a few GCs in the Harris catalogue have been discovered in the infrared, and while this selection function does not include that possibility, the number of such GCs is so small that the effect is insignificant. Further, our results are not noticeably affected by any choice of threshold value greater than $4$. In fact, excluding a selection function all together does not dramatically alter our results. Below we keep $\mathrm{E(B-V)}_{\textup{max}}=4$.

Since we require values for $E(B-V)$ over the whole sky and it is essential to use a consistent set of values, we do not use the values for individual clusters given by Harris.

3.3.2 Evaluating the likelihood

We now use Monte Carlo methods to approximate the integrals in equations (22) and (25). We introduce a sampling density $f_{\mathrm{S}}(\bm{x},\bm{v})$ that approximates the phase space distribution of the GC data. This ensures evaluations are concentrated where the integrand is largest. The normalisation constant in equation (25) can now be evaluated as

[TABLE]

where we draw $N_{\textup{S}}$ points $(\bm{x}_{k},\bm{v}_{k})$ from the sampling density $f_{\textup{S}}$. The integral in equation (22) for the $\eta^{\textup{th}}$ star becomes

[TABLE]

where we draw $N_{\eta}$ points $(\bm{x}_{k},\bm{v}_{k})$ from the sampling density

[TABLE]

The normalisation constant $C_{\eta}$ depends only on the data and not on the model $f(\bm{J})$, hence need not be computed in what follows.

3.3.3 Choice of sampling density

A good choice for the sampling density $f_{\textup{S}}(\bm{x},\bm{v})$ is one that approximates a typical model $f(\bm{J}(\bm{x},\bm{v}))$. We have chosen $f_{\textup{S}}$ to be a product of the density profile $\rho_{\textup{halo}}(R,z)$ given in equation (2) with $1/r_{\textup{cut}}=0$ and a triaxial Gaussian velocity distribution with principal axes aligned along the $v_{R}$, $v_{z}$ and $v_{\phi}$ directions.

A maximum likelihood fit to the GC data yields $a=0$, $b=4.43$, $r_{0}=2.49\,\mathrm{kpc}$ and $q=0.83$ for the density profile. We have not attempted to attach any confidence intervals to these numbers since all we require is a crude first guess at the phase space distribution. For this same reason, we have chosen to approximate the density profile using only $\rho_{\textup{halo}}(R,z)$, rather than with an appropriate combination of $\rho_{\textup{halo}}(R,z)$ and $\rho_{\textup{disc}}(R,z)$.

Choosing parameters for the velocity distribution of $f_{\textup{S}}$ is a more nuanced task because only 64 of the 157 GCs have complete velocity data. Naturally we take the means of $v_{R}$ and $v_{z}$ to be zero, while we set $\langle v_{\phi}\rangle=v_{\textup{rot}}$, where $v_{\textup{rot}}=79\,\mathrm{km}\,\mathrm{s}^{-1}$ is the rotation velocity returned by the algorithm of Frenk & White (1980) when applied to our data set with the circular speed at the Sun set to $v_{\textup{c}}=240\,\mathrm{km}\,\mathrm{s}^{-1}$. When the Frenk-White algorithm is applied to the metal-rich subpopulation, we find $v_{\textup{rot}}=137\,\mathrm{km}\,\mathrm{s}^{-1}$, while the metal-poor population yields $v_{\textup{rot}}=48\,\mathrm{km}\,\mathrm{s}^{-1}$. These values will be useful later. Note that these results are larger than those obtained by Thomas (1989) (see Section 1). This is because we have a larger data set and have also used a larger value for the local circular speed.

As for the velocity dispersions of $f_{\textup{S}}$, we have chosen $(\sigma_{R},\sigma_{\phi},\sigma_{z})=(140,140,100)\,\mathrm{km}\,\mathrm{s}^{-1}$, based on the dispersions calculated from the sample of 64 GCs with proper motion data. It is not necessary that these dispersion values describe the GC system as a whole. To ensure that we have made a reasonable choice, we have repeated our analysis choosing different values of the dispersions in the range $100$–$250\,\mathrm{km}\,\mathrm{s}^{-1}$. We find our results are not significantly affected for dispersions within this range.

3.4 Posterior distribution

We MCMC sample the posterior distribution as follows:

Sample $N_{\textup{S}}$ points from $f_{\textup{S}}(\bm{x},\bm{v})$ and $N_{\eta}$ points for each GC from $\xi(\bm{u}|\bm{u}^{\eta})$. These are independent of $f(\bm{J})$ so need only be sampled once at the beginning. 2. 2.

Use the Schlegel et al. (1998) map to determine the selection function $S(\bm{u})$ for each of the above $N_{\textup{S}}$ points. 3. 3.

Use the Stäckel fudge (Binney, 2012, 2014) with the gravitational potential $\Phi$ in Section 2.1 to calculate $\bm{J}$ for each of these points. 4. 4.

Pick a point $\Pi$ in the space of model parameters at random. 5. 5.

Calculate the prior $\textup{Pr}(\Pi)$ using equations (17) and (18). 6. 6.

If $\textup{Pr}(\Pi)\neq 0$, calculate the likelihood $\mathcal{L}$ for the model $f(\bm{J}|\Pi)$ using equations (27) and (28). 7. 7.

Calculate the posterior $\textup{Pr}(\Pi|\textup{Data})$ using equation (16). 8. 8.

We use a robust adaptive Metropolis algorithm (Vihola, 2012) to decide on the next point $\Pi$ in the Markov chain. 9. 9.

Repeat steps 5 to 8 until the desired number of MCMC points have been sampled.

We have sampled $50\,000$ MCMC points, and we have used $N_{\eta}=100$ and $N_{\textup{S}}=157\times N_{\eta}$, the latter choice made such that an equal number of points are evaluated for the numerator and denominator of $\mathcal{L}$.

3.5 Pseudo-catalogues

We can probe the impact of noise, and the extent to which a model can match the data by using a model to generate pseudo-catalogues of clusters. We generate a pseudo-catalogue as follows:

Sample a phase-space point $(\bm{x},\bm{v})$ for a ‘pseudo-GC’ from $f(\bm{J}(\bm{x},\bm{v})|\Pi)$ and compute the observables $\bm{u}$. 2. 2.

Accept the pseudo-GC with probability given by the selection function $S(\bm{u})$ and return to the previous step until 157 pseudo-GCs have been accepted. 3. 3.

The observables of each pseudo-GC are ascribed errors $\bm{\sigma}$ according to the error model described below. 4. 4.

The observables $\bm{u}$ are scattered by their errors.

The line-of-sight velocities $v_{\mathrm{los}}$ have errors that approximately follow a lognormal distribution (see Fig. 2). A priori, we might expect that these errors depend on the heliocentric distance $s$, but by plotting errors in $v_{\mathrm{los}}$ versus distance $s$ in the Harris catalogue, the lower plot of Fig. 2 demonstrates that errors in $v_{\mathrm{los}}$ do not increase significantly with $s$. Consequently, each GC in a pseudo-catalogue is given an error randomly drawn from the same lognormal distribution.

As more than half of the GCs do not have proper motion data, quantities calculated using the proper motions will come with very large uncertainties, and we do not consider these useful to decide if a model agrees with data. We therefore have not concerned ourselves with modelling the error distribution of the proper motions, and will restrict ourselves to analysing pseudo-catalogues using only the observables $(l,b,s,v_{\mathrm{los}})$.

4 Results

4.1 Parameters to fix

Equation (10) applied to the halo and disc DFs defines four parameters $k_{\textup{h}}$, $L_{\textup{h}}$, $k_{\textup{d}}$ and $L_{\textup{d}}$ that control the amplitude and central steepness of each system’s rotation curve. Along early MCMC chains the posterior distributions of $L_{\textup{h}}$ and $L_{\textup{d}}$ were indistinguishable from their priors (uniform in $\log L_{i}$). Moreover, no correlation was apparent between $L_{i}$ and any other parameter of the DF. Evidently, the data do not usefully constrain the steepness of the central rotation curves of the components, so in the following we simply adopt $L_{i}=100\,\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$, which we shall find causes the rotation curves of the disc and halo populations to have central slopes $\sim 100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{kpc}^{-1}$.

Fig. 3 shows the posterior distribution of $k_{\textup{d}}$, which is crowded near the upper limit of its permitted range. The crowding of the posterior distribution of $k_{\textup{d}}$ near unity implies a lack of evidence for a counter-rotating disc of clusters. Consequently, in the following we set $k_{\textup{d}}=1$.

In early MCMC exploration of parameter space we adopted a uniform prior on the parameter $\gamma$ that controls the radial variation of the disc’s velocity dispersions (eqn. 13). The posterior distribution of $\gamma$ then extended both sides of zero. Since we consider increases in velocity dispersion with radius implausible, we then ran chains with the prior on $\gamma$ taken to be uniform in $\log\gamma$. The resulting posterior distribution of $\gamma$ was essentially uniform in $\log\gamma$ for $\gamma<10^{-6}\,\mathrm{kpc}^{-1}$. That is, the data only require that $\gamma$ is so small as to have negligible impact in the region $r<100\,\mathrm{kpc}$ to which the data are confined. This being so, we subsequently set $\gamma=0$, i.e., we made the velocity-dispersion parameters independent of $J_{\phi}$.

In summary, we fix the values of four parameters:

[TABLE]

leaving the posterior distributions of eight parameters to be explored by MCMC chains. Of these four describe the halo DF ($\alpha$, $\beta$, $k_{\textup{h}}$ and $J_{0}$), three describe the disc DF ($\mathcal{R}_{\textup{d}}$, $\sigma_{r}$ and $\sigma_{z}$), and the eighth and final adjustable parameter is $F_{\textup{disc}}$, which is the fraction of the probability associated with the disc.

4.2 Favoured models

Fig. 4 illustrates the structure of the posterior distributions of all adjustable parameters by showing on the diagonal histograms for each parameter after marginalising over all other parameters, and in the off-diagonal panels the correlations between each pair of model parameters, again after marginalising over all other parameters. The white contours in the off-diagonal panels enclose 68 per cent of the probability.

Parameters that are quite well determined are the extent to which the halo cluster distribution rotates $k_{\rm h}$, the fraction of disc clusters $F_{\textup{disc}}$, the velocity-dispersion parameters of the disc clusters $\sigma_{r}$ and $\sigma_{z}$, and the scale length of the system of disc clusters $\mathcal{R}_{\textup{d}}$. For each parameter, Table 2 gives the expectation value and the standard deviation along the MCMC chain.

Three parameters for the halo clusters are strongly correlated: the inner and outer slope parameters $\alpha,\beta$ and the scale action $J_{0}$. Large values of $J_{0}$ are associated with large values of $\alpha$ and $\beta$. This correlation makes perfect sense physically when one recalls that $\alpha$ and $\beta$ control the inner and outer slopes of the real-space density profile of the population of halo clusters, while $J_{0}$ controls the break radius that divides the two regimes. In so far as the density profile of the halo cluster distribution steepens smoothly with increasing radius, such a correlation between $\alpha$, $\beta$ and $J_{0}$ and is inevitable. Naively, one expects $J_{0}$ to be roughly the product of the circular speed and the radius $\sim 2.5\,\mathrm{kpc}$ of the break in the density profile, so $\sim 600\,\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$. The probable values of $J_{0}$ are significantly larger than this. However, we show below that the recovered parameter values do reproduce the expected break radius.

4.3 What do we learn from the MCMC sample?

Fig. 5 shows the distribution of the likelihoods of the data given the models sampled by MCMC. This distribution – essentially the $\chi^{2}$ distribution of the models – is more than two orders of magnitude wide, so significant probability is associated with models that make the data hundreds of times less probable than does the most probable model. This state of affairs is commonplace when models with significant numbers of parameters are employed in Bayesian inference.

Why must we consider models that make the data so much less probable than the maximum-likelihood model? Because the maximum-likelihood model achieves its high likelihood in large measure by fitting not only the signal but also the the noise in the data. In Fig. 6 we demonstrate this by comparing the observables predicted by samples of 157 clusters drawn from (a) the maximum-likelihood model (top row) with (b) models drawn at random from the MCMC chain (lower row). We see that the scatter around the real observables (red curves) of the observables predicted by the maximum-likelihood model and the models drawn at random are indeed similar. This result confirms that the excess likelihood of the maximum likelihood model over typical models in the MCMC chain indeed reflects its ability to fit the high level of noise inherent in there being only 157 clusters.

In Fig. 6 the only panels in which the red line of the data lies outside the region explored by the 20 realisations are those for the distance, $s$: in these panels the full red line for the data rises steeply with increasing $s$ at values of $s$ that are $\sim 10$ per cent smaller than the broken curves of the pseudo-data. The number of clusters increases rapidly as $s$ becomes comparable to $R_{0}$ and the clusters gathered around the Galactic Centre enter the sample.

The red histogram in Fig. 7 shows the distribution of the Galactocentric azimuths of the real clusters: there is a clear excess of clusters with $|\phi|<20$ degrees, indicating that they lie in front of the Galactic Centre. By contrast, the black dashed histogram, which shows the azimuthal distribution of the clusters in twenty models from the MCMC chain, is very uniform, so we have not succeeded in modelling the distribution of clusters in the central region. The cause of this discrepancy could be either (a) that distances to clusters are systematically too short, or (b) that our value $R_{0}=8.3\,\mathrm{kpc}$ is too long, or (c) that observers have failed to identify significant numbers of clusters that lie in or behind the bulge. The latter could be due either to more extinction than we have adopted, or confusion in crowded fields. However, following the Vista Variables in the Via Lactea (VVV) survey (Minniti et al., 2010) of the bulge region in the near IR and with excellent spatial resolution, it seems unlikely that item (c) is a significant problem.

Apart from this question surrounding the distribution of cluster in the bulge region, we consider that the plots of Fig. 6 are consistent with the actual clusters being drawn from a model in the MCMC chain.

4.4 The expected distribution of clusters

An MCMC chain encodes the probability density of each part of model space, and the DF quantifies the probability density of GCs in phase space for a given model. Consequently, if we average $f({\bm{x}},{\bm{v}})$ over models in an MCMC chain, we obtain our best estimate of the probability of finding a GC at $({\bm{x}},{\bm{v}})$. We now present plots obtained by averaging $f$ over 50 models drawn from an MCMC chain.

Fig. 8 shows the probability density of disc and halo clusters in action space. Specifically, the logarithm to base 10 of the density is shown on three slices at constant radial action, from $J_{r}=0$ (circular orbits) at the top to larger values of $J_{r}$ lower down. The two types of cluster are seen to have very different phase-space distributions. The disc clusters are most dense along the $J_{\phi}$ axis and their density declines steeply with increasing $J_{r}$ or $J_{z}$, whereas the halo clusters are dense only at the origin of action space but their density declines relatively slowly with increasing $|{\bf J}|$.

The model provides probabilities for a cluster to be a disc rather than a halo cluster based solely on the cluster’s phase-space position. It is natural to ask how these probabilities relate to the cluster’s metallicities: the latter have a clearly bimodal distribution with a minimum at $\hbox{[Fe/H]}\simeq-0.8$ (Zinn, 1985; Harris et al., 2016) dividing the metal-poor halo cluster from the metal-rich disc clusters.

Fig. 9 shows for each cluster in the Harris catalogue the ratio $\langle P_{\rm disc}\rangle/\langle P_{\rm halo}\rangle$ versus [Fe/H], where, for example $P_{\rm disc}$ is the integral in one model of the disc DF over the cluster’s error ellipsoid (eqn 22). For 73 clusters this ratio exceeds unity, so they are deemed more likely to be disc than halo clusters. Clusters with $\hbox{[Fe/H]}\ga-1.3$ are never considered much more likely to be halo clusters than disc cluster, while only a few clusters with $\hbox{[Fe/H]}<-2$ are much more likely to be disc than halo clusters. The only clusters that will be considered much more likely to be halo than disc clusters are those that counter-rotate or lie at very large radii. Inevitably, at least as many halo clusters at moderate $r$ co- rather than counter-rotate, and these clusters will have significant probabilities to be disc clusters. Hence it is to be expected that more clusters have high probabilities to belong to the disc than have $\hbox{[Fe/H]}>-0.8$. Hence Fig. 9 is fully consistent with the conjecture that all metal-rich clusters belong to the disc population, while the metal-poor clusters all belong to the halo population.

Fig. 10 shows the real-space density of disc clusters (left) and halo clusters (right) estimated from 50 models drawn at random from the MCMC chain. The colour scale shows the logarithm of the density in a slice that includes the $z$ axis. For semi-major axes of length $a\la 5\,\mathrm{kpc}$, the system of disc clusters has the isodensity contours of an oblate body with axis ratio $\sim 1:3$. Larger isodensity surfaces have a deep depression around the minor axis. Any system in which all particles rotate in the same sense about the $z$ axis can be expected to have a low density on that axis. The right panel of Fig. 10 indicates that the system of halo clusters forms a simpler oblate structure of axis ratio $q\sim 0.7$ out to semi-major axes $a\sim 20\,\mathrm{kpc}$, but a similar plot on a bigger scale shows that around $r\sim 30\,\mathrm{kpc}$ the system becomes quite spherical because the potential is then dominated by the spherical dark halo rather than the disc.

Fig. 11 shows the density of clusters as a function of radius in the equatorial plane (left panel) and along the $z$ axis (right panel). In each panel we show the contributions of the disc and halo clusters (blue and red lines, respectively) and the total density in black. Straight black lines at the lower right of each panel have slope $-5$, revealing that at the largest radii the cluster density is declining somewhat faster than as $r^{-5}$. From results in Posti et al. (2015) it is easy to show that if the rotation curve were flat, $\rho\sim r^{-5}$ would imply $f\sim|{\bm{J}}|^{-5}$. Table 2 indicates a markedly steeper decline in the DF, $f\sim|{\bm{J}}|^{-9}$, because the MCMC chain favours large values of the scale action $J_{0}$, so the asymptotic regime has not been reached even at $r\sim 80\,\mathrm{kpc}$.

The profiles of the system of halo clusters show wiggles at $R,z\ga 20\,\mathrm{kpc}$. The wiggle on the $z$ axis is most pronounced. This feature disappears when the disc’s potential is removed and the dark halo’s potential is enhanced in compensation, so it arises from the response of the cluster system to the fading of the disc’s gravitational field. It coincides with the decrease in flattening noted above in connection with Fig. 10.

Fig. 12 shows the kinematics computed from 50 models as a function of radius in the equatorial plane. As the centre is approached, the three broken red curves, which show $\sigma_{R}$, $\sigma_{\phi}$ and $\sigma_{z}$ in the halo system, approach one another at quite a large value: $\sim 170\,\mathrm{km}\,\mathrm{s}^{-1}$. While the centre of the halo system is isotropic, with increasing $R$ the curves for $\sigma_{\phi}$ and $\sigma_{z}$ fall much more steeply than the curve for $\sigma_{R}$, so the body of the system is quite radially biased: the conventional anisotropy parameter $\beta_{\rm s}$ reaches a peak value $0.68$ at $R=12\,\mathrm{kpc}$ and then gently falls to $\beta_{\rm s}=0.53$ at $R=30\,\mathrm{kpc}$. The full red curve shows that the mean rotation rate of the halo system is small – it peaks at $\langle v_{\phi}\rangle\sim 40\,\mathrm{km}\,\mathrm{s}^{-1}$ at $R\sim 1\,\mathrm{kpc}$ and from there falls to a minimum of $12\,\mathrm{km}\,\mathrm{s}^{-1}$ around $R=15\,\mathrm{kpc}$.

The disc system rotates fast: its rotation curve reaches $\langle v_{\phi}\rangle=185\,\mathrm{km}\,\mathrm{s}^{-1}$ at $R=5\,\mathrm{km}\,\mathrm{s}^{-1}$ and from there $\langle v_{\phi}\rangle$ falls almost linearly to $118\,\mathrm{km}\,\mathrm{s}^{-1}$ at $R=30\,\mathrm{kpc}$. The velocity dispersion tensor is quite anisotropic, with $\sigma_{z}\sim 100\,\mathrm{kpc}$ being largest and $\sigma_{\phi}$ falling from $\sim 50\,\mathrm{km}\,\mathrm{s}^{-1}$ at $R=3\,\mathrm{kpc}$ to $\sim 30\,\mathrm{km}\,\mathrm{s}^{-1}$ at $R=30\,\mathrm{kpc}$.

As we reported in Section 3.3.3, the algorithm of Frenk & White (1980) applied to the complete Harris catalogue yields a rotation velocity $v_{\textup{rot}}=79\,\mathrm{km}\,\mathrm{s}^{-1}$. When we assemble a catalogue of 157 clusters from each of 20 models from the MCMC chain and apply the Frenk & White algorithm to each catalogue, the recovered values of $v_{\textup{rot}}$ have mean $69\,\mathrm{km}\,\mathrm{s}^{-1}$ and standard deviation $29\,\mathrm{km}\,\mathrm{s}^{-1}$. Thus our models are consistent with the data from the perspective of rotation rate.

5 Conclusions

We have constructed the first fully dynamical model of the Galaxy’s GC system that includes a realistic Galactic potential, that determined by Piffl et al. (2014). The GC system is treated as a set of 157 identical, non-interacting point particles in dynamic equilibrium orbiting in this static potential. Motivated by the bimodal distribution of cluster metallicities, the DF consists of two components: a disc and halo. After fixing a number of parameters that were either barely constrained by the data or essentially fixed by the data, the final model had eight parameters to be fitted to the data: inner and outer slope parameters, a scale action and a rotation rate for the halo, a scale length, in-plane and vertical velocity dispersions for the disc, and the fraction of the probability provided by the disc.

Although 157 clusters prove too few to constrain tightly any of the parameters, the disc fraction, the three disc parameters and the halo’s rotation rate all produce well-defined peaks in the likelihood. The remaining three halo parameters have a natural degeneracy, in which an increase in the scale action (and corresponding scale radius) can be compensated by changes to the inner and outer slopes of the density profile.

Given the similarities between the metallicity distribution functions of the disc clusters and the thick disc (Wyse & Gilmore, 1995), the DF of the disc component invites comparison with the Extended Distribution Function (EDF) Sanders & Binney 2015, hereafter SB15 fitted to the stellar thick disc using data from the Geneva-Copenhagen survey (Holmberg et al., 2007, 2009; Casagrande et al., 2011). We found that the cluster data did not provide useful constraints on any radial variation in the velocity-dispersion parameters $\sigma_{r}$ and $\sigma_{z}$ whereas SN15 found that these decreased outwards with a scale length $\sim 6.2\,\mathrm{kpc}$. If we compare our values, $\sigma_{r}\sim 94\,\mathrm{km}\,\mathrm{s}^{-1}$ $\sigma_{z}\sim 130\,\mathrm{km}\,\mathrm{s}^{-1}$, with those, $\sigma_{r}\sim\sigma_{z}\sim 100\,\mathrm{km}\,\mathrm{s}^{-1}$, predicted by SB15 at $R\simeq 4\,\mathrm{kpc}$, where most of the disc cluster lie, the agreement is good except that we find $\sigma_{r}<\sigma_{z}$. Our scale length $\mathcal{R}_{\textup{d}}\sim 1.9\,\mathrm{kpc}$ is a bit shorter than that, $\mathcal{R}_{\textup{d}}=2.3\,\mathrm{kpc}$ obtained by SB15, but entirely in line with the findings of (Bovy et al., 2016) from an analysis of APOGEE data Majewski et al. (2016). Even though we have used no metallicity information, our value for the fraction of disc clusters, $F_{\rm disc}\sim 0.32\pm 0.07$, is consistent with the $\sim 31\pm 9$ per cent of clusters that are metal-rich (Harris et al., 2016).

It is natural to compare our halo DF with the EDF fitted to halo K giants by Das & Binney (2016). For definiteness, we compare our metallicity-blind halo DF with the EDF evaluated at $\hbox{[Fe/H]}=-2$. The K giants required essentially the same scale action ($\sim 5000\,\mathrm{kpc}\,\mathrm{km}\,\mathrm{s}^{-1}$) as the clusters, but our inner slope parameter ($\alpha\sim 0.8$) is slightly smaller than that ($\alpha\sim 1.3$) fitted to the K giants. Our outer slope parameter, $\beta\sim 8.8$, is definitely larger than that $\sim 5$ fitted to the K giants, with the consequence that in real space we predict the cluster density falls off as $\sim r^{-5}$ rather than $\sim r^{-4}$ for the K giants. Given the degeneracy between $\alpha,\beta$ and $J_{0}$, it is not clear that these differences are significant, but we have reason to expect we might expect the distribution of GCs to be less centrally peaked than that of halo giants: dynamical friction can drag globular clusters in to the densest part of the bulge, where they will be tidally destroyed (Tremaine et al., 1975; Gnedin et al., 2014). The axis ratio of the distribution of K giants seems to increase from $q\sim 0.7$ at small radii to unity at large radii just as does that of the cluster distribution.

We recover a probability density of clusters in real space that is consistent with previous work (Bica et al., 2006). We note, however, that a DF such as our halo DF, which is featureless in action space, gives rise to an interesting feature at $r\sim 10\,\mathrm{kpc}$ in real space as a consequence of the rapid decay of the quadrupole in the disc’s gravitational field. Specifically, as $r$ increases through $10\,\mathrm{kpc}$ the cluster system becomes spherical quite rapidly. Since it is natural for the DF to have a simple form in action space, a rapid reduction in flattening in all halo components around $r\sim 10\,\mathrm{kpc}$ is a robust prediction.

While neither the K giants (Das & Binney, 2016) nor Blue Horizontal Branch stars show clear rotation (Das et al., 2016), the great majority of the halo DFs in our MCMC chain have parts odd in $J_{\phi}$ that cause the system to rotate in the same sense as the disc. The rate of rotation is, however, slow, typically peaking at $\langle v_{\phi}\rangle<40\,\mathrm{km}\,\mathrm{s}^{-1}$, which is close to the upper limit of the rotation of the system of K giants halo. Samples of clusters drawn from favoured models of the entire cluster system yield values of the measure of rotation $v_{\rm rot}$ defined by Frenk & White (1980) that are consistent with the value obtained from the real clusters.

The only respect in which mock catalogues extracted from favoured models materially deviate from the data is the distribution of clusters close to the Galactic centre. With our adopted distance to the Galactic centre, $R_{0}=8.3\,\mathrm{kpc}$, the cluster distances in the Harris catalogue place significantly more clusters in front of the Galactic centre than behind it. Given that a small fractional error in the distance $s$ to a cluster near the Galactic centre gives rise to a large change in the cluster’s Galactocentric distance $r$, uncertainty in the distribution of clusters around the centre is inevitable. To make progress with this issue one would need to reconsider the distance to every cluster within, say $r=5\,\mathrm{kpc}$ to ensure that it is consistent with data that point to $R_{0}=8.3\,\mathrm{kpc}$. An alternative explanation of the excess of clusters in front of the Galactic centre is that we have under-estimated the bias arising from dust and confusion against discovering GCs located behind the bulge.

In this paper we have developed a robust framework within which theories about the GC system can be formulated and tested. With only 157 objects, models cannot be tightly constrained, but as the currently very sparse proper-motion data grow in volume and precision, it should be rewarding to revisit the present analysis. In the next few years our knowledge of the thick disc will increase markedly and this understanding should be encapsulated in an EDF. In view of the preliminary results we have obtained here, a promising line of enquiry would be to require the DF of the system of disc clusters to coincide with that of the stellar thick disc and see what halo DF is required to complement it.

As the data become more precise, it will be interesting to fit DFs separately to the high- and low-metallicity clusters: will the parameters of the disc and halo DFs that emerge from this exercise be similar to those found here?

Several interesting lines of enquiry are made possible by possession of a DF for the GCs. For example, we expect halo GCs to be clustered in action space as a consequence of more than one GC being stripped from a single satellite galaxy. One could seek evidence for clustering of GCs by comparing the frequency of pairs of objects at separation $\Delta$ in action space when (a) both objects are real GCs, (b) one object is a real GC and one is a pseudo GC drawn from the DF, and (c) both objects are pseudo GCs. Another interesting investigation would consider the speed with which a plausible GC DF would evolve through the action of dynamical friction and tidal destruction since the rates of both processes are fully specified by ${\bm{J}}$.

Acknowledgements

We thank the referee and members of the Oxford dynamics group for valuable comments on drafts of this paper. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 321067.