# Two major accretion epochs in M31 from two distinct populations of globular clusters

Dougal Mackey, Geraint Lewis, Brendon Brewer, Annette Ferguson, Jovan Veljanoski, Avon Huxor, Michelle Collins, Patrick Côté, Rodrigo A Ibata, Mike Irwin, et al.

► To cite this version:

Dougal Mackey, Geraint Lewis, Brendon Brewer, Annette Ferguson, Jovan Veljanoski, et al. Two major accretion epochs in M31 from two distinct populations of globular clusters. Nature, 2019, 574 (7776), pp.69-71. 10.1038/s41586-019-1597-1. hal-03095182

HAL Id: hal-03095182

https://hal.science/hal-03095182v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Two major accretion epochs in M31 from two distinct populations of globular clusters

Dougal Mackey<fnref n="1" />*, Geraint F. Lewis², Brendon J. Brewer³, Annette M. N. Ferguson⁴, Jovan Veljanoski⁵, Avon P. Huxor⁶, Michelle L. M. Collins⁷, Patrick Côté⁸, Rodrigo A. Ibata⁹, Mike J. Irwin¹⁰, Nicolas Martin⁹,¹¹, Alan W. McConnachie⁸, Jorge Peñarrubia⁴, Nial Tanvir¹² & Zhen Wan²

Large galaxies grow through the accumulation of dwarf galaxies<fnref n="1" />,². In principle it is possible to trace this growth history via the properties of a galaxy's stellar halo³⁻⁵. Previous investigations of the galaxy Messier 31 (M31, Andromeda) have shown that outside a galactocentric radius of 25 kiloparsecs the population of halo globular clusters is rotating in alignment with the stellar disk⁶,⁷, as are more centrally located clusters⁸,⁹. The M31 halo also contains coherent stellar substructures, along with a smoothly distributed stellar component¹⁰⁻¹². Many of the globular clusters outside a radius of 25 kiloparsecs are associated with the most prominent substructures, but some are part of the smooth halo¹³. Here we report an analysis of the kinematics of these globular clusters. We find two distinct populations rotating perpendicular to each other. The rotation axis for the population associated with the smooth halo is aligned with the rotation axis for the plane of dwarf galaxies¹⁴ that encircles M31. We interpret these separate cluster populations as arising from two major accretion epochs, probably separated by billions of years. Stellar substructures from the first epoch are gone, but those from the more recent second epoch still remain.

We conducted the Pan-Andromeda Archaeological Survey to obtain a multi-band panoramic view of the M31 stellar halo out to a radius of about 150 kpc in projection¹⁵,¹⁶. As well as revealing extensive halo substructure from tidally disrupted dwarfs, the Pan-Andromeda Archaeological Survey uncovered a substantial number of globular clusters at large galactocentric radii¹⁷⁻¹⁹ (Fig. 1). The focus of the present study is the population of globular clusters in the outer halo of M31, which we define as the 92 objects known at projected radii larger than 25 kpc (ref. ¹³). We have previously undertaken a series of spectroscopic observations to obtain line-of-sight velocities for the majority of this sample (77 of 92 clusters) with a typical measurement uncertainty per target of about 10 km s⁻<fnref n="1" /> (refs. ⁶,⁷,¹³). Treated as a single population, these objects exhibit net rotation about the centre of M31 with amplitude approximately 85 km s⁻<fnref n="1" /> (refs. ⁶,⁷).

Our recent work¹³ has allowed us to classify many of the clusters in our sample as being physically associated with underlying halo substructures (a subgroup hereafter referred to as 'GC-Sub'), while other clusters exhibit no association with substructure ('GC-Non') and are plausibly part of the smooth halo. The remaining clusters are labelled as 'ambiguous'. With this division in place, we revisit our understanding of the kinematics of M31's outer globular cluster population. Our analysis is restricted to the region outside 25 kpc because at smaller radii our classification scheme is compromised by the high degree of overlapping substructure, plus extensive stellar material stripped from Andromeda's disk; moreover, the shorter dynamical times tend to erase any association between clusters and substructures more quickly.

We consider a family of three kinematic models to explain the observed cluster velocities, each consisting of an underlying rotational component and a superimposed velocity dispersion (see Extended Data Tables 1 and 2). For a given model, we examine two scenarios: first, a single global rotation plus dispersion to describe the entire globular cluster population, and then a solution comprising two distinct components corresponding to the GC-Sub and GC-Non cluster subgroups and with the ambiguous objects treated in a probabilistic fashion. We employ Bayesian inference to compute the probabilities of the models given the data, and determine their veracity using Bayes factors. The results are summarized in Extended Data Table 3.

Examining these, it is immediately apparent that the globular cluster population of the outer halo of Andromeda exhibits substantial global rotation, in close agreement with previous results⁶,⁷. However, it is also clear that, for each of the three models considered, the statistical evidence overwhelmingly favours the case in which two rotational components are present over that with a single rotational component. That is, dividing the globular cluster population into two physically motivated subgroups, each of which exhibits its own distinct rotational signature, yields a statistically significant (see Extended Data Table 3) improvement in explaining the observed data (as indicated by the >100-fold increase in Bayes factor in all three cases).

We turn to a brief description of the favoured model, $V_2$ in Extended Data Table 1, for which the inferred parameter values are listed in Extended Data Table 4 and the best-fit rotational kinematics are illustrated in Fig. 2. In general, similar conclusions hold for all of the models considered in our analysis.

