Chemo-dynamical substructure in the M31 inner halo globular clusters:
Further evidence for a recent accretion event

The data underpinning our analysis is the publicly-available catalogue of Caldwell & Romanowsky (2016), originally obtained through spectroscopic observations for nearly all of the known globular clusters in the inner regions of Andromeda using the Hectospec multi-fibre spectrograph (Fabricant et al., 2005) on the 6.5m MMT Observatory telescope in Arizona. Caldwell & Romanowsky (2016) used these spectra to provide a set of uniformly-derived line-of-sight velocities and metallicity estimates for approximately $94\%$ of globular clusters with projected galactocentric radii inside $21$ kpc. Our previous analysis of globular clusters in Andromeda’s outer halo (Mackey et al., 2019a, b) had an inner cut-off of $25$ kpc, so we searched the literature, and in particular the cluster catalogue from the Pan-Andromeda Archaeological Survey (PAndAS; Veljanoski et al., 2014), to ensure we were not missing any objects in the range $21-25$ kpc. We also updated the velocities from Caldwell & Romanowsky (2016) with more precise PAndAS measurements for a few clusters (G268, H14, and SK255B). In summary, the total base catalogue comprises 345 globular clusters extending to a projected galactocentric radius of $25$ kpc.

Caldwell & Romanowsky (2016) convincingly established that the most metal-rich globular clusters in the central regions of M31 are kinematically and spatially associated with its disk. Since our analysis is focused on Andromeda’s inner halo, we removed this component by excluding all clusters with metallicity $FeH\geq-0.4$²²2Throughout we use $FeH$ to represent $[\textrm{Fe}/\textrm{H}]$.. We further removed all globular clusters for which the observations were of insufficient quality to determine a reliable metallicity estimate.

The final sample for consideration contains $n=278$ globular clusters. Each cluster possesses the following data: $(x_{i},y_{i},v_{i},\sigma_{v_{i}},\widehat{FeH_{i}},\sigma_{FeH_{i}})$ – namely the Cartesian tangent-plane position centred on Andromeda, $(x_{i},y_{i})$, the Andromeda-centric line-of-sight velocity and associated uncertainty, $(v_{i},\sigma_{v_{i}})$, and the reported metallicity and associated uncertainty, $(\widehat{FeH_{i}},\sigma_{FeH_{i}})$; the metallicity and velocities of the globular cluster population, coupled with their uncertainties, is presented graphically in the Appendix as Figure 7. Cartesian positions correspond to tangent-plane distances centred on M31; these coordinates, and the Andromeda-centric velocities, were calculated according to the prescription outlined in Mackey et al. (2019a). Note that the transformation of Cartesian coordinates into regular polar coordinates is defined as:

Here, $\theta_{i}$ is not the position angle $(PA)$ on the sky that is usually quoted in astronomy, which is measured from North through East, but instead is defined from East through North. For convenience, the results of this study are converted into astronomical position angles at the completion of the analysis.

The typical individual $1\sigma$ uncertainty on $\widehat{FeH_{i}}$ is $\sigma_{FeH_{i}}=0.2$ dex, while that on the Andromeda-centric velocity $v_{i}$ is $\sigma_{v_{i}}=7.5$ km/s. This latter value is dominated by the measurement uncertainty (typically $\approx 6$ km/s), with a small contribution due to the assumed systemic motion of M31 (Mackey et al., 2019a).

As with previous studies (c.f. Mackey et al., 2019b), we explored three different kinematic models. Each of these is a nonlinear regression model that predicts line-of-sight velocity $v_{i}$ as a function of position $(x_{i},y_{i})$ or $(r_{i},\theta_{i})$, and potentially (when there are two populations) the metallicity $FeH_{i}$. The three models are denoted $V,S,F$ and possess different forms for the rotational component as presented in Table 1; note that these models will be subscripted with ${\rm I}$ and ${\rm II}$ if they represent models with a single or double kinematic components. $V$ is a simple sinusoidal function previously used in a variety of studies of the rotation of the globular cluster population in the halo of Andromeda (e.g., Veljanoski et al., 2013, 2014; Caldwell & Romanowsky, 2016), whereas $S$ corresponds to a solid-body rotation. Finally, $F$, represents an asymptotically-flat rotation curve which smoothly transitions from one side of the galaxy to the other. For a single-population model (e.g. $V_{\rm I}$), the index $k$ is fixed at $0$. When the metallicity cut is applied and the model as two kinematic components (e.g. $V_{\rm II}$), the globular cluster population below the metallicity cut is represented by subscript $0$, whereas the globular cluster population above the cut is represented by subscript $1$.

