Since the Gaia Data Release 2 (DR2) in April 2018, a number of studies have presented new estimates or lower bounds for the mass of the Milky Way (MW) Galaxy (e.g. Helmi et al., 2018; Malhan & Ibata, 2018; Posti & Helmi, 2018; Sohn et al., 2018b; Watkins et al., 2018; Vasiliev, 2018). Approaches have included maximum likelihood and Bayesian analyses, and have used data for kinematic tracers such as globular clusters (GCs) and stellar streams provided by DR2.
These studies usually interpret results from the inferred gravitational potential profile or the circular velocity profile given an assumed model and the data. Results are often reported as a mass within a specific distance from the Galactic center, or as or the virial mass. While it is useful to have visualizations of the potential and circular velocity profiles in addition to individual estimates of the MW mass within specific distances, it would also be beneficial for studies to present a cumulative mass profile (CMP) of the MW with credible regions. A CMP makes it straightforward to compare mass results reported within different Galactocentric distances, and it can also be compared to the CMP of MW-type galaxies from cosmological, hydrodynamical simulations. Moreover, CMPs resulting from different model assumptions can be compared and used to assess differences in model fits at different distances.
In Eadie et al. (2015) (Paper 1), Eadie & Harris (2016) (Paper 2), and Eadie et al. (2017) (Paper 3), a hierarchical Bayesian method was developed for estimating the total mass and CMP of the MW. This Bayesian model was applied to the MW’s GC population data, using the Harris (2010) catalog and supplemented with proper motion measurements made by many other studies (e.g. Majewski & Cudworth, 1993; Zoccali et al., 2001; Feltzing & Johnson, 2002; Casetti-Dinescu et al., 2010, 2013; Fritz & Kallivayalil, 2015; Rossi et al., 2015, see Paper 2 for a complete table).
In this paper, we provide the method to calculate a CMP given any model for the total gravitational potential and the posterior distribution of the model parameters acquired from a Bayesian analysis. As a motivating example, we employ a simple model for the total gravitational potential and CMP, and confront the model with Gaia DR2 data of GC kinematics to arrive at a new mass estimate and CMP for the MW. For reasons outlined in Section 2, we only use GCs that reside beyond 15kpc. This limits the sample considerably, and also increases the percentage of GCs without proper motion measurements.
The kinematic measurements of GCs released by the Gaia DR2 collaboration do not vastly improve the completeness of tracer data beyond 15kpc, but the HSTPROMO project has helped contribute valuable measurements for GCs at large distances (Sohn et al., 2018b). Additionally, Vasiliev (2018) estimated the mean proper motion of 150 GCs using the Gaia DR2 dataset, thereby increasing the number of proper motion measurements at larger distances considerably. The catalog provided by Vasiliev is appealing not only because the data is more complete, but also because (1) all estimates of GC proper motions are measured in a consistent manner, and (2) the measurements agree very well with the available Gaia and HSTPROMO observations of 20 distant GCs by Sohn et al. (2018b).
Therefore, we estimate the mass and CMP of the MW using first the Gaia collaboration and the HSTPROMO project data, and then again using the catalog presented in Vasiliev (2018). To obtain the posterior distribution of model parameters that will be used to determine the CMP, we use the hierarchical Bayesian method presented in Papers 1–3, with small adjustments in light of the results in Eadie et al. (2018) (hereafter EKH).
The paper is organized as follows:
Section 3 discusses the details and differences between the two data sets (Gaia DR2 + HSTPROMO, and Vasiliev).
Section 4 presents the Bayesian estimates of the MW’s CMP, given each data set, and some discussion and interpretation of results.
Section 5 summarizes our findings and the impact on future studies.
Specific details about our hierarchical Bayesian model and sampling methods for the posterior distribution can be found in Papers 1-3. Here, we provide a brief review for completeness.
The method allows the user to supply (1) kinematic and position data of tracers (such as GCs or halo stars), (2) an analytic distribution function (DF) for the specific energies and angular momenta
of these tracers (based on a gravitational potential and tracer density profile), and (3) hyperprior distributions for the model parameters. The hierarchical Bayesian model includes incomplete and complete velocity data simulataneously, and also accounts for observational error.
The observations of the kinematic tracers are the Galactocentric distance , line-of-sight velocity , and proper motions (,
), and their uncertainties. These measurements are assumed to be samples drawn from normal distributions centered on the tracers’ true but unknown values with standard deviations equal to the measurement uncertainties. In other words, the true tracer velocity components and distances are treated as parameters. This measurement model defines the likelihood in the hierarchical Bayesian model.
The prior distribution on the parameters for the true values of , , and is a distribution function (DF) of the specific energies and angular momenta , given a model for the total gravitational potential and density profile of the tracer population . In Papers 2 and 3, we used a simple power law profile for both the gravitational potential,
and the density profile of the tracer population,
where are parameters. Together, equations 1 and 2, determine the analytic DF first presented in Evans et al. (1997). The DF also includes the constant velocity anisotropy parameter for the tracer population as a free parameter. The derivation of this DF is shown in Evans et al. (1997), with a condensed version shown in Paper 2 (although notations differ).
Hyperprior distributions are also defined for the four model parameters: . In Papers 2–3, the parameters were given uniform prior distributions with upper and lower bounds of , , and respectively, while the prior distribution on
was a Gamma distribution determined by data not used in the analysis.
The posterior distribution is the probability distribution of model parameters given the likelihood, prior, hyperpriors, and the data. Instead of calculating the normalization constant in Bayes’ theorem, samples from a distribution proportional to the posterior distribution are acquired. This is done by running multiple Markov chains in parallel, until they have reached a common, stationary distribution assumed to be proportional to the posterior distribution. Both visual inspection and statistical diagnostics are used to assess the mutual convergence of the chains. The code, calledGME, was written in the R Statistical Software environment (see Papers 1–3).
The hierarchical Bayesian model outlined above was used to estimate the total mass and CMP of the MW, given kinematic data for the MW GC population before Gaia DR2 became available. More recently, EKH tested this method on tracer data from simulated galaxies made in the McMaster Unbiased Galaxy Simulations 2 (MUGS2).
The study of MUGS2 simulated data was done as a strictly blind test. There were two main caveats to the study: (1) the tracer data from MUGS2 were globular cluster analogs (i.e. star tracer particles instead of actual GCs), and (2) half of the galaxies created in the MUGS2 simulations did not follow the standard stellar-mass-to-halo-mass relation. Despite these caveats, the galaxies had the basic components of a bulge, disk, and dark matter halo, and it was therefore instructive to observe how well the single power law gravitational potential model could describe the overall system.
The results of the blind tests were mixed; the total mass () of some galaxies were estimated well within the 95% Bayesian credible regions, while others were not. In particular, there was a tendency for the mass to be underestimated. EKH suggested a number of compounding reasons for the bias (e.g. the physical location of incomplete data and the spatial distribution of tracers in the simulation), but a main factor seemed to be the model for the total gravitational potential. The model had difficulty predicting both the inner and outer regions of the true CMPs for the simulated galaxies when tracers at all Galactocentric radii were included.
Of the galaxies whose masses were poorly estimated, the posterior distribution was pushed towards and reached the edge of allowable parameter space. Posterior distributions that are at an extreme end of parameter space allowed by the prior distribution usually indicate a poor fit to the data, and can be a red flag that the results should be interpreted with caution. Moving forward, we keep this in mind as we apply the method of Paper 3 to the new Gaia data.
EKH also showed that using only the outermost tracers from the simulated galaxies improved the mass estimates and increased the chance of containing the true mass profile within the Bayesian credible regions. These results echoed the sensitivity test in Paper 3 which also used a power law potential. In the latter case, the removal of inner GC data from the analysis caused a slight increase in the mass estimate for the MW. Thus, EKH recommended that in future studies, only outer kinematic tracers be used in the analysis when a single power law model for the gravitational potential is employed. For this study, we use GCs at .
EKH also found that using a uniform distribution as a prior onmay have been too broad. The uniform distibution was centered on 0.5, but also allowed values as low and as high as 0.3 and 0.7. When , Equation 1 approaches a Navarro-Frenk-White-type halo at large distance, and this is the prior information we intended to include. Thus, for this study, we replace the uniform distribution for with a normal distribution with mean and standard deviation (Figure 1).
In this study, we continue to use a Gamma distribution for the
with hyperparameters defined by GC data not used in the analysis because of their lack ofand measurements (Figure 2). Because we are excluding GCs within 15kpc, there are only three GCs rith that are lacking both line-of-sight and proper motion measurements: AM 4, Ko 1, and Ko 2.
2.1 Calculating the CMP with Bayesian credible regions
To calculate Bayesian credible regions for a CMP, given samples from the posterior distribution, perform the following steps:
Create a sequence of closely spaced values for the desired range from the Galactic center.
At each for in , calculate for every sample from the posterior distribution. For each value, you should have values of , where is the number of samples from the posterior distribution (e.g. the length of the Markov chain).
At each , sort the values of
, and calculate the credible intervals of interest (e.g. inner 50%, 75%, and 95%).
In our study, the total gravitational potential assumed in Equation 1 leads to a CMP of the form
(Deason et al., 2012b). Thus, for a Markov chain of length , each pair in the chain is used to calculate the mass within a given radius
. This results in a vector of masses within, which are sorted and used to calculate the marginal distribution of the mass within . We use values logarithmically spaced from to . Our chains are of length 90000, with effective sample sizes of for all parameters.
As part of Gaia DR2, Gaia Collaboration et al. (2018) published kinematic data for 75 GCs in the Milky Way, including proper motions, positions, distances, and line-of-sight velocities111https://www.astro.rug.nl/~ahelmi/research/dr2-dggc/. We use these data to replace the relevant GC measurements in the Harris (2010) catalog, with some exceptions.
The distance measurements to GCs as measured by Gaia show a systematic difference that should improve with the next data release (Gaia Collaboration et al., 2018). In the meantime, we follow the Gaia Collaboration et al. (2018) guidelines and continue to use the distances provided in Harris (2010). The Gaia data are also not entirely complete— only 57 of the 75 GCs have line-of-sight velocity measurements. For those that are missing measurements, we use the ones listed in Harris (2010) when available.
The HSTPROMO team recently provided proper motions for 20 GCs at larger distances (Sohn et al., 2018b). Four GCs measured by HSTPROMO are also measured by Gaia (NGC 362, NGC 2298, NGC 2808, NGC 3201), and in these cases we use the DR2 data. In all other cases, we let the HSTPROMO measurements of GCs supersede any previous measurements in the literature. We then supplement other missing proper motions measurements with those reported in other studies (Majewski & Cudworth, 1993; Zoccali et al., 2001; Casetti-Dinescu et al., 2010, 2013; de la Fuente Marcos et al., 2015; Rossi et al., 2015).
The total number of GCs in our compiled data set is 157, but four GCs — Arp 2, Pal 12, Terzan 7, and Terzan 8 — were recently confirmed by Sohn et al. (2018b) to be associated with the Sagittarius (Sgr) dwarf galaxy. Thus, Sohn et al. (2018b) used only one of these GCs (Arp 2) in their study of the MW’s mass. They also report no significant change in the mass estimate if Pal 12, Terzan 7, or Terzan 8 is used instead. We follow their lead and exclude all but Arp 2 from our analysis. Law & Majewski (2010) listed a few other GCs that may be associated with Sgr, but these have yet to be confirmed and so we leave these GCs in our analysis. Thus, our sample consists of 154 GCs, with a breakdown of complete and incomplete data shown in Figure 3.
Together, the Gaia collaboration and the HSTPROMO team have greatly increased the number of proper motion measurements of GCs in the MW. Only 52 out of 154 GCs in our updated Gaia HSTPROMO catalog have incomplete measurements. This is a significant improvement from the data set used previously (Papers 2–3), in which approximately 50% of the data were incomplete.
Even more recently, Vasiliev (2018) used the Gaia DR2 data to calculate the mean proper motions of 150 GCs in the MW and substantially increase the percentage of complete data. The new Vasiliev catalog is an appealing data set not only because of its completeness, but also because the proper motions are calculated in a consistent manner. Thus, we use this new data set to form a second catalog of GCs and estimate the MW mass and CMP for comparison.
We generate the second catalog from the Vasiliev data such that it contains the same GCs as our GaiaHSTPROMO catalog. We again start with the Harris catalog and then replace all GC measurements of proper motion and line-of-sight velocity with the new measurements found by Vasiliev. The Vasiliev catalog contains some new line-of-sight velocity measurements that were previously missing in the Harris Catalog, and so we adopt these values.
The Vasiliev catalog also contains information for the GC Leavens 1 (Crater), but we exclude this data point because it is not in the GaiaHSTPROMO catalog. Crater is also the most distant GC known to date (Laevens et al., 2014), and its proper motion measurement is quite uncertain (Vasiliev, 2018). Although we do not include the cluster in this study, its distance could make it an important GC for future MW mass estimates once its proper motion is known with more certainty.
Overall, the main difference between the Gaia HSTPROMO and Vasiliev catalogs is completeness, with the latter containing complete measurements for 143 GCs, and incomplete measurements for 11 GCs (Figure 4).
As mentioned in Section 2, we use GCs with Galactocentric distances , thus limiting our sample size to 35 in both our Gaia HSTPROMO and Vasiliev catalogs. In this reduced Gaia HSTPROMO catalog, 16 GCs are missing measurements, three of which are also missing measurements. In contrast, the latter three GCs (AM 4, Ko 1, and Ko 2) are the only incomplete data in the Vasiliev catalog. Because only the distances and positions of AM 4, Ko 1, and Ko 2 are known, these GCs are removed from both catalogs and used to define the prior distribution’s hyperparamters (see Section 2, Figure 2, and Papers 2-3).
Figure 5 shows the proper motion components and line-of-sight velocities as a function of Galactocentric distance for GCs at . The Gaia
HSTPROMO measurements are shown with open blue circles and the Vasiliev’s measurements are grey diamonds. There is excellent agreement between the catalogs’ measurements. The only outlier is themeasurements for Pal 3, which can be seen in the second panel at 95.7kpc.
4 Results and Discussion
Figure 6 shows the CMP of the MW with Bayesian credible regions (50%, 75%, and 95%) given the Gaia + HSTPROMO catalog (left) and the Vasiliev catalog (center). These CMPs were calculated using the posterior distribution of model parameters (Section 2.1), and GCs at . The points with error bars in Figure 6 are mass estimates reported in other studies, which are identified in the legend (Kochanek, 1996; Wilkinson & Evans, 1999; Sakamoto et al., 2003; Battaglia et al., 2005; Xue et al., 2008; Gnedin et al., 2010; Watkins et al., 2010; McMillan, 2011; Deason et al., 2012a, c; Kafle et al., 2012; Gibbons et al., 2014; Eadie et al., 2015; Küpper et al., 2015), including five which also use Gaia DR2 data (Malhan & Ibata, 2018; Posti & Helmi, 2018; Sohn et al., 2018a; Watkins et al., 2018; Vasiliev, 2018).
The two CMPs in Figure 6 are nearly identical, but the Vasiliev catalog provides a better constraint on the mass because the data are more complete. The median estimates of given each catalog are also extremely similar at
where numbers in brackets are the lower and upper bounds of the 95% Bayesian credible region. The quantiles for the model parameters calculated from the posterior distributions are also extremely similar between the two data sets. Henceforth, we refer to the Vasiliev results in the discussion.
A large number of GCs in our catalogs lie between 15 and 20kpc of the Galactic centre (Figure 5), and yet we extrapolate the results out to larger distances with the CMP. To test how sensitive our method is to these inner GCs, we perform the analysis again on the Vasiliev catalog, using a new cut of . Figure 8 compares the marginal posterior distributions and medians for , given values of 15kpc and 20kpc. The median estimate for is for . The new CMP is also virtually unchanged, indicating that the method in insensitive to the tracers between 15kpc and 20kpc.
The main advantage of our CMP for the MW is the ability to easily compare our mass estimate to many other studies both visually (i.e. Figure 6) and quantitatively, within any distance from the Galactic center. To illustrate the quantitative comparisons, we present slices of our CMP taken at through , by steps of 25kpc, to show the marginal posterior distributions for at those distances (Figure 7). Marginal distributions for the mass within any are also easily calculated.
Our CMP agrees with many of the results reported in the literature within their uncertainties, except for three higher-mass estimates in Figure 6: Watkins et al. (2010) (with Draco) at 100kpc, Eadie et al. (2015) at 125kpc, and Vasiliev (2018) at 80kpc. Notably, the first two of these studies included dwarf galaxies in their samples, whereas we do not. There is evidence suggesting that some dwarf galaxies may inflate the mass estimate of the MW. Case in point, the result of Watkins et al. (2010) without Draco is in good agreement with our estimate at 100kpc.
Other studies that use dwarf galaxies, or even high/extreme velocity stars, also tend to infer “heavier” mass estimates of the MW. Many of these studies’ results cannot be included on the CMP (Figure 6) because they are reported as or the virial mass . Instead, we compare three of these studies’ results to our marginal distribution for given the Vasiliev catalog (Figure 8).
Patel et al. (2017) performed a Bayesian analysis in conjunction with cosmological simulations of MW-type galaxies to infer the Galaxy’s mass. They used the kinematics of one massive satellite, and inferred and a 95% Bayesian credible of . Hattori et al. (2018) also deduced a high mass using newly discovered extreme velocity stars in the Gara DR2 data. If these stars are bound to our Galaxy, then they imply . As a final example, Monari et al. (2018) used counter-rotating halo stars in Gaia to estimate the escape velocity of the MW between 5-10kpc. Their results also implied a heavy MW with a dark matter virial mass of . All three of these mass estimates are almost completely ruled out by our estimate of (Figure 8), which did not include high/extreme velocity stars or dwarf galaxies.
It is less clear why the estimate at 80kpc by Vasiliev (2018) is significantly higher than our estimate. Two major differences exist between our studies. The first difference is that they used all GCs in their catalog, whereas we only used GCs beyond 15kpc. The second difference is the choice of gravitational potential model. While we chose a simple power-law model and excluded inner GCs in our analysis, they used a three-component model and included GCs at all distances.
We also compare our results to studies that use GC data from Gaia and/or HSTPROMO. For example, Sohn et al. (2018b) estimated the MW mass within 39.5kpc using 16 GCs with proper motions measurements by HSTPROMO. Utilizing the estimator from Watkins et al. (2010) and Monte Carlo simulations to estimate the uncertainty, they find . Our estimate for the mass within 39.5kpc is
which does not agree with their results within the uncertainties. However, more recently Watkins et al. (2018) combined Gaia DR2 data with HSTPROMO data, and used the same estimator as Sohn et al. (2018b). They found . Reassuringly, the latter agrees with our estimate at this distance within the uncertainties, despite their use of a different model for .
Posti & Helmi (2018) also use the kinematic data of 75 GCs provided by Gaia DR2 and HSTPROMO to estimate the mass of the MW in a Bayesian analysis. In contrast to this study, they use a two-component DF for the GC system of the MW. The assumption is that the GC population dynamics can be decomposed into disk-like and halo-like components, with the latter containing parameters for the mean rotational velocity and velocity anisotropy of the GC system. They also include a shape parameter for the dark matter halo. Because their study uses a Bayesian paradigm, like ours, their method includes measurements uncertainties and incomplete data directly in the analysis. Despite using a more complicated DF than this study, their mass estimate within 20kpc () is within the bounds of our 95% Bayesian credible regions:
Vasiliev (2018) also estimates the mass of the MW using their GC catalog in a maximum likelihood analysis. Their analysis assumes a single-component DF for the GC population, but a three-component gravitational potential to account for the bulge, disc, and halo components. They find The uncertainties almost overlap with the upper bound of the 95% Bayesian credible for our mass estimate within 50kpc:
Finally, Malhan & Ibata (2018) used Gaia data for the GD-1 stellar stream and found , which is in agreement with
given by our method.
Overall, it appears that studies are reaching a consensus on the mass of the MW within 40kpc since the release of Gaia DR2 and HSTPROMO. Even when different analysis methods and models are confronted with the same data, similar conclusions are reached.
Mass estimates of the Galaxy out to larger distances, however, are still rather uncertain. For example, Sohn et al. (2018b) and Watkins et al. (2018) report virial mass estimates of and respectively. Posti & Helmi (2018) provide a tighter constraint on the virial mass at , as does our Bayesian estimate of (Equation 5). Although there is some tension between results, our marginal credible regions for (Figure 8) do overlap with those of Posti & Helmi (2018).
5 Conclusions & Future Work
The state of knowledge about the MW’s total mass is still uncertain, but the state of data for tracers has vastly improved. Multiple studies have now estimated the mass within 40kpc, and in general they are in good agreement with one another. Our CMP of the MW within is also in agreement with many recent results.
Our CMP overlaps more than half of the mass estimates at different radii, suggesting a less massive MW. However, mass estimates that disagree with our results tend to be at larger distances. These discrepancies may be attributed to the inclusion/exclusion of tracers such as dwarf galaxies and/or high velocity stars in the analyses, but it could also be due to differing model assumptions. Thus, the Galaxy’s mass within larger distances is still up for debate. Until more complete measurements for tracers at larger distances are acquired, the total mass of the MW will remain uncertain.
Comparing and verifying mass estimates at large distances will remain challenging without more data, because of the variety of models and data analysis techniques. As a community, we compare results by looking at the uncertainties presented in each others’ studies, with the caveat that everyone’s methods and data are different. We have seen how the Gaia DR2 and HSTPROMO data have greatly improved our understanding of the mass within 40kpc of the MW within this context. Thus, it is imperative to obtain complete velocity measurements of tracers at large distances to reach a similar consensus about the Galaxy’s total mass.
At the same time, studies continue to report mass estimates within different distances from the Galactic center, making comparisons difficult. The CMP presented in this work instead provides a holistic view of our understanding and estimate of the MW’s mass at all distances. Thus, we hope CMPs with Bayesian credible regions become a standard way to present and compare estimates of the MW’s mass in future studies, and openly share instructions to create these profiles.
In summary, future work on the mass of the MW will benefit greatly from complete velocity measurements of tracers at greater distances, and by comparing CMPs from different studies.
This research was funded by an eScience Institute Postdoctoral Fellowship to G. Eadie, made possible from generous donations from the three Moore, Sloan, and Washington Research Foundations at the eScience Institute, University of Washington (UW). G. Eadie and M. Juric acknowledge support from the University of Washington Institute for Data Intensive Research in Astrophysics and Cosmology (DIRAC). The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. M. Juric acknowledges the support of the Washington Research Foundation Data Science Team Chair fund, and the UW Provost’s Initiative in Data-Intensive Discovery. The author would like to thank members of the DIRAC Institute for insightful conversations regarding this work. The author would like to extend graitude to Dr. Vasiliev for sharing their catalog of new proper motions measurements. R Statistical Software Environment(R Development Core Team, 2012), including the following packages: stats (part of base R), CODA: Convergence Diagnosis and Output Analysis for MCMC (Plummer et al., 2006a, b), emdbook: Ecological Models and Data in R (Bolker, 2018, 2008), ggplot2 (Wickham, 2016), limma (vennDiagram function) (Ritchie et al., 2015), MASS: Modern Applied Statistics with S Venables & Ripley (2002), moments: Moments, cumulants, skewness, kurtosis and related tests (Komsta & Novomestky, 2015), pracma: Practical Numerical Math Functions (Borchers, 2017), RColorBrewer: ColorBrewer Palettes, and SNOW: Simple Network of Workstations (Tierney et al., 2013). (Neuwirth, 2014).
- Battaglia et al. (2005) Battaglia, G., Helmi, A., Morrison, H., et al. 2005, MNRAS, 364, 433
- Bolker (2018) Bolker, B. 2018, emdbook: Ecological Models and Data in R, r package version 1.3.10. https://cran.r-project.org/web/packages/emdbook/index.html
- Bolker (2008) Bolker, B. M. 2008, Ecological Models and Data in R (Princeton University Press)
- Borchers (2017) Borchers, H. W. 2017, pracma: Practical Numerical Math Functions, r package version 2.0.7. https://CRAN.R-project.org/package=pracma
- Casetti-Dinescu et al. (2013) Casetti-Dinescu, D. I., Girard, T. M., Jílková, L., et al. 2013, AJ, 146, 33
- Casetti-Dinescu et al. (2010) Casetti-Dinescu, D. I., Girard, T. M., Korchagin, V. I., van Altena, W. F., & López, C. E. 2010, AJ, 140, 1282
- de la Fuente Marcos et al. (2015) de la Fuente Marcos, R., de la Fuente Marcos, C., Moni Bidin, C., Ortolani, S., & Carraro, G. 2015, A&A, 581, A13
- Deason et al. (2012a) Deason, A. J., Belokurov, V., Evans, N. W., & An, J. 2012a, MNRAS, 424, L44
- Deason et al. (2012b) Deason, A. J., Belokurov, V., Evans, N. W., & McCarthy, I. G. 2012b, ApJ, 748, 2
- Deason et al. (2012c) Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012c, MNRAS, 425, 2840
- Eadie & Harris (2016) Eadie, G., & Harris, W. 2016, ApJ, 829, 108
- Eadie et al. (2015) Eadie, G., Harris, W., & Widrow, L. 2015, ApJ, 806, 54
- Eadie et al. (2018) Eadie, G. M., Keller, B. W., & Harris, W. E. 2018, ApJ, 865, 72
- Eadie et al. (2017) Eadie, G. M., Springford, A., & Harris, W. E. 2017, ApJ, 835, 167
- Evans et al. (1997) Evans, N. W., Hafner, R. M., & de Zeeuw, P. T. 1997, MNRAS, 286, 315
- Feltzing & Johnson (2002) Feltzing, S., & Johnson, R. A. 2002, A&A, 385, 67
- Fritz & Kallivayalil (2015) Fritz, T. K., & Kallivayalil, N. 2015, ApJ, 811, 123
- Gaia Collaboration et al. (2018) Gaia Collaboration, Helmi, A., van Leeuwen, F., et al. 2018, ArXiv e-prints, arXiv:1804.09381
- Gibbons et al. (2014) Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2014, MNRAS, 445, 3788
- Gnedin et al. (2010) Gnedin, O. Y., Brown, W. R., Geller, M. J., & Kenyon, S. J. 2010, ApJ, 720, L108
- Harris (2010) Harris, W. E. 2010, astro-ph, arXiv:1012.3224
- Hattori et al. (2018) Hattori, K., Valluri, M., Bell, E. F., & Roederer, I. U. 2018, ArXiv e-prints, arXiv:1805.03194
- Helmi et al. (2018) Helmi, A., van Leeuwen, F., McMillan, P. J., et al. 2018, A&A, 616, A12
- Kafle et al. (2012) Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2012, ApJ, 761, 98
- Kochanek (1996) Kochanek, C. S. 1996, ApJ, 457, 228
- Komsta & Novomestky (2015) Komsta, L., & Novomestky, F. 2015, moments: Moments, cumulants, skewness, kurtosis and related tests, r package version 0.14. https://CRAN.R-project.org/package=moments
- Küpper et al. (2015) Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80
- Laevens et al. (2014) Laevens, B. P. M., Martin, N. F., Sesar, B., et al. 2014, The Astrophysical Journal Letters, 786, L3. http://stacks.iop.org/2041-8205/786/i=1/a=L3
- Law & Majewski (2010) Law, D. R., & Majewski, S. R. 2010, ApJ, 718, 1128
- Majewski & Cudworth (1993) Majewski, S. R., & Cudworth, K. M. 1993, PASP, 105, 987
- Malhan & Ibata (2018) Malhan, K., & Ibata, R. A. 2018, ArXiv e-prints, arXiv:1807.05994
- McMillan (2011) McMillan, P. J. 2011, MNRAS, 414, 2446
- Monari et al. (2018) Monari, G., Famaey, B., Carrillo, I., et al. 2018, A&A, 616, L9
- Neuwirth (2014) Neuwirth, E. 2014, RColorBrewer: ColorBrewer Palettes, r package version 1.1-2. https://CRAN.R-project.org/package=RColorBrewer
- Patel et al. (2017) Patel, E., Besla, G., & Mandel, K. 2017, MNRAS, 468, 3428
- Plummer et al. (2006a) Plummer, M., Best, N., Cowles, K., & Vines, K. 2006a, R News, 6, 7. http://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf
- Plummer et al. (2006b) —. 2006b, CODA: Convergence Diagnosis and Output Analysis for MCMC. https://cran.r-project.org/web/packages/coda/index.html
- Posti & Helmi (2018) Posti, L., & Helmi, A. 2018, ArXiv e-prints, arXiv:1805.01408
- R Development Core Team (2012) R Development Core Team. 2012, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0. http://www.R-project.org/
- Ritchie et al. (2015) Ritchie, M. E., Phipson, B., Wu, D., et al. 2015, Nucleic Acids Research, 43, e47
- Rossi et al. (2015) Rossi, L. J., Ortolani, S., Barbuy, B., Bica, E., & Bonfanti, A. 2015, MNRAS, 450, 3270
- Sakamoto et al. (2003) Sakamoto, T., Chiba, M., & Beers, T. C. 2003, A&A, 397, 899
- Sohn et al. (2018a) Sohn, S. T., van der Marel, R. P., Deason, A., et al. 2018a, in IAU Symposium, Vol. 334, Rediscovering Our Galaxy, ed. C. Chiappini, I. Minchev, E. Starkenburg, & M. Valentini, 47–50
- Sohn et al. (2018b) Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018b, ApJ, 862, 52
- Tierney et al. (2013) Tierney, L., Rossini, A. J., Li, N., & Sevcikova, H. 2013, snow: Simple Network of Workstations, r package version 0.3-13. http://CRAN.R-project.org/package=snow
- Vasiliev (2018) Vasiliev, E. 2018, ArXiv e-prints, arXiv:1807.09775
- Venables & Ripley (2002) Venables, W. N., & Ripley, B. D. 2002, Modern Applied Statistics with S, 4th edn. (New York: Springer), iSBN 0-387-95457-0. http://www.stats.ox.ac.uk/pub/MASS4
- Watkins et al. (2010) Watkins, L., Evans, N., & An, J. 2010, MNRAS, 406, 264
- Watkins et al. (2018) Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2018, submitted to AAS Journals, arXiv:1804.11348v2
- Wickham (2016) Wickham, H. 2016, ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag New York). http://ggplot2.org
- Wilkinson & Evans (1999) Wilkinson, M. I., & Evans, N. W. 1999, MNRAS, 310, 645
- Xue et al. (2008) Xue, X., Rix, H., Zhao, G., et al. 2008, ApJ, 684, 1143
- Zoccali et al. (2001) Zoccali, M., Renzini, A., Ortolani, S., Bica, E., & Barbuy, B. 2001, AJ, 121, 2638