In $V_2$, the inferred amplitude of rotation for each of the two cluster subgroups is similar, and non-zero at >3.5σ. In both cases it is comparable in magnitude (about 100 km s⁻<fnref n="1" />) to the velocity dispersion. In the Cartesian coordinate system adopted for our fits (see Fig. 1), the axis of rotation for GC-Non is oriented at -6.2° ± 13.4°, whereas that for GC-Sub is at 80.0° ± 10.0°, where the unceratinties are ±1 standard deviation (see Extended Data Table 4). It is notable that these rotational axes are orthogonal to each other, within the uncertainties. The projected angular momentum vectors for the two cluster subgroups are represented as arrows in Fig. 2 (adopting the right-hand convention), and are oriented at about 84° for GC-Non and about 170° for GC-Sub in the standard astronomical coordinate system in which angles are taken from North through East. For comparison, the rotation axis of Andromeda's stellar disk is oriented at 128° in the standard astronomical system, with a disk rotation amplitude of about 250 km s⁻<fnref n="1" /> and dispersion of about 60 km s⁻<fnref n="1" /> (refs. ²⁰,²¹). Within a radius of 20 kpc, the population of metal-poor globular clusters rotates with a similar amplitude to that inferred for our cluster subgroups (at about 100 km s⁻<fnref n="1" />), but exhibits a higher dispersion of about 150 km s⁻<fnref n="1" /> (refs. ⁸,⁹). The

<footnote n="1">¹Research School of Astronomy and Astrophysics, Australian National University, Canberra, Australian Capital Territory, Australia. ²Sydney Institute for Astronomy, School of Physics, The University of Sydney, Sydney, New South Wales, Australia. ³Department of Statistics, The University of Auckland, Auckland, New Zealand. ⁴Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, UK. ⁵Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands. ⁶HH Wills Physics Laboratory, University of Bristol, Bristol, UK. ⁷Department of Physics, University of Surrey, Guildford, UK. ⁸NRC Herzberg Astronomy and Astrophysics, Victoria, British Columbia, Canada. ⁹Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, Strasbourg, France. ¹⁰Institute of Astronomy, University of Cambridge, Cambridge, UK. ¹¹Max-Planck-Institut für Astronomie, Heidelberg, Germany. ¹²Department of Physics and Astronomy, University of Leicester, Leicester, UK. *e-mail: dougal.mackey@anu.edu.au</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1 |** Map showing the distribution of metal-poor red giant stars in Andromeda's halo. The density has been statistically corrected for contamination, and smoothed with a Gaussian kernel. Members of Andromeda's globular cluster system are marked with filled symbols. Those outside 25 kpc are coloured according to their line-of-sight velocities (where available) and have symbols according to their classification (circles for GC-Non, stars for GC-Sub, squares for ambiguous and diamonds for no data). Dashed circles represent projected galactocentric radii of 25 kpc and 150 kpc, respectively. The axes define the Cartesian coordinate system used in our models; following astronomical convention, East is to the left and North is up.</figcaption></figure>

rotation axis for the inner metal-poor clusters is in the same sense as Andromeda's stellar disk, but may plausibly be misaligned by about 20° to 30° in the direction of the GC-Non axis.

How are we to interpret the existence of rotational signatures with approximately perpendicular orientation amongst the outer globular cluster populations in M31? The GC-Sub sample consists of clusters that are members of well defined stellar streams or small dynamically cold subgroups in the halo. These clearly represent debris from one or more accretions that must have occurred relatively recently in order for the underlying structures to be still coherent$^{7,13,19}$. It is well established that Andromeda has undergone at least one substantial late merger, accreting the fourth or fifth largest galaxy in the Local Group approximately one billion years ago and producing the Giant Stellar Stream in the process$^{22,23}$. Currently available models follow only the late stages of this system's infall, and suggest that at the point of destruction the progenitor's orbit was probably highly radial$^{23}$. This makes it difficult to assess the degree to which the observed rotation of the GC-Sub sample could be due to this event.

Some recent analyses have advocated for a scenario in which M31 experienced a much larger (approximately 1:4 mass ratio) merger event within the past two to three billion years$^{24,25}$, based on more indirect lines of evidence. The high-amplitude rotational signal observed here for GC-Sub is certainly consistent with this overall picture, but we are at present unable to distinguish between the accretion of one large progenitor, or several smaller satellite systems with correlated angular momenta. One hypothesis arising from this scenario is that the compact elliptical satellite galaxy M32 may be the stripped core of the accreted progenitor$^{25}$. If true, we predict that future astrometric measurements of M32's transverse motion will imply a projected angular momentum vector aligned with that observed for GC-Sub; however, we also note that no globular clusters are known to be bound to this satellite.

The GC-Non sample consists of clusters that exhibit no evidence for association with halo substructure. As such, they have previously been linked with the smooth halo component in Andromeda$^{13}$. That these objects exhibit a substantial net angular momentum again suggests that they are plausibly due to either a single large merger or the accretion of several lesser progenitors with correlated orbits. However, this assembly epoch must have occurred billions of years earlier than the GC-Sub accretion(s) for phase mixing to erase coherent substructure over

multiple orbital periods; some memory of the original large-scale angular momentum is retained during this process, provided the shape of Andromeda's gravitational potential does not vary strongly with time.

It is particularly intriguing that the GC-Non sample rotates in the same sense, and is very closely aligned with, Andromeda's prominent plane of satellites$^{14,26}$. The rotational axis of this plane, comprising nearly half of the dwarf galaxy population orbiting Andromeda, is oriented at an angle of approximately -9.7° in the adopted Cartesian coordinates$^{26}$. This is within a few degrees of the GC-Non axis, well within the uncertainty on both quantities. This close alignment is surprising because satellite planes are fragile, short-lived structures, destroyed by precession and interaction with dark-matter halo substructure unless they are finely aligned within a spherical gravitational potential$^{27,28}$. Assuming that the properties of the GC-Non sample do indeed represent an ancient accretion epoch, the progenitor systems cannot have been members of the present-day plane. If the strong rotational alignment between the two structures is not coincidental, it might instead point to the larger-scale accretion of material onto Andromeda from a preferred direction in the local environment (for example, along a dominant dark-matter filament$^{29}$) over cosmic history.