As well as rotation, it is also assumed that each globular cluster population possesses a velocity dispersion. We adopt the functional form of the velocity dispersion over the populations to be

where $\Sigma_{k0}$ and $s_{k}$ are parameters to be determined; clearly if $s_{k}=0$ then the velocity dispersion is spatially constant over the population.

The models considering two populations assume that there exists a metallicity threshold $FeH_{\rm cut}$ that cleaves the globular cluster sample into two dependent subgroups based upon whether an individual globular cluster’s metallicity, ${FeH_{i}}$, is above or below this limit. However, as the metallicity uncertainties, $\sigma_{FeH_{i}}$, can be significant, this is analogous to a curve fitting problem with error bars on the $x$-values (Hogg et al., 2010). Thus, we included all $n=278$ true metallicities as extra parameters to be inferred. We take the prior for each true metallicity to be a normal distribution centered at the measured value and with the appropriate standard deviation:

where $\widehat{FeH_{i}}$ is the given measurement and $\sigma_{FeH_{i}}$ the reported uncertainty.

For simplicity we did not use a hierarchical model, so the metallicities do not ’borrow strength’ from each other (i.e., there is no tendency for an unknown metallicity value to be dragged towards the preponderance of the other metallicities).

For a given rotational model $M_{j}=\{V_{j},\,S_{k},\,{\rm or}\ F_{j}\}$ (where $j={\rm I\ or\ II}$) which are dependent upon a set of parameters $\Theta_{k}$, and assuming a normal distribution for the velocity uncertainties on the globular clusters, we define the corresponding likelihood to be

where the velocity dispersion and the measurement uncertainty are added in quadrature such that

When considering a single-population model, such that $k=0$, the total likelihood is given by

However, when considering a model with two populations split by metallicity such that $k=0,1$, it is appropriate to select the corresponding kinematic component for the likelihood term in the product.

As typical in numerical calculations involving likelihood, practical considerations require that the log of the likelihood is computed.

For each model, we determine the posterior distribution for the parameters given the data, and the marginal likelihood of the model as a whole, using Diffusive Nested Sampling (Brewer et al., 2011; Brewer & Foreman-Mackey, 2018), a variant of Nested Sampling (Skilling, 2004, 2006).

For simplicity, we employed uniform priors for all parameters, with astronomically informed limits; these are given in Table 2. In Bayesian model comparison, the priors can have a significant effect on the marginal likelihoods. As a consequence, the marginal likelihoods we compute in this study are sensitive to the prior widths. However, this is not a major concern for two reasons. Firstly, we used the same priors for all parameters shared in common between models so they have similar predictive ability, and secondly, the main caveat to Bayesian model comparison with uniform priors is that very wide limits might overly penalise the more complex model. Since we found that the two-component model was favoured despite the wide uniform priors, a more sophisticated choice of priors would likely strengthen this result, rather than weakening it.

The Bayesian evidences from our nested sampling exploration of the posterior probability spaces are presented in Table 3. One immediate result is that, irrespective of the rotational model under consideration, the Bayesian evidence explicitly favours the two component rotational model over the single component, with a separation of the populations at metallicity $FeH\sim-2.2$ consistent across the models. For the best-fitting model, $\{F_{\rm II}\}$, the ratio of the Bayesian evidence between the single and two components fits, known as the Bayes factor, is 12.06, with it taken as "strong" evidence that the two component model is the better description of the data (Kass & Raftery, 1995). We summarise the properties of the best-fit model in the corner plot presented in Figure 1 and tabulated in Table 4; again, we emphasise that the characteristics of this best-fit are seen across all of the two component fits. Essentially, the the metal-poor component, which consists of $\sim 10-20$ globular clusters, has a rotational velocity amplitude of $\sim 340\textrm{km/s}$, substantially higher than the more metal-rich population (with a rotational amplitude of $\sim 55\textrm{km/s}$), \textcolorblackalthough we note that the large rotational amplitude is potentially driven by a small number of globular clusters at the velocity extremes. Furthermore, the velocity dispersion of the metal-poor component, at $\sim 200\textrm{km/s}$, is higher than that for the more metal-rich component $(\sim 170\textrm{km/s})$. Whilst these values refer to the $50^{\rm th}$ percentile of the posterior probability distributions, it should be noted that some of these show substantial skewness, and we note that with this the velocity dispersions are roughly consistent within the uncertainties. \textcolorblackWe comment on this potentially large velocity dispersion in the conclusions to this paper.

