Global Properties of M31's Stellar Halo from the SPLASH Survey. III. Measuring the Stellar Velocity Dispersion Profile^∗

The stellar spectra were obtained with 124 multi-object spectroscopic slitmasks in 50 fields spread throughout M31's stellar halo (Figure 1). The masks span a large range in azimuth and projected radii from the center of M31. The fields were chosen to target relatively smooth areas of M31's halo, individual tidal debris features, and dwarf satellites (Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The majority of the multi-object spectroscopic slitmasks were designed from images in the Washington system M and T₂ filters and the intermediate-width DDO51 filter, obtained with the Mosaic camera on the 4 m Mayall telescope at Kitt Peak National Observatory (KPNO).¹³ (Ostheimer 2003; Beaton 2014). The innermost spectroscopic slitmasks were designed using g'- and i'-band imaging obtained with MegaCam instrument on the 3.6 m Canada–France–Hawaii Telescope (CFHT¹⁴ ; Kalirai et al. 2006b). A small number of spectroscopic masks in the outer halo of M31 were designed from Johnson–Cousins V- and I-band imaging obtained with the Suprime-Cam instrument on the Subaru Telescope (fields "streamE" and "streamF"; Tanaka et al. 2010) and the William Herschel Telescope (field "and10"; Zucker et al. 2007).

Stars were prioritized for inclusion on the spectroscopic masks based on their colors and magnitudes. Stars with colors and magnitudes consistent with RGB stars at the distance of M31 were assigned high priority for inclusion on the spectroscopic masks, with brighter RGB stars (within ∼1–1.5 magnitudes of the tip of the RGB) given highest priority. When available, the surface-gravity sensitive intermediate-band DDO51 photometry was also used to prioritize stars with a high probability of being RGB stars (based on their location in the M–T₂, M–DDO51 color–color diagram; Majewski et al. 2000).

The spectra were obtained with the DEIMOS spectrograph on the Keck II 10 m telescope over 10 observing seasons (Fall 2002–2011). All of the spectra were obtained using the 1200 line mm⁻¹ grating, which produces a dispersion of 0.33 Å pix⁻¹. The survey used a slit width of 1'', which yields a resolution of 1.6 Å FWHM. The wavelength range of the observed spectra is λλ ∼ 6450–9150 Å, which includes the Ca ii triplet absorption feature (∼8500 Å) and the Na i absorption feature (8190 Å). The Keck/DEIMOS spectra were reduced using the spec2d (flat-fielding, night-sky emission line removal, and extraction of one-dimensional spectra) and spec1d (redshift measurement) software developed at the University of California, Berkeley (Cooper et al. 2012; Newman et al. 2013).

Only stellar spectra with secure velocity measurements are included in the final data set (Gilbert et al. 2006). A heliocentric correction is applied to the measured line-of-sight velocities as well as a correction for imperfect centering of the star within the slit (using the observed position of the atmospheric A-band absorption feature relative to night-sky emission lines; Simon & Geha 2007; Sohn et al. 2007).

The mean (median) signal-to-noise ratio (S/N) per Angstrom of the stellar spectra with secure velocity measurements is 11.7 (8.2). The mean (median) S/N per Angstrom of the spectra of stars that are more probable to be M31 stars than MW stars (Section 3.2) is 7.8 (4.4). The mean (median) velocity measurement uncertainty of all stellar spectra with secure velocity measurements is 5.7 km s⁻¹ (4.8 km s⁻¹), while the mean (median) velocity uncertainty of stars more probable to be M31 stars than MW stars is 6.5 km s⁻¹ (5.6 km s⁻¹). The uncertainties on the velocity measurements are calculated by adding in quadrature the random velocity measurement uncertainty, estimated from the cross-correlation routine, and a systematic uncertainty of 2.2 km s⁻¹, estimated from repeat observations of stars (Simon & Geha 2007).

Almost one-third of the fields that are farther than 4° from M31's center targeted dwarf satellite galaxies (Figure 1; Majewski et al. 2007; Kalirai et al. 2009, 2010; Tollerud et al. 2012). In these fields, a significant number of the stars observed on each spectroscopic mask are dSph members, rather than M31 halo stars or MW stars along the line of sight.

Stars that are likely to be gravitationally bound to a dwarf satellite galaxy are identified following the method outlined by Gilbert et al. (2009b). M31's dSphs are spatially compact, have small velocity dispersions, and span a limited range of [Fe/H]. The spatial extents of the majority of the M31 dSphs are small enough that they cover only a portion of a DEIMOS spectroscopic slitmask. Only stars within the King limiting radius are considered potential dSph members. Stars that are outside the King limiting radius are included in our final data set. Thus we explicitly classify any extra-tidal dSph stars as M31 halo stars. This number is small: only ∼5% of the stars beyond the King limiting radius in dSph fields have velocities and [Fe/H] values consistent with the dSph. In addition, any stars within the King limiting radius of the dSph but well removed from the distribution of dSph stars in line-of-sight velocity or [Fe/H] are classified as M31 halo stars. The interested reader can find examples for And I and And III in Figures 3 and 4 of Gilbert et al. (2009b).

After we removed stars classified as dSph members, the final data set contained 5299 stars.

We use the method established by Gilbert et al. (2006) to determine the likelihood that an individual star is an RGB star in M31 or a foreground MW dwarf star along the line of sight to M31. The Gilbert et al. method determines the probability that a star is an M31 RGB or MW dwarf star from multiple photometric- and spectroscopic-based diagnostics. The full set of diagnostics includes (1) line-of-sight velocity, (2) location in the (M–T₂, M–DDO51) color–color diagram (when available), (3) strength of the Na i doublet absorpt-ion line as a function of (V − I) color¹⁵ , (4) location in the (I, V − I) color–magnitude diagram, and (5) comparison of spectroscopic and photometric metallicity estimates. Each diagnostic provides a (log-)likelihood that the star is an RGB or dwarf star: L = log₁₀(P_RGB/P_dwarf). The overall likelihood, $\langle {L}_{i}\rangle$, that a star is an M31 red giant or MW dwarf is defined as the sum of the individual log-likelihoods for each diagnostic (i.e., the product of the likelihoods).

An additional factor not included in the Gilbert et al. method is the projected radius from M31's center: the relative stellar density of M31 and MW populations changes dramatically as a function of M31-centric radius. For a given set of values of the stellar properties included in $\langle {L}_{i}\rangle$, a star is in fact more likely to be an M31 star in the inner regions compared to the outer regions. This is described in more detail in Section 4.2.3, where we explicitly include this in the model.

Since our aim is to model the stellar velocities, we do not include the velocity diagnostic in the computation of $\langle {L}_{i}\rangle$. Figure 3 shows the velocity of stars as a function of projected radius from M31, color-coded by $\langle {L}_{i}\rangle$ values computed with (left) and without (right) the velocity diagnostic. When velocity is included in the $\langle {L}_{i}\rangle$ computation, there is a strong trend of increasing $\langle {L}_{i}\rangle$ with decreasing line-of-sight velocity, which is a direct result of using the velocity diagnostic in the $\langle {L}_{i}\rangle$ calculation. Although weaker, a correlation of $\langle {L}_{i}\rangle$ with line-of-sight velocity is still evident when velocity is not included in the $\langle {L}_{i}\rangle$ computation, reflecting the velocity distributions of the M31 and MW populations.

Zoom In Zoom Out Reset image size

In previous papers, we defined samples of M31 and MW stars based on $\langle {L}_{i}\rangle$ thresholds, and stars significantly bluer than the most metal-poor, 10 Gyr RGB isochrone have been classified as securely identified MW stars regardless of their $\langle {L}_{i}\rangle$ values (Gilbert et al. 2006). This acknowledged that many of the empirical diagnostics have little discriminating power for stars with blue colors, as well as the fact that these stars are much more likely to be MW stars than RGB stars in M31. In this paper we use the $\langle {L}_{i}\rangle$ values directly, rather than using subsets of stars classified as belonging to MW or M31. Thus we must choose how to treat stars bluer than the most metal-poor RGB isochrone in the following analysis. Rather than removing them from the sample altogether, we set their $\langle {L}_{i}\rangle$ values to −5. This places these blue stars in the tail of the $\langle {L}_{i}\rangle$ distribution of MW stars, acknowledging that it is very unlikely for them to be M31 RGB stars. This is equivalent to what has been done in our previous papers. The sensitivity of the final results to this choice is explored in Section 5.2.

Kinematically cold tidal debris features have been identified in a significant fraction of the spectroscopic masks (Figures 1 and 2). Many of the tidal debris features identified in individual fields are related to a single accretion event: the Giant Southern Stream and its associated shell features (e.g., Ibata et al. 2001; Fardal et al. 2007; Gilbert et al. 2007, 2009b).

Each of the tidal debris features in the data set have been identified and characterized via maximum-likelihood, multi-Gaussian fits to the velocity distribution of M31 stars in the field, and presented in previous papers (Guhathakurta et al. 2006; Kalirai et al. 2006b; Gilbert et al. 2007, 2009b, 2012). Details of the fitting technique can be found in the papers by Gilbert et al. (2007, 2012). Each maximum-likelihood fit includes a Gaussian with a large velocity dispersion, representing the underlying, kinematically hot stellar halo, and additional Gaussian components with small velocity dispersions, representing the kinematically cold stellar streams.

The velocity distribution of the M31 halo was held fixed in all the published fits, with a mean heliocentric velocity of $\langle {v}_{\mathrm{hel}}\rangle =-300\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (M31's systemic velocity) and a velocity dispersion of σ_v = 129 km s⁻¹ (Gilbert et al. 2007). The maximum-likelihood fits provide an estimate of the mean velocity, velocity dispersion, and fraction of M31 halo stars in each of the kinematically cold tidal debris features in a field. Using the maximum-likelihood fits, the probability that any individual M31 star belongs to the kinematically hot halo or a kinematically cold tidal debris feature can be computed. Figure 4 shows the line-of-sight velocity as a function of projected radius of all M31 and MW stars in the data set, color-coded by the probability that the star belongs to tidal debris.

Zoom In Zoom Out Reset image size

In previous studies, we successfully identified M31 and MW samples using cuts on the various membership criteria. This work has included analyses of substructure (Gilbert et al. 2007, 2009a, 2009b), the surface brightness profile of M31's halo (Gilbert et al. 2012), and the metallicity profile of M31's halo (Kalirai et al. 2006a; Gilbert et al. 2014). Depending on the analysis, the sample selection was tuned to favor either a complete sample with some amount of contamination, or a clean but incomplete sample with minimal contamination. To produce the cleanest sample of M31 stars, many of the analyses required that a star be three times as likely to be an M31 red giant than a foreground MW dwarf star to be included in the M31 sample (Section 3.2).

However, simple threshold cuts on membership criteria are insufficient for the current analysis of the M31 halo velocity distribution, with the exception of stars belonging to the dwarf satellites of M31. We maintain the simple cut described above (Section 3.1) to identify and remove dSph stars from the sample. The small velocity dispersions, limited range of color (due to the limited range in [Fe/H] of the dSph stars) and relatively well-defined spatial boundaries of the M31 dSphs allow us to identify and remove probable dwarf satellite members while minimally affecting the observed velocity distribution of M31 halo and MW stars in those fields. Even in dSph fields with relatively high densities of halo stars, at most two or three M31 halo stars are removed from the sample. Moreover, the dSphs in the sample have a broad range of systemic line-of-sight velocities. Thus, the small amount of incompleteness introduced into the M31 and MW data set by removing dSph stars will not introduce a bias in the resulting velocity distribution.

In contrast, stars belonging to tidal debris features span the full spatial range of the spectroscopic masks. Moreover, while the stars belonging to some tidal debris features have a limited range of colors and velocities, this is not universally true. The tidal debris features are also significantly less dominant in a field than a dSph is, making the overlap of halo stars and tidal debris features in color and velocity space more significant. Finally, the probability that a given star belongs to the kinematically hot halo or a kinematically cold tidal debris feature, based on the fits described in Section 3.3, is calculated assuming a fixed mean velocity and velocity dispersion of the kinematically hot M31 halo. Thus, using these probabilities directly in our fits would introduce an internal inconsistency in the velocity measurement of the M31 halo as a function of radius, and may result in a systematic bias in the measurements.

MW stars and M31 stars have complete overlap in their spatial distribution, and large overlap in their velocity and color distributions. However, the velocity diagnostic cannot be included in the likelihoods, as justified in Section 3.2. Since the velocity diagnostic provides strong discriminating power between the M31 and MW populations, there is substantially more overlap of the $\langle {L}_{i}\rangle$ distributions of the M31 and MW populations when using $\langle {L}_{i}\rangle$ computed without the velocity diagnostic. This results in an unavoidable increase in contamination and decrease in completeness in any sample defined using only an $\langle {L}_{i}\rangle$ threshold. Furthermore, since the relative densities of MW and M31 stars vary with projected radius, the level of contamination and completeness of any given sample will also vary with projected radius. Both contamination and completeness will affect the resulting velocity distribution of the final M31 sample (Gilbert et al. 2007).

The biases in the velocity distribution of any sample defined via an $\langle {L}_{i}\rangle$ threshold occur primarily in the region of signi-ficant overlap between the MW and M31 velocity distributions (∼−200 < v_hel < −125 km s⁻¹). There are expected to be very few MW stars with heliocentric velocities v_hel < −300 km s⁻¹ in our sample (Robin et al. 2003). It is thus conceivable to consider using only stars with line-of-sight velocities more negative than M31's systemic velocity to identify an M31 sample using the Gilbert et al. (2006) likelihood technique (Section 3.2). However, there are several arguments against adopting this approach. First, it increases the sensitivity of the results to substructure, since many of the tidal debris features observed in SPLASH have mean line-of-sight heliocentric velocities $\langle {v}_{\mathrm{hel}}\rangle \lt -300\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Figure 4; many of these individual features are part of the Giant Southern Stream). It also reduces the size of the M31 halo sample by half, and in the sparse outer halo, the number of observed M31 stars is already small. Finally, it requires assuming a mean velocity for M31's halo, removing any sensitivity to departures in the mean line-of-sight velocity of M31's halo from M31's systemic velocity. Hence, we do not include velocity information in the membership likelihoods, but rather explicitly model the velocity distributions of all the known populations.

The full set of fields spans a significant area on the sky, with the largest angular separations between fields surpassing 20°. To eliminate the effects of perspective motion, all line-of-sight velocities are transformed to the Galactocentric frame, and the bulk motion of M31 along the line of sight to each star is removed. After this transformation, a star with no peculiar velocity relative to M31's bulk motion will have v_pec = 0 km s⁻¹, regardless of its position on the sky.

To facilitate comparison with measurements of the velocity dispersion of M31's globular cluster population, we follow Veljanoski et al. (2014) in using the Courteau & van den Bergh (1999) relation, with the McMillan (2011) estimate of the circular speed of the Galaxy's disk at the Sun (239 km s⁻¹) and the Schönrich et al. (2010) values for the solar peculiar motion [(U, V, W)_⊙ = (11.1, 12.24, 7.25) km s⁻¹]. Thus, the transformation of the observed heliocentric line-of-sight stellar velocities, v_helio,los to the Galactocentric frame, v_Gal is given by

$$\begin{eqnarray}{v}_{\mathrm{Gal}} & = & {v}_{\mathrm{helio}}+251.24\sin (l)\cos (b)\\ & & +11.1\cos (l)\cos (b)+7.25\sin (b),\end{eqnarray} \tag{ 1 }$$

where l and b are the Galactic coordinates (longitude and latitude) of the star.

Performing the transformation to the Galactocentric frame is vital, as it removes perspective effects in the outer halo fields that are on the same order as the dispersion we are trying to measure (several tens of km s⁻¹). However, the analysis presented here is not sensitive to the exact values assumed in Equation (1) (for alternate values, see Bland-Hawthorn & Gerhard 2016). With the SPLASH data set, the difference in Galactocentric transformations assuming alternate values, such as 218 km s⁻¹ for the circular speed of the disk (Bovy & Rix 2013), versus the transformations found using the nominal values above, is minimal: the typical effect on the measured halo dispersion would be on the order 2–3 km s⁻¹ or smaller. This is less (by a factor of two) than the typical velocity measurement error for the sample, and is more than an order of magnitude smaller than the expected halo dispersion.

To remove the bulk motion of M31 along the line of sight to each star, we assumed a heliocentric velocity for M31 of v_M31,helio = −301 km s⁻¹, corresponding to a Galactocentric radial velocity of v_M31,r = −109 km s⁻¹ (e.g., van der Marel & Guhathakurta 2008), and a transverse velocity (in the Galactocentric frame) of v_M31,t = 17 km s⁻¹, with a position angle of θ_M31,t = 287° (van der Marel et al. 2012). The removal of M31's motion from the line-of-sight velocities transformed to the Galactocentric frame, resulting in peculiar line-of-sight velocities for each star, v_pec, is calculated following van der Marel & Guhathakurta (2008):

$$\begin{eqnarray}{v}_{\mathrm{pec}} & = & {v}_{\mathrm{Gal}}-{v}_{{\rm{M}}31,{\rm{r}}}\cos (\rho )\\ & & +{v}_{{\rm{M}}31,{\rm{t}}}\sin (\rho )\cos (\phi -{\theta }_{{\rm{M}}31,{\rm{t}}}),\end{eqnarray} \tag{ 2 }$$

where ρ is the angular separation between M31's center and the star and ϕ is the position angle of the star with respect to M31's center. The uncertainties in M31's tangential motion are rather large (the 1σ confidence interval on v_M31,t is ≤34.3 km s⁻¹). However, as with the uncertainties in the transformation to the Galactic reference frame, we calculate the typical effect of these uncertainties on the measured dispersion to be small compared to the expected halo dispersion, on the order of ∼3 km s⁻¹ or less.

In the analysis that follows, all velocities have been transformed to v_pec, using Equations (1) and 2. Readers are referred to Veljanoski et al. (2014) for a broader discussion of the above transformations in the context of M31 kinematical analyses.

Our goal is to determine the most likely model parameters that describe the observed velocity distribution of stars along the line of sight to M31 as a function of projected radius from M31's center. We accomplish this by inferring the probability distributions for the parameters of a probabilistic generative model for the data, using Bayes' theorem.

The primary challenge in constructing the model is in determining the likelihood function. We construct the likelihood function by making the simplifying assumption that the line-of-sight velocity distribution of each of the individual stellar populations present in any line of sight can be adequately modeled by a Gaussian distribution. We then describe the likelihood function as a combination of M31 (Section 4.2.1) and MW (Section 4.2.2) models, each of which is composed of a combination of normalized Gaussians. Each individual stellar population in M31 and the MW thus contributes a set of parameters to the model, namely the mean velocity (μ), velocity dispersion (σ), and a relative fraction at which it contributes to the M31 or MW populations (f). In the likelihood function, each Gaussian component in the model is evaluated at the velocity of each data point, which is notated in Sections 4.2.1 and 4.2.2 by ${\mathscr{N}}({v}_{i}| \mu ,\sigma )$.

If MW or M31 membership were known precisely for each star, the likelihood of the observed line-of-sight velocity for an individual star given the model parameters Θ could be written as

$$\begin{eqnarray}&&{{\mathscr{L}}}_{i}={\eta }_{i}\,{{\mathscr{L}}}_{i}^{{\rm{M}}31}+(1-{\eta }_{i}){{\mathscr{L}}}_{i}^{\mathrm{MW}},\end{eqnarray} \tag{ 3 }$$

where η_i = 1 if the star is an RGB star in M31 and η_i = 0 if the star is an MW star along the line of sight, and ${{\mathscr{L}}}_{i}^{{\rm{M}}31}$ and ${{\mathscr{L}}}_{i}^{\mathrm{MW}}$ are the likelihoods for the M31 and MW models. We do not know a priori which stars are M31 and MW stars. However, we do have prior information on the probability of M31 membership for each star, obtained by evaluating the Gilbert et al. (2006) photometric and spectroscopic diagnostics (Section 3.2).

Thus, we combine the M31 and MW likelihood functions using a mixture model. For an individual star, the likelihood of the observed line-of-sight velocity v_i given the model parameters Θ is

$$\begin{eqnarray}&&{{\mathscr{L}}}_{i}={p}_{{\rm{M}}31,{\rm{i}}}\,{{\mathscr{L}}}_{i}^{{\rm{M}}31}+(1-{p}_{{\rm{M}}31,{\rm{i}}}){{\mathscr{L}}}_{i}^{\mathrm{MW}},\end{eqnarray} \tag{ 4 }$$

with p_M31,i describing the prior probability that the star is an M31 RGB star (Section 4.2.3). The likelihood of the observed data set given the model parameters Θ is simply the product of the individual likelihoods:

$$\begin{eqnarray}&&{{\mathscr{L}}}_{\theta }=\displaystyle \prod _{i=1}^{{N}_{\mathrm{stars}}}{{\mathscr{L}}}_{i,\theta }.\end{eqnarray} \tag{ 5 }$$

Finally, for a set of model parameters Θ, we compute

$$\begin{eqnarray}&&p({\rm{\Theta }}| {{v}_{i}}_{i=1}^{N},I)\propto p({{v}_{i}}_{i=1}^{N}| {\rm{\Theta }},I)p({\rm{\Theta }}| I),\end{eqnarray} \tag{ 6 }$$

where ${{v}_{i}}_{i=1}^{N}$ is the set of observed velocities, I represents our prior knowledge, and $p({{v}_{i}}_{i=1}^{N}| {\rm{\Theta }},I)$ is the likelihood term (Equation (5)). Equation (6) simply asserts that the probability of the set of model parameters Θ, given the observed data and all prior information, is proportional to the probability of the observed data given the model and all prior information, multiplied by the probability of the model given all prior information.

4.2.1. M31 Model

We characterize the population of M31 stars as a mixture of all known stellar components in our spectroscopic fields: a kinematically hot halo with multiple distinct, kinematically cold tidal debris features. This means that the M31 model, and thus the M31 likelihood function that is evaluated in Equation (4), is dependent on the spectroscopic field, sf, in which the star i is observed.

The number of kinematically cold tidal debris features is based on our previous analyses of individual spectroscopic fields (Section 3.3; Gilbert et al. 2007, 2009b, 2012). For the purposes of building the likelihood function, each kinematically cold component (KCC) corresponding to an observed tidal debris feature is assumed to contribute only to the field in which it was observed. The data set does contain observations of large tidal debris features that span multiple spectroscopic fields, most notably the Giant Southern Stream and the Southeast Shelf (both of which are related to a single accretion event; Fardal et al. 2007; Gilbert et al. 2007). However, in all cases, the mean velocity, velocity dispersion, and/or surface density of the feature is spatially dependent, resulting in different measured values in each spectroscopic field. Therefore, each kinematically cold tidal debris feature observed in a spectroscopic field is included as a separate Gaussian component in the likelihood function, with parameters that are fit independently of detections of the same tidal debris feature in other fields.

However, the properties of the underlying dynamically hot stellar halo change much more slowly and smoothly with spatial position than do the presence and properties of tidal debris features. We therefore fit a single underlying M31 halo component across all spectroscopic fields included in a given fit.

This results in field-independent and field-dependent (denoted with the subscript sf) M31 model parameters. Field-independent model parameters include the mean velocity and velocity dispersion of the M31 halo components. These model parameters are present in the M31 likelihood function for every spectroscopic field.

Field-dependent model parameters include the mean velocity and velocity dispersion of each of the M31 tidal debris features present in a field and the relative fractions of each of the M31 components (tidal debris features and halo). These model parameters are present only in the M31 likelihood function for a single spectroscopic field. Our spectroscopic fields have at most three kinematically distinct M31 components: the kinematically hot halo, and two kinematically cold tidal debris features. Therefore, the likelihood function, ${{\mathscr{L}}}_{i,{sf}}^{{\rm{M}}31}$, for a given star i in a spectroscopic field sf takes one of three forms.

All spectroscopic fields without tidal debris features, and thus without any field-dependent model parameters, are described by the M31 likelihood function:

$$\begin{eqnarray}&&{{\mathscr{L}}}_{i}^{{\rm{M}}31}={\mathscr{N}}({v}_{i}| {\mu }_{{\rm{M}}31\mathrm{halo}},{\sigma }_{{\rm{M}}31\mathrm{halo}}).\end{eqnarray} \tag{ 7 }$$

If one KCC has been identified in a field,

$$\begin{eqnarray}{{\mathscr{L}}}_{i}^{{\rm{M}}31} & = & {f}_{\mathrm{KCC}1}{\mathscr{N}}({v}_{i}| {\mu }_{\mathrm{KCC}1},{\sigma }_{\mathrm{KCC}1})\\ & & +{f}_{{\rm{M}}31\mathrm{halo}}{\mathscr{N}}({v}_{i}| {\mu }_{{\rm{M}}31\mathrm{halo}},{\sigma }_{{\rm{M}}31\mathrm{halo}}).\end{eqnarray} \tag{ 8 }$$

If there are two KCCs in a field,

$$\begin{eqnarray}{{\mathscr{L}}}_{i}^{{\rm{M}}31} & = & Z\ {f}_{\mathrm{KCC}1}{\mathscr{N}}({v}_{i}| {\mu }_{\mathrm{KCC}1},{\sigma }_{\mathrm{KCC}1})\\ & & +{f}_{\mathrm{KCC}2}{\mathscr{N}}({v}_{i}| {\mu }_{\mathrm{KCC}2},{\sigma }_{\mathrm{KCC}2})\\ & & +{f}_{{\rm{M}}31\mathrm{halo}}{\mathscr{N}}({v}_{i}| {\mu }_{{\rm{M}}31\mathrm{halo}},{\sigma }_{{\rm{M}}31\mathrm{halo}}).\end{eqnarray} \tag{ 9 }$$

The relative fractions of the M31 components in a spectroscopic field, sf, are normalized such that

$$\begin{eqnarray}&&\displaystyle \sum _{k=1}^{{N}_{{\rm{M}}31,\mathrm{sf}}}{f}_{k}=1.\end{eqnarray} \tag{ 10 }$$

This is enforced by calculating the M31 halo fraction (the only M31 component present in every field) as

$$\begin{eqnarray}&&{f}_{{\rm{M}}31\mathrm{halo}}=1-\displaystyle \sum _{k=1}^{{N}_{\mathrm{KCC}}}{f}_{k}.\end{eqnarray} \tag{ 11 }$$

In fields without any KCCs, the M31 halo fraction (f_{M31 halo}) is equal to one.

4.2.2. MW Model

We characterize the population of MW stars in the sample with three Gaussian components. These include a component with a mean heliocentric line-of-sight velocity near zero and a relatively small velocity dispersion (corresponding to the MW thin disk in the direction of M31), a second component with a slightly more negative mean velocity and a slightly larger velocity dispersion (corresponding to the MW thick disk), and a third component with a significantly more negative mean velocity and a large velocity dispersion (the MW halo). In our data set, the MW halo appears as a population of stars with blue (V − I)₀ colors extending to large negative heliocentric velocities (Figure 3); these stars are significantly bluer than the most metal-poor 10 Gyr RGB isochrone. The MW disk populations span a range of color roughly similar to that spanned by M31 RGB stars, with MW stars that have heliocentric velocities closest to 0 km s⁻¹ in general also having the reddest (I, V − I) colors (e.g., Figure 3 of Gilbert et al. 2012).

Figure 5 shows the velocity histogram of a sample of stars identified as MW stars using a simple cut of $\langle {L}_{i}\rangle$ < −0.5 (Section 3.2), along with the preferred Gaussian mixture model to the velocity distribution. Based on both the Aikake information criterion (AIC) and the Bayesian information criterion (BIC), the preferred number of Gaussian components for this data set is three. Adding additional model components beyond three improves neither the AIC nor the BIC. The best-fit values of the mean velocity, dispersion, and mixture fraction for each component depends on the MW sample selection method (e.g., inclusion of the velocity diagnostic in computing $\langle {L}_{i}\rangle$, and/or inclusion of marginally identified MW stars, with −0.5 ≤ $\langle {L}_{i}\rangle$ < 0). However, the statistical preference for three Gaussian components, with decreasing mean velocity and increasing dispersion, is robust regardless of the MW sample selection: the minimum value in both the AIC and BIC metrics occurs with three Gaussian components for all MW samples tested.

Zoom In Zoom Out Reset image size

The range in Galactic latitude spanned by the spectroscopic fields is substantial (∼20°, Figure 1). The relative fraction of the three MW components is known to change over the full field of view of the survey (e.g., Martin et al. 2013), and it is reasonable to expect that the mean velocity and dispersion of the disk components may change detectably as well. Ideally, the changing mixture of MW stellar populations over the field of the survey (e.g., with galactic latitude) could be included in our model, either empirically or through use of a physical model for the Galaxy.

A fundamental impediment to using a physical MW model, or a parameterized empirical model, is that the MW stars displayed in Figure 5 are not a simple and representative sampling of the MW components along the line of sight to M31. The MW stars with measured velocities are instead a complicated function of the target selection process, which was optimized to preferentially target M31 stars based on the available photometric data in each field, and the success of the cross-correlation routine used to measure velocities (dependent on both the observing conditions for each mask (affecting the S/N of the spectra) and the properties of the star itself (affecting the strength of absorption lines). We nevertheless explored drawing from the Besançon Model of the Galaxy (Robin et al. 2003) at the location of each field, limiting the model results only to those stars in the same apparent magnitude and color ranges as the targets in each field of the spectroscopic survey. We found that the velocity distribution of stars drawn from the Besançon Model does not match the observed velocity distribution of MW stars in our survey.

Given the above considerations, we do not fix the MW component parameters based on fits to a selected sample of MW stars nor attempt a hierarchical fit including Galactic latitude or longitude. Rather, we fit all the parameters for the three MW components simultaneously with the parameters for the M31 components. As was assumed for the properties of the M31 halo, the properties of each of the three MW components is assumed to change relatively slowly and smoothly with spatial position within M31's stellar halo (this is consistent with the Besançon Model of the Galaxy; Robin et al. 2003). Thus, the MW likelihood function is not field dependent and is simply

$$\begin{eqnarray}&&{{\mathscr{L}}}_{i}^{\mathrm{MW}}=\displaystyle \sum _{j=1}^{{N}_{\mathrm{MW}}=3}{f}_{j}{\mathscr{N}}({v}_{i}| {\mu }_{j},{\sigma }_{j}).\end{eqnarray} \tag{ 12 }$$

To ensure a normalized ${{\mathscr{L}}}_{i}^{\mathrm{MW}}$, the sum of the fractions of the MW components is constrained to equal one:

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{{N}_{\mathrm{MW}}=3}{f}_{j}=1.\end{eqnarray} \tag{ 13 }$$

Enforcing this normalization reduces the number of model parameters by one: the fraction of the third MW component is set to

$$\begin{eqnarray}&&{f}_{\mathrm{MW}3}=1-{f}_{\mathrm{MW}1}-{f}_{\mathrm{MW}2}.\end{eqnarray} \tag{ 14 }$$

4.2.3. Probability of M31 Membership

The probability of a given star i being an M31 star (p_M31, Equation (4)) is primarily derived from the overall likelihood $\langle {L}_{i}\rangle$ that the star is an M31 star, which is based on multiple spectroscopic and photometric measurements (excluding the velocity diagnostic, Section 3.2). For each star, $\langle {L}_{i}\rangle$ is the loga-rithm of the odds ratio that the star is an RGB star at the distance of M31 or an MW dwarf star. Thus, the overall likelihood can be converted into a probability of M31 membership via

$$\begin{eqnarray}&&\displaystyle \frac{{10}^{\langle {L}_{i}\rangle }}{1+{10}^{\langle {L}_{i}\rangle }}.\end{eqnarray} \tag{ 15 }$$

The distribution of $\langle {L}_{i}\rangle$ values of stars in a field or radial region is in reality the superposition of two independent and overlapping distributions: that of MW stars, and that of M31 stars. The MW (M31) $\langle {L}_{i}\rangle$ distribution has a tail that extends to positive (negative) $\langle {L}_{i}\rangle$ values (Gilbert et al. 2006, 2007, 2012). When one population (either the MW or M31) is strongly dominant, contamination (in $\langle {L}_{i}\rangle$ space) of the minority population by the tail of the distribution of the dominant population can be significant. For the present analysis, this is most relevant when considering M31's outer halo, where MW stars outnumber M31 stars in the spectroscopic sample by 50:1 or more. If this effect is not accounted for, the tail of the MW star distribution, which has M31-like $\langle {L}_{i}\rangle$ values (greater than zero), results in the fit being driven to include MW stars as part of the M31 model. This in turn will drive the M31 mean velocity to higher values with increasing projected radius.

We therefore introduce a hyperparameter C to the prior probability of M31 membership to account for the possibility that a significant fraction of stars in a field may have $\langle {L}_{i}\rangle$-based probabilities that are not in line with their actual M31 or MW membership. Furthermore, since the density of M31 stars is the single largest driver in the amount of contamination in M31 samples, we add a field-dependent parameter, α_sf, to the fit in fields with identified tidal debris, which typically have higher stellar densities than nearby halo fields without tidal debris features in M31's halo. Thus,

$$\begin{eqnarray}&&{p}_{{\rm{M}}31}={\alpha }_{{sf}}C\displaystyle \frac{{10}^{\langle {L}_{i}\rangle }}{1+{10}^{\langle {L}_{i}\rangle }},\end{eqnarray} \tag{ 16 }$$

and the probability that a star is an MW star is simply 1 − p_M31. In fields without tidal debris features, α_sf ≡ 1.

We restrict α_sf C to be less than or equal to one. This means that it can only reduce the probability of M31 membership based on the empirical photometric and spectroscopic diagnostics. The greatest contrast between M31 and MW stellar density is in the outer M31 halo fields, where MW stars greatly outnumber M31 stars, and improperly categorized MW stars can have a significant effect on the M31 model parameters. In the innermost M31 halo fields included here, the converse is true: M31 stars outnumber MW stars, and improperly categorized M31 stars may affect the MW model parameters. However, the MW model parameters are nuisance parameters that we marginalize over: the MW model parameters have no direct physical interpretation due to our experimental setup (e.g., a varying spectroscopic selection function, averaging over Galactic latitude and longitude when binning the data based on projected radius from M31, and removing M31's bulk motion from the stellar velocities). Moreover, under the assumption that the kinematically hot halo is a well-mixed population, the fact that the tail of M31 stars with $\langle {L}_{i}\rangle$ < 0 contributes little to the M31 model will not introduce a systematic bias in the measurement of the M31 halo parameters.

As discussed above (Section 4.2), there are both field-independent and field-dependent model parameters. We discuss our choice of priors for each of these sets of parameters in turn.

4.3.1. Field-independent Parameters

4.3.2. Field-dependent Parameters

Following Bayes' theorem, we multiplied the likelihood of the data given the model (Section 4.2) with the prior distribution for each parameter in the model (Section 4.3) to compute the likelihood of the model given the data. We used MCMC methods to efficiently sample the parameter space, and marginalized over all model parameters to obtain posterior probability distributions for each parameter of interest.

We sampled the posterior probability distributions for all of the model parameters described above using the open-source emcee python package (Foreman-Mackey et al. 2013a, 2013b), which provides an efficient implementation of the Goodman & Weare (2010) Affine Invariant MCMC Ensemble sampler. The number of parameters in the model is dependent on the number of tidal debris features (Section 4.2), and thus varies based on the spectroscopic fields included in the fit. In order to balance the required computational resources with the final number of independent samples, the number of MCMC chains was set to be at least 10 times, and no more than 20 times, the number of model parameters included in the fit. However, if that value was lower than 300, we instead ran the MCMC analysis with a minimum threshold of 300 chains (Table 3).

Table 3.  Parameters Describing the Velocity Distribution of M31 Halo Stars

Rmin	Rmax	No.	No. M31	No.	No.	Mean	Velocity	Fields with
(kpc)	(kpc)	Stars	Starsa	Model	MCMC	Velocityb	Dispersionb	Substructure
				Parameters	Chains	μM31	σM31
						(km s−1)	(km s−1)
8.2	14.1	617	525	23	460	$-{16.5}_{-4.5}^{+4.3}$	${108.2}_{-6.6}^{+6.8}$	f115, f116, H11
14.1	24	896	697	40	400	$-{20.6}_{-4.2}^{+4.4}$	${98.1}_{-5.0}^{+5.3}$	f207, f123, H13s, f115, f135
24	40	382	240	19	380	$-{18.9}_{-4.4}^{+4.4}$	${98.0}_{-6.6}^{+7.2}$	a3, and9
40	63	589	202	29	436	$-{17.1}_{-4.8}^{+4.9}$	${93}_{-10}^{+11}$	and1, a13, m4
63	90	1068	247	15	300	$-{11.8}_{-4.7}^{+4.7}$	${76.0}_{-7.6}^{+8.6}$	R06A220
90	130	1013	202	11	300	$-{17.4}_{-4.8}^{+4.8}$	${88}_{-10}^{+13}$
130	200	684	104	11	300	$-{17.2}_{-5.0}^{+5.0}$	${92}_{-31}^{+40}$

Notes.

^aNumber of stars more likely to be M31 stars than MW stars based on their $\langle {L}_{i}\rangle$ values, without inclusion of velocity in the calculation of the M31/MW probabilities (Section 3.2). ^bResults are the 50th percentile of the one-dimensional posterior probability distribution, marginalized over all other model parameters. The quoted errors are the 16th and 84th percentiles of the posterior probability distribution. Mean velocities are in the Galactocentric frame, with the bulk motion of M31 removed (Section 4).

Download table as:  ASCIITypeset image

Chains were initialized using a random uniform distribution over the valid parameter space implemented in the prior for each parameter. The chains were run for a minimum of 12,000 steps. The marginalized one- and two-dimensional projections of the posterior probability distributions were computed by drawing the values of the model parameters from the last half of each MCMC chain. The choice to use the first half of the chains as the burn-in period was made as a conservative choice that could be applied uniformly to all chains, for all radial bins. The autocorrelation time of the chains stabilized by half way through the chain for all the parameters considered in the following analysis. In most cases, the autocorrelation time stabilized well before the half-way point.

After the burn-in, the chains for all parameters were inspected to ensure that they had settled to an equilibrium. The parameter values and confidence limits as a function of the number of steps in the chain were also inspected to ensure that stability had been achieved. The chains were run long enough to supply a sufficient number of independent samples (estimated using the autocorrelation length) to estimate the uncertainties on the parameters of interest to a level of ∼1%.

The results of the MCMC analysis presented below define the best estimate for each parameter as the 50th percentile of the one-dimensional posterior probability distribution, marginalized over all other model parameters. The associated uncertainties on the best estimates are the 16th and 84th percentiles of the posterior probability distributions.

To parameterize the change in the velocity dispersion of M31's halo with projected radius, we leverage the one-dimensional posterior probability distributions, i.e., those that are marginalized over all other model parameters. A visualization of the posterior probability distributions for the M31 halo velocity dispersion as a function of radius is shown in Figure 8. In this representation, it is very clear that the halo dispersion in outermost radial bins is poorly constrained: the probability distribution extends over a large range of values and does not have a prominent peak.

Zoom In Zoom Out Reset image size

We randomly sample the marginalized one-dimensional posterior probability distributions of the M31 halo dispersion in each radial bin to parameterize the change in the velocity dispersion of M31's halo with radius. We fit a power law of the form

$$\begin{eqnarray}&&\sigma ={\sigma }_{0}{({R}_{\mathrm{med}}/{R}_{0})}^{-\gamma },\end{eqnarray} \tag{ 17 }$$

where the power law is normalized at the projected radius R₀. We set R₀ = 30 kpc, which is the scale radius chosen for the analysis of the M31 GC dispersion profile in Veljanoski et al. (2014). R_med is the median value of the projected radii of all stars in each radial bin.

We make 10,000 random draws (with replacement) from the one-dimensional posterior probability distributions for the M31 halo velocity dispersion in the first six radial bins and perform a least-squares fit on each draw to determine the power-law parameters σ₀ and γ. The resulting distribution functions of σ₀ and γ provide estimates of the most likely values and uncertainty on each parameter (Figure 9). The 50th percentile values, with uncertainties estimated from the 16th and 84th percentiles of the distributions, are γ = −0.12 ± 0.05 and ${\sigma }_{0}={96.5}_{-3.2}^{+3.3}\,\mathrm{km}\,{{\rm{s}}}^{-1}$. A random drawing of 1000 of these power-law fits is shown in Figure 10. This parameterization confirms the likely existence of a weak gradient in the velocity dispersion with projected radius: only a small percentage of the fits are consistent with a flat or increasing dispersion with radius.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The uncertainties in the power-law slope are estimated based on the method and do not include formal propagation of uncertainties such as those in the velocity measurements or transformations. However, as discussed in Section 5, these uncertainties are estimated to have a minimal effect on the measured velocity dispersions, and this will propagate to a minimal effect on the power-law slope uncertainties. We discuss the effect of modeling choices on the power-law slope below (Section 5.2).

We made several choices regarding both the data and the modeling with the potential to influence the results. We performed alternative fits to explore the impact of these choices on the measured dispersion profile.

We instituted a strong prior probability that stars significantly bluer than the most metal-poor RGB isochrone are MW stars along the line of sight, rather than M31 RGB stars, by assigning them very negative $\langle {L}_{i}\rangle$ values of −5. This places them on the tail end of the distribution of MW stars in $\langle {L}_{i}\rangle$ space, and is consistent with the approach taken in previous papers. However, we also ran the above analysis after assigning less negative $\langle {L}_{i}\rangle$ values to these stars, based on the $\langle {L}_{i}\rangle$ distribution of stars classified as MW stars: they were assigned an $\langle {L}_{i}\rangle$ value 1σ below the mean $\langle {L}_{i}\rangle$ of all MW stars. This still places them firmly as having a high probability of being MW stars, but is a less draconian $\langle {L}_{i}\rangle$ assignment. This resulted in no appreciable difference in the power-law parameters found for the stellar halo dispersion profile (${\sigma }_{0}={96.5}_{-3.2}^{+3.3}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, γ = −0.12 ± 0.05).

We also explored fits with variations on the M31 halo mean velocity parameter, including a delta prior on the mean velocity at M31's systemic velocity (0 km s⁻¹). This also made little difference in the power-law parameters, resulting in values consistent with those found in the nominal fit (${\sigma }_{0}={100.0}_{-3.3}^{+3.4}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, γ = −0.09 ± 0.05). Using a normal prior on the M31 halo mean velocity parameter, but centered at M31's systemic velocity and with a width of σ = 10 km s⁻¹, resulted in a slightly flatter halo velocity dispersion profile, but with power-law parameters statistically consistent (within one sigma) with the nominal fit (${\sigma }_{0}={99.3}_{-3.3}^{+3.5}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, γ = −0.07 ± 0.05).

Relaxing the boundary constraints on the substructure prior (allowing values within ±10 times the uncertainties on the published best-fit values) also did not result in a significant change in the power-law parameters (${\sigma }_{0}={96.3}_{-3.2}^{+3.4}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, γ = −0.11 ± 0.05).

Finally, we have measured the power law based on projected radii calculated using a distance modulus to M31 of 2447 (Section 1). If instead we use a distance modulus of 2438, this results in an ∼4% difference in the estimated projected radii in kpc. The primary effect of assuming a different distance modulus on the estimated power-law slope is easy to compute by simply recalculating the power law after recomputing the median projected radius for each bin. The resulting power-law parameters (σ₀ = 96.1 ± 3.3 km s⁻¹, γ = −0.12 ± 0.05) are fully consistent with the nominal fit.

Only one other measurement of the dispersion profile of stars in M31's halo over a large range of radii has been attempted; this is reported in the paper by Chapman et al. (2006, Section 1). The Chapman et al. fields were primarily within 40 kpc of M31's center and almost exclusively located on M31's disk. This required the application of strict window functions to the velocity distribution to remove M31 disk stars and MW stars. This early analysis could not account for recent discoveries of rotation in M31's inner spheroid (v_rot ∼ 50 km s⁻¹; Dorman et al. 2012) and the relatively large line-of-sight velocity dispersion of disk RGB stars (σ_v = 90 km s⁻¹; Dorman et al. 2015), both of which will affect the interpretation of the relative mix of stellar populations assumed to fall within the Chapman et al. velocity window functions. Nevertheless, the results are consistent with the recent measurements of Veljanoski et al. (2014), who measured the dispersion profile of M31 halo globular clusters (GCs) over a large radial range. The study by Veljanoski et al. (2014) sampled projected radii of 25–145 kpc, and while it includes far fewer tracer objects (72 beyond 30 kpc), there are fewer populations from which these objects are drawn: all are firmly in the M31 system. After accounting for the rotation observed in the GC system, Veljanoski et al. measured a power-law dispersion profile with both a higher value of σ₀ and a steeper power-law slope, −0.45, with only a 1% posterior probability of γ = 0. While Veljanoski et al. find clear evidence of coherent velocity patterns among groups of GCs that are spatially correlated with tidal debris features, the analysis did not account for this. The ability to do so is an advantage of using much more abundant halo RGB stars as tracers of the halo.

Measurements of the velocity dispersion of MW halo stars have found a sharp decrease in the velocity dispersion in the inner regions of the MW halo from ∼150 to 100 km s⁻¹ over the inner 20 kpc, settling to a relatively flat dispersion profile at large radii, with measurements of σ_r ∼ 100 km s⁻¹ from ∼20–80 kpc (Battaglia et al. 2005; Xue et al. 2008; Bond et al. 2010; Brown et al. 2010; Deason et al. 2012; Kafle et al. 2012, 2014). A graphical summary of the MW's velocity dispersion profile can be found in Bland-Hawthorn & Gerhard (2016).

M31 also appears to have a sharply decreasing velocity dispersion in the inner regions (Dorman et al. (2012) measured a velocity dispersion of 140 km s⁻¹ at R_proj = 7 kpc in M31), followed by a relatively flat dispersion to large radii. However, the reader should note that the MW profiles measure primarily the radial velocity of MW halo stars. Given the large spread of the SPLASH spectroscopic fields on the sky, the M31 line-of-sight velocity dispersion profile measures a combination of the stars' tangential and radial velocities in the M31 coordinate frame, with the relative contributions changing with field position.

To date, there have been few analyses of the velocity dispersion profiles of MW- or M31-like stellar halos in ΛCDM simulations (one example is Abadi et al. 2006). The stellar density profiles, substructure characteristics, and metallicity profiles of the M31 and MW halos have proven to be useful constraints and checks on ΛCDM simulations of stellar halo formation, and comparisons of the simulations to observations have provided insight into the physical origins of the stellar halos of M31 and the MW (e.g., Font et al. 2006, 2008; Zolotov et al. 2010; Font et al. 2011; Gilbert et al. 2012, 2014; McCarthy et al. 2012). We expect that future comparisons of the observed MW and M31 velocity dispersion profiles with simulated halos will yield further insights into the origins of stellar halos.