Can we deduce any properties of the destroyed GC-Sub and GC-Non progenitors from their now-assimilated globular cluster populations? Empirically, the overall mass of globular clusters in a galaxy is a fraction (about $2.9 \times 10^{-5}$) of the total mass of the system that holds constant across a range of five orders of magnitude in host galaxy mass$^{30}$. In our favoured model $V_2$ the total number of clusters in GC-Sub is 43, and in GC-Non it is 34. These include probabilistically assigned members of the ambiguous set. Assuming a median cluster mass of $1.3 \times 10^5 M_\odot$ (where $M_\odot$ is the solar mass; refs. $^{13,30}$), the implied progenitor masses are then $1.9 \times 10^{11} M_\odot$ and $1.5 \times 10^{11} M_\odot$ respectively. These are lower limits, as the accretion events almost certainly deposited additional clusters into Andromeda's halo at radii inside 25 kpc. As an individual accretion, the GC-Sub event would correspond to a merger with M31 of a satellite galaxy having mass ratio 1:10, comparable to the Large Magellanic Cloud or the inferred Giant Stream progenitor$^{23}$. The GC-Non event is similar, but as it occurred much earlier in Andromeda's history, it would probably represent a more substantial relative merger.

The distinct rotational signatures we have uncovered for globular cluster subsystems in the outskirts of M31 offer an insight into the
<figure><img src="image_3.png" /><figcaption>**Fig. 2** | Best-fit rotational kinematics for the favoured model $V_2$. The best-fit parameters are listed in Extended Data Tables 1 and 4. **a**, Results for GC-Non; **b**, results for GC-Sub. The underlying colour illustrates the rotational kinematic signatures, while globular clusters are coloured according to their measured Andromeda-centric velocity. Clusters not in the applicable subsample are black symbols, and ambiguous clusters are</figcaption></figure>

not plotted but are apparent in Fig. 1. Dashed lines indicate the inferred rotational axes, while arrows represent the normalized angular momentum vectors. The fainter underlying arrows show the previously measured⁷ rotation axis for the entire outer cluster sample (matching model $V_1$), which is closely aligned with that for Andromeda's disk.

assembly history of our nearest large neighbour. Our results provide empirical evidence that Andromeda has undergone at least two major epochs of satellite accretion—one relatively recent in cosmological terms, and one that probably occurred many billions of years ago. The apparent kinematic alignment between one or more long-destroyed ancient progenitor systems and the present-day plane of satellites may indicate long-term coherent accretion from the cosmic web. These represent important steps towards a detailed reconstruction of the main events that led to the Local Group of galaxies that we see today.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-019-1597-1.

Received: 22 March 2019; Accepted: 16 July 2019;

Published online 2 October 2019.

1. Searle, L. & Zinn, R. Compositions of halo clusters and the formation of the galactic halo. Astrophys. J. 225, 357–379 (1978).

2. Springel, V. et al. Simulations of the formation, evolution and clustering of galaxies and quasars. Nature 435, 629–636 (2005).

3. Bullock, J. S. & Johnston, K. V. Tracing galaxy formation with stellar halos. I. Methods. Astrophys. J. 635, 931–949 (2005).

4. Cooper, A. P. et al. Galactic stellar haloes in the CDM model. Mon. Not. R. Astron. Soc. 406, 744–766 (2010).

5. Johnston, K. V. et al. Tracing galaxy formation with stellar halos. II. Relating substructure in phase and abundance space to accretion histories. Astrophys. J. 689, 936–957 (2008).

6. Veljanoski, J. et al. Kinematics of outer halo globular clusters in M31. Astrophys. J. Lett. 768, L33 (2013).

7. Veljanoski, J. et al. The outer halo globular cluster system of M31—II. Kinematics. Mon. Not. R. Astron. Soc. 442, 2929–2950 (2014).

8. Perrett, K. M. et al. The kinematics and metallicity of the M31 globular cluster system. Astron. J. 123, 2490–2510 (2002).

9. Caldwell, N. & Romanowsky, A. J. Star clusters in M31. VII. Global kinematics and metallicity subpopulations of the globular clusters. Astrophys. J. 824, 42 (2016).

10. Ferguson, A. M. N., Irwin, M. J., Ibata, R. A., Lewis, G. F. & Tanvir, N. R. Evidence for stellar substructure in the halo and outer disk of M31. Astron. J. 124, 1452–1463 (2002).

11. Ibata, R. A. et al. The haunted halos of Andromeda and Triangulum: a panorama of galaxy formation in action. Astrophys. J. 671, 1591–1623 (2007).

12. Ibata, R. A. et al. The large-scale structure of the halo of the Andromeda galaxy. I. Global stellar density, morphology and metallicity properties. Astrophys. J. 780, 128 (2014).

13. Mackey, A. D. et al. The outer halo globular cluster system of M31—III. Relationship to the stellar halo. Mon. Not. R. Astron. Soc. 484, 1756–1789 (2019).

14. Ibata, R. A. et al. A vast, thin plane of corotating dwarf galaxies orbiting the Andromeda galaxy. Nature 493, 62–65 (2013).

15. McConnachie, A. W. et al. The remnants of galaxy formation from a panoramic survey of the region around M31. Nature 461, 66–69 (2009).

16. McConnachie, A. W. et al. The large-scale structure of the halo of the Andromeda galaxy. II. Hierarchical structure in the Pan-Andromeda Archaeological Survey. Astrophys. J. 868, 55 (2018).