blackBefore closing this section, we return to the fact that metallicity cut was made to the initial sample at $FeH\geq-0.4$ to remove the galactic component identified by Caldwell & Romanowsky (2016, see Section 2.1). To ensure this selection has not biased the findings in this paper, several additional tests were undertaken which are summarised here. Essentially, single and multiple rotational components, using model $F$ (see Table 1), were fit for each of the components, and the Bayesian evidence calculated. In fitting a single rotational component, the resultant Bayesian evidence $Z_{1}=−2165.62$. The second model comprised of two components, but the prior on the metallicity cut was expanded to cover the entire metallicity range. The resultant evidence for this was $Z_{2}=−2162.58$, demonstrating that this is favoured over the single component model. However, the posterior for the metallicity cut, presented in the left-hand panel of Figure 3, possesses two distinct peaks, one at $FeH\sim-2$ and the other at $FeH\sim-0.4$, indicating that the modelling suggests that the globular cluster populations could be split into three. Fixing the metallicity cut at $FeH=-2.0$ recovers the results presented in this paper, namely that there are two kinematically distinct populations with significantly different rotational axes. However, fixing the metallicity cut at $FeH=-0.4$ shows that both populations share a common axis of rotation, with those at $FeH>-0.4$ possessing a rotation amplitude of $\sim 55$ km/s, whilst the lower metallicity population possesses a rotational amplitude of $\sim 160$ km/s. Finally, a three component model was fit to the entire population, with fixed metallicity cuts at $FeH=-2.0$ and $FeH=-0.4$. The resulting Bayesian evidence was $Z_{3}=−2156.5$, demonstrating that this is the favoured model. In the right-hand panel, we present the rotational axes of the three components from the posterior exploration, with $FeH<-2.0$ in blue, $-2.0\leq FeH\leq-0.4$ in red and $FeH>-0.4$ in green. The distribution of red and green axes overlap and are aligned with rotational axis of Andromeda, whilst the most metal-poor component in blue is distinctly offset from these. We conclude that disregarding globular clusters with $FeH>-0.4$ has not biased the findings in this paper.

If we assume that the inner globular clusters of Andromeda represent a mix of those formed in-situ coupled with phase-mixed accretion events, a natural explanation for the properties of the metal-poor globular clusters is that it is simply a kinematically hotter component of the overall population. However, if this was the correct interpretation, we would expect the orientation of the rotation of the globular cluster population to reflect that of the overall Andromeda Galaxy. To explore this, we present the orbital orientation of the two components for the best-fit model is displayed in Figure 2. The globular clusters are shown shown as filled circles colour-coded with their metallicity (left-hand panel) and velocity (right-hand panel). In both panels, more metal-poor clusters being larger. In the left-hand panel, the more metal-rich globular clusters are emphasised through their opacity, whereas it is the more metal-poor clusters emphasised in the right-hand panel. Underlying the globular clusters are solid lines denoting the major and minor axes of the Andromeda Galaxy, emphasised with overlying ellipses.

In the left-hand panel of Figure 2 the axes of rotation of the metal-poor component (blue) and metal-rich component (red) are presented underlying the globular cluster population; again these are taken from the $50^{\rm th}$ percentile of the posterior probability distribution, with one $\sigma$ uncertainties corresponding to the $16^{\rm th}$ and $84^{\rm th}$ percentiles (see Figure 1). The rotational axis of the metal-rich component is reasonably aligned with the minor axis of the Andromeda Galaxy, and so we conclude that their kinematic properties are related to the larger-scale galactic rotation. However, the rotational axis of the metal-poor component is substantially offset from the minor axis of the Andromeda Galaxy by $\sim 75\degr$. This is emphasised in the angular momentum vectors of the two components, with the metal-rich population at $PA\sim 115\degr$ whereas the metal-poor component is at $PA\sim 190\degr$³³3 The best-fit values for the orientation of the rotation are $\phi_{0}=1.41$ and $\phi_{1}=2.72$ radians for the Dulais Structure and dominant galactic population respectively (Table 4). However, given the form of the underlying model (Table 1), the direction of the rotational vector (given the right-hand rule) is these values $\pm\pi$. This corresponds to $\theta_{0}\sim-100\degr$ and $\theta_{1}\sim-25\degr$ in the model coordinates (Equation 1). Given that $\theta$ and $PA$ are related through $\theta+PA=90\degr$, this gives the quoted $PA$s for the orientations of the spin vectors of the two components. Hence, we concluded that the metal-poor globular clusters reveal the presence of a kinematically distinct population from the overall inner globular clusters; it is this that we identify as the Dulais Structure. We do, however, add a word of caution in interpreting these findings on the metallicity of the Dulais Structure as it is unlikely that globular cluster populations can cleanly be separated by metallicity alone and some of the low metallicity clusters are likely to be drawn from the galactic population, whereas the Dulais Structure  will have globular clusters above $FeH_{\rm cut}$. What the results of this study show is that the kinematics of the Dulais Structure  dominate at low metallicity, whereas it is lost against the kinematics of the galactic population at higher metallicities. This does indicate that the Dulais Structure  is, overall, more metal-poor than the galactic population of globular clusters. We explore this in more detail below.

