The outer halo globular cluster system of M31 – II. Kinematics
We present a detailed kinematic analysis of the outer halo globular cluster system of the Andromeda galaxy (M31). Our basis for this is a set of new spectroscopic observations for 78 clusters lying at projected distances between - kpc from the M31 centre. These are largely drawn from the recent PAndAS globular cluster catalogue; 63 of our targets have no previous velocity data. Via a Bayesian maximum likelihood analysis we find that globular clusters with kpc exhibit coherent rotation around the minor optical axis of M31, in the same direction as more centrally-located globular clusters, but with a smaller amplitude of kms. There is also evidence that the velocity dispersion of the outer halo globular cluster system decreases as a function of projected distance from the M31 centre, and that this relation can be well described by a power law of index . The velocity dispersion profile of the outer halo globular clusters is quite similar to that of the halo stars, at least out to the radius up to which there is available information on the stellar kinematics. We detect and discuss various velocity correlations amongst subgroups of globular clusters that lie on stellar debris streams in the M31 halo. Many of these subgroups are dynamically cold, exhibiting internal velocity dispersions consistent with zero. Simple Monte Carlo experiments imply that such configurations are unlikely to form by chance, adding weight to the notion that a significant fraction of the outer halo globular clusters in M31 have been accreted alongside their parent dwarf galaxies. We also estimate the M31 mass within 200 kpc via the Tracer Mass Estimator (TME), finding . This quantity is subject to additional systematic effects due to various limitations of the data, and assumptions built in into the TME. Finally, we discuss our results in the context of formation scenarios for the M31 halo.
keywords:Local Group — galaxies: individual (M31) — galaxies: kinematics and dynamics — galaxies: halos — globular clusters: general
Even after decades of active research, our understanding of the various processes governing galaxy formation and evolution remains incomplete. In the currently favoured CDM cosmological models and their semi-analytic extensions, the extended dark matter and stellar haloes of galaxies are at least partly formed through hierarchical build-up of smaller fragments, akin to the dwarf galaxies we observe today (e.g. Abadi et al., 2003a, b; Purcell et al., 2007; Guo et al., 2011). Observations of various stellar streams in the halo of our Galaxy (e.g. Belokurov et al., 2006; Grillmair & Dionatos, 2006; Grillmair, 2006; Martin et al., 2014), many of which are thought to be tidally disrupted dwarf galaxies, support this idea. Indeed, observations of the Sagittarius dwarf galaxy, which is currently being accreted onto the Milky Way, provide direct evidence that this process is still on-going (e.g. Ibata et al., 1994; Majewski et al., 2003; Koposov et al., 2012; Slater et al., 2013). However, searching for and studying such features is challenging. Halo streams are typically very faint and therefore difficult to observe even in the Milky Way. Seeking the progenitor systems of tidal streams is also complicated, because satellite galaxies that are losing stellar mass to tides tend to become faint and cold (Peñarrubia et al., 2008), and because their survival time depends strongly on both their (typically unknown) mass distribution (Peñarrubia et al., 2010) and the shape of the dark matter halo of their host (Peñarrubia et al., 2002).
Another way to probe galaxy haloes is through their globular cluster (GC) systems. Due to their high luminosities, GCs are observed much more easily than the underlying stellar field components in these remote parts of galaxies. Various studies have found correlations between GC systems and their host galaxy properties that can shed light on galaxy formation mechanisms (Brodie & Strader, 2006). Since GCs are frequently found at large radii, their kinematics constitute a particularly useful tool. GC motions provide information about the assembly history of the host galaxy, its total mass, the shape of the gravitational potential and the dark matter distribution (e.g. Schuberth et al., 2010, 2012; Strader et al., 2011). Of relevance to the present paper is that GCs also provide an alternative way to look for and study past accretion events, by searching for spatially and dynamically linked GC groups that can serve as tracer populations for their (now disrupted) parent systems.
In the Milky Way, various properties of the halo GCs such as their ages, horizontal branch morphologies, luminosities, sizes and kinematics, are consistent with them having an external origin, in line with expectations from hierarchical formation models (e.g. Searle & Zinn, 1978; Mackey & Gilmore, 2004; Mackey & van den Bergh, 2005; Marín-Franch et al., 2009; Forbes & Bridges, 2010; Keller et al., 2012). Indeed, it has been clearly demonstrated that a number of outer halo GCs are kinematically associated with the Sagittarius stream (e.g., Ibata et al., 1995; Da Costa & Armandroff, 1995). However, despite this, and despite the fact that GCs provided the first clues that the Milky Way halo formed at least to some degree through the coalescence of smaller fragments (Searle & Zinn, 1978), it has proven difficult to establish the existence of dynamically linked Galactic GC groups (e.g., Palma et al., 2002), or locate additional examples where Galactic GCs are clearly kinematically (or even spatially) associated with stellar streams from the disruption of their host dwarfs.
The Andromeda galaxy (M31) provides a unique opportunity for detailed study of galaxy assembly processes. Its close proximity of kpc (Conn et al., 2012) makes M31 the only massive spiral galaxy other than our own in which both star clusters and the diffuse stellar field can be resolved into individual stars. Subtending a large area on the sky, M31 provides a much clearer view of a typical spiral galaxy halo than our own Galaxy, where one must observe vast angular regions and battle with projection effects and hugely variable extinction. M31 hosts a rich GC system, with over 450 confirmed members listed in the Revised Bologna Catalogue111http://www.bo.astro.it/M31/ (Galleti et al., 2004), most of which lie within 30 kpc in projection from the galactic centre.
Kinematics have been of particular interest for studies of the M31 GC system. The first radial velocities for M31 GCs were obtained by van den Bergh (1969). Later Hartwick & Sargent (1974) used the available data to estimate the mass of M31 using its GCs as kinematic tracers. A number of works followed (Huchra et al., 1982; Federici et al., 1990; Huchra et al., 1991; Federici et al., 1993), increasing the GC velocity data set, updating the kinematic mass estimate and providing velocity dispersions for different subsamples. Using higher precision data for over 200 GCs, all having projected radii smaller than 22 kpc, Perrett et al. (2002) showed that unlike in the Milky Way, the GC system of M31 exhibits a strong rotation around the minor optical axis of the galaxy. More recent investigations have further enlarged the total number of radial velocity measurements, and presented updated kinematic analyses and mass estimates (Galleti et al., 2006, 2007; Lee et al., 2008; Caldwell et al., 2011). It is important to emphasise that all these past investigations have focused on GCs at projected distances less than 30 kpc from the M31 centre.
A number of studies have used the available kinematic data to search for and detect possible velocity sub-clustering amongst the M31 GCs, which in turn may reflect past mergers or accretions in the system. The first significant attempt was by Ashman & Bird (1993), who used a technique in which groups of GCs were isolated based on deviations of the global mean velocity and velocity dispersion between each cluster and its nearest companions. This technique yielded a number of groups, but Ashman & Bird (1993) warned that their method may produce false positives if the GC system were to exhibit significant rotation, which was later found to indeed be the case (Perrett et al., 2002). Using an improved and enlarged data set, Perrett et al. (2003) searched for sub-clustering in the inner M31 GC system as evidence of past merger remnants. These authors employed a modified friends-of-friends algorithm which can detect the elongated groups that are expected to be found along tidal debris streams. Perrett et al. (2003) detected 10 unique groups of at least 4 GCs in each. They performed additional tests and found that even though the majority of these might be chance groupings, there was a high probability that at least some might be genuine dynamically linked units. More recently, Perina et al. (2009) attempted to identify clusters in the inner parts of M31 sitting away from the global trend in metallicity with position. They located three such GCs, sitting at similar projected radii and possessing matching velocities quite distinct from the kinematics of the M31 disk, altogether suggestive of being physically part of a coherent structure.
We again emphasise that these previous attempts at finding coherent GC groups in M31 were undertaken for objects lying at projected radii kpc, where such searches are extremely challenging. It is difficult, if not impossible, to trace a single stellar debris stream reliably due to the presence of many intertwined stellar substructures, as well as the comparatively high stellar densities of the M31 spheroid. In addition, the high number density and wide range of GC properties makes it difficult for distinct kinematic groups to be robustly detected. Because the dynamical timescales in the inner parts of M31 are also comparatively short, it is also likely that any accreted objects presently found in these regions are now well mixed with the host population of stars and clusters and no longer retain their initial kinematic relationships.
In recent years, our international collaboration has initiated a sequence of wide-area imaging surveys in order to explore in detail the far outer halo of M31, culminating in the recent Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al., 2009). Our high quality data has enabled us to search for, and on account of their partially resolved nature, unambiguously detect GCs out to a projected radius () of kpc (Huxor et al., 2005, 2008; Huxor et al., 2014), and, in at least one case, a 3D radius of kpc (Martin et al., 2006; Mackey et al., 2010a). In total, our group has discovered 80 GCs in M31 with galactocentric distances larger than 30 kpc in projection. For comparison, prior to our work only 3 globular clusters were known at such large radii. One particularly interesting conclusion drawn from the results of the PAndAS survey is that a large fraction (50-80%) of these remote GCs preferentially lie projected on top of various stellar streams and other tidal debris features uncovered in the M31 halo (Mackey et al., 2010b). Monte Carlo simulations suggest that the probability of such alignment arising by chance is lower than 1%. This finding is a major step forward in understanding how the M31 outer halo GC population formed, and supports the idea first put forward by Côté et al. (2000) that accretion was the dominant mode of formation for this GC system.
In a recent Letter (Veljanoski et al., 2013), we presented an initial survey of the kinematics of the M31 outer halo GC system. We showed for the first time that the outer halo GC system of M31 appears to be rotating and that there is a hint of a decreasing velocity dispersion amongst the population as a function of increasing projected radius.
In the present contribution, we significantly increase the GC sample for which radial velocities have been measured, and present the first detailed kinematic analysis of the outer halo GC system of M31. This paper is structured in the following manner. Section 2 contains a complete description of the data, the data reduction, and our methodology for determining radial velocities. In Section 3 we use a Bayesian framework to derive the global kinematic properties of the M31 outer halo GCs, while in Section 4 we focus on various velocity correlations and sub-clustering observed for GCs that lie along particular stellar debris features. Finally, we discuss the implications of our results in Section 5, followed by the summary in Section 6.
2 Observations and data reduction
2.1 The sample
Our spectroscopy campaign involves eight separate observing runs conducted with three different facilities: the ISIS spectrograph mounted on the 4.2m William Herschel Telescope (WHT), the RC spectrograph on the 4.0m Mayall Telescope at the Kitt Peak National Observatory (KPNO), and the GMOS instrument installed on the 8.1m Gemini-North telescope. The targets were selected from a catalogue of outer halo M31 GCs comprised mainly of objects discovered in the PAndAS survey (Huxor et al., 2014, Paper I in this series), but also including clusters found in previous searches by our group (Martin et al., 2006; Huxor et al., 2008), as well as objects listed in the Revised Bologna Catalogue (RBC, Galleti et al., 2004). In total our list of possible targets contained 83 GCs situated at larger than 30 kpc, plus those at smaller radii presented in Paper I. Table 1 shows the log of observations for all eight observing runs. Altogether we made 90 separate GC observations, acquiring spectra for 78 different clusters, of which 63 had no previous spectroscopic information. Repeated observations of some GCs were made primarily to facilitate consistency checks but also to supplement lower quality data in a few instances.
Throughout the campaign, priority was given to clusters lying on top of stellar substructures, and to those having larger . Our final observed sample consists of GCs with between 18 and 141 kpc. Most, however, lie beyond 30 kpc in projection from the centre of M31 – a region which, throughout the remainder of this paper, we will refer to as the “outer halo”. This radius corresponds to the clear break in the GC radial number density profile observed by Huxor et al. (2011, see also Mackey et al. 2014, in prep.). In total we acquired spectra for 71 clusters with kpc, corresponding to 85.5% of the known globular clusters in the M31 outer halo. Of these, there are 20 in our sample beyond 80 kpc including 10 beyond 100 kpc. The full radial distribution of our observed clusters is shown in Figure 1.
In the following subsections we describe the data and the data reduction process. Even though the reduction procedure is standard and similar in the case of the 4m and 8m class telescopes, we discuss it separately in order to highlight any differences.
|Instrument||Date of obs.||Program Number||PI||Obs. mode||No. GCs|
|WHT/ISIS||29/09-02/10 2005||.||A. P. Huxor||Visitor||19|
|WHT/ISIS||16/08-18/08 2009||.||A. P. Huxor||Visitor||12|
|WHT/ISIS||09/09-11/09 2010||.||A. P. Huxor||Visitor||13|
|KPNO/RC||13/08-17/08 2009||.||A. M. N. Ferguson||Visitor||17|
|Gemini/GMOS-N||20/07-02/09 2010||GN-2010B-Q-19||A. D. Mackey||Service||4|
|Gemini/GMOS-N||02/08-05/09 2011||GN-2011B-Q-61||A. D. Mackey||Service||11|
|Gemini/GMOS-N||29/07-13/09 2012||GN-2012B-Q-77||A. D. Mackey||Service||7|
|Gemini/GMOS-N||02/08-31/08 2013||GN-2013B-Q-66||A. D. Mackey||Service||7|
2.2 WHT and KPNO data
We used the ISIS instrument mounted on the 4.2m WHT for three observing runs, performing longslit spectroscopic observations of 41 different GCs in our sample. ISIS has two detectors (“arms”), that independently sample two separate wavelength ranges, a bluer and a redder one. In all runs we set the slit width at 1.5 - 2. For the blue arm we used the R600B grating to cover the wavelength range between 3500 and 5100 Å, and the EEV12 detector with a dispersion of 0.45 Å pixel. The spectral resolution was R . For the red arm we used the R600B grating covering the wavelength range between 7500 - 9200 Å, and the REDPLUS camera, achieving a dispersion of 0.49 Å pixel. The spectral resolution was R . The only exception to this set-up was for the observing run conducted in 2005, when only the blue arm of ISIS was used. We observed each GC as a series of short exposures, with the total integration time varying between 600 and 7200 seconds depending on the cluster brightness. The data are unbinned in both the spatial and wavelength direction. The typical signal-to-noise of the spectra is per Å, while reaching per Å for the brightest targets.
We used the RC spectrograph mounted on the KPNO-4m telescope in single slit mode to obtain spectra of 17 GCs during a single observing run in 2009. For this purpose, we used a slit width of 2 and the T2KB detector, along with the KPC007 grating to select the wavelength range 3500 - 6500 Å with a dispersion of 0.139 Å pixel and a spectral resolution of R . A similar observing strategy to that which we used for the WHT observations was adopted. Each cluster was observed with multiple exposures, and the total integration time ranged between 600 and 6400 seconds depending on the brightness of the target. There is no binning of the data in either the wavelength or the spatial directions. The typical signal to noise of the spectra is 25 - 50 per Å.
The data obtained with the WHT and the KPNO telescopes were reduced using standard iraf222 iraf is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. procedures. The basic reduction of the spectra (bias and overscan subtraction, flat-fielding, illumination correction) was done with dedicated standard tasks, which are part of the ccdred package. The apall task in the kpnoslit package was used to extract one dimensional spectra from the two-dimensional frames. The extraction apertures had radii of 2 - 2.5. The same task was also used interactively to select background sky regions and to find the trace. The sky in the selected regions of the target spectra was then fit with a 2 order Chebyshev polynomial and subtracted. The spectra were traced using a 3 order cubic spline function, and were extracted using the optimal variance weighting option in apall. An advantage of apall is that it also produces an error spectrum based on the Poisson noise of the spectrum that is being extracted.
Wavelength calibration was based on Cu-Ne-Ar and He-Ne-Ar lamps for WHT and KPNO spectra respectively. Comparison “arcs” were obtained before and after each program target exposure. The arc spectra were extracted using the same apall parameters as the target GCs they were used to calibrate. The identify task was used to identify RC, ISIS blue and ISIS red lines in the arc spectra, and the dispersion solution was fit with a 3 order cubic-spline function. The RMS residuals of the fits were 0.08 0.01 Å 0.05 0.01 Å and 0.02 0.01 Å for the data obtained with the RC, ISIS blue arm and ISIS red arm instruments respectively. Since two wavelength solutions were found for each target from the ‘before’ and ‘after’ arcs, they were averaged and assigned to the GC spectrum via the dispcor task. To check whether the wavelength calibration is reliable, we measured the positions of sky lines in separately extracted sky spectra. We found that the wavelength calibration is accurate to 0.08 Å with no systematic shifts for all data observed with the 4m class telescopes.
The multiple exposures of each target GC from a given observing run were stacked in the following manner. First, all exposures were shifted into the heliocentric frame and interpolated onto a common wavelength scale. They were then combined together as described by Equation 1:
where S represents a spectrum, while is the corresponding error spectrum. The index corresponds to a particular pixel in an exposure . Finally the spectra were continuum subtracted, for the purpose of measuring radial velocities. Several examples of fully reduced spectra are shown in Figure 2. The displayed spectra are continuum normalized rather than continuum subtracted in order to preserve the relative strengths of the absorption features for better visualization.
2.3 Gemini data
To observe the fainter as well as the more diffuse and extended GCs, which typically have lower surface brightness than classical compact globular clusters, we used the GMOS instrument mounted on the Gemini-North telescope. Spectra were taken for a total of 29 objects over the course of four separate observing runs executed in service mode between 2010 and 2013. The observations were conducted using a longslit mask with a slit width of 0.75. The grating of choice was R831, which was used to cover the wavelength range between 7450 and 9500 Å. To account for the gaps between the chips of the GMOS detector, two or three sets of three exposures were taken with slightly different grating angles and hence slightly different central wavelengths. For each science exposure in each set, we dithered the telescope by 15 along the spatial direction of the slit. This helps to minimise any effects coming from imperfections in the detector as well as any systematic effects of the background sky subtraction. The typical total integration time was 5700 seconds. The data were binned by a factor of two in the wavelength direction, obtaining a resolution of 0.68 Å pixel, as well as in the spatial direction to achieve a spatial resolution of 0.147 per pixel. The spectral resolution achieved with this setup was R . The average signal to noise of the data is per Å.
The data reduction was undertaken using iraf, employing tasks from the dedicated gemini/gmos package provided and maintained by the Gemini staff. The reduction was carried out separately for data sets observed with different central wavelengths. A master bias frame was created with the gbias task from raw bias frames acquired near to the time of the program observations. The standard overscan and bias subtraction, flat-fielding and mosaicing of the three chips of the detector into a single frame was done with the gsreduce task. Unlike the data taken with the 4m class telescopes, here we wavelength calibrated the two dimensional frames before one-dimensional spectra were extracted. The wavelength calibration is based on Cu-Ar arcs, and such frames were taken before and after each set of program target exposures, with the central wavelength of the arcs matching the central wavelength of the observed GCs in a specific set. The wavelength task was employed to identify lines in the arc spectrum and to fit a dispersion solution using a fourth order Chebyshev polynomial. The typical RMS of the fit is 0.02 Å. The gstransform task was employed to assign a wavelength solution to each GC frame.
|HD 4388||00 46 27.0||+30 57 05.6||K3III||-27.5 0.3||4||KPNO/RC||2009||Udry et al. (1999)|
|HD 12029||01 58 41.9||+29 22 47.7||K2III||38.5 0.3||4||KPNO/RC||2009||Udry et al. (1999)|
|HD 145001||16 08 04.5||+17 02 49.1||G8III||-10.3 0.3||4||KPNO/RC||2009||Udry et al. (1999)|
|HD 149803||16 35 54.3||+29 44 43.3||F7V||-7.5 0.7||1||KPNO/RC||2009||Udry et al. (1999)|
|HD 154417||17 05 16.8||+00 42 09.2||F9V||-18.6 0.3||7||KPNO/RC||2009||Udry et al. (1999)|
|HD 171391||18 35 02.4||-10 58 37.9||G8III||7.4 0.2||5||KPNO/RC||2009||Udry et al. (1999)|
|G1||00 32 46.8||+39 34 42.6||-332 3||1||KPNO/RC||2009||Galleti et al. (2004)|
|MGC1||00 50 42.4||+32 54 58.7||-355 2||2||WHT/ISIS||2009||Alves-Brito et al. (2009)|
One dimensional spectra were extracted with the standard iraf apall package. As the Gemini observations primarily targeted faint and diffuse clusters, often multiple apertures were used to extract the light coming from individual bright stars within the cluster. Typical effective aperture radii range from 0.7 to 2.5. The internal velocity dispersion of extended and low luminosity clusters is smaller than the measured radial velocity uncertainty of each star in such a cluster, which makes this approach an appropriate one. The apall task was also used to subtract the background sky, and to find the trace, in a similar way as for the 4m data. The extraction was done with the variance weighting option on. Finally, all exposures for a given GC observation were stacked together in the same manner as for the 4m telescope data. If multiple extraction apertures were used, they were also stacked together following the prescription in Equation 1, producing a final one dimensional science spectrum. These science spectra were then continuum subtracted as required for measuring radial velocities. An example of a representative, fully reduced spectrum is shown in Figure 2.
2.4 Radial velocity measurements
For the purpose of determining the radial velocities of the GCs, throughout our observing campaign we also performed multiple observations of 6 different radial velocity standard stars. The standard stars were chosen to have a stable and accurately known radial velocity, to be of a certain spectral type so their spectra would be similar to the GC spectra, and to be sufficiently bright so a high S/N spectrum could be obtained with a very short exposure.
In addition we also used two M31 GCs, G1 and MGC1, as radial velocity templates. These clusters have high precision radial velocity information obtained from high resolution spectra (Galleti et al., 2004; Alves-Brito et al., 2009), comprise some of the brightest GCs in the outskirts of M31, and possess metallicities spanning the expected range for outer halo GCs FeH. These properties allow for high S/N spectra to be observed with relatively short exposures and makes the two GCs appropriate radial velocity templates. Table 2 displays the relevant data regarding the radial velocity standard stars and template GCs.
Heliocentric radial velocities were determined via a customized routine that performs a minimisation between the target and template spectra. First, the template is adjusted to the wavelength scale of the target spectrum. The template spectrum is then Doppler shifted by an input velocity, which is systematically varied between -1000 to 500 km s in increments of 1 km s in the heliocentric frame. The chosen velocity search range is large enough to comfortably encompass the expected velocities of all GCs that belong to the M31 system. The match between the target and template spectra is then calculated via:
where is the pixel index, d is the spectrum of the target GC (‘data’ in the statistical sense), and M is the template spectrum (the ‘model’ against which the data are tested). The uncertainties in the target and template spectra are and respectively. The model M is a function of two free parameters. The first is the input velocity . The second parameter, , is due to the different width of the absorption lines in the target cluster and the template star spectrum, caused by the internal velocity dispersion of the stars that comprise a certain GC. However, as the resolution of the spectrographs we have employed is not sufficient to probe the internal velocity dispersions of the GCs, this parameter can be ignored. The parameter accounts for the flux difference between the target and the template spectra. It is not independent, and can be calculated via:
where the symbols are as in Equation 2. The input velocity corresponding to the minimum of the function is the measured velocity of the GC. In order to remove the large telluric features, any regions of higher sky subtraction residuals, and the edges of the spectra where the S/N is low, we selected certain wavelength ranges over which the function was calculated. For the data observed with GMOS-N and the red arm of ISIS, the window was selected just around the Caii triplet (CaT) lines with range of 8400 - 8750 Å. For the KPNO data the window was in the range of 3831 - 6000 Å, and for the data observed with the blue arm of ISIS the corresponding window was 3900 - 4900 Å.
This technique yielded smaller velocity uncertainties by 23% on average compared to the results coming from the standard cross correlation method. Usage of the uncertainties in both the template and target spectra helps to eliminate spurious peaks in the function that might arise due to imperfectly subtracted sky lines, which are especially strong for the faint GCs near the CaT.
|ID||Alternative ID||Position (J2000.0)||Reference|
|G1||Mayall-II||00 32 46.5||+39 34 40||-332 3||average, see Galleti et al. (2004)|
|G2||Mayall-III||00 33 33.7||+39 31 18||-313 17||average, see Galleti et al. (2006)|
|B514||MCGC4||00 31 09.8||+37 54 00||-458 23||Galleti et al. (2007)|
|B517||.||00 59 59.9||+41 54 06||-272 54||Galleti et al. (2009)|
|G339||BA30||00 47 50.2||+43 09 16||33 30||Federici et al. (1993)|
|EXT8||.||00 53 14.5||+41 33 24||-152 30||Federici et al. (1993)|
|H1||MCGC1/B520||00 26 47.7||+39 44 46||-219 15||Galleti et al. (2007)|
|H10||MCGC5||00 35 59.7||+35 41 03||-358 2||Alves-Brito et al. (2009)|
|H14||MCGC7||00 38 49.4||+42 22 47||-248 24||Caldwell et al. (2011)|
|H23||MCGC8||00 54 24.9||+39 42 55||-381 15||Galleti et al, in preparation|
|H24||MCGC9||00 55 43.9||+42 46 15||-147 20||Galleti et al, in preparation|
|H27||MCGC10||01 07 26.3||+35 46 48||-291 2||Alves-Brito et al. (2009)|
|MGC1||.||00 50 42.4||+32 54 58||-355 2||Alves-Brito et al. (2009)|
|PAndAS-07||PA-7||00 10 51.3||+39 35 58||-433 8||Mackey et al. (2013)|
|PAndAS-08||PA-8||00 12 52.4||+38 17 47||-411 4||Mackey et al. (2013)|
For GCs observed in a single observing run, the final adopted radial velocity and its corresponding uncertainty are given by the mean and standard deviation from all individually obtained velocities resulting from the minimisation between the spectrum of that cluster and all available template spectra, respectively. Regarding the GCs which were observed with WHT in 2009 and 2010, two independent radial velocities were measured from the blue and the red arms of ISIS. It is important to note that these measurements are consistent with each other, and there is no systematic offset between them. A comparison between these two independent sets of measurements is shown in Figure 3. The mean offset is km s, which is much smaller than the RMS deviation, found to be 20 km s. Because of the excellent agreement between them, to obtain final velocities for objects in the 2009 and 2010 WHT runs we simply computed the error-weighted average of the blue and red arm measurements.
There are 7 clusters that were repeatedly observed in different observing runs. For the radial velocity of those objects we adopt the error weighted mean of the velocities measured in each run individually, including the measurements conducted with the two ISIS arms if available.
A total of 15 objects in our sample, listed in Table 3, have measured radial velocities in the literature. We show a comparison between these values and the velocity measurements from our present study in Figure 4. The mean offset is km s, while the RMS deviation is found to be 39 km s. Removing the two velocity measurements coming from Federici et al. (1993), which is the set of points that are the most deviant, results in a mean offset of km s with a RMS scatter of 18 km s. The origin of the discrepancy between our velocity measurements and those from Federici et al. (1993) is unclear, but may stem from the very different measurement and calibration techniques employed in the two different studies. Apart from this, we find excellent agreement between the velocities derived in this work and those present in the literature, with no signficant systematic offset between our radial velocity measurements and those collated in the RBC.
2.5 Corrections for perspective
Because our GC sample is spread over a large area of sky spanning , we converted our measured radial velocities from the heliocentric to the Galactocentric frame in order to remove any effects the solar motion could have on the kinematic analysis. The conversion was computed via the relation found in Courteau & van den Bergh (1999), with updated values of the solar motion from McMillan (2011) and Schönrich et al. (2010):
where and are the Galactic latitude and longitude.
The wide angular span of our GC sample on the sky introduces additional factors that must be considered. As per van der Marel & Guhathakurta (2008), the observed (Galactocentric) line-of-sight velocity for a target that is part of the extended M31 system, but separated from its centre by an angle on the sky, can be decomposed as:
Here, we take M31 to have a systemic radial velocity (measured along the line-of-sight to its centre) of , and a systemic transverse velocity in a direction on the sky given by the position angle . The position angle of the target with respect to the centre of M31 is , while is its peculiar line-of-sight velocity.
The first two terms in Equation 5 tell us that with increasing separation , a decreasing fraction of the systemic M31 radial velocity is observed along the line-of-sight to the target, but an increasing fraction of the transverse motion is carried on this vector. This induces the appearance of a solid body rotation for targets at wide separations from the centre of the system, around an axis sitting away from the direction of the transverse velocity on the sky.
One of the main things we wish to test in this paper is whether the outer halo GC system of M31 exhibits coherent rotation as suggested by Veljanoski et al. (2013). It is therefore important to consider whether we need to make a correction for the “perspective rotation” described above. The most precise measurement of the M31 transverse velocity to date comes from van der Marel et al. (2012) who found km s with respect to the Milky Way, at a position angle east of north. Their confidence region is km s, consistent with M31 being on a completely radial orbit towards our Galaxy. The small transverse velocity of M31 means that the induced perspective rotation for our GC sample is negligible – at most a few km s even for the most remote objects (which have ). This is smaller than our typical measurement uncertainties. In principle we could, for completeness, still use Equation 5 to correct for the rotation; however the formal uncertainties on the individual components of the van der Marel et al. (2012) transverse velocity (i.e., the components in the north and west directions on the sky) are km s each. Hence making the correction would introduce significantly larger random uncertainties into our final velocity measurements than ignoring this effect entirely.
We use Equation 5 with the second term set to zero to obtain the peculiar line-of-sight velocity of each GC in our sample. That is, we remove the component due to the radial systemic motion of M31 by solving for . In this study we adopt a heliocentric velocity of -301 1 km s for M31 (van der Marel & Guhathakurta, 2008), which corresponds to a Galactocentric radial velocity -109 4 km s (see also van der Marel et al., 2012).
One remaining effect to consider is that each of our final corrected velocities lies along a slightly different vector due to the different lines-of-sight to the individual GCs in our sample. In principle we ought to correct these to lie parallel – adopting the direction of the line-of-sight to the centre of the galaxy might be a logical choice – before assessing, for example, the how the one-dimensional velocity dispersion of the system varies with projected radius. However to make this correction for a given GC requires knowledge of its peculiar proper motion, because, in analogy with Equation 5, a small component of this transverse velocity is carried onto the new vector. In the absence of this information we choose to leave our measurements unaltered; in any case the expected magnitude of the corrections is, for most targets, smaller than the random uncertainties on our velocities.
Table 4 lists the radial velocity measurements for all GCs in our sample. For each object the raw heliocentric velocity is reported, followed by the Galactocentric velocity from Equation 4, and the peculiar line-of-sight velocity in the M31 frame obtained by solving Equation 5 as described towards the end of Section 2.5. We hereafter refer to this latter velocity as .
The results of our measurements are also shown in Figure 5, where the observed GCs from this study are overlaid as coloured points on top of the most recent PAndAS metal-poor ([Fe/H] ) red giant branch stellar density map. The colour of each marked GC corresponds to .
As described in the next Section, we will concentrate our kinematic analysis on the region outside kpc. For completeness we searched through the RBC and found that there is only one cluster with kpc that we have not observed but which has a radial velocity in the literature. To improve statistics, we add this cluster – dubbed HEC12, or alternatively EC4 (Collins et al., 2009) – to our sample, and list its relevant data in Table 4. Therefore, our final sample of outer halo GCs ( kpc) contains 72 objects in total.
|Cluster ID||Position (J2000.0)||PA||Prominent||Instrument|
|RA||Dec||[kpc]||[deg]||[deg]||[km s]||[km s]||[km s]||[km s]||substructure|
|B514||00 31 09.8||+37 54 00||55.2||214.4||4.04||-471 8||-279 8||-169 8||-84||.||W|
|B517||00 59 59.9||+41 54 06||44.9||77.5||3.40||-277 13||-93 13||16 13||-57||Stream C/D (overlap)||K|
|EXT8||00 53 14.5||+41 33 24||27.2||80.8||1.99||-194 6||-7 7||102 7||.||W|
|G001||00 32 46.5||+39 34 40||34.7||229.1||2.54||-335 5||-141 6||-31 6||58||Association 2||W, K, G|
|G002||00 33 33.7||+39 31 18||33.8||225.7||2.47||-352 19||-158 19||-49 19||77||Association 2||G|
|G268||00 44 10.0||+42 46 57||21.0||9.8||1.54||-277 8||-84 8||25 8||.||W|
|G339||00 47 50.2||+43 09 16||28.8||26.2||2.11||-97 6||95 6||204 6||.||W|
|H1||00 26 47.7||+39 44 46||46.3||244.6||3.39||-245 7||-48 7||61 7||143||.||W|
|H2||00 28 03.2||+40 02 55||41.6||247.5||3.04||-519 16||-322 16||-212 16||-133||Association 2||W|
|H3||00 29 30.1||+41 50 31||34.7||284.1||2.54||-86 9||113 9||221 9||266||.||W|
|H4||00 29 44.9||+41 13 09||33.4||269.9||2.44||-368 8||-170 8||-61 8||1||.||W|
|H5||00 30 27.2||+41 36 19||31.8||279.3||2.33||-392 12||-194 12||-85 12||-34||.||W|
|H7||00 31 54.5||+40 06 47||32.2||241.5||2.24||-426 23||-231 23||-121 23||-38||Association 2||K|
|H8||00 34 15.4||+39 52 53||29.1||229.9||2.13||-463 3||-269 4||-160 4||-73||Association 2||G|
|H9||00 34 17.2||+37 30 43||56.1||204.2||4.10||-374 5||-184 6||-74 6||7||.||W|
|H10||00 35 59.7||+35 41 03||78.4||193.8||5.47||-352 9||-165 9||-56 9||12||.||W|
|H11||00 37 28.0||+44 11 26||42.1||342.1||3.08||-213 7||-15 7||93 7||54||.||W|
|H12||00 38 03.8||+37 44 00||49.9||194.7||3.65||-396 10||-207 10||-98 10||-23||.||W|
|H14||00 38 49.4||+42 22 47||18.2||327.0||1.32||-271 15||-76 15||33 15||.||K|
|H15||00 40 13.2||+35 52 36||74.0||185.4||5.42||-367 10||-182 10||-73 10||-6||.||W|
|H17||00 42 23.6||+37 14 34||55.0||181.0||3.97||-246 16||-60 16||48 16||122||.||K|
|H18||00 43 36.0||+44 58 59||50.8||2.4||3.72||-206 21||-10 21||99 21||35||.||W|
|H19||00 44 14.8||+38 25 42||39.0||174.1||2.85||-272 18||-85 18||24 18||79||.||W|
|H22||00 49 44.6||+38 18 37||44.4||155.0||3.25||-311 6||-127 6||-17 6||12||.||W|
|H23||00 54 24.9||+39 42 55||37.0||124.0||2.71||-377 11||-193 11||-84 11||-100||Stream D||W|
|H24||00 55 43.9||+42 46 15||38.8||57.0||2.96||-121 15||66 15||175 15||91||Stream C/D (overlap)||K|
|H25||00 59 34.5||+44 05 38||57.2||46.2||4.19||-204 14||-16 14||93 14||6||.||W|
|H26||00 59 27.4||+37 41 30||65.8||136.6||4.81||-411 7||-233 7||-124 7||-121||Stream C||W|
|H27||01 07 26.3||+35 46 48||99.9||136.7||7.31||-291 5||-121 6||-12 6||-9||.||W|
|HEC1||00 25 33.8||+40 43 38||44.9||261.9||3.28||-233 9||-34 9||74 9||145||.||K, G|
|HEC2||00 28 31.5||+37 31 23||63.5||217.4||4.64||-341 9||-148 9||-39 9||48||.||G|
|HEC6||00 38 35.4||+44 16 51||42.5||346.2||3.11||-132 12||65 12||174 12||130||.||G|
|HEC10||00 54 36.4||+44 58 44||58.7||29.3||4.30||-98 5||93 6||202 6||119||.||G|
|HEC11||00 55 17.4||+38 51 01||46.6||134.2||3.41||-215 5||-33 6||76 6||75||Stream D||G|
|HEC12||00 58 15.4||+38 03 01||60.0||135.9||4.39||-288 2||-109 4||0 4||1||Stream C||.|
|HEC13||00 58 17.1||+37 13 49||68.8||142.1||5.04||-366 5||-188 6||-79 6||-68||Stream C||G|
|MGC1||00 50 42.4||+32 54 58||116.2||168.6||8.50||-355 7||-181 7||-73 7||-23||.||W, G|
|PAndAS-01||23 57 12.0||+43 33 08||118.9||289.0||8.62||-333 21||-119 21||-11 21||28||.||K|
|PAndAS-02||23 57 55.6||+41 46 49||114.7||277.2||8.40||-266 4||-54 4||53 4||108||.||G|
|PAndAS-04||00 04 42.9||+47 21 42||124.6||315.1||9.12||-397 7||-183 7||-75 7||-74||NW stream||W|
|PAndAS-05||00 05 24.1||+43 55 35||100.6||294.3||7.36||-183 7||28 7||136 7||168||.||G|
|PAndAS-06||00 06 11.9||+41 41 20||93.7||276.5||6.75||-327 15||-119 15||-10 15||45||.||K|
|PAndAS-07||00 10 51.3||+39 35 58||86.0||257.2||6.29||-452 18||-248 18||-139 18||-46||SW cloud||G|
|PAndAS-08||00 12 52.4||+38 17 47||88.3||245.0||6.46||-416 8||-215 8||-106 8||-19||SW cloud||G|
|PAndAS-09||00 12 54.6||+45 05 55||90.8||307.7||6.60||-444 21||-235 21||-126 21||-115||NW stream||K|
|PAndAS-10||00 13 38.6||+45 11 11||90.0||308.9||6.59||-435 10||-226 10||-117 10||-108||NW stream||G|
|PAndAS-11||00 14 55.6||+44 37 16||83.2||305.7||6.09||-447 13||-239 13||-130 13||-116||NW stream||W|
|PAndAS-12||00 17 40.0||+43 18 39||69.2||295.9||5.06||-472 5||-267 5||-157 5||-129||NW stream||G|
|PAndAS-13||00 17 42.7||+43 04 31||68.0||293.4||4.90||-570 45||-365 45||-256 45||-224||NW stream||K|
|PAndAS-14||00 20 33.8||+36 39 34||86.2||224.9||6.31||-363 9||-167 9||-58 9||29||SW cloud||W|
|PAndAS-15||00 22 44.0||+41 56 14||51.9||281.8||3.80||-385 6||-183 6||-74 6||-26||NW stream||G|
|PAndAS-16||00 24 59.9||+39 42 13||50.8||246.6||3.60||-490 15||-292 15||-183 15||-102||.||K|
|PAndAS-17||00 26 52.2||+38 44 58||53.9||231.6||3.83||-279 15||-84 15||25 15||112||.||K|
|PAndAS-18||00 28 23.2||+39 55 04||41.6||244.8||3.08||-551 18||-354 18||-245 18||-163||Association 2||G|
|PAndAS-19||00 30 12.2||+39 50 59||37.9||240.2||2.77||-544 6||-348 6||-239 6||-155||Association 2||G|
|PAndAS-21||00 31 27.5||+39 32 21||37.7||232.1||2.76||-600 7||-405 7||-296 7||-210||Association 2||W|
|PAndAS-22||00 32 08.3||+40 37 31||28.7||253.0||2.10||-437 1||-241 3||-132 3||-55||Association 2||G|
|PAndAS-23||00 33 14.1||+39 35 15||33.7||227.9||2.47||-476 5||-282 6||-172 6||-86||Association 2||G|
|PAndAS-27||00 35 13.5||+45 10 37||56.6||341.3||4.14||-46 8||154 8||262 8||225||.||W|
|PAndAS-36||00 44 45.5||+43 26 34||30.1||9.6||2.21||-399 7||-205 7||-96 7||-167||.||W|
|PAndAS-37||00 48 26.5||+37 55 42||48.1||161.3||3.50||-404 15||-220 15||-111 15||-72||.||K|
|PAndAS-41||00 53 39.5||+42 35 14||33.1||56.1||2.42||-94 8||94 8||203 8||118||Stream C/D (overlap)||W|
|PAndAS-42||00 56 38.0||+39 40 25||42.2||120.0||3.09||-176 4||7 5||115 5||93||Stream D||G|
|PAndAS-43||00 56 38.8||+42 27 17||38.9||64.2||2.85||-135 6||52 7||160 7||79||Stream C/D (overlap)||G|
|PAndAS-44||00 57 55.8||+41 42 57||39.4||79.8||2.99||-349 11||-164 11||-54 11||-126||Stream C/D (overlap)||K|
|PAndAS-45||00 58 37.9||+41 57 11||41.7||75.7||3.05||-135 16||50 16||159 16||85||Stream C/D (overlap)||G|
|PAndAS-46||00 58 56.3||+42 27 38||44.3||67.1||3.36||-132 16||54 16||162 16||82||Stream C/D (overlap)||K|
|PAndAS-47||00 59 04.7||+42 22 35||44.3||68.7||3.35||-359 16||-174 16||-64 16||-144||Stream C/D (overlap)||K|
|PAndAS-48||00 59 28.2||+31 29 10||141.3||159.7||10.34||-250 5||-83 6||25 6||62||.||G|
|PAndAS-49||01 00 50.0||+42 18 13||48.2||71.5||3.53||-240 7||-55 7||53 7||-24||Stream C/D (overlap)||G|
|PAndAS-50||01 01 50.6||+48 18 19||106.7||24.1||7.81||-323 7||-131 7||-22 7||103||.||G|
|PAndAS-51||01 02 06.6||+42 48 06||53.4||65.3||3.91||-226 5||-41 6||67 6||-14||.||G|
|PAndAS-52||01 12 47.0||+42 25 24||78.1||75.9||5.71||-297 9||-118 9||-9 9||-84||.||W|
|PAndAS-53||01 17 58.4||+39 14 53||95.9||103.9||7.01||-253 10||-82 10||26 10||-18||.||W|
|Cluster ID||Position (J2000.0)||PA||Prominent||Instrument|
|RA||Dec||[kpc]||[deg]||[deg]||[km s]||[km s]||[km s]||[km s]||substructure|
|PAndAS-54||01 18 00.1||+39 16 59||95.8||103.6||7.01||-336 8||-165 8||-56 8||-101||.||bf|
|PAndAS-56||01 23 03.5||+41 55 11||103.3||81.7||7.56||-239 8||-66 8||42 8||-26||.||bf|
|PAndAS-57||01 27 47.5||+40 40 47||116.4||90.3||8.52||-186 6||-18 7||90 7||30||Eastern Cloud||bf|
|PAndAS-58||01 29 02.1||+40 47 08||119.4||89.4||8.74||-167 10||1 10||109 10||48||Eastern Cloud||bf|
|SK255B||00 49 03.0||+41 54 57||18.4||60.8||1.34||-191 10||-1 10||107 10||.||bf|
3 Kinematic Analysis
3.1 The tools: Bayesian inference
Bayesian inference is a powerful statistical tool for estimating a set of parameters in a model , as well as discriminating between different models. Given some data and certain prior information , the mathematical form of Bayes theorem is:
where is the posterior probability distribution function (pdf) of the model parameters, is the prior, and is the likelihood function. The expression is the Bayesian evidence which is the average value of the likelihood weighted by the prior and integrated over the entire parameter space. It is mathematically expressed as:
When more of a model’s parameter space has high likelihood values, the evidence is large. However, the evidence is very small for models with large parameter spaces having low likelihood values, even if the likelihood function itself is highly peaked. This Bayesian quantity is key when one wants to discriminate between two different models, and . The typical question that one needs to answer in this case is which model is a better fit to the data. This can be done via the Bayes factor , which is essentially the ratio between the evidence of the models that are considered, and formally defined as:
The Bayes factor is a powerful tool for model selection, especially since it does not depend on any single set of parameters as the integrations are over the entire parameter space in each model. This allows for significantly different models to be compared. In addition, the Bayesian model comparison implicitly guards against overfitting (Kass & Raftery, 1995).
3.2 Kinematic models
One of the main goals of this paper is to constrain the overall kinematic properties of the M31 outer halo GC population. Working in the Bayesian framework provides the ability to discriminate between different kinematic models, while simultaneously deriving probability distribution functions for the free parameters in each model.
We construct two kinematic models, and . Model comprises two components: an overall rotation of the M31 outer halo GC system, and the velocity dispersion of the GCs. Model contains only the velocity dispersion of the GC population. By considering both a rotating and a non rotating model, we can quantify the statistical significance of any detected rotation of the M31 outer halo GCs.
The rotation component in is modelled as prescribed in Côté et al. (2001):
where is the rotational velocity of the GC system at position angle , measured east of north, and is the position angle of the rotation axis of the GC system. The rotation amplitude is labelled as , while denotes the systemic radial motion of the GC system333Note that in practice we set this term to be zero, having already removed the fixed systemic motion of M31 from our GC velocities.. As detailed in Côté et al. (2001), this method assumes that the GC system being investigated is spherically symmetric and that the rotation axis is perpendicular to the line of sight – i.e., it lies in the plane of the sky. The model also assumes that the three-dimensional angular velocity is a function of radial distance only (constant on a sphere). Together these assumptions imply that the projected angular velocity is a function of projected radius only, justifying the use of a sinusoid to describe the rotation of the system.
The velocity dispersion is assumed to have a Gaussian form and to decrease as a function of projected radius from the M31 centre in a power-law manner. The observed dispersion () is comprised of two components – the intrinsic dispersion of the GC system, and the effect of the measurement uncertainties in the GC radial velocities. This is mathematically described in Equation 10, where is the aggregate uncertainty of the GC velocities, is the projected radius, and is the power-law index describing how the velocity dispersion changes as a function of :
The scale radius is fixed at 30 kpc, as this is the point at which the M31 halo begins to dominate (c.f. Geehan et al., 2006); is the corresponding intrinsic velocity dispersion at kpc.
|Strength of evidence|
|Negative (supports )|
|0 to 0.5||1 to 3.2||Barely worth mentioning|
|0.5 to 1||3.2 to 10||Positive|
|1 to 1.5||10 to 32||Strong|
|1.5 to 2||32 to 100||Very strong|
Similarly, the model which does not contain an overall rotation component is simply constructed as:
Note that model is clearly a member of the family of models – it is the special case where the amplitude of rotation is zero. In principle, therefore, we could assess the likelihood of this model relative to the favoured model in the family simply by considering the marginalised probability distributions for the latter. However, for clarity we prefer to make an explicit comparison between the two models and using the Bayesian evidence.
Having defined our models, the likelihood function for each of them is:
in which are the observed properties of the GCs, and are the free parameters of the models we are trying to determine. The index loops over each individual data point. In all our models we assume flat priors. Previous studies (e.g, Lee et al., 2008; Veljanoski et al., 2013) have found the velocity dispersion and the overall rotation of the M31 GCs in both the halo and the disk to be roughly equivalent in magnitude. Thus, it is important to note that in our proposed model we are attempting to describe the rotation and velocity dispersion simultaneously rather than separately as has been the case in the majority of past studies. This is done in order to avoid any possible bias that can arise from measuring these quantities in succession, because in such cases the latter measurement depends on the first.
As a reminder, our input sample of GCs is defined by the 72 objects in Table 4 with kpc (our 71 observed targets plus HEC12). The vast majority of velocity measurements for this sample come from our observations as defined in Section 2, except in a handful of cases where previous measurements from the literature are more precise. The spatial coverage of the input sample is high but non-uniform, being slightly biased towards GCs that (i) project onto visible substructures in the field halo, and (ii) lie at larger .
Calculating the likelihood function, the evidence and the posterior probability distributions as described in Equations 6, 7, 13 and 14 can be numerically challenging. Various Monte Carlo algorithms (e.g. Lewis & Bridle, 2002; Skilling, 2004; Feroz & Hobson, 2008; Feroz et al., 2009) have been introduced to make the calculation of these quantities more efficient. Even though these methods greatly reduce the computation time, and have been thoroughly tested and widely applied, they do not fully sample the entire parameter space and there is always danger that a secondary peak in a posterior distribution might remain undetected, or that the algorithm might get stuck in a local minimum. Because our models contain a low number of free parameters, we choose to fully sample the parameter space via a brute-force exploration method. The likelihood function is systematically calculated for each combination of the free parameters stated in Equations 13 and 14. In this calculation the amplitude ranges between 0 and 200 kms with a step size of 3 kms, ranges between 0 and radians with an interval of 0.1, ranges between 50 and 600 kms with a 5 kms increment, and ranges between -1.5 and 0.5 with a step size of 0.025. Careful testing has shown that this combination of parameters and sampling gives excellent balance between computational speed and resolution of the likelihood function and the posterior probability distributions. Finally, the integral in Equation 7 is evaluated by applying the Simpson rule in multiple dimensions.
3.3 Overall halo kinematics
In Veljanoski et al. (2013) we presented the first kinematic analysis for a significant number of outer halo M31 GCs. We discovered that (i) these clusters exhibit substantial net rotation; (ii) they share the same rotation direction and a similar rotation axis to centrally-located GCs; and (iii) this rotation axis is approximately coincident with the optical minor axis of M31. We also observed a hint of decreasing velocity dispersion with increasing galactocentric radius. In our present study, we want to determine the statistical significance of these phenomena by employing the models and methodology presented in the previous subsection.
This allows us to derive posterior probability distribution functions for the free parameters of each model. Since these distribution functions are not necessarily Gaussian (or even symmetric), we report both their peak and mean values accompanied by their 68% confidence limits in Table 6. This table also displays the logarithm of the Bayesian evidence for each model, which is used to discriminate between them. We find that the rotating model is decisively preferred over the non rotating one, with log(. The inferred amplitude of the rotation is kms.
The posterior probability distribution functions for the free parameters of model are shown in Figure 6, along with the 68% confidence limits around the peak and mean of each distribution. The strong preference for the favoured model over the non rotating model (for which ) can clearly be seen from the upper left panel. The position angle of the M31 optical minor axis is 135 east of north, matching the inferred rotation axis of the M31 outer halo GC system. As expected, the rotation of the outer GCs is in the same direction as their inner region counterparts albeit with a smaller amplitude. This is best seen in Figure 7, which shows the Galactocentric radial velocities of the GCs in our sample, corrected for the systemic radial motion of M31, versus their projected radii along the major (left panel) and minor (right panel) optical axes. The left panel of Figure 7 clearly shows that the rotation is observed even for the GCs with the largest projected distances, and is not driven solely by clusters projected onto major halo substructures or by clusters not lying on any visible substructure.
When modelling the rotation of the outer halo GC population, we assumed that the rotation axis lies in the plane of the sky – i.e., perpendicular to the line of sight. Thus, so far we have determined the projected rotation amplitude, and there is an additional unknown factor to account for, where is the inclination angle of the rotation axis to the plane of the sky. As we have found the rotation axis of the M31 outer halo GC population to coincide with the minor optical axis of this galaxy, it is possible that the rotation axis lies perpendicular to the disk of M31. In this case, taking the inclination of M31 with respect to our line of sight to be (Ferguson et al., 2002), the peak and mean of the deprojected rotation amplitude posterior probability distribution function would both be kms, barely different from the projected values.
We also find substantial evidence for decreasing velocity dispersion with increasing . Looking at the bottom right panel of Figure 6, it is seen that the peak and the mean of the posterior probability distribution function are inconsistent with . In fact, the posterior probability to measure is less than 1%. This is shown in more detail in Figure 8, which shows the 1, 2 and 3 levels of the likelihood on the - plane. It can easily be seen that a constant velocity dispersion as a function of may be discarded at approximately the 2 level.
Figure 9 shows the Galactocentric radial velocities of the outer halo GCs, corrected for both their bulk rotation as per model as well as for the systemic radial motion of M31, as a function of their projected radii (for convenience, we list the rotation-corrected velocities in Table 4). Different groups of GCs that lie along specific stellar streams are marked (see Section 4). The GC halo dispersion profile is displayed as a solid line as described by Equation 10 using the best fit parameters from Table 6. We also plotted the stellar velocity dispersion profile determined by Chapman et al. (2006) for metal-poor giant stars in the range between 10 and 70 kpc in projection, with the majority of the data points lying between 10 and 50 kpc. The stellar profile was assumed to be linear in shape. Note that beyond 70 kpc we have used a simple linear extrapolation. Figure 9 shows a close similarity between the velocity dispersions of the M31 halo stars and GCs, despite being fitted by different models, at least out to 80 kpc in projection. This similarity might imply that the spatial density profiles of the M31 halo stars and globulars are also similar, and indeed Huxor et al. (2011) have shown this to be the case by comparing the radial number density profile of the M31 GCs to that of the metal poor () stars (see their Figure 9).
For consistency, we note that the best-fit parameters for model are in very good agreement with the results presented in Veljanoski et al. (2013), where we determined the rotation and velocity dispersion separately and using only a subset of our currently-available radial velocity sample.
4 GC groups on streams
In the previous section we treated the M31 outer halo GCs as a single system and attempted to describe its global continuous properties. However, the M31 halo is rich with various field substructures in the form of stellar streams, loops and filaments. Many remote GCs appear spatially associated in projection with prominent features that are visible in the stellar maps (Mackey et al., 2010b), including a significant number of the clusters in our spectroscopic sample. Examination of Figure 5 reveals that objects projected onto a given feature tend to exhibit correlated velocities. It might naively be expected that coherent velocity patterns amongst GC groups would be quite unlikely to arise in the case where the GCs are randomly-selected members of a pressure-supported halo (even if a substantial rotation component is also present), but would be unsurprising in the case where they are associated with an underlying kinematically cold stellar debris feature. In the following sub-sections we consider several GC groups that project onto the main stellar substructures seen in the M31 halo, as marked on Figure 10, and attempt to indicatively assess the significance of any observed velocity patterns.
We proceed by employing simple Monte Carlo experiments similar to the one devised by Mackey et al. (2013). Although these tests are tailored to each specific instance, they all share a common basis. Our most-likely global kinematic model from the previous Section tells us what halo velocity dispersion and systemic rotation to expect at each GC position. For the GCs in a given group, we first subtract the global rotation signal from the observed velocities (which have, of course, already been corrected for the M31 systemic motion)444We remind the reader that these rotation-corrected GC velocities are listed in Table 4 for easy reference., and then generate 10 sets of GCs with positions matching those of the real set, but with each individual velocity randomly drawn from a Gaussian distribution centred on zero and with a width set by the dispersion model described by Equation 10 at the appropriate projected radius. This shows us what typical velocity configuration(s) to expect for the null hypothesis that all GCs are independent, uncorrelated members of the M31 halo, and thus allows us to broadly quantify how unusual any observed velocity pattern might be in this context.
Note that, ideally, we would simply match our GC velocities to kinematic measurements derived directly from a given stellar substructure to establish or refute any association between them. However, determining velocities for these extremely low surface-brightness features is a challenging observational problem and at present few such measurements have been published. In what follows we highlight only a couple of cases where it is possible to make such a direct comparison using extant data.
4.1 The North-West Stream
The North-West Stream is a narrow stellar debris feature extending radially over a range kpc. Projected on top lie 7 GCs for which we have measured velocities. Because of the radial nature of the stream, it is interesting to examine how the velocities of these GCs behave as a function of . This is shown in Figure 11, where we plot Galactocentric velocity, corrected for the measured rotation and systemic motion of M31, against projected radius.
Six of the seven NW Stream GCs share a clear trend in corrected radial velocity as a function of , in that the velocity becomes more strongly negative the closer a GC is to the centre of M31. However the innermost GC, PA-15 (the spectrum for which has S/N per Å), deviates substantially from this trend and, assuming its measured velocity is correct, it is difficult to see how this object could be associated with the NW Stream despite the fact that its position projects precisely onto the feature.
As marked in Figure 11, the relationship between the outermost five GCs on the stream is very close to linear, with a gradient of kms per kpc, a zero-point of kms, and a Pearson correlation coefficient . This is remarkable – we do not know of any compelling reason to expect a highly linear correlation between velocity and radial distance. Indeed we ascribe no important physical insight into this specific form of the relationship – fitting a straight line to the data is merely the simplest means of quantifying the observed trend. It is also notable that the NW Stream clusters lie substantially displaced by a magnitude kms from zero velocity, which is where the mean of the distribution of corrected halo velocities should sit. This is larger than the dispersion of the GC system at commensurate radii (see Figure 9).
Although the sixth NW Stream cluster, PA-13, does share in the trend of increasingly negative velocity with smaller , it is quite displaced from the linear relationship described above. However, this object has the lowest S/N spectrum in our entire GC sample (S/N per Å), and its velocity is thus accompanied by a large uncertainty such that its relationship to the stream is ambiguous. A more precise measurement for PA-13 would clearly be valuable.
We conducted a Monte Carlo experiment to consider the outermost five GCs on the NW Stream. We (conservatively) counted what fraction of our mock GC sets satisfied or , which is considered an indicator of high (anti)correlation. Around 3% of the simulated sets satisfy this criterion. If we only consider cases of infall, meaning the sets only need to satisfy , the probability of finding such a pattern falls to 2%. We also counted how many times all five GCs fell outside either from the mean (i.e., outside the measured dispersion at given radius). This is a very unusual configuration, occurring only % of the time. These two simple tests show that the kinematic properties of the five NW Stream GCs are almost certainly not due to a chance occurrence among independent halo GCs, providing additional convincing evidence, beyond their spatial alignment, for an association with each other and the underlying stellar stream. The observed velocity gradient amongst the GCs likely represents the infall trajectory of the progenitor satellite. The scatter of the five GC velocities about the best-fit line is very small, suggestive of a dynamically cold system with a dispersion consistent with zero.
4.2 The South-West Cloud
The South-West Cloud is a large diffuse overdensity near the major axis of M31 at kpc. It has been studied in detail by Bate et al. (2014) (see also Lewis et al., 2013). There are three GCs projected onto the Cloud, two of which (PA-7, PA-8) were shown by Mackey et al. (2013) to have velocities consistent with being members of this substructure.
Here we have measured a velocity for the third possible member of the sub-group, PA-14, as well as confirming the velocities for PA-7 and PA-8 from Mackey et al. (2013). Since the SW Cloud closely resembles an arc tracing roughly constant , we consider the rotation-corrected velocities of these three GCs as a function of position angle relative to the M31 centre (Figure 12). There is a clear velocity gradient along the arc, in that the corrected velocities become more negative with increasing position angle (i.e., in the anti-clockwise direction on Figure 10, or from south to north along the arc of the stream).
Once again, a linear fit does an excellent job of describing this trend. The best-fit line has a gradient of kms per degree, a zero-point of kms, and a correlation coefficient . As before, we do not ascribe any particular significance to this assumed form for the relationship – a linear fit is just the simplest means of quantifying the observed trend in velocity with position angle.
We ran a Monte Carlo experiment to consider the three SW Cloud GCs, and assessed what fraction of our mock sets exhibit a linear correlation with or . This is moderately common, arising % of the time. Our calculation supersedes that of Mackey et al. (2013), as their model did not include any correction for the systemic rotation because its existence was not known at that time. However, since the SW Cloud lies near the M31 major axis, it is important to account for the rotation signal.
This result is, alone, insufficient to robustly associate this sub-group of three GCs with each other; in addition, these objects do not have rotation-corrected velocities offset far from the expected mean of zero. However, we recall that the chance of all three GCs being projected directly onto the SW Cloud in the first place is quite small (%, Mackey et al., 2010b). Furthermore, Bate et al. (2014) have noted that spectroscopic measurements of the M31 field halo near to PA-14 by Gilbert et al. (2012) revealed a cold kinematic peak at a very similar velocity to that of the cluster: kms in the corrected frame, with an intrinsic dispersion kms (see Figure 12)555Here we have assumed that the M31 field halo is a subject to the same rotation effects as the GCs projected onto the SW Cloud, which is not unreasonable if the clusters trace the motions of the underlying substructure.. We conclude that, on balance, PA-7, PA-8, and PA-14 are all likely associated with each other and the SW Cloud. A radial velocity gradient along the arc, as suggested by the GCs, may imply substantial motion in the plane of the sky and hence a significant line-of-sight depth to the feature – as tentatively detected by both Mackey et al. (2013) and Bate et al. (2014).
4.3 Streams C and D
Streams C and D are two well-defined arc-like substructures due east of M31 (see Ibata et al., 2007; Richardson et al., 2011). This is a complex part of the M31 halo – the streams appear to overlap, in projection, at their northern end; in addition, Stream C is known to split into two kinematically distinct constituents (Chapman et al., 2008) – a metal-poor component, Cp, which is the narrow stream visible in Figures 5 and 10, and a metal-rich component, Cr, which overlaps Cp in projection but is spatially much broader. This latter feature is not evident in Figures 5 and 10 because its member stars fall above the metal-rich cut-off used in the construction of these plots; see however maps in, e.g., Martin et al. (2013); Ibata et al. (2014).
There is a total of 15 GCs projected on top of Streams C and D, all of which we have measured velocities for. The northern area of overlap is particularly rich in clusters, with 9 contained inside a relatively small region on the sky. Our observations suggest that these GCs split into two distinct kinematic subgroups. We employed the biweight location and scale estimators (Beers et al., 1990) to determine the mean velocity and velocity dispersion of each. The first contains five GCs (H24, PA-41, PA-43, PA-45, PA-46), has a mean rotation-corrected velocity of kms and dispersion of kms, while the second features three GCs (B517, PA-44, PA-47) with a mean rotation-corrected velocity of kms and a dispersion of kms. The ninth cluster in the region, PA-49 cannot be robustly identified with either kinematic subgroup. These measurements supersede those from our earlier work (Veljanoski et al., 2013), where we identified two similar kinematic groupings, but with an incomplete data set.
We conducted a Monte Carlo experiment for each GC subgroup to assess the likelihood that they are chance assemblies of unrelated objects. At kpc the expected velocity dispersion of halo GCs is kms. For the first subgroup, the fraction of mock sets where all five members lie outside either but with an internal velocity dispersion less than 10 kms, is extremely small at 0.02%. For the second subgroup, the configuration where three GCs all sit outside either while having an internal dispersion below 40 kms is somewhat more frequent, occurring of the time. Combined with the low probability of so many GCs clustering spatially (see Mackey et al., 2010b), we conclude that our GC groups are very likely associated with two of the underlying substructures. It is notable that the mean velocities of the subgroups fall either side of zero. Thus they, and, in all likelihood, two of the three overlapping streams, are in counter-rotating orbits about M31.
Following Streams C and D anti-clockwise in Figure 10, both have three GCs projected onto their southern regions. The three clusters projecting onto Stream D have velocities encompassing a range of kms; thus these objects do not form a kinematic subgroup. At present we are unable to robustly assess which, if any, of these GCs might be associated with Stream D.
The three Stream C clusters also have quite disparate velocities and do not constitute a kinematic subgroup. One of these objects, HEC12 (also known as EC4), is quite well studied. It lies precisely on the narrow metal-poor component Cp, and shares a common velocity with this stream (Chapman et al., 2008; Collins et al., 2009). Of the other two GCs, HEC13 lies less than a degree away and has a velocity very similar to that measured for the metal-rich component of Stream C by Chapman et al. (2008) ( kms and kms, respectively) – i.e., HEC13 is plausibly associated with Cr. The velocity of the third GC, H26, indicates that it is most likely a chance alignment that is not a member of either component of Stream C.
4.4 The Eastern Cloud
The Eastern Cloud is a small arc-like stellar debris feature located at kpc due east of M31. There are two GCs projected onto this overdensity – PA-57 and PA-58. These have a velocity difference of less than 20 kms and a mean rotation-corrected velocity that sits away from the expected systemic mean of zero (as before, is the dispersion of halo GCs at the appropriate radius).
The explore how commonly such an arrangement might occur we conducted a Monte Carlo experiment for these two GCs. The fraction of mock pairs that, conservatively, have both members sitting outside but with an absolute velocity difference between them of smaller than 30 kms, is around . However, the Eastern Cloud is a comparatively small overdensity in terms of its projected surface area. Mackey et al. (2010b) did not consider this feature as it had not been discovered at the time of their work. The arc subtends a position angle of and spans, generously, the radial range kpc (). Its projected surface area is hence deg. We also know that the surface density of GCs at this radius in the M31 halo is very close to deg (Mackey et al. 2014, in prep.). Applying simple Poisson statistics, we infer that the probability of observing two or more GCs falling in this region by chance is %. The chance that they also have very similar velocities, as per the calculation above, is then %. We conclude that the two GCs projected onto the Eastern Cloud are almost certainly associated with this substructure.
4.5 Association 2
Mackey et al. (2010b) identified a spatial overdensity of GCs, dubbed “Association 2”, close to the western major axis of M31 at a distance of kpc. It consists of 10 clusters666Two additional GCs were discovered in this region since the analysis by Mackey et al. (2010b). sitting within a small area, constituting the single highest local density enhancement of GCs, relative to the azimuthal average, seen in the M31 halo. This is a complex region where the outer disk and stellar halo of M31 overlap, and it is difficult to assess whether there might be one or more distinct underlying stellar features (see Figure 5).
We have measured velocities for all 10 GCs lying within the Association 2 region. The ensemble splits naturally into two kinematic sub-groups: (i) a set of four clusters (H2, PA-18, PA-19, PA-21) for which the biweight estimators give a mean rotation-corrected velocity of kms and a dispersion of kms; and (ii) a second set of four objects (H7, H8, PA-22, PA-23) for which the biweight indicators suggest a mean rotation-corrected velocity of kms and a dispersion of kms. The two remaining clusters, G1 and G2, cannot be robustly associated with either ensemble. Intriguingly, our first GC subgroup consists only of objects with in the range - kpc, while the second consists only of GCs with kpc. That is, splitting the overall ensemble by velocity, as we have done, naturally also results in clustering by position.
To assess the plausibility of these two apparently coherent sub-units forming by chance we conducted a Monte Carlo experiment for each. We find that the occurrence of 4 random GCs all lying outside either but having an internal dispersion less than kms, as per our first observed subgroup, is very rare at %. Our second subgroup, for which all 4 members sit outside either but have an inter-group dispersion smaller than kms is also very rare and arises % of the time in our model.
In summary, Association 2 is not a single kinematically coherent ensemble of clusters. Instead it is primarily comprised of two clearly distinct sub-groups, and thus may possibly represent the projected convergence of two relic systems. In this context it is interesting that Association 2 sits very close to the expected base of the North-West Stream. Our first GC sub-group, which has kpc and a rotation-corrected velocity of kms, matches very closely to an extrapolation of the linear velocity gradient along the NW Stream that we described above (see Figure 11), and it is tempting to speculate that this sub-group may be linked to that feature. Similarly, Ibata et al. (2005) studied a discrete fragment of the M31 outer disk very nearby to this region, and found its velocity to sit near kms (see also Reitzel et al., 2004; Faria et al., 2007). This matches well with the heliocentric velocities of the GCs in our second sub-group (see Table 4) – the weighted mean for which is kms. It is therefore possible that these GCs may be associated with the outer disk of M31.
It it also worth nothing that the GC G1 (also known as Mayall II) lies in the Association 2 region. This cluster is probably the brightest and most massive member of the M31 GC halo system, and it exhibits a number of peculiar properties analogous to the Galactic GC -Cen (e.g. Ma et al., 2007, 2009), which make it a likely galaxy remnant candidate (e.g. Meylan et al., 2001). It is therefore perhaps surprising that this cluster does not belong in either of the two kinematically identified groups within Association 2; in particular it does not fit with the trend set by the GCs associated with the NW Stream.
4.6 And XVII
Irwin et al. (2008), in their discovery paper for the M31 dwarf spheroidal (dSph) satellite And XVII, noted that three GCs lie very close to this system – H11, HEC6, and HEC3, with projected distances of 2.0 kpc, 3.7 kpc, and 5.9 kpc respectively. This is the only known instance of a possible association between GCs and an M31 dSph satellite. We have obtained velocity measurements for both H11 and HEC6, allowing us to assess whether either of these GCs might plausibly be gravitationally bound to And XVII.
The heliocentric radial velocity of And XVII is kms (Collins et al., 2013) while that for HEC6 is kms. This GC, the most distant of the three candidates, is clearly not associated with the dwarf. The situation is more complex for H11, which is the closest candidate. Its velocity is only separated from that of And XVII by kms in the heliocentric frame, and hence there is a higher chance it might be gravitationally bound. We investigate this via a simple Newtonian escape velocity argument:
where is the escape velocity, is the gravitational constant, is the total mass of the system and is the distance between the two objects. The only available mass estimate of And XVII comes from Collins et al. (2014), who find within the half-light radius of the galaxy. Under the assumption that kpc and the escape velocity is found to be just kms. Applying the same principles, And XVII would be required to have a total mass of at least in order for H11 to be gravitationally bound.
5.1 GC Kinematics
Our high quality PAndAS data has made exploring the true outer halo of M31 ( kpc) in a continuous and complete fashion possible for the first time. This region is seen to be dominated by various stellar debris features, thought to be the remnants of accreted dwarf galaxies (e.g. McConnachie et al., 2009; Ibata et al., 2014). In addition, a significant portion of the GCs discovered in the outer halo appear to preferentially lie projected on top of these debris features (Mackey et al., 2010b).
Using the measurements presented in the current paper, we have demonstrated that various discrete groups of such GCs – specifically those projecting onto the most luminous halo streams and overdensities – exhibit clear kinematic patterns. In Section 4 we used our global kinematic measurements from Section 3, in conjunction with simple Monte Carlo experiments, to indicatively assess how frequently these velocity trends or correlations might occur in the case where all the GCs in a given group are independent members of the M31 halo. Each instance we examined (apart from the ostensible And XVII association) showed clear evidence for non-random behaviour, indicative of a dynamical link between the GCs and the streams or overdensities that they project onto. Together these results strongly reinforce the notion from Mackey et al. (2010b) that a substantial fraction of the outer halo GC population of M31 consists of objects accreted along with their now-defunct host galaxies. A striking feature of many of the ensembles we considered is the coldness of their kinematics – cluster groups on the NW Stream, SW Cloud, and Eastern Cloud, as well as sub-groups in the Stream C/D overlap area and in Association 2, all exhibit velocity dispersions consistent with zero. At present, measurements of stream velocities directly from the constituent stars are available only in a handful of cases; however we have demonstrated that these few instances largely support our assertions.
In light of these results, our discovery of the high overall degree of coherent rotation exhibited by the outer halo GC population in M31 is rather surprising. It is relevant to note that this rotation is not predominantly driven by either the subset of remote clusters clearly associated with underlying stellar streams, or the complement of this subset. This is clearly evident from inspection of the left-hand panel of Figure 7. We have also found that the outer halo GCs share the same rotation axis as the inner halo GCs; indeed, the rotation of these two subsystems is virtually indistinguishable save for the difference in amplitude. This is in contrast with observations in the Milky Way where the halo GC population appears to exhibit at most only a mild net rotation (Harris, 2001; Brodie & Strader, 2006; Deason et al., 2011).
It is unfortunate that to date there is little opportunity to compare the kinematics of the M31 outer halo GCs to those of other massive spiral galaxies apart from the Milky Way. Olsen et al. (2004) derived kinematic properties for 6 spiral galaxies in the Sculptor group. However, the GCs in these galaxies mainly lie in the inner regions of their hosts, and the results of their study are likely to be affected by small number statistics. Nantais & Huchra (2010) presented a discovery of rapid rotation in the GC system of M81. However, all of the GCs with available radial velocity data in that galaxy lie at projected radii of less than 20 kpc. In elliptical galaxies it is common to find rotating GC populations towards more central regions, but finding significant rotation beyond a few tens of kpc appears to be a rare occurrence (e.g. Woodley et al., 2010; Strader et al., 2011; Blom et al., 2012; Pota et al., 2013).
It is interesting to consider how to reconcile our discovery of significant rotation in the outer halo of M31 with the chaotic accretion of dwarf galaxies implied by hierarchical models. One way this phenomenon might arise is through the major merger of two spiral galaxies. For example, the numerical simulations of Bekki (2010) suggest that a major merger between M31 and a similar spiral galaxy could give rise to the rapid rotation observed in the inner GC system of M31, including the rotation of the halo populations. More generally, a large fraction of the halo GC subsystem might have been brought into the potential well of M31 via a single moderate-mass satellite. In this event, the satellite seeds its GCs in the halo as it spirals in towards the centre of M31.
This idea is supported by both the thick disk of M31, which is found to rotate in the same sense (although more rapidly than) the outer halo GC population (Collins et al., 2011), and the kinematic properties of the inner spheroid, which also exhibits substantial rotation (Dorman et al., 2012). However, in order for such a satellite to be able to deposit several tens of GCs it would need to have a considerable mass – perhaps akin to the Large Magellanic Cloud, which possesses 16 old GCs (e.g., Mackey et al., 2006). If an encounter between M31 and such a massive satellite did occur, the question must arise as to how disruptive such an event would have been on the M31 disk. In addition, in this scenario (and indeed that involving the merger of two spirals) it may also be difficult to explain the observed presence of distinct dynamically cold subgroups of GCs as well as the typically narrow stellar debris streams in the halo. Detailed modelling is required to address these uncertainties.
Another possibility is that the outer halo GC system of M31 is indeed the product of the assimilation of multiple dwarf galaxies, but that these were accreted into the M31 potential well from a preferred direction on the sky. This would be consistent with the observation that many of the surviving dwarf galaxies associated with M31 lie in a thin rotating planar structure, as reported by Ibata et al. (2013). It is interesting to note that this plane of dwarfs also appears to be rotating in the same direction as the outer halo GC system, although its rotation axis is inclined by with respect to the minor axis of M31. A similar plane of dwarfs is observed in the Milky Way (e.g., Metz et al., 2007), and it has been shown that some of the GCs in the outer Galactic halo share this planar alignment (Keller et al., 2012; Pawlowski et al., 2012). It has been hypothesised that the formation of these planes of dwarfs, and by extension the possibility of accretion from a preferred direction, could occur as dwarf galaxies move along large scale dark matter filaments or sheets, in which case they are expected to have aligned angular momenta as seen in some recent cosmological simulations (Libeskind et al., 2005; Libeskind et al., 2011; Lovell et al., 2011). An alternative hypothesis is that these are tidal dwarfs formed during an early gas-rich merger (e.g., Pawlowski et al., 2013).
It is important to note that we are not suggesting that the properties of the outer halo GC system of M31 stem directly from the rotating plane of satellites reported by Ibata et al. (2013); indeed, there are almost no GCs with kpc lying close to this plane. Nonetheless, the observed fact that a substantial number of dwarf satellites of M31 possess correlated angular momenta raises the possibility that the parent galaxies of the accreted outer halo GC population may once have shared a similar, but now disrupted, configuration. In this context it is relevant that almost all of the dwarf galaxies thought to be members of the planes presently observed in both M31 and the Milky Way are insufficiently massive to host GCs. Hence the outer halo GC population in M31 might still have been assembled from only a few larger host systems, even if their accretion was related to a previous planar structure.
5.2 Application: the mass of M31
An accurate measurement of the mass of M31 is important for constraining the dynamics of Local Group galaxies, and for testing various cosmological models and predictions. Even though M31 is the closest massive galaxy to our own, it is striking that we have yet been unable to measure its mass to high precision; indeed there is even still debate as to whether M31 or the Milky Way is more massive. A number of studies have employed a variety of methods in order to estimate the mass of M31. Some of the more recent such estimates are displayed in Table 7. For summaries of older mass estimates we refer the reader to Federici et al. (1990); Federici et al. (1993); Evans & Wilkinson (2000).
|Fardal et al. (2013)||200 times the critical density||Inferred from the Giant Stream|
|Veljanoski et al. (2013)||within 200 kpc||Dynamical tracers - 50 outer halo GCs|
|van der Marel et al. (2012)||within the virial radius||Timing argument + literature + M33|
|Tollerud et al. (2012)||within the virial radius||Dynamical tracers - 19 dwarf galaxies|
|Watkins et al. (2010)||within 300 kpc||Dynamical tracers - 23 dwarf galaxies|
|Lee et al. (2008)||within 100 kpc||Dynamical tracers - 504 inner regions GCs|
|Evans et al. (2003)||within 100 kpc||Dynamical tracers - 89 inner regions GCs|
The M31 outer halo GCs can serve as kinematic mass tracers and therefore be used to provide an estimate of the total mass of their host galaxy. One way to do this is to solve the Jeans equations (Binney & Tremaine, 1987). Since we have found the M31 outer halo GC system to exhibit a significant degree of rotation, in is necessary to separate the solution into two parts: a rotating and a non-rotating, pressure supported component, the sum of which comprises the total mass of M31.
Since we have assumed the rotation of the halo GCs to occur only on simple circular orbits, the rotating mass component, , is simply determined via:
where is the rotation amplitude of the GC system, is the projected radius of the outermost GCs in the considered sample, and is the gravitational constant.
In order to determine the pressure supported mass contribution , we use the solution to the non-rotating Jeans equations proposed by Evans et al. (2003), dubbed the Tracer Mass Estimator (TME):
where is the projected radius from the M31 centre for a given GC, is the radial velocity of that GC with the rotational component removed, and is the total number of clusters in the sample under consideration. The index loops over each GC in the sample. The constant is dependent on the shape of the underlying gravitational potential, the radial distribution of the tracers and the anisotropy of the system. Here we assume that the M31 outer halo system is spherical and isotropic, and therefore takes the form of:
In the above definition of C, and are, respectively, the smallest and largest deprojected radii of the GCs in the ensemble being studied. For our present sample, the value of is taken to be the distance at which the halo begins to dominate, i.e. 30 kpc. The value of is chosen to be 200 kpc to reflect the measured radius of MGC1, which is the most remote known M31 cluster (Mackey et al., 2010a).
As required by the TME, the GC radial volume density distribution is approximated by a power-law, the index of which in the case of spherical symmetry is found to be (Mackey et al. 2014, in prep.). Hence, in the case of an isothermal halo potential, for which the parameter in Equation 18 is zero, we find the total mass enclosed within 200 kpc from the centre of M31 to be . If we assume a NFW profile for the M31 halo (Navarro et al., 1996), for which , the estimated mass is somewhat smaller, at .
In our earlier publication (Veljanoski et al., 2013), we used the same method with a sample of 50 GCs spanning 3D radii between 30 and 200 kpc in order to estimate the total mass of M31. Here we update this estimate by by enlarging our halo GC sample to contain a total of 72 GCs covering the same radial range. Our updated estimate of the dynamical mass of M31 is, perhaps not surprisingly, consistent with our previous value. Although these are also comparable to the majority of the dynamical mass estimates found in the literature that sample a similar spatial range (Table 7), it is important to state that there are a number of possible caveats related to using the TME in the present situation.
One important issue is that many of the M31 halo GCs appear to be spatially associated with cold stellar debris features, and indeed in this paper we have demonstrated that various groups of such GCs have correlated velocities as described in Section 4. This has two consequences. First, the TME assumes that the tracer population is in a steady state equilibrium, which is not necessarily true as many of the M31 halo GCs are likely to be relatively new arrivals. Second, because groups of GCs have correlated velocities, it is almost certainly not appropriate to treat all 72 objects that we used for the mass estimate as independent data points. In this case we are effectively weighting some data points more heavily than others, introducing a bias which is not understood for this specific case. That said, studies that have explored the presence of substructure in tracer populations found their results to be biased only by (Yencho et al., 2006; Deason et al., 2012).
We have also assumed, due to a lack of information to the contrary, that the velocity anisotropy of the GC orbits is zero, while in reality this is unlikely to be true for the whole halo population. Nonetheless, Di Cintio et al. (2012) found the anisotropy parameter to have a negligible effect on mass estimates derived from tracers for which only radial velocity information is available. Finally, there is no theoretical motivation to assume that the entire dark matter halo of a massive galaxy follows a single power law, and thus fixing to a single number might introduce additional biases. These caveats give rise to systematic uncertainties in our M31 mass estimate that are not incorporated in the quoted errors, which only contain the statistical uncertainty. Given the complex nature of the M31 halo GC population, a more reliable mass estimation may well require a method specifically tailored to this system.
In this contribution we presented radial velocity measurements for 78 GCs around M31, 63 of which have no previous spectroscopic information. Our sample extends from kpc out to kpc in projection, and at least up to 200 kpc in 3D, which enables us to explore the kinematic properties of the GCs located in the true outer halo of M31.
Our global kinematic analysis detected a significant degree of net rotation exhibited by the outer halo GC population of M31. Interestingly, this population shares the same rotation axis and direction as the GCs located in the inner regions of M31, as well as the M31 disk. We also find evidence for decreasing velocity dispersion as a function of projected distance from the centre of M31. The dispersion profile for the halo GC population is similar to that of the stellar halo, consistent with previous observations that the GCs and stars share similar spatial density profiles.
Our measurements further revealed a variety of velocity correlations for the various groups of GCs that lie projected on top of distinct stellar debris features in the field halo. In particular, several such GC groups appear to be kinematically cold, possessing velocity dispersions consistent with zero. Simple Monte Carlo experiments showed that it is highly unlikely that these velocity correlations are due to chance arrangements, but rather are most likely due to a common origin for the GCs and the stellar substructures. This further supports the idea that a significant fraction of the M31 halo GC system has an external origin (e.g., Mackey et al., 2010b). Definitive proof of the association between GCs and their underlying streams will require matching velocities between the GCs and the stellar members of the underlying substructure to be robustly established in each case. We highlighted the few cases where such a correlation can already tentatively be shown to exist.
In light of the clear association between many groups of halo GCs and underlying field substructures, our finding that the GC population as a whole possesses a substantial degree of coherent rotation out to very large radii is quite puzzling. It is difficult to reconcile this property with the chaotic accretion of parent dwarf galaxies into the halo as implied by our kinematic observations and hierarchical galaxy formation models. We speculate that the solution to this problem may be related to the recent discovery that many dwarf galaxies, both in the Milky Way and M31, appear to lie in thin rotating planar configurations such that their angular momenta are correlated. Alternatively, it may be that most of the outer halo GCs in M31 arrived with one or two relatively large host galaxies; however, detailed modelling is required to assess whether this scenario is compatible with the observed narrow stellar debris streams and the main features of the M31 disk.
Finally, using the halo GCs as kinematic tracers, we estimated the total mass of M31 enclosed within a deprojected radius of 200 kpc via the Tracer Mass Estimator. Even though our value of is in agreement with other recent mass estimates which employed kinematic tracers extending to similar radii as our sample, it is likely to be subject to several poorly-understood biases due to the various assumptions that are built into the TME.
We would like to thank the referee, Flavio Fusi Pecci, for a detailed and constructive report that helped improve the paper. We also thank Luca Ciotti for carefully reading the manuscript.
ADM is grateful for support by an Australian Research Fellowship (Discovery Projects Grant DP1093431) from the Australian Research Council. APH was partially supported by Sonderforschungsbereich SFB 881 “The Milky Way System” of the German Research Foundation.
The WHT is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.
This paper is 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 Australian Research Council (Australia), Ministério da Ciência, Tecnologia e Inovação (Brazil) and Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina). These observations were obtained under programmes GN-2010B-Q-19, GN-2011B-Q-61, GN-2012B-Q-77, and GN-2013B-Q-66.
- Abadi et al. (2003a) Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003a, ApJ, 591, 499
- Abadi et al. (2003b) Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003b, ApJ, 597, 21
- Alves-Brito et al. (2009) Alves-Brito A., Forbes D. A., Mendel J. T., Hau G. K. T., Murphy M. T., 2009, MNRAS, 395, L34
- Ashman & Bird (1993) Ashman K. M., Bird C. M., 1993, AJ, 106, 2281
- Bate et al. (2014) Bate N. F. et al., 2014, MNRAS, 437, 3362
- Beers et al. (1990) Beers T. C., Flynn K., Gebhardt K., 1990, AJ, 100, 32
- Bekki (2010) Bekki K., 2010, MNRAS, 401, L58
- Belokurov et al. (2006) Belokurov V. et al., 2006, ApJL, 642, L137
- Binney & Tremaine (1987) Binney J., Tremaine S., 1987, Galactic dynamics. Princeton, NJ, Princeton University Press
- Blom et al. (2012) Blom C., Forbes D. A., Brodie J. P., Foster C., Romanowsky A. J., Spitler L. R., Strader J., 2012, MNRAS, 426, 1959
- Brodie & Strader (2006) Brodie J. P., Strader J., 2006, ARAA, 44, 193
- Caldwell et al. (2011) Caldwell N., Schiavon R., Morrison H., Rose J. A., Harding P., 2011, AJ, 141, 61
- Chapman et al. (2008) Chapman S. C. et al., 2008, MNRAS, 390, 1437
- Chapman et al. (2006) Chapman S. C., Ibata R., Lewis G. F., Ferguson A. M. N., Irwin M., McConnachie A., Tanvir N., 2006, ApJ, 653, 255
- Collins et al. (2011) Collins M. L. M. et al., 2011, MNRAS, 413, 1548
- Collins et al. (2009) Collins M. L. M. et al., 2009, MNRAS, 396, 1619
- Collins et al. (2013) Collins M. L. M. et al., 2013, ApJ, 768, 172
- Collins et al. (2014) Collins M. L. M. et al., 2014, ApJ, 783, 7
- Conn et al. (2012) Conn A. R. et al., 2012, ApJ, 758, 11
- Côté et al. (2000) Côté P., Marzke R. O., West M. J., Minniti D., 2000, ApJ, 533, 869
- Côté et al. (2001) Côté P. et al., 2001, ApJ, 559, 828
- Courteau & van den Bergh (1999) Courteau S., van den Bergh S., 1999, AJ, 118, 337
- Da Costa & Armandroff (1995) Da Costa G. S., Armandroff T. E., 1995, AJ, 109, 2533
- Deason et al. (2011) Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 411, 1480
- Deason et al. (2012) Deason A. J. et al., 2012, MNRAS, 425, 2840
- Di Cintio et al. (2012) Di Cintio A., Knebe A., Libeskind N. I., Hoffman Y., Yepes G., Gottlöber S., 2012, MNRAS, 423, 1883
- Dorman et al. (2012) Dorman C. E. et al., 2012, ApJ, 752, 147
- Dubath & Grillmair (1997) Dubath P., Grillmair C. J., 1997, AAp, 321, 379
- Evans & Wilkinson (2000) Evans N. W., Wilkinson M. I., 2000, MNRAS, 316, 929
- Evans et al. (2003) Evans N. W., Wilkinson M. I., Perrett K. M., Bridges T. J., 2003, ApJ, 583, 752
- Fardal et al. (2013) Fardal M. A. et al., 2013, MNRAS, 434, 2779
- Faria et al. (2007) Faria D., Johnson R. A., Ferguson A. M. N., Irwin M. J., Ibata R. A., Johnston K. V., Lewis G. F., Tanvir N. R., 2007, AJ, 133, 1275
- Federici et al. (1993) Federici L., Bonoli F., Ciotti L., Fusi-Pecci F., Marano B., Lipovetsky V. A., Niezvestny S. I., Spassova N., 1993, AAp, 274, 87
- Federici et al. (1990) Federici L., Marano B., Fusi Pecci F., 1990, AAp, 236, 99
- Ferguson et al. (2002) Ferguson A. M. N., Irwin M. J., Ibata R. A., Lewis G. F., Tanvir N. R., 2002, AJ, 124, 1452
- Feroz & Hobson (2008) Feroz F., Hobson M. P., 2008, MNRAS, 384, 449
- Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
- Forbes & Bridges (2010) Forbes D. A., Bridges T., 2010, MNRAS, 404, 1203
- Galleti et al. (2009) Galleti S., Bellazzini M., Buzzoni A., Federici L., Fusi Pecci F., 2009, AAp, 508, 1285
- Galleti et al. (2007) Galleti S., Bellazzini M., Federici L., Buzzoni A., Fusi Pecci F., 2007, AAp, 471, 127
- Galleti et al. (2006) Galleti S., Federici L., Bellazzini M., Buzzoni A., Fusi Pecci F., 2006, AAp, 456, 985
- Galleti et al. (2004) Galleti S., Federici L., Bellazzini M., Fusi Pecci F., Macrina S., 2004, AAp, 416, 917
- Geehan et al. (2006) Geehan J. J., Fardal M. A., Babul A., Guhathakurta P., 2006, MNRAS, 366, 996
- Gilbert et al. (2012) Gilbert K. M. et al., 2012, ApJ, 760, 76
- Grillmair (2006) Grillmair C. J., 2006, ApJL, 645, L37
- Grillmair & Dionatos (2006) Grillmair C. J., Dionatos O., 2006, ApJL, 643, L17
- Guo et al. (2011) Guo Q. et al., 2011, MNRAS, 413, 101
- Harris (2001) Harris W. E., 2001, in L. Labhardt & B. Binggeli ed., Saas-Fee Advanced Course 28: Star Clusters. p. 223
- Hartwick & Sargent (1974) Hartwick F. D. A., Sargent W. L. W., 1974, ApJ, 190, 283
- Huchra et al. (1982) Huchra J., Stauffer J., van Speybroeck L., 1982, ApJL, 259, L57
- Huchra et al. (1991) Huchra J. P., Brodie J. P., Kent S. M., 1991, ApJ, 370, 495
- Huxor et al. (2011) Huxor A. P. et al., 2011, MNRAS, p. 452
- Huxor et al. (2014) Huxor A. P. et al., 2014, MNRAS, submitted
- Huxor et al. (2008) Huxor A. P., Tanvir N. R., Ferguson A. M. N., Irwin M. J., Ibata R., Bridges T., Lewis G. F., 2008, MNRAS, 385, 1989
- Huxor et al. (2005) Huxor A. P., Tanvir N. R., Irwin M. J., Ibata R., Collett J. L., Ferguson A. M. N., Bridges T., Lewis G. F., 2005, MNRAS, 360, 1007
- Ibata et al. (2005) Ibata R., Chapman S., Ferguson A. M. N., Lewis G., Irwin M., Tanvir N., 2005, ApJ, 634, 287
- Ibata et al. (2007) Ibata R., Martin N. F., Irwin M., Chapman S., Ferguson A. M. N., Lewis G. F., McConnachie A. W., 2007, ApJ, 671, 1591
- Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
- Ibata et al. (1995) Ibata R. A., Gilmore G., Irwin M. J., 1995, MNRAS, 277, 781
- Ibata et al. (2013) Ibata R. A. et al., 2013, Nature, 493, 62
- Ibata et al. (2014) Ibata R. A. et al., 2014, ApJ, 780, 128
- Irwin et al. (2008) Irwin M. J., Ferguson A. M. N., Huxor A. P., Tanvir N. R., Ibata R. A., Lewis G. F., 2008, ApJL, 676, L17
- Jeffreys (1961) Jeffreys H., 1961, The Clarendon Press, Oxford University Press
- Kass & Raftery (1995) Kass R. E., Raftery A. E., 1995, JSTOR, 90, 773
- Keller et al. (2012) Keller S. C., Mackey D., Da Costa G. S., 2012, ApJ, 744, 57
- Koposov et al. (2012) Koposov S. E. et al., 2012, ApJ, 750, 80
- Lee et al. (2008) Lee M. G., Hwang H. S., Kim S. C., Park H. S., Geisler D., Sarajedini A., Harris W. E., 2008, ApJ, 674, 886
- Lewis & Bridle (2002) Lewis A., Bridle S., 2002, PhRev, 66, 103511
- Lewis et al. (2013) Lewis G. F. et al., 2013, ApJ, 763, 4
- Libeskind et al. (2005) Libeskind N. I., Frenk C. S., Cole S., Helly J. C., Jenkins A., Navarro J. F., Power C., 2005, MNRAS, 363, 146
- Libeskind et al. (2011) Libeskind N. I., Knebe A., Hoffman Y., Gottlöber S., Yepes G., Steinmetz M., 2011, MNRAS, 411, 1525
- Lovell et al. (2011) Lovell M. R., Eke V. R., Frenk C. S., Jenkins A., 2011, MNRAS, 413, 3013
- Ma et al. (2007) Ma J. et al., 2007, MNRAS, 376, 1621
- Ma et al. (2009) Ma J. et al., 2009, Research in Astronomy and Astrophysics, 9, 641
- Mackey et al. (2010a) Mackey A. D. et al., 2010a, MNRAS, 401, 533
- Mackey & Gilmore (2004) Mackey A. D., Gilmore G. F., 2004, MNRAS, 355, 504
- Mackey et al. (2010b) Mackey A. D. et al., 2010b, ApJL, 717, L11
- Mackey et al. (2013) Mackey A. D. et al., 2013, MNRAS, 429, 281
- Mackey et al. (2006) Mackey A. D., Payne M. J., Gilmore G. F., 2006, MNRAS, 369, 921
- Mackey & van den Bergh (2005) Mackey A. D., van den Bergh S., 2005, MNRAS, 360, 631
- Majewski et al. (2003) Majewski S. R., Skrutskie M. F., Weinberg M. D., Ostheimer J. C., 2003, ApJ, 599, 1082
- Marín-Franch et al. (2009) Marín-Franch A. et al., 2009, ApJ, 694, 1498
- Martin et al. (2006) Martin N. F., Ibata R. A., Irwin M. J., Chapman S., Lewis G. F., Ferguson A. M. N., Tanvir N., McConnachie A. W., 2006, MNRAS, 371, 1983
- Martin et al. (2013) Martin N. F., Ibata R. A., McConnachie A. W., Dougal Mackey A., Ferguson A. M. N., Irwin M. J., Lewis G. F., Fardal M. A., 2013, ApJ, 776, 80
- Martin et al. (2014) Martin N. F. et al., 2014, ApJ, submitted
- McConnachie et al. (2009) McConnachie A. W., Irwin M. J., Ibata R. A. e., 2009, Nature, 461, 66
- McMillan (2011) McMillan P. J., 2011, MNRAS, 414, 2446
- Metz et al. (2007) Metz M., Kroupa P., Jerjen H., 2007, MNRAS, 374, 1125
- Meylan et al. (2001) Meylan G., Sarajedini A., Jablonka P., Djorgovski S. G., Bridges T., Rich R. M., 2001, AJ, 122, 830
- Nantais & Huchra (2010) Nantais J. B., Huchra J. P., 2010, AJ, 139, 2620
- Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
- Olsen et al. (2004) Olsen K. A. G., Miller B. W., Suntzeff N. B., Schommer R. A., Bright J., 2004, AJ, 127, 2674
- Palma et al. (2002) Palma C., Majewski S. R., Johnston K. V., 2002, ApJ, 564, 736
- Pawlowski et al. (2013) Pawlowski M. S., Kroupa P., Jerjen H., 2013, MNRAS, 435, 1928
- Pawlowski et al. (2012) Pawlowski M. S., Pflamm-Altenburg J., Kroupa P., 2012, MNRAS, 423, 1109
- Peñarrubia et al. (2010) Peñarrubia J., Benson A. J., Walker M. G., Gilmore G., McConnachie A. W., Mayer L., 2010, MNRAS, 406, 1290
- Peñarrubia et al. (2002) Peñarrubia J., Kroupa P., Boily C. M., 2002, MNRAS, 333, 779
- Peñarrubia et al. (2008) Peñarrubia J., Navarro J. F., McConnachie A. W., 2008, ApJ, 673, 226
- Perina et al. (2009) Perina S., Federici L., Bellazzini M., Cacciari C., Fusi Pecci F., Galleti S., 2009, AAp, 507, 1375
- Perrett et al. (2002) Perrett K. M., Bridges T. J., Hanes D. A., Irwin M. J., Brodie J. P., Carter D., Huchra J. P., Watson F. G., 2002, AJ, 123, 2490
- Perrett et al. (2003) Perrett K. M., Stiff D. A., Hanes D. A., Bridges T. J., 2003, ApJ, 589, 790
- Peterson (1989) Peterson R., 1989, in Merritt D., ed., Dynamics of Dense Stellar Systems. pp 161–166
- Pota et al. (2013) Pota V. et al., 2013, MNRAS, 428, 389
- Purcell et al. (2007) Purcell C. W., Bullock J. S., Zentner A. R., 2007, ApJ, 666, 20
- Reitzel et al. (2004) Reitzel D. B., Guhathakurta P., Rich R. M., 2004, AJ, 127, 2133
- Richardson et al. (2011) Richardson J. C. et al., 2011, ApJ, 732, 76
- Schönrich et al. (2010) Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829
- Schuberth et al. (2010) Schuberth Y., Richtler T., Hilker M., Dirsch B., Bassino L. P., Romanowsky A. J., Infante L., 2010, AAp, 513, A52
- Schuberth et al. (2012) Schuberth Y., Richtler T., Hilker M., Salinas R., Dirsch B., Larsen S. S., 2012, AAp, 544, A115
- Searle & Zinn (1978) Searle L., Zinn R., 1978, ApJ, 225, 357
- Skilling (2004) Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American Institute of Physics Conference Series Vol. 735, American Institute of Physics Conference Series. pp 395–405
- Slater et al. (2013) Slater C. T. et al., 2013, ApJ, 762, 6
- Strader et al. (2011) Strader J. et al., 2011, ApJS, 197, 33
- Tollerud et al. (2012) Tollerud E. J. et al., 2012, ApJ, 752, 45
- Udry et al. (1999) Udry S. et al., 1999, in Hearnshaw J. B., Scarfe C. D., eds, Astronomical Society of the Pacific Conference Series Vol. 185, IAU Colloq. 170: Precise Stellar Radial Velocities. p. 383
- van den Bergh (1969) van den Bergh S., 1969, ApJS, 19, 145
- van der Marel et al. (2012) van der Marel R. P., Fardal M., Besla G., Beaton R. L., Sohn S. T., Anderson J., Brown T., Guhathakurta P., 2012, ApJ, 753, 8
- van der Marel & Guhathakurta (2008) van der Marel R. P., Guhathakurta P., 2008, ApJ, 678, 187
- Veljanoski et al. (2013) Veljanoski J. et al., 2013, ApJL, 768, L33
- Watkins et al. (2010) Watkins L. L., Evans N. W., An J. H., 2010, MNRAS, 406, 264
- Woodley et al. (2010) Woodley K. A., Gómez M., Harris W. E., Geisler D., Harris G. L. H., 2010, AJ, 139, 1871
- Yencho et al. (2006) Yencho B. M., Johnston K. V., Bullock J. S., Rhode K. L., 2006, ApJ, 643, 154