17. Huxor, A. P. et al. Globular clusters in the outer halo of M31: the survey. Mon. Not. R. Astron. Soc. 385, 1989–1997 (2008).

18. Huxor, A. P. et al. The outer halo globular cluster system of M31—I. The final PAndAS catalogue. Mon. Not. R. Astron. Soc. 442, 2165–2187 (2014).

19. Mackey, A. D. et al. Evidence for an accretion origin for the outer halo globular cluster system of M31. Astrophys. J. Lett. 717, L11–L16 (2010).

20. Ibata, R. A. et al. On the accretion origin of a vast extended stellar disk around the Andromeda galaxy. Astrophys. J. 634, 287–313 (2005).

21. Dorman, C. E. et al. A clear age-velocity dispersion correlation in Andromeda's stellar disk. Astrophys. J. 803, 24 (2015).

22. Ibata, R. A., Irwin, M., Lewis, G., Ferguson, A. M. N. & Tanvir, N. A giant stream of metal-rich stars in the halo of the galaxy M31. Nature 412, 49–52 (2001).

23. Fardal, M. A. et al. Inferring the Andromeda galaxy's mass from its giant southern stream with Bayesian simulation sampling. Mon. Not. R. Astron. Soc. 434, 2779–2802 (2013).

24. Hammer, F. et al. A 2–3 billion year old major merger paradigm for the Andromeda galaxy and its outskirts. Mon. Not. R. Astron. Soc. 475, 2754–2767 (2018).

25. D'Souza, R. & Bell, E. F. The Andromeda galaxy's most important merger about 2 billion years ago as M32's likely progenitor. Nature Astron. 2, 737–743 (2018).

26. Conn, A. R. et al. The three-dimensional structure of the M31 satellite system; strong evidence for an inhomogeneous distribution of satellites. Astrophys. J. 766, 120 (2013).

27. Fernando, N. et al. On the stability of satellite planes—I. Effects of mass, velocity, halo shape and alignment. Mon. Not. R. Astron. Soc. 465, 641–652 (2017).

28. Fernando, N., Arias, V., Lewis, G. F., Ibata, R. A. & Power, C. Stability of satellite planes in M31 II: effects of the dark subhalo population. Mon. Not. R. Astron. Soc. 473, 2212–2221 (2018).

29. Libeskind, N. I. et al. Planes of satellite galaxies and the cosmic web. Mon. Not. R. Astron. Soc. 452, 1052–1059 (2015).

30. Harris, W. E., Blakeslee, J. P. & Harris, G. L. H. Galactic dark matter halos and globular cluster populations. III. Extension to extreme environments. Astrophys. J. 836, 67 (2017).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019
## METHODS

**Observational data.** The globular clusters used in this study were discovered as part of the Pan-Andromeda Archaeological Survey, a multi-band panoramic imaging programme targeting the Andromeda (M31) and Triangulum (M33) galaxies out to physical distances (in projection) of ~150 kpc and ~50 kpc respectively, which used the MegaCam wide-field camera on the 3.6-m Canada-France-Hawaii Telescope<sup>15,16,31</sup>. We focus on the 92 objects known at projected radii larger than 25 kpc (ref. <sup>13</sup>). Our motivation for selecting this as the inner radial limit for our analysis is described in more detail below.

Line-of-sight velocities for globular clusters in this sample were obtained through a series of spectroscopic observations, using the Intermediate-dispersion Spectrograph and Imaging System (ISIS) on the 4.2-m William Herschel Telescope at the Observatorio del Roque de los Muchachos on the island of La Palma, the Ritchey-Chrétien (RC) spectrograph at the 4-m Mayall Telescope at the Kitt Peak National Observatory in Arizona, the Gemini Multi-Object Spectrograph (GMOS) on the 8-m Gemini North Telescope in Hawaii, and the Deep Imaging Multi-Object Spectrograph (DEIMOS) on the 10-m Keck II Telescope, also in Hawaii. The full set of outer halo clusters with measured velocities numbers 77 objects of the 92 known outside 25 kpc (refs. <sup>6,7,13</sup>). For the present analysis, the line-of-sight velocities were transformed into the Galactocentric frame and then an M31-centred frame using the procedure outlined in these previous contributions. This accounts for projection effects due to both the solar motion and the systemic velocity of M31, induced by the large angular span of the M31 halo on the sky. The typical measurement uncertainty per cluster is of the order of 10 km s<sup>-1</sup>.

Overall, some 500 globular clusters have been catalogued in M31<sup>9,13,32–35</sup>, the majority of which (~80%) sit at radii smaller than our inner limit of 25 kpc. Another ~50 clusters are known to be members of large Andromedan satellites such as M33, NGC 147, NGC 185 and NGC 205 (ref. <sup>15</sup>). Although these various cluster populations were excluded from our analysis, we used these previous comparisons to plot them in Fig. 1 for context.

**Relationship to stellar substructure.** Details of the relationship between the outer globular cluster population and the underlying stellar halo in Andromeda were developed in a recent publication<sup>13</sup>. In summary, the clusters can be assigned to one of three groups: those that exhibit strong spatial and/or kinematic evidence for a link with a substructure in the field halo, those that show no such evidence, and those for which the evidence is weak or conflicting. Throughout this paper, these three sets are labelled as 'substructure globular clusters' (GC-Sub), 'non-substructure globular clusters' (GC-Non) and 'ambiguous globular clusters', respectively.