The right-hand panel of Figure 2 presents the best-fit rotational model for the Dulais Structure component, colour-coded for velocity; clearly the direction of rotation is distinctly different to the rotation of Andromeda, which is aligned along the minor axis. Also in this figure, the blue arrow denotes the direction of the angular momentum vector for this rotation of the Dulais Structure component. The red arrow, however, presents the direction of the angular momentum for the outer globular cluster population that was associated with underlying stellar substructure, and therefore indicative of a relatively recent accretion event into the halo of the Andromeda Galaxy (Mackey et al., 2019b). Clearly, given that the direction and rotational sense of these inner and outer globular cluster populations are extremely similar, it is natural to suggest that the two components are related. One obvious conclusion is that the Dulais Structure  represents the remnants of a stellar system accreted as part of the larger-scale accretion seen in the outer halo of Andromeda.

As discussed previously, given the limitations of the data the sub-populations of globular clusters were defined solely on a metallicity cut. However, care must be taken in interpreting the metallicity distribution of the Dulais Structure based upon the metallicity cut as the dominant population is likely to extend below the cut, and conversely, there are likely to be globular clusters in the Dulais Structure which are more metal-rich than the metallicity cut. Hence, the metallicity cut is more likely to be indicative of where the Dulais Structure population of globular clusters comes to dominate over the larger galactic population. Given the best-fit model presented in Section 3.2, we here undertake an analysis of the posterior distributions to provide further insights into metallicity distribution of the Dulais Structure, although we readily note that this is indicative only, and more detailed modelling with more complete data is required to make robust statements.

For each globular cluster in our sample, we assesses which population they are more likely drawn from by considering the ratio of the likelihoods, $R$, of the two models. We present this in Figure 4, noting that small $R$ values correspond to globular clusters being more likely members of the more dominant galactic population rather than the Dulais Structure. For the sake of clarity in the coming analysis, we take an (arbitrary) line at $R=2$ and colour-code this sub-population of globular clusters in red. As well as a low metallicity tail, we see this sub-population has members that extend to higher metallicity.

With this separation of the populations, it is interesting to explore whether those more likely to be members of Dulais Structure possess clustering across any other properties. Of course, without full 3D positions and velocities, we are working within a limited representation of phase space. However, we define two quantities representing the projected angular momentum of the form;

where ${x,r}$ are the tangent plane coordinates and $v$ is line-of-sight velocity as the described in Section 2.1. These two quantities are presented in Figure 5 as a function of the metallicity, incorporating the same colour-coding with respect to $R$ as shown in Figure 4. In the upper panel, it is immediately apparent that the globular clusters are more likely members of the Dulais Structure are found at higher total angular momentum than the dominant galactic population. This is reinforced in the lower panel of Figure 5 which presents $L_{x}$ as a function of metallicity, which delineates the potential members of the Dulais Structure from the dominant population by having large positive values of this quantity. We note that this clustering is significant in this quantity due to the alignment of the rotational axis of the Dulais Structure with the $y$-coordinate of the tangent plane. It is again worth emphasising that whilst this posterior analysis cannot be considered a robust method of separating the dominant and Dulais Structure populations, such clustering in terms of angular momentum adds credence to there being two distinct sub-populations of globular clusters present.

Finally, in Figure 6, we present the spatial distribution of the globular cluster population (c.f. Figure 2) colour-coded with the relative likelihood parameter, $R$, with blue representing those most likely drawn from the dominant population, whilst those from the Dulais Structure are in red. Here, there is no obvious division between the two populations, and more in-depth chemical and kinematic analysis will be needed to cleanly separate the dominant population from the Dulais Structure.