‘Spatial’ evidence is quantified by measuring the density of halo stars locally around each cluster, and finding those objects with high local densities relative to the observed distribution at commensurate galactocentric radii. ‘Kinematic’ evidence is limited to identifying objects that are members of small-scale dynamically cold groupings; no information about global velocity patterns is used. For the 77 clusters considered here, this classification scheme yields subgroups of 32 objects in GC-Sub and 26 in GC-Non, with the remaining 19 identified as ambiguous. In the Methods section ‘Handling of ambiguous globular clusters’ we detail how these ambiguous globular clusters are incorporated into our models.

The inner halo of M31 is highly complex. At radii smaller than 25 kpc, the stellar substructures become so pervasive, and the degree of overlap so great, that the classification scheme outlined above breaks down<sup>13</sup>. Moreover, the dynamical times become sufficiently short that the association between stellar substructures and globular clusters from accreted satellites is far less persistent than at larger radii. This inner region is also complicated by the presence of extensive stellar material disrupted from the disk of Andromeda<sup>20,21,36</sup>, as well as a substantial in situ globular cluster population (that is, non-accreted objects that formed in Andromeda’s disk or bulge)<sup>9</sup>. These issues motivate our choice of 25 kpc as the inner radius defining the sample for analysis; our key results are not sensitive to mild variations in the adopted limit (20–40 kpc).

**Kinematic models.** To explore the velocity characteristics of M31’s outer globular clusters, we consider models in which the population possesses an overall velocity dispersion superimposed upon an underlying net rotation. We employ three physically motivated models for the rotational component, which are detailed in Extended Data Table 1. The V models, which were used in a previous study on the bulk rotation of globular cluster systems<sup>6,7</sup>, consider a constant amplitude of rotational velocity modulated by an angular dependence. The S models represent solid-body rotation, where the rotational velocity increases linearly from the axis of rotation. Finally, the F models consider an asymptotically flat velocity profile with a smooth transition through the axis of rotation, akin to the rotation curves of spiral galaxies. These models represent three broad families of potential rotation and, as will be seen in detail below, the choice of a specific model does not affect the need for two rotational components to explain the observed velocity properties of Andromeda’s globular cluster population.

For each rotational model, we consider two distinct scenarios. Initially, we treat the entire globular cluster population as a whole, considering a single underlying

global rotation<sup>6,7</sup>. Following this, we treat GC-Sub and GC-Non as distinct populations, each possessing their own characteristic rotational component.

In addition to the rotational component, the models also include a velocity dispersion as a function of position given by:

$$ \Sigma(r, \theta; \Lambda) = \sigma \left( \frac{r}{R_0} \right)^{\gamma} \quad (1) $$

where r is the projected galactocentric radius, θ is the angular coordinate in the chosen system (detailed below), γ is a slope parameter and σ is the velocity dispersion at projected radius R₀; to compare with previous work, and as this parameter is degenerate with σ and γ, we fix R₀ = 30 kpc (refs. <sup>6,7</sup>). The model parameters to be estimated, including those for the rotational components of the model, are denoted collectively by Λ. This functional form for Σ was used for all of the rotational models. For the models where we consider two distinct rotational components, we also consider each to have their own distinct dispersion profiles, as described by equation (1). **Bayesian inference.** To test these models for the rotational pattern exhibited by the globular clusters, we employ Bayesian inference to compute the probabilities of the models given the data. Clearly, the models have differing numbers of parameters, but we can employ the marginal likelihood (also known as the evidence) to test the various propositions.

Throughout this paper we adopt a Cartesian coordinate system in the tangent plane centred on M31, with the x axis aligned with West and the y axis aligned with North in the equatorial astronomical system (see Fig. 1). In detailing our approach, we let (rᵢ, θᵢ) be the position of the ith globular cluster on the sky, in plane polar coordinates with respect to the above coordinate system, and let vᵢ be its velocity along the line of sight, including rotation, velocity dispersion and measurement error. Our model attempts to predict vᵢ as a function of (rᵢ, θᵢ), in a manner analogous to regression. The rotational component is specified in terms of a functional form, f(r, θ; Λ), which gives the model-predicted velocity as a function of position. This is combined with the model-predicted line-of-sight velocity dispersion, Σ(r, θ; Λ) (see equation (1)), to specify the full kinematic model.

The probability distribution for the data (that is, the measured velocities) vᵢ, given the parameters and one of the models from Extended Data Table 1, is:

$$ v_i | \Lambda \sim \text{Normal}(f(r, \theta; \Lambda), \Sigma(r, \theta; \Lambda)^2 + s_i^2) \quad (2) $$

where ~ means ‘is distributed as’ and the {sᵢ} quantities are the uncertainties on the velocity measurements {vᵢ}. With this choice, the likelihood function can be written explicitly as:

$$ L(\Lambda) = \prod_{i=1}^{N} \frac{1}{\sqrt{2\pi(\Sigma(r_i, \theta_i; \Lambda)^2 + s_i^2)}} \exp\left(-\frac{(v_i - f(r_i, \theta_i; \Lambda))^2}{2(\Sigma(r_i, \theta_i; \Lambda)^2 + s_i^2)}\right) \quad (3) $$

where N is the total number of globular clusters in the sample. In this paper, ‘parameter estimation’ involves computing the posterior distributions for the parameters Λ given the data assuming a model from Extended Data Table 1 (that is, given a particular choice for f), and ‘model comparison’ involves the assessment of different choices for the form of the function f.

**Handling of ambiguous globular clusters.** For 19 out of the 77 clusters, it is unclear whether or not a substructure is present. For each of these ambiguous objects, we assign a flag fᵢ which takes the value of 1 if the cluster really is associated with a substructure and 0 if it is not. We set the prior for the {fᵢ} parameters by introducing a hyperparameter p<sub>subs</sub>, which sets the prior expected proportion of the ambiguous globular clusters that actually are associated with a substructure. The prior for p<sub>subs</sub> and the {fᵢ} flags is then:

$$ p_{\text{subs}} \sim \text{Uniform}(0, 1) \quad (4) $$

$$ f_j \sim \text{Bernoulli}(p_{\text{subs}}) \quad (5) $$

We have also tested alternative approaches for treating the ambiguous globular clusters, such as omitting them completely from the sample, and found this to have only minor effects on the results. Hence our conclusions are not driven by the treatment of the ambiguous clusters.

**Prior distributions.** The joint prior distribution for all parameter and data values reflects all the assumed information used in the analysis apart from the data itself. This is specified in Extended Data Table 2. Single-component models may be obtained by ignoring all parameters referring to the substructure component. For positive parameters of unknown order of magnitude, we have assigned truncated Cauchy priors to the logarithm of the parameter. The truncation is for numerical convenience, and the Cauchy distribution weakly emphasizes values within a few orders of unity. For two-component models, the prior for the second component implies a high probability that it is within an order of magnitude of the first component.
This avoids any fatal problems caused by the Jeffreys–Lindley paradox when using naive wide priors for additional parameters in Bayesian model comparison problems. **Computation and results.** We used DNest4³⁷,³⁸ to sample the posterior distributions for the parameters of each model, and to compute the marginal likelihood of each model. The marginal likelihood of a model *M* is the average value of the likelihood function with respect to the prior distribution for the parameters, and plays the role of a likelihood when comparing models.

The estimated marginal likelihoods are given in Extended Data Table 3. The data are far more probable under models $V_2$ and $F_2$ than the others. Within each class of models (S, V, F), the two-component model makes the data far more likely than the one-component model. By Bayes' theorem, the two-component models are thus much more probable than the one-component models, assuming they have nontrivial prior probabilities. If prior probabilities of 1/6 are assigned to each model, the posterior probability for two substructures (that is, for the proposition $V_2 \lor S_2 \lor F_2$) is 0.9994. Crucially, this is conditional on $\top \equiv S_1 \lor S_2 \lor V_1 \lor V_2 \lor F_1 \lor F_2$ (the proposition that it is one of these models), and other sets of assumptions may always be considered. The marginal likelihood of $\top$ is -468.51, which can be used to compare our findings with any future analysis³⁹.

Estimates for the parameters in the overall favoured model, $V_2$, are reported in Extended Data Table 4 and were used in the construction of Fig. 2. A corner plot showing the marginal distributions for the parameters is provided in Extended Data Fig. 1. We note that the second-best model by marginal likelihood, $F_2$, while still favouring a two-rotational-component model over a single-rotational-component model, possesses a slightly different angle for the orientation of GC-Sub. that is, $\phi_1 \approx 60^\circ$. Under $F_2$, the orientation of GC-Non is $\phi_0 = 0^\circ$, consistent with the other models. Investigating the posterior distributions derived for $F_2$ reveals that the preferred scale parameter, $L$, becomes unphysically small, resulting in a step function in velocity across the rotational axis. This results in step functions in the likelihood, with the few individual globular clusters close to the rotation axis strongly constraining $\phi_1$.

The analysis in this paper considered the 77 globular clusters with measured velocities, from an overall sample of 92. We note that the sample is therefore incomplete. However, examining the marginal distributions presented in Extended Data Fig. 1, we see that the velocity amplitudes are strongly constrained to be nonzero, and similarly the orientations are precisely constrained. It is unlikely that the velocities for clusters in the remaining sample are extreme enough to counter these outcomes.

## Data availability

All data analysed for this study are publicly available. Pan-Andromeda Archaeological Survey data products, including the stellar photometry catalogue, reduced individual images, and image stacks, may be downloaded from the Canadian Astronomical Data Center (CADC; http://www.cADC-cda.cda.hia-iha.nrc-cnrc.gc.ca/en/community/pandas/query.html)¹⁶. Globular cluster locations, radial velocities, and classifications are published online⁷,¹³,¹⁸ and are also found in the code repository (https://github.com/eggplantbren/AndromedaMixture).

## Code availability

Readers may access the code used for the inference, with the data included, at https://github.com/eggplantbren/AndromedaMixture. The code has been released under the permissive MIT licence. The README file describes how to reproduce the results of the paper.

31. Martin, N. F. et al. The PAndAS view of the Andromeda satellite system. I. A Bayesian search for dwarf galaxies using spatial and color-magnitude information. *Astrophys. J.* **776**, 80 (2013).

32. Galleti, S., Federici, L., Bellazzini, M., Fusi Pecci, F. & Macrina, S. 2MASS NIR photometry for 693 candidate globular clusters in M 31 and the Revised Bologna Catalogue. *Astron. Astrophys.* **416**, 917–924 (2004).

33. Caldwell, N. et al. Star clusters in M31. I. A catalog and a study of the young clusters. *Astron. J.* **137**, 94–110 (2009).

34. Peacock, M. B. et al. The M31 globular cluster system: ugriz and K-band photometry and structural parameters. *Mon. Not. R. Astron. Soc.* **402**, 803–818 (2010).

35. Caldwell, N., Schiavon, R., Morrison, H., Rose, J. A. & Harding, P. Star clusters in M31. II. Old cluster metallicities and ages from Hectospec data. *Astron. J.* **141**, 61 (2011).

36. Dorman, C. E. et al. A new approach to detailed structural decomposition from the SPLASH and PHAT surveys: kicked-up disk stars in the Andromeda galaxy? *Astrophys. J.* **779**, 103 (2013).

37. Brewer, B. J., Párta, L. B. & Csányi, G. Diffusive nested sampling. *Stat. Comput.* **21**, 649–656 (2011).

38. Brewer, B. J. & Foreman-Mackey, D. DNest4: diffusive nested sampling in C++ and Python. *J. Stat. Softw.* **86**, 1–33 (2018).

39. Skilling, J. Nested sampling for general Bayesian computation. *Bayesian Anal.* **1**, 833–859 (2006).

40. Foreman-Mackey, D. corner.py: scatterplot matrices in Python. *J. Open Source Softw.* **1**, 24 (2016).

**Acknowledgements** This work is based in part on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. This work is further based in part on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), the Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), the Ministério da Ciência, Tecnologia e Inovação (Brazil) and the Korea Astronomy and Space Science Institute (Republic of Korea). Some of the data presented here were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. We wish to recognise and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. The William Herschel Telescope (WHT) is operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This work is based in part on observations at Kitt Peak National Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation. We are honoured to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham. D.M. is supported by an Australian Research Council (ARC) Future Fellowship (FT160100206). G.F.L. acknowledges support from a Partnership Collaboration Award between the University of Sydney and the University of Edinburgh. D.M. and G.F.L. appreciate the hospitality of the Royal Observatory, Edinburgh, where the final stages of the preparation of this paper were undertaken. B.J.B. thanks the Marsden Fund of the Royal Society of New Zealand. This work has been published under the framework of the IdEx Unistra and benefits from funding from the state managed by the French National Research Agency as part of the investments for the future program. Z.W. is supported by a Dean's International Postgraduate Research Scholarship at the University of Sydney.

**Author contributions** A.W.M., R.A.I., M.J.I., A.M.N.F., G.F.L. and N.T. initiated the Pan-Andromeda Archaeological Survey, with extensive data analysis and interpretation undertaken with N.M., M.L.M.C., P.C. and J.P. D.M., A.M.N.F., J.V. and A.P.H. were responsible for the discovery and characterization of Pan-Andromeda Archaeological Survey globular clusters, for measuring their line-of-sight velocities, and for conducting detailed earlier analyses of this population. G.F.L. was responsible for the development of the kinematic models employed in this study. G.F.L. and B.J.B. undertook the statistical analysis, with B.J.B. responsible for implementing and running the kinematic models in DNest4. Z.W. was responsible for undertaking geometric transformations into the Andromeda frame to enable comparisons with previous studies. All authors assisted in the interpretation of the results and writing of the paper.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Correspondence and requests for materials** should be addressed to D.M.

**Reprints and permissions information** is available at http://www.nature.com/ reprints.
<figure><img src="image_4.png" /><figcaption>**Extended Data Fig. 1** | Corner plot of the posterior distribution of the parameters for the V₂ model. This model is defined in Extended Data Table 1. Subscripts 0 and 1 correspond to the parameters for the non-</figcaption></figure>

substructure (GC-Non) and substructure (GC-Sub) samples respectively, and a summary of the marginal distributions is given in Extended Data Table 4. Corner plots were generated using Python⁴⁰.
<table><caption>Extended Data Table 1 | Functional form of the rotational component for the models under consideration</caption><thead><tr><th>Model</th><th>Description</th><th>Functional form of f()</th></tr></thead><tbody><tr><td>V<sub>1</sub></td><td>Veljanoski et al., one component</td><td>A<sub>0</sub> sin(θ - φ<sub>0</sub>)</td></tr><tr><td>V<sub>2</sub></td><td>Veljanoski et al., two components</td><td>A<sub>k</sub> sin(θ - φ<sub>k</sub>)</td></tr><tr><td>S<sub>1</sub></td><td>Solid-body, one component</td><td>A<sub>0</sub> (x sin φ<sub>0</sub> - y cos φ<sub>0</sub>)</td></tr><tr><td>S<sub>2</sub></td><td>Solid-body, two components</td><td>A<sub>k</sub> (x sin φ<sub>k</sub> - y cos φ<sub>k</sub>)</td></tr><tr><td>F<sub>1</sub></td><td>Asymptotically-flat rotation curve</td><td>A<sub>0</sub> tanh((x sin φ<sub>0</sub> - y cos φ<sub>0</sub>)/L<sub>0</sub>)</td></tr><tr><td>F<sub>2</sub></td><td>Asymptotically-flat, two components</td><td>A<sub>k</sub> tanh((x sin φ<sub>k</sub> - y cos φ<sub>k</sub>)/L<sub>k</sub>)</td></tr><tr><td>⊤</td><td>Disjunction of all of the above</td><td></td></tr></tbody></table>

The third column gives the model-predicted line-of-sight velocity as a function of position on the sky, expressed in a combination of Cartesian coordinates (x, y) and plane polar coordinates (r, φ) defined as specified in the text. For two-component models the subscript k takes values of 0 for GC-Non and 1 for GC-Sub. For one-component models, only subscript 0 applies.
<table><caption>Extended Data Table 2 | The joint prior distribution for all parameters and data</caption><thead><tr><th>Quantity</th><th>Meaning</th><th>Distribution</th></tr></thead><tbody><tr><td colspan="3"><strong>Parameters</strong></td></tr><tr><td>A<sub>0</sub></td><td>Amplitude of rotation</td><td>ln(A<sub>0</sub>) ~ Cauchy(0, 5)T(-100, 100)</td></tr><tr><td>A<sub>1</sub></td><td>Amplitude of rotation</td><td>ln(A<sub>1</sub>) ~ Normal(ln A<sub>0</sub>, 1)</td></tr><tr><td>φ<sub>0</sub></td><td>Orientation of rotation axis</td><td>Uniform(-180°, 180°)</td></tr><tr><td>φ<sub>1</sub></td><td>Orientation of rotation axis</td><td>Uniform(-180°, 180°)</td></tr><tr><td>σ<sub>0</sub></td><td>Velocity dispersion at distance R<sub>0</sub></td><td>ln σ<sub>0</sub> ~ Cauchy(0, 5)T(-100, 100)</td></tr><tr><td>σ<sub>1</sub></td><td>Velocity dispersion at distance R<sub>0</sub></td><td>ln σ<sub>0</sub> ~ Normal(ln σ<sub>0</sub>, 1)</td></tr><tr><td>γ<sub>0</sub></td><td>Velocity dispersion slope</td><td>Uniform(-3, 0)</td></tr><tr><td>γ<sub>1</sub></td><td>Velocity dispersion slope</td><td>Uniform(-3, 0)</td></tr><tr><td>p<sub>subs</sub></td><td>Hyperparameter for {f<sub>j</sub>}</td><td>Uniform(0, 1)</td></tr><tr><td>{f<sub>j</sub>}</td><td>Substruct. flags for ambiguous GCs</td><td>Bernoulli(p<sub>subs</sub>)</td></tr><tr><td>L<sub>0</sub></td><td>Length scale of transition</td><td>ln L<sub>0</sub> ~ Uniform(ln(10<sup>-3</sup>R<sub>0</sub>), ln R<sub>0</sub>)</td></tr><tr><td>L<sub>1</sub></td><td>Length scale of transition</td><td>ln L<sub>1</sub> ~ Uniform(ln(10<sup>-3</sup>R<sub>1</sub>), ln R<sub>1</sub>)</td></tr><tr><td colspan="3"><strong>Data</strong></td></tr><tr><td>{v<sub>i</sub>}</td><td>Measured line-of-sight GC velocities</td><td>Equation 2</td></tr><tr><td colspan="3"><strong>Prior info.</strong></td></tr><tr><td>N</td><td>Number of GCs</td><td>Given (77)</td></tr><tr><td>{(r<sub>i</sub>, θ<sub>i</sub>)}</td><td>Positions of GCs on sky</td><td>Given</td></tr><tr><td>{g<sub>i</sub>}</td><td>Classification flags for all GCs</td><td>Given (GC-Non / GC-Sub / Ambiguous)</td></tr></tbody></table>

For two-component models, subscript 0 corresponds to GC-Non and subscript 1 applies to GC-Sub. For one-component models, only subscript 0 applies. The length scale parameters L₀ and L₁ apply only to models F₁ and F₂ and have no effect in the other models. We note that the orientation angles φ are defined according to the Cartesian coordinate system specified in the text, which differs from the astronomical convention of angles measured North through East by 90°. GC, globular cluster.
<table><caption>Extended Data Table 3 | Marginal likelihoods for each of the models, along with the resulting Bayes factor compared to the most favoured model $V_2$</caption><thead><tr><th>Model</th><th>Seed</th><th>ln(Z)</th><th>Z/Z<sub>max</sub></th></tr></thead><tbody><tr><td>V<sub>1</sub></td><td>0</td><td>-474.41</td><td>7.85 × 10<sup>-4</sup></td></tr><tr><td>V<sub>2</sub></td><td>1</td><td>-467.26</td><td>1</td></tr><tr><td>S<sub>1</sub></td><td>2</td><td>-479.64</td><td>4.20 × 10<sup>-6</sup></td></tr><tr><td>S<sub>2</sub></td><td>3</td><td>-475.29</td><td>3.26 × 10<sup>-4</sup></td></tr><tr><td>F<sub>1</sub></td><td>4</td><td>-475.45</td><td>2.77 × 10<sup>-4</sup></td></tr><tr><td>F<sub>2</sub></td><td>5</td><td>-467.59</td><td>0.719</td></tr><tr><td>T</td><td></td><td>-468.51</td><td></td></tr></tbody></table>

Models $V_2$ and $F_2$ are strongly favoured by the data. 'Seed' refers to the random number generator seed employed in each run of the Markov chain Monte Carlo exploration, and is provided for reproducibility reasons. Z is the marginal likelihood.
<table><caption>Extended Data Table 4 | Estimates of the parameters for the favoured model $V_2$</caption><thead><tr><th>Parameter</th><th>Post. Mean ± Post. SD</th><th>Max. Likelihood</th></tr></thead><tbody><tr><td>A<sub>0</sub> (km/s)</td><td>92 ± 27</td><td>120</td></tr><tr><td>A<sub>1</sub> (km/s)</td><td>110 ± 15</td><td>110</td></tr><tr><td>φ<sub>0</sub> (°)</td><td>-6.2 ± 13.4</td><td>-3.9</td></tr><tr><td>φ<sub>1</sub> (°)</td><td>80.0 ± 10.0</td><td>86.0</td></tr><tr><td>σ<sub>0</sub> (km/s)</td><td>134 ± 36</td><td>101</td></tr><tr><td>σ<sub>1</sub> (km/s)</td><td>133 ± 30</td><td>130</td></tr><tr><td>γ<sub>0</sub></td><td>-0.66 ± 0.32</td><td>-0.52</td></tr><tr><td>γ<sub>1</sub></td><td>-0.84 ± 0.33</td><td>-0.95</td></tr><tr><td>p<sub>subs</sub></td><td>0.67 ± 0.20</td><td>0.60</td></tr></tbody></table>

These are reported as the posterior mean ± the posterior standard deviation. Since none of the marginal posterior distributions are extremely asymmetric or non-Gaussian (Extended Data Fig. 1), these are approximately 68% credible intervals. As noted previously, a subscript of 0 corresponds to GC-Non and 1 to GC-Sub. The maximum likelihood estimates, given in the final column, were used in the construction of Fig. 2.