Estimating the Mass of the Milky Way Using the Ensemble of Classical Satellite Galaxies

Satellite galaxies around the Milky Way (MW) are often used to study the structure and assembly history of the Galaxy's dark matter and stellar halo. In particular, the kinematics of halo tracers (satellites, globular clusters, and stellar streams) have been used to constrain the Galaxy's gravitational potential and total mass. Many efforts have been made to constrain the current factor of two spread in the total mass range of the MW, but a high precision estimate has yet to be made.

Leo I, as one of the highest speed MW satellites, is often used to place lower limits on the mass of the MW. However, because its relative velocity hovers around the MW's escape speed at its current separation of ∼260 kpc, any MW mass constraint requires the assumption that Leo I is bound. Boylan-Kolchin et al. (2013, hereafter BK13) illustrate that unbound satellites within the virial radius are rare. By assuming that Leo I is bound at its current position and velocity, BK13 used the Aquarius cosmological zoom simulations (Springel et al. 2008) to infer that the virial mass of the MW must be >10¹² ${M}_{\odot }$, with a median mass of 1.6 × 10¹² ${M}_{\odot }$. Other studies relying on the boundedness of Leo I, such as Li and White (2008), used a radial timing argument analysis similar to that of Zaritsky et al. (1989) to determine that the MW's mass is 2.43 × 10¹² ${M}_{\odot }$ with a lower limit of 0.8 × 10¹² ${M}_{\odot }$. Non-radial timing argument studies predict a MW mass as high as >3 × 10¹² ${M}_{\odot }$ if Leo I is bound (Sohn et al. 2013).

If Leo I is unbound, however, these methods would likely overestimate the true mass of the MW. The case of Leo I strongly motivates a study to calibrate the impact of using a single satellite to constrain the mass of its host. We argue in this study that a single satellite can result in a significantly biased mass estimate, especially when satellites are on extreme orbits (see also Sales et al. 2007; Patel et al. 2017b). Instead, we have developed a novel method of estimating the mass of the MW using an ensemble of observed satellites and analogs in cosmological simulations.

The HST Proper Motion collaboration (Kallivayalil et al. 2006a, 2006b, 2013; Massari et al. 2013; Sohn et al. 2013, 2017, e.g.,) and other authors (Scholz & Irwin 1994; Piatek et al. 2016, 2008, 2007, 2006, 2005, 2003, 2002; Walker et al. 2008) have now measured the proper motions of several low mass dwarf galaxies in the MW's halo. With the combined 6D phase space information derived from the proper motions for the classical MW satellites (LMC, SMC, Carina, Draco, Fornax, Sculptor, Leo I, Leo II, Ursa Minor, Sextans, and the Sagittarius dSph), we can now use the 3D dynamics of a population of halo tracers to further our understanding of the MW's dark matter halo and its global properties. In this work, we utilize the magnitude of specific orbital angular momentum of each satellite, rather than their instantaneous position and velocity, as motivated in Patel et al. (2017b, hereafter P17B). An important benefit is that we make no assumption of whether a satellite is bound—simply whether it currently resides within the virial radius of the MW.

As the two most massive MW satellites, the Magellanic Clouds (MCs) are also often used to characterize properties of the MW, specifically via analogs selected from cosmological simulations. Both Boylan-Kolchin et al. (2011a) and Busha et al. (2011) have used the instantaneous characteristics of the MCs, such as their current position and velocity, to make predictions for the MW's mass in a frequentist and Bayesian fashion, respectively (see also González et al. 2013). Using the Kallivayalil et al. (2013) proper motions, the Bayesian posterior mean mass estimate for the MW using the Large Magellanic Cloud (LMC) alone is 1.70 × 10¹² ${M}_{\odot }$ (P17B) and approximately 1.15 × 10¹² ${M}_{\odot }$ using both MCs (González et al. 2013).

While the MCs are well-studied members of the MW's halo, their current properties, including the spatial proximity of the MCs to the MW (∼50 kpc) and to each other (∼23 kpc), and their unique orbital configuration (just past pericenter), make them rare in a cosmological context (Besla et al. 2007, 2012; Boylan-Kolchin et al. 2011a). Less than 5% of simulated MW mass halos host two massive satellites that have made a recent first pericentric approach as close as 50 kpc. We showed in P17B that MW mass estimates are highly susceptible to the orbital phase of the LMC, and thus conclusions based on the current properties of only the MCs should be taken with caution. As a result, we must turn to other satellite properties to estimate the mass of the MW in an unbiased fashion.

The momentum method developed in P17B relies on orbital angular momentum and thus affords a larger simulated data set from which MW mass estimates can be formed, since it does not limit satellite analogs to a narrow range of position and velocity combinations as in the instantaneous method. For the MW-LMC system, a nearly tenfold increase in the number of statistically significant satellite analogs from the dark-matter-only Illustris simulation (or Illustris-1-Dark) enables us to combine inferences for the MW's mass even for rare systems (i.e., unique orbital configurations), thus providing a powerful method moving forward as more high precision data becomes available for halo tracers.

The impact of only one satellite on the mass estimate of the MW is clearly evident with regard to Leo I. This motivates the need to use several satellites simultaneously to form a MW mass estimate. Here, we extend the P17B method to all classical low mass MW satellites (satellites less massive than the MCs) to constrain the mass of the MW. This will test whether satellites like Leo I are still outliers in these MW mass estimation techniques, or if other satellites will become more critical players. For example, recent studies by Gibbons et al. (2014), Belokurov et al. (2014) have suggested a lower mass bound for the MW based on the properties of the Sagittarius dSph stellar stream, indicating a lower limit of about 7 × 10¹¹ ${M}_{\odot }$. With this statistical framework, such assertions can be tested in a cosmological context for the first time.

In Section 2, we compile properties of the nine low mass satellites used in this study and describe the details of the Illustris dark-matter-only cosmological simulation. Section 3 provides the details of the Bayesian framework that has been extended from P17B to accommodate a population of lower mass satellites. The results of these methods are presented in Section 4. Section 5 contains a discussion of the global trend between the distribution of satellite specific orbital angular momenta and host halo mass. Finally, Section 6 summarizes our results.

In P17B, we combined the properties of the most massive satellite galaxies in the halos of the MW and M31(the LMC and M33) with a state-of-the-art cosmological simulation in a Bayesian statistical framework to infer the masses of their host halos. In this work, we extend our previous analysis to lower mass satellites and additionally implement a method that uses multiple satellites simultaneously to refine MW mass estimates. This section describes the observational data that are currently available for the low mass MW satellites considered and the specifications of the Illustris-1-Dark simulation, which we will use to select low mass satellite analogs.

2.1. Observed Properties of Nine Low Mass MW Satellites

To date, the proper motions of the nine brightest MW dwarf satellite galaxies (after the MCs) have been measured with high astrometric precision in the last decade. In principle, any halo tracer (e.g., satellite galaxies, globular clusters) with proper motion information can be used in this type of statistical analysis to estimate the mass of its host if a sufficiently large set of simulated analogs are also available.

The classical⁵ low mass satellites considered in this work include Carina, Draco, Fornax, Sculptor, Leo I, Leo II, Ursa Minor, Sextans, and the Sagittarius dSph. The stellar masses of these dwarf spheroidal (dSph) galaxies are between 10⁵ and 10⁷ ${M}_{\odot }$ (see McConnachie 2012), and their masses at infall are predicted to be 10⁸–10¹⁰ ${M}_{\odot }$ from cosmological expectations (Boylan-Kolchin et al. 2012; Moster et al. 2013). For reference, the LMC's stellar mass is about 3 × 10⁹ ${M}_{\odot }$ at present, and its total mass at infall could be as high as a few times 10¹¹ ${M}_{\odot }$ (Kim et al. 1998; Wang et al. 2006; Besla et al. 2010).

The observed data (${\boldsymbol{d}}$) for each low mass satellite and the LMC are provided in Table 1. These data consist of: (i) the observed position relative to the MW (r^obs), (ii) the total velocity relative to the MW (${v}_{\mathrm{tot}}^{\mathrm{obs}}$), and (iii) the magnitude of the specific orbital angular momenta about the MW (j^obs), where $j=| {\boldsymbol{r}}\times {\boldsymbol{v}}|$. The mean values of these quantities and the uncertainties associated with them have been calculated from a set of 10,000 Monte Carlo samples drawn from the 4σ error space of distance, radial velocity, and proper motion of each satellite. For the low mass satellites, all adopted distances and radial velocities are taken from the compilation presented in McConnachie (2012) and the references therein. The proper motions of the galaxies come from a variety of groups and programs (most measurements come from the Hubble Space Telescope), as indicated in the final column of Table 1. In this analysis we do not use the properties of the SMC explicitly, due to the low frequency of LMC and SMC analogs in cosmological simulations (see Boylan-Kolchin et al. 2011a), but SMC analogs may exist about the halos included in the prior sample (see Section 3.1).

Table 1.  Observational MW Satellite Data

	robs	1σ	${\text{}}{v}_{\mathrm{tot}}^{\mathrm{obs}}$	1σ	jobs	1σ	References
	(kpc)	(kpc)	(km s−1)	(km s−1)	(kpc km s−1)	(kpc km s−1)	Dist. + vrad	PM
Leo II	236	14	138	42	32,105	9,998	1	4
Fornax	149	12	193	45	28,471	8,712	1	5
Leo I	258	15	202	19	27,918	9,146	1	3
Sculptor	86	6	213	10	17,232	1,267	1	2
Sextans	89	4	206	22	17,097	2319	1	8
LMC	50	5	321	24	15,688	1,788	1,10	10
Ursa Minor	78	3	173	49	12,545	4,375	1	7
Draco	76	6	183	4	12,241	1,229	1	2
Carina	107	6	97	41	10,322	4,532	1	6
Sagittarius dSph	18	2	320	23	5,249	1,000	1	9

Note. Observational data (${\boldsymbol{d}}$) for the LMC and the nine MW dSph galaxies with proper motions (PMs). The satellites are listed in order of descending j^obs. Data for the distance, radial velocity, and PMs of each satellite were taken from the following references.

References (1) McConnachie (2012, and references therein), (2) Sohn et al. (2017), (3) Sohn et al. (2013), (4) Piatek et al. (2016), (5) Piatek et al. (2007), (6) Piatek et al. (2003), (7) Piatek et al. (2005), (8) Casetti-Dinescu et al. (2018), (9) Sohn et al. (2015), (10) Kallivayalil et al. (2013).

Download table as:  ASCIITypeset image

2.2. The Illustris-Dark Simulation

In this section, we identify the subset of host halos from P17B's prior that have exactly one massive satellite analog and at least one satisfactory low mass subhalo analogous to the classical satellites. Note that while the MW actually hosts several low mass satellites, we only require a minimum of one low mass satellite analog per prior sample host halo due to the simulation's resolution limits. We also outline the modified selection criteria $({\boldsymbol{C}}^{\prime} )$ for the prior sample, likelihood functions tailored for low mass satellites, and the statistical approximation used to infer the mass of the MW using the ensemble of low mass satellites simultaneously.

3.1. Prior Samples

3.1.1. Massive Satellite Analogs

While this work focuses on the low mass MW satellites, it is important to include the presence of a massive satellite analog (weighing approximately 10% of their host's total mass), since they could alter the gravitational potential of their host halos and subsequently affect the kinematics of other nearby satellites (see Section 5). The existence of a massive satellite analog also ensures that MW analogs with similar mass assembly histories are chosen. We also note that the subhalo abundances around host halos with and without a massive satellite analog differ, which we discuss in detail in Appendix B.

To build a prior sample, we first require that exactly one massive satellite analog satisfying the following selection criteria (${\boldsymbol{C}}$) from P17B exists in all halos from the final 20 snapshots (z ≈ 0.26 to z = 0) of the Illustris-Dark simulation. We consider only the final 20 snapshots to be consistent with P17B and previous work (e.g., Busha et al. 2011; González et al. 2013).

• C₁: A subhalo is considered a massive satellite analog only if v_max > 70 km s⁻¹.
• C₂: The massive satellite analog must reside within its host's virial radius (R${}_{\mathrm{vir}}$) at z ≈ 0.
• C₃: The massive satellite analog must have a minimal subhalo mass of 10¹⁰ ${M}_{\odot }$ at z ≈ 0.

Note that subhalo mass is provided in units of ${M}_{\odot }$ h⁻¹, where h = 0.704, in the Illustris-Dark halo catalogs, and is used as given. v_max is the maximum circular velocity of a subhalo. These selection criteria return a total of 19,653 host halos, each with a companion massive satellite analog. This data constitutes the prior sample used in P17B.

3.1.2. Low Mass Satellite Analogs

Analogs of the low mass classical MW satellites must belong to the first 15 subhalos in each FoF group, which are ranked by decreasing subhalo mass. This ensures that systems with multiple massive satellite analogs (see previous section) are omitted. Ideally, the prior sample would contain only those halos that host a massive satellite analog in addition to about 10 low mass satellite analogs. However, simulations with large volumes and sufficient resolution are currently unable to resolve a statistically significant number of halos that represent true analogs of the MW in this way, so these constraints are relaxed. Small scale cosmological volumes with higher mass resolution and perhaps finer redshift spacing are likely the most ideal simulations to build such a prior sample in the future.

We require that all low mass satellite analogs have v_max < 45 km s⁻¹ to ensure that we are truly selecting just the low mass subhalos (i.e., the MCs have v_max > 50 km s⁻¹). To avoid the "Too Big to Fail" (TBTF) problem (Boylan-Kolchin et al. 2011b), only an upper v_max bound is used in these criteria. TBTF is a discrepancy that arises in dissipation-less ΛCDM simulations where subhalos with masses analogous to the classical MW dwarfs fail to host bright satellites. This discrepancy leads to a mismatch between simulated v_max values and their observed counterparts such that the observed values are much lower than that of their simulated analogs.

Our goal in this work is to choose dynamical analogs, not analogs based purely on mass or energetics. Since dynamical friction plays an insignificant role in the orbital evolution of these low mass satellites, the kinematics (positions, velocities, specific angular momenta) are not expected to decrease drastically between infall and today. The selection criteria for low mass satellite analogs (denoted by C') are summarized below.

• ${{\rm{C}}}_{1}^{{\prime} }$: A subhalo is considered a low mass satellite analog only if v_max < 45 km s⁻¹.
• ${{\rm{C}}}_{2}^{{\prime} }$: The satellite analog must reside within its host's virial radius (R${}_{\mathrm{vir}}$) at z = 0.
• ${{\rm{C}}}_{3}^{{\prime} }$: The satellite analog must have a subhalo mass ≥10⁹ ${M}_{\odot }$ at z = 0.

A minimum subhalo mass of 10⁹ ${M}_{\odot }$ corresponds to ∼133 dark matter particles in Illustris-Dark.

3.1.3. Prior 1

Prior 1 is the sample of host halos that host exactly one massive satellite analog (criteria ${\boldsymbol{C}}$, Section 3.1.1) and one or more low mass satellite analogs (criteria ${\boldsymbol{C}}^{\prime}$, Section 3.1.2). Approximately 92.8% of host halos from the prior sample used in P17B satisfy these combined criteria. This subset of 18,236 host halos and the properties of their associated massive satellite analogs will be referred to as Prior 1 from here on. Note that this prior is only trivially different from that selected in P17B. The LMC's properties will be used with Prior 1 to compute the mass of the MW.

3.1.4. Prior 2

We find 87,598 subhalos that satisfy the low mass satellite analog criteria in a total of 18,236 unique host halos. Thus, many host halos harbor more than one low mass satellite analog, as expected. This set of host halos and the properties of their associated low mass satellite analogs will be referred to as Prior 2. Prior 2 will be used to compute the mass of the MW using the properties of all classical MW satellites less massive than the MCs.

Figure 1 shows the distribution of latent properties ${\boldsymbol{x}}$ for a select set of subhalos in Prior 2. We have plotted the properties of subhalos residing in the bottom, middle, and top 5% of the host halo mass distribution. Each panel indicates a pair of two parameters in ${\boldsymbol{x}}=[r,{v}_{\mathrm{tot}},j]$ colored by the corresponding host halo mass. Overplotted are the observed properties ${\boldsymbol{d}}=[{r}^{\mathrm{obs}},{v}_{\mathrm{tot}}^{\mathrm{obs}},{j}^{\mathrm{obs}}]$ for all MW satellites listed in Table 1. The overlap between the observed properties and prior sample properties shows that the low mass satellite analogs accurately represent the true MW satellite properties. Without implementing any statistical techniques, the current properties of the classical MW satellites are most similar to those subhalos residing in host halos with masses <1.4 × 10¹² ${M}_{\odot }$.

Zoom In Zoom Out Reset image size

The greatest number of low mass satellite analogs (by definition) hosted by any individual host halo in Prior 2 is 14 subhalos. Only 0.05% of all halos considered in the prior sample host exactly 14 low mass satellite analogs. About 6% of all host halos have 10 or more satisfactory subhalos within their virial radii. Typically halos host between 2 and 5 low mass satellite analogs. See Appendix B for more information on the abundance of subhalos within the virial radii of host halos in Prior 2.

The covariance between measurement errors of r^obs and ${v}_{\mathrm{tot}}^{\mathrm{obs}}$ is small, as shown in Busha et al. (2011), so we treat these two measurements independently in the instantaneous method, which uses these two satellite properties to compute likelihoods for all halos in the prior sample. We consider the magnitude of specific orbital angular momenta separately in a parallel analysis, which we refer to as the momentum method. The most ideal candidates for our analysis are those residing at 100–300 kpc from the center of mass of their host halos, since they are least affected by strong tides from the host (see Appendix A for more details).

3.1.5. Multiplicity of Host Halos

In Prior 2, all host halos with more than one low mass satellite analog are counted toward the prior distribution of host halo masses multiple times. For example, if a host halo has a virial mass of ${\text{}}{M}_{\mathrm{vir}}$ = 1.23 × 10¹² ${M}_{\odot }$ and hosts four low mass satellite analogs, its virial mass will appear four times in the list of halo masses that correspond to low mass satellite analogs. In P17B (and subsequently Prior 1), there is a one to one relation between host halos and massive satellite analogs since that prior was limited to hosts with just one massive satellite analog.

3.2. Likelihood Functions for Low Mass Satellites

For those host halos and low mass satellite analogs included in Prior 2, the physical parameters of interest are ${\boldsymbol{\theta }}=\{{\boldsymbol{x}},{M}_{\mathrm{vir}}\}$, where ${\text{}}{M}_{\mathrm{vir}}$ is the virial mass of the corresponding host halo for any given subhalo. The parameters ${\boldsymbol{x}}$ are the latent, observable properties for satellite subhalos in Illustris-Dark. Figure 1 illustrates these for a fraction of the subhalos in the prior. The observational data $({\boldsymbol{d}})$ listed in Table 1 are measurements of the parameters ${\boldsymbol{x}}$, so if the measurement errors are zero, then ${\boldsymbol{d}}={\boldsymbol{x}}$. Subsets of the physical parameters ${\boldsymbol{x}}$ are used in the different likelihood functions to estimate the mass of the MW.

The likelihood functions from P17B are altered such that they no longer rely on ${v}_{\max }^{\mathrm{obs}}$ due to the TBTF problem. We implement two methods to compute the likelihood of a given MW mass, as given below.

Instantaneous

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{x}}| \,{\boldsymbol{d}})=N({r}^{\mathrm{obs}}| \,r,{\sigma }_{r}^{2})\times N({v}_{\mathrm{tot}}^{\mathrm{obs}}| \,{v}_{\mathrm{tot}},{\sigma }_{v}^{2})\end{eqnarray} \tag{ 1 }$$

Momentum

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{x}}| \,{\boldsymbol{d}})=N(\,{j}^{\mathrm{obs}}| \,j,{\sigma }_{j}^{2}).\end{eqnarray} \tag{ 2 }$$

All measured satellite properties are assumed to have Gaussian error. The posterior distribution for the virial halo mass of the MW using the properties of each low mass satellite is computed using Prior 2 and these likelihood functions via importance sampling.

3.2.1. Importance Sampling

Bayes's theorem is written as

$$\begin{eqnarray}&&P({\boldsymbol{x}},{M}_{\mathrm{vir}}| \,{\boldsymbol{d}},{\boldsymbol{C}}^{\prime} )\propto P({\boldsymbol{d}}| \,{\boldsymbol{x}})\times P({\boldsymbol{x}},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} ),\end{eqnarray} \tag{ 3 }$$

where we denote the dependence on the prior selection criteria ${\boldsymbol{C}}^{\prime}$. The left-hand side is the posterior probability distribution. $P({\boldsymbol{x}},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )$ represents the prior probability distribution, and $P({\boldsymbol{d}}| \,{\boldsymbol{x}})$ is the likelihood (equivalent to ${ \mathcal L }({\boldsymbol{x}}| {\boldsymbol{d}})$).

A posterior probability density function (PDF) is then calculated by drawing a set of samples (n) from an importance sampling function. Here, we have chosen the importance sampling function to be the prior PDF (as in P17B). For each sample in n, an importance sampling weight proportional to the likelihood (using either Equations (1) or (2)) is assigned. Using these weights, integrals that summarize the posterior PDF for halo mass are calculated as follows, where the denominator is a normalization constant.

$$\begin{eqnarray}&&\displaystyle \int f({\boldsymbol{\theta }})P({\boldsymbol{x}},{M}_{\mathrm{vir}}| \,{\boldsymbol{d}},C^{\prime} )\,d{\boldsymbol{\theta }}\\ &&\quad =\,\displaystyle \frac{\displaystyle \int f({\boldsymbol{\theta }})P({\boldsymbol{d}}| \,{\boldsymbol{x}})P({\boldsymbol{x}},{M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )\,d{\boldsymbol{\theta }}}{\displaystyle \int P({\boldsymbol{d}}| \,{\boldsymbol{x}})P({\boldsymbol{x}},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )\,d{\boldsymbol{\theta }}}\\ &&\quad \approx \,\displaystyle \frac{{\displaystyle \sum }_{j}^{n}f({{\boldsymbol{\theta }}}_{j})P({\boldsymbol{d}}| \,{{\boldsymbol{x}}}_{j})}{{\displaystyle \sum }_{j}^{n}P({\boldsymbol{d}}| \,{{\boldsymbol{x}}}_{j})}.\end{eqnarray} \tag{ 4 }$$

If f(θ) only depends on ${\text{}}{M}_{\mathrm{vir}}$, a representation of the marginal posterior PDF for ${\text{}}{M}_{\mathrm{vir}}$ can be derived by computing Equation (5) over a grid of potential halo mass values (i.e., using kernel density estimation).

$$\begin{eqnarray}&&\displaystyle \int f({M}_{\mathrm{vir}})P({M}_{\mathrm{vir}}| \,{\boldsymbol{d}},{\boldsymbol{C}}^{\prime} )\,{{dM}}_{\mathrm{vir}}\\ &&\quad =\,\displaystyle \int \displaystyle \int f({M}_{\mathrm{vir}})P({\boldsymbol{x}},{M}_{\mathrm{vir}}| \,{\boldsymbol{d}},{\boldsymbol{C}}^{\prime} )\,d{\boldsymbol{x}}\,{{dM}}_{\mathrm{vir}}\\ &&\quad \approx \,\displaystyle \frac{{\displaystyle \sum }_{j}^{n}f({M}_{\mathrm{vir}}^{j})\,P({\boldsymbol{d}}| \,{{\boldsymbol{x}}}_{j})}{\displaystyle \sum _{j}^{n}P({\boldsymbol{d}}| \,{{\boldsymbol{x}}}_{j})}\\ &&\quad =\,{\displaystyle \sum }_{j}^{n}f({M}_{\mathrm{vir}}^{j})\,{w}_{j}\end{eqnarray} \tag{ 5 }$$

where ${w}_{i}=P({\boldsymbol{d}}| \,{{\boldsymbol{x}}}_{i})/{\sum }_{j}^{n}P({\boldsymbol{d}}| \,{{\boldsymbol{x}}}_{j})$ are importance weights. Setting $f({\boldsymbol{\theta }})={\text{}}{M}_{\mathrm{vir}}$ gives the posterior mean value for virial halo mass of the MW. For more details on the importance sampling technique, see Section 3.2 of P17B.

Note that in practice all calculations are carried out in log₁₀${\text{}}{M}_{\mathrm{vir}}$ (by directly replacing ${\text{}}{M}_{\mathrm{vir}}$ with log₁₀${\text{}}{M}_{\mathrm{vir}}$ in all equations) because the posterior distribution of log₁₀${\text{}}{M}_{\mathrm{vir}}$ is more roughly Gaussian than the posterior distribution of ${\text{}}{M}_{\mathrm{vir}}$, and therefore more easily summarized by a central value. Consequently, all results reported on a physical scale as ${\text{}}{M}_{\mathrm{vir}}={X}_{-L}^{+U}$ ${M}_{\odot }$ should be interpreted on a log scale. For example, log₁₀X is the posterior mean of log₁₀${\text{}}{M}_{\mathrm{vir}}$ and [log₁₀(X − L), log₁₀(X + U)] is the 68% credible interval in log₁₀${\text{}}{M}_{\mathrm{vir}}$. These summaries should not be naively translated to constraints on a linear scale, as probability densities do not trivially transform under a nonlinear change of variables (see Jensen's inequality). The MW mass inferred by each satellite is discussed in Section 4.1.

3.3. A Statistical Approximation to Include Several Low Mass Satellites Simultaneously

Thus far, we have only outlined how to infer the virial mass of the MW using the properties of any individual satellite galaxy for which the proper motion has been measured. While these individual estimates are interesting on their own, such an analysis leads to the question, What if the phase space information of all satellites is used simultaneously to infer the mass of the MW? One might expect that additional information from multiple satellites will yield a more precise MW mass.

Below we outline a statistical approximation to simultaneously infer the MW's mass using the properties of all nine low mass satellites. An approximation is necessary to make a combined MW mass estimate, since Prior 2 does not exclusively include only halos with approximately 10 low mass satellites each (see Section 3.1). To consider an arbitrary number of satellites with measurements $\{{{\boldsymbol{d}}}_{s}\}$, Bayes's theorem yields the joint posterior of host halo virial mass ${{\boldsymbol{M}}}_{\mathrm{vir}}$:

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \,\{{{\boldsymbol{d}}}_{s}\},{\boldsymbol{C}}^{\prime} )\\ &&\qquad =\,\displaystyle \frac{P(\{{{\boldsymbol{d}}}_{s}\}| \,\{{{\boldsymbol{x}}}_{s}\})\times P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )}{P(\{{{\boldsymbol{d}}}_{s}\}| \,{\boldsymbol{C}}^{\prime} )},\end{eqnarray} \tag{ 6 }$$

where s = 1... N_sat and N_sat is the total number of satellites considered. Because the measurements of each satellite are independent from the others, the likelihood can be written

$$\begin{eqnarray}&&P(\{{{\boldsymbol{d}}}_{s}\}| \,\{{{\boldsymbol{x}}}_{s}\})=\displaystyle \prod _{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{d}}}_{s}| \,{{\boldsymbol{x}}}_{s}).\end{eqnarray} \tag{ 7 }$$

Next, we note that the prior factor can be written, using the definition of conditional probability,

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )=P({M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )\times P(\{{{\boldsymbol{x}}}_{s}\}| \,{M}_{\mathrm{vir}},{\boldsymbol{C}}^{\prime} ).\end{eqnarray} \tag{ 8 }$$

If we make the naive Bayes assumption that, given the mass M_vir, the satellites properties $\{{{\boldsymbol{x}}}_{s}\}$ are conditionally independent, we have

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\}| \,{M}_{\mathrm{vir}},{\boldsymbol{C}}^{\prime} )=\displaystyle \prod _{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{x}}}_{s}| \,{M}_{\mathrm{vir}},{\boldsymbol{C}}^{\prime} ),\end{eqnarray} \tag{ 9 }$$

and therefore

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )=P({M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )\times \displaystyle \prod _{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{x}}}_{s}| \,{M}_{\mathrm{vir}},{\boldsymbol{C}}^{\prime} ).\end{eqnarray} \tag{ 10 }$$

Putting this all together, we have

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \,\{{{\boldsymbol{d}}}_{s}\},{\boldsymbol{C}}^{\prime} )\\ &&=\,\displaystyle \frac{[{\displaystyle \prod }_{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{d}}}_{s}| \,{{\boldsymbol{x}}}_{s})P({{\boldsymbol{x}}}_{s}| \,{M}_{\mathrm{vir}},{\boldsymbol{C}}^{\prime} )]\times P({M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )}{P(\{{{\boldsymbol{d}}}_{s}\}| \,{\boldsymbol{C}}^{\prime} )}.\end{eqnarray} \tag{ 11 }$$

The individual prior factors can be written as

$$\begin{eqnarray}&&P({{\boldsymbol{x}}}_{s}| \,{M}_{\mathrm{vir}},{\boldsymbol{C}}^{\prime} )=\displaystyle \frac{P({{\boldsymbol{x}}}_{s},{M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )}{P({M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )},\end{eqnarray} \tag{ 12 }$$

so the joint posterior becomes

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \{{{\boldsymbol{d}}}_{s}\},{\boldsymbol{C}}^{\prime} )\\ &&=\,\displaystyle \frac{[{\displaystyle \prod }_{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{d}}}_{s}| {{\boldsymbol{x}}}_{s})P({{\boldsymbol{x}}}_{s},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )]\times P{({M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )}^{1-{N}_{\mathrm{sat}}}}{P(\{{{\boldsymbol{d}}}_{s}\}| {\boldsymbol{C}}^{\prime} )}.\end{eqnarray} \tag{ 13 }$$

Next we notice that the posterior of given the data for a single satellite s is

$$\begin{eqnarray}&&P({{\boldsymbol{x}}}_{s},{M}_{\mathrm{vir}}| \,{{\boldsymbol{d}}}_{s},{\boldsymbol{C}}^{\prime} )=\displaystyle \frac{P({{\boldsymbol{d}}}_{s}| \,{{\boldsymbol{x}}}_{s})P({{\boldsymbol{x}}}_{s},{M}_{\mathrm{vir}}| \,{\boldsymbol{C}}^{\prime} )}{P({{\boldsymbol{d}}}_{s}| \,{\boldsymbol{C}}^{\prime} )},\end{eqnarray} \tag{ 14 }$$

so the joint posterior is

$$\begin{eqnarray}&&P(\{{{\boldsymbol{x}}}_{s}\},{M}_{\mathrm{vir}}| \,\{{{\boldsymbol{d}}}_{s}\},{\boldsymbol{C}}^{\prime} )=\left[\displaystyle \prod _{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{x}}}_{s},{M}_{\mathrm{vir}}| \,{{\boldsymbol{d}}}_{s},{\boldsymbol{C}}^{\prime} )\right]\\ &&\quad \times \,P{({M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )}^{1-{N}_{\mathrm{sat}}}\times \displaystyle \frac{{\displaystyle \prod }_{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{d}}}_{s}| \,{\boldsymbol{C}}^{\prime} )}{P(\{{{\boldsymbol{d}}}_{s}\}| \,{\boldsymbol{C}}^{\prime} )}.\end{eqnarray} \tag{ 15 }$$

Integrating out $\{{{\boldsymbol{x}}}_{s}\}$, we find the marginal posterior for the mass, given all of the satellites' data.

$$\begin{eqnarray}&&P({M}_{\mathrm{vir}}| \,\{{{\boldsymbol{d}}}_{s}\},{\boldsymbol{C}}^{\prime} )=\left[\displaystyle \prod _{s=1}^{{N}_{\mathrm{sat}}}P({M}_{\mathrm{vir}}| \,{{\boldsymbol{d}}}_{s},{\boldsymbol{C}}^{\prime} )\right]\\ &&\times \,P{({M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )}^{1-{N}_{\mathrm{sat}}}\times \displaystyle \frac{{\displaystyle \prod }_{s=1}^{{N}_{\mathrm{sat}}}P({{\boldsymbol{d}}}_{s}| \,{\boldsymbol{C}}^{\prime} )}{P(\{{{\boldsymbol{d}}}_{s}\}| \,{\boldsymbol{C}}^{\prime} )}\end{eqnarray} \tag{ 16 }$$

The last factor is a normalization constant that does not depend on the parameters. Since we already know how to calculate $P({M}_{\mathrm{vir}}| \,{{\boldsymbol{d}}}_{s},{\boldsymbol{C}}^{\prime} )$ for one satellite at a time using Equation (5), a useful expression for the marginal posterior is

$$\begin{eqnarray}&&P({M}_{\mathrm{vir}}| \,\{{{\boldsymbol{d}}}_{s}\},{\boldsymbol{C}}^{\prime} )\propto \left[\displaystyle \prod _{s=1}^{{N}_{\mathrm{sat}}}P({M}_{\mathrm{vir}}| \,{{\boldsymbol{d}}}_{s},{\boldsymbol{C}}^{\prime} )\right]\\ &&\quad \times \,P{({M}_{\mathrm{vir}}| {\boldsymbol{C}}^{\prime} )}^{1-{N}_{\mathrm{sat}}}.\end{eqnarray} \tag{ 17 }$$

By writing the marginal posterior distribution in this form, the multiplicity from including Prior 2 N_sat times is eliminated. Again, all calculations are computed in log₁₀${\text{}}{M}_{\mathrm{vir}}$ as noted in Section 3.2. Results for the combined MW mass estimates using this statistical approximation are discussed in Section 4.2.

3.4. The Conditional Independence Assumption and Computing the Joint Posterior Distribution

The conditional independence assumption states that the properties of a subhalo in a given host halo are independent from those of another subhalo in the same host halo for each given value of host halo mass. Here, we demonstrate the validity of this assumption for our statistical framework.

The conditional independence assumption requires that the correlation based on $P({j}_{\mathrm{sat},1},{j}_{\mathrm{sat},2}| {M}_{\mathrm{vir}})$, for example, is zero for all ${\text{}}{M}_{\mathrm{vir}}$, where j_{sat, 1} and j_sat,2 are the values of total specific orbital angular momentum for the first two low mass satellite analogs in each host halo (where at least two analogs exist). In Figure 2, we plot j_{sat, 1} versus j_sat,2 for the fraction of Prior 2 as four subsamples split by host halo mass. The assumption should hold for all values of host halo mass. The Pearson correlation coefficient (r) between j_sat,1 and j_sat,2 for each subsample yields values between 0.12 and 0.27 with uncertainties <0.07.

Zoom In Zoom Out Reset image size

Physically, this weak correlation shows that there is reasonable scatter among satellite specific orbital angular momenta for the same host halo. This suggests that satellites' angular momenta vectors are not set by the strength of the large scale tidal field but rather by more complex processes such as varying accretion and orbital histories, even for satellites orbiting the same host. Such values suggest there is only a weak 2-point correlation and therefore that the conditional independence is a useful and reasonable approximation. The conditional independence assumption actually requires the N-dimensional joint distribution factors as

$$\begin{eqnarray}&&P({j}_{\mathrm{sat},1},{j}_{\mathrm{sat},2}...,\,{j}_{\mathrm{sat},{\text{}}{\rm{N}}}| {M}_{\mathrm{vir}})\\ &&=\,P({j}_{\mathrm{sat},1}| {M}_{\mathrm{vir}})\times P({j}_{\mathrm{sat},2}| {M}_{\mathrm{vir}})\times ...\,\times \,P({j}_{\mathrm{sat},{\text{}}{\rm{N}}}| {M}_{\mathrm{vir}}).\end{eqnarray} \tag{ 18 }$$

However, we have only demonstrated this for N = 2 in Figure 2, as it is difficult to rigorously show this for N > 2.

Calculating Equation (17) directly can sometimes lead to numerical underflow or overflow errors. To prevent this, we compute the logarithm of Equation (17), apply it to our data, and exponentiate to retrieve the final results. This strategy is successful in a majority of cases unless the product of all satellite kernel density estimate (KDE) posteriors and the KDE estimate of the prior in Equation (17) both approach zero in a numerically unstable way. In such cases, the limiting edge effects must be carefully considered to produce a stable result. Such numerical caveats will be unnecessary when more advanced high resolution simulations are available.

We now infer the mass of the MW using both the instantaneous and momentum likelihood methods with the properties of each individual dwarf satellite and the Illustris-Dark cosmological simulation. While we already demonstrated in P17B that the momentum method is more reliable as a function of time and satellite orbital phase, we will report results from both likelihood functions for comparison. We also include the results for the ensemble MW mass estimates using all classical satellites.

In what follows, we provide mass estimates of the MW using the LMC's orbital properties and Prior 1. However, this mass estimate cannot be combined with those resulting from the lower mass satellite galaxies, as the two sets of results are computed from two different prior samples (Prior 1 versus Prior 2).

4.1. MW Mass Estimates from Individual Low Mass Dwarf Satellites

4.1.1. Instantaneous Likelihood

The top left panel of Figure 3 shows the MW mass estimates and associated uncertainties for each low mass satellite using the instantaneous likelihood method. The inferred MW masses are plotted against the total velocity of each satellite relative to the MW today (Table 1, Column 4). Error bars indicate the MW's posterior mean mass included in the 68% credible interval of the posterior distribution for host halo mass. The colored lines in the bottom left panel of Figure 3 show the corresponding posterior distributions. All results are also listed in Table 2. The LMC is indicated by a black square in Figure 3.

Zoom In Zoom Out Reset image size

Table 2.  Summary Statistics for MW Mass Estimates

	Instantaneous	Momentum
	${\text{}}{M}_{\mathrm{vir}}$ [1012 M⊙]	${\text{}}{M}_{\mathrm{vir}}$ $[{10}^{12}\,{M}_{\odot }]$
LMC	1.45+1.75−0.41	1.0+0.79−0.56
Leo II	1.52+0.58−0.59	1.54+1.14−0.7
Fornax	1.11+0.84−0.45	1.44+1.11−0.63
Leo I	2.02+0.60−0.69	1.41+1.13−0.62
Sculptor	1.12+0.55−0.39	1.14+0.94−0.51
Sextans	1.03+0.64−0.35	1.12+0.93−0.50
Ursa Minor	0.80+0.65−0.33	0.96+0.76−0.48
Draco	0.91+0.57−0.34	0.92+0.70−0.45
Carina	0.66+0.46−0.33	0.91+0.76−0.45
Ensemble (8 sats)	1.19+0.19−0.21	0.96+0.29−0.28
Sagittarius dSph	⋯	0.78+0.76−0.39
Ensemble (9 sats)	⋯	0.85+0.23−0.26

Note. The posterior mean and 68% credible interval in halo mass for the MW using the properties of each of the nine low mass satellites and the LMC. The MW halo mass computed using host halos in Prior 2 and the properties of eight (nine) low mass satellites simultaneously are given in the third to last row. The MW mass derived using the LMC was calculated with Prior 1 (first row). The LMC is excluded from the ensemble mass estimates, since these calculations use fundamentally different prior samples (Prior 1 versus Prior 2). The SMC's properties are not considered due to the low frequency of LMC-SMC analog pairs in Illustris-Dark.

Download table as:  ASCIITypeset image

The results of the instantaneous method display a fairly large scatter (standard deviation of ∼0.14 dex) among the posterior mean mass estimates for the MW using each individual satellite. The lowest value for the MW's mass comes from Carina, which suggests a posterior mean MW mass of only 0.66 × 10¹² ${M}_{\odot }$. As expected, the highest MW mass estimate (2.02 × 10¹² ${M}_{\odot }$) with the instantaneous method comes from Leo I. This mass is consistent with the bound mass arguments given by BK13 and the mass inferred by the radial timing argument, but lower than that given by the non-radial timing argument (Li & White 2008; Sohn et al. 2013). Carina is on a low energy, fairly circular orbit about the MW (Pasetto et al. 2011), whereas Leo I is on a high energy orbit and has recently completed its first pericentric approach (Sohn et al. 2013). The difference in orbital energies suggests a strong correlation between energy and halo mass.

As any one individual satellite can produce an MW mass estimate that could be misleading given its orbital energy, this warrants combining as many satellites as possible to infer the mass of the MW in tandem. The instantaneous method essentially reproduces expectations for the MW's mass from traditional methods, such that Leo I pushes the MW's mass to higher values and satellites on more circular orbits that experience little to no dynamical friction, preferring lower MW halo masses.

The Sagittarius dSph has been omitted from the results for the instantaneous method, as the effective sample size reduces to just a few during the importance sampling step, given its unique combination of position and velocity. Upon examining the cumulative number density of subhalos in Prior 2 at any given distance relative to their respective host, we find that the minimum host-satellite separation is ∼30 kpc. Therefore, very few subhalos are found at the separation of Sagittarius (20 kpc) in Illustris-Dark. This problem may be twofold, as the gravitational smoothing length of Illustris-Dark may not be able to sufficiently resolve distinct halos at these small separations and there may also be a depletion of subhalos in the inner regions of halos due to tidal disruption, especially for subhalos on radial orbits (Garrison-Kimmel et al. 2017).

While only ∼20% of the known MW satellites (e.g., Drlica-Wagner et al. 2015) are considered in this analysis, the individual estimates using the instantaneous method already demonstrate a factor of three scatter in MW halo mass, even larger than that in the literature to date. The instantaneous method is not recommended if unbiased MW mass estimates are desired.

4.1.2. Momentum Likelihood

MW mass estimates for the momentum method are plotted in the top right panel of Figure 3 against the current magnitude of specific orbital angular momentum for each satellite.The bottom right panel shows the corresponding posterior distributions. In general, we see a similar trend among the low mass satellites as we did for the LMC and M33 analysis in P17B—satellites that have a higher total specific orbital angular momentum value (see Table 1) infer higher MW masses. The same trend is not strictly true for ${v}_{\mathrm{tot}}^{\mathrm{obs}}$ (Figure 3, top left panel) or r^obs. Overall, the scatter among individual MW mass estimates using the momentum method is far better constrained to a range of 0.78 × 10¹² ${M}_{\odot }$ (Sagittarius dSph)—1.54 × 10¹² ${M}_{\odot }$ (Leo II), in agreement with the current factor of two spread in MW mass. Excluding Sagittarius dSph, the scatter of our results is narrowed even further to 0.91–1.54 × 10¹² ${M}_{\odot }$. Comparing the standard deviation across individual estimates from the instantaneous method and the momentum method gives ∼0.14 dex versus ∼0.09 dex, clearly demonstrating the improvement that the momentum method provides if satellites are used individually.

In P17B, we demonstrated that using instantaneous properties like position and velocity skew the resulting MW mass estimates because they change significantly with time. The specific orbital angular momentum vector of Sagittarius dSph has likely undergone significant changes since its infall into the halo of the MW, due to dynamical friction and tidal stripping which lead to the subsequent formation of the Sagittarius stellar stream (see Belokurov et al. 2014). It is therefore unsurprising that Gibbons et al. (2014) and our analysis yield such low MW masses using the Sagittarius dSph (or its associated stream). The MW mass inferred by the Sagittarius dSph using our statistical framework should therefore be taken with caution.

Contrary to previous studies, Leo II, rather than Leo I, puts limits on the MW's mass. Leo II has the highest specific orbital angular momentum of the satellites considered in this analysis, so it brackets the upper end of the MW's plausible mass range.

Satellites on the most extreme orbits, where they have only made one pericentric passage about the MW (Leo I or the LMC), also show the biggest deviations between the posterior mean mass estimates using the instantaneous versus the momentum method. The fact that the MW mass estimates vary so drastically for these satellites provides further evidence that the combination of position and velocity is not reliable for recovering the mass of the MW accurately, and cautions against using one satellite alone to infer the MW's mass.

We conclude that the momentum method not only is a more consistent estimator of host halo mass as a function of time and orbital evolution but is also a more consistent method for determining the mass of a host halo given the 6D phase space information for a population of satellite galaxies. The momentum method directly correlates inferred host halo mass and satellite specific orbital angular momentum, thereby distinguishing between halos that can host a population of satellites exhibiting a given observationally constrained distribution of specific orbital angular momenta and those that cannot. The correlation between the distribution of specific orbital angular momentum for a population of satellites and host halo mass will be discussed further in Section 5.

While the momentum method for low mass satellites only considers one physical property associated with each satellite—the magnitude of the specific orbital angular momentum—it encompasses some 6D phase space information of each satellite too as $j=| {\boldsymbol{r}}\times {\boldsymbol{v}}|$, where ${\boldsymbol{r}}$ and ${\boldsymbol{v}}$ are 3D vectors. The direction of the specific orbital angular momentum method is not utilized in this method, as given. In principle, such information could be incorporated into the momentum likelihood function, but this is beyond the scope of this work.

4.1.3. Caveats for the Individual MW Mass Estimates

In P17B, we carried out an additional bootstrap analysis to address the sampling noise associated with our technique. This is especially important for satellites in unique orbital configurations, like the LMC and the Sagittarius dSph. Since Prior 2 is significantly larger than that considered in P17B, the effective sample size does not rapidly decrease to zero for most of the classical MW satellites (except Sagittarius dSph) considered. We also showed in Section 5.2 of P17B that smaller measurement errors can improve host halo mass estimates, but when the precision is already as low as a few percent (as is the case for the LMC, Sculptor, Draco), the change in results is insignificant. The Gaia mission will be able to reduce the measurement errors on proper motions and derived quantities for the nearest classical satellite galaxies (van der Marel & Sahlmann 2016), further narrowing MW mass estimates using those satellites.

The barrier to achieving very high precision MW mass estimates for the individual satellites is the irreducible uncertainty owing to cosmic variance (Section 5.3, P17B), or the intrinsic correlation between host halo mass and satellite dynamics. We expect that the most significant improvements to the precision of MW mass estimates will therefore arise from using the properties of several satellites simultaneously.

4.2. MW Mass Estimates from the Ensemble of Low Mass Dwarf Satellites

In Section 4.1.2, we note the strong correlation between the specific orbital angular momenta distribution of the observed satellite population and the MW mass inferred by those satellites. Naturally this leads to the question, Is there a correlation between host halo mass and the distribution of specific orbital angular momentum for a population of subhalos? We use simulated analogs from Illustris-Dark to explore if this intrinsic relationship exists and how it might inform our knowledge of the MW's true mass.

5.1. The Halo Mass–Specific Angular Momenta Distribution Relation in the Presence of a Massive Satellite

Figure 4 shows the distribution of total specific orbital angular momenta for subhalos in Prior 2 (Section 3.1.4) binned by their corresponding host halo mass. The specific orbital angular momenta shown are the median values for each host halo mass bin, and the error bars indicate the extents of the 25th and 75th quartiles. The gray shaded area shows the median and quartiles of the specific orbital angular momentum distribution for the population of classical MW satellites considered in this work.

Zoom In Zoom Out Reset image size

The high concentration of subhalos residing in halos with masses in the range of ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\approx 11.8\mbox{--}12.3$ suggests that this is the most typical MW halo mass independent of our results from Section 4. This mass range is in agreement with the frequentist MW mass predictions from Patel et al. (2017a) using the energetics of LMC analogs in Illustris-Dark. Note that this region is also coincident with the gray shaded area, which represents the true distribution of orbital angular momentum for eight of the classical MW satellites. The green dashed line and shaded region represent the posterior mean and 68% credible interval included in the ensemble MW mass estimate from the momentum method (see Table 2, last row). The black solid line is the line of best fit with a slope of m = 0.779 ± 0.0280 and an intercept of b = −5.116 ± 0.339 in units of dex (i.e., the fits are calculated using log₁₀(j_sat) and ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })$).

When we further explore if the specific orbital angular momenta-host halo mass trend is correlated with the orbital history of the massive satellite analogs, we find that the linear relationship shown in Figure 4 is generally unaffected. The subhalos residing in halos whose massive satellite analogs have crossing times⁶ less than 4 Gyr ago (first infall scenarios) and those hosting massive satellite analogs with crossing times more than 4 Gyr ago (multiple passages about their hosts) exhibit slopes in agreement with the best fit line in Figure 4 (black solid line).

Using Figure 4, we can also make predictions for the specific orbital angular momenta of ultra-faint dwarf galaxies. The red box, which encompasses the region of highest subhalo abundance (darkest blue squares), our MW mass estimate (green shaded region), and the observed distribution in specific orbital angular momentum of satellites (gray shaded region) warrants the prediction that ultra-faint dwarfs should generally exhibit 5 × 10³ kpc km s⁻¹ < j^obs < 5 × 10⁴ kpc km s⁻¹ if the MW's true mass is ∼10¹² ${M}_{\odot }$. Upcoming observations of ultra-faint dwarf galaxies will be some of the most interesting in this regard, given their expected abundance.

5.2. The Halo Mass–Specific Angular Momenta Distribution Relation in the Absence of a Massive Satellite

We now explore whether the absence of a massive satellite analog changes the satellite specific angular momentum distribution as a function of host halo mass. To do so, we choose a new sample including all host halos in the final five snapshots of Illustris-Dark that pass the low mass satellite analog selection criteria (C') given in Section 3.1.2 without requiring that a massive satellite analog (Section 3.1.1) also resides in the halo. Thus, the host halos included in Prior 1 are strictly excluded from this sample. The new sample contains 179,381 low mass satellite analogs residing in 104,362 unique host halos.

Figure 5 shows the distribution of subhalo specific orbital angular momenta binned by halo mass for this alternate set of Illustris-Dark subhalos. These subhalos reside in halos with masses extending down to ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\approx 9.3$ due to the absence of a massive satellite analog, roughly two orders of magnitude lower than in Figure 4. The gray shaded region is identical to that shown in Figure 4 and represents the range of specific orbital angular momenta for eight low mass classical MW satellites. The green shaded region represents the posterior mean and 68% credible interval included in the ensemble MW mass estimate from the momentum method calculated using this new sample (i.e., no massive satellite analogs).

Zoom In Zoom Out Reset image size

The black solid line is the line of best fit from Figure 4, and it indicates that the overall relationship between halo mass and specific orbital angular momentum distribution (i.e., the slope of the distribution) still exists. However, at very low host halo masses, the corresponding orbital angular momentum distribution begins to deviate from this trend such that the median specific orbital angular momentum is higher than if a massive satellite analog is included. This is important because it means that massive satellites (i.e., ∼10% of their host's mass) may affect the orbital histories and kinematics of lower mass satellites, but it is harder to see explicitly in the higher host halo masses because a wide range of massive satellite analogs was included in the analysis. We conclude that there is a relationship between these two quantities only in host halos with masses ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\gt 11.5$, whether or not a massive satellite analog (i.e., LMC-mass companion) is present.

A few additional points of interest are also worth noting upon comparing Figures 4 and 5. The abundance of low mass satellite analogs as a function of host halo mass (depicted by the color bars) is highest in halos with ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\approx 11.8\mbox{--}12.3$ in Figure 4 when massive satellite analogs are present, compared to a range of log₁₀(${\text{}}{M}_{\mathrm{vir}}$/${M}_{\odot }$) ≈ 10.6–12.1 in Figure 5 when they are not present. Such differences are expected from the hierarchical assembly of galaxy halos and are probably closely linked to the subhalo abundance functions for halos with and without a massive satellite companion (see Appendix B). Despite these differences, the analysis still yields narrow MW mass ranges that are in agreement with each other (green shaded regions). This suggests there is no bias associated with the presence or absence of a massive satellite analog and that our prior selection criteria are justified.

The existence of a strong correlation between halo mass and the distribution of orbital angular momenta for subhalos provides a generalized trend that can be applied to halos and their subsequent satellite populations. While the MW system is often considered common in a cosmological context, recent work by Geha et al. (2017) cautions against thinking about the MW as a typical halo in their study of observational MW analogs and their satellite populations' abundance and properties. Our notes on this trend are therefore widely applicable.

Building on the Bayesian framework used to estimate the MW's virial mass from P17B, we provide a new methodology to combine the 6D phase space information of multiple low mass dwarf satellites. By doing so, we are able to estimate the mass of the MW using each of the following satellites individually and as an ensemble: Fornax, Sculptor, Carina, Draco, Leo I, Ursa Minor, Sextans, Leo II, Sagittarius dSph. This method is straightforwardly generalizable to include more dwarf satellites as new data is obtained. The main conclusions of this work are summarized below.

• 1.  
  When the mass of the MW is inferred using each low mass satellite individually, we find a much larger scatter in the mass of the MW using the instantaneous method. This method yields a MW mass range of 0.6–2 × 10¹² ${M}_{\odot }$ across eight low mass satellites (Sagittarius dSph excluded), while the momentum method spans a much narrower range from 0.9 to 1.5 × 10¹² ${M}_{\odot }$ for eight satellites. The latter results are more constrained than the current values in the literature.
• 2.  
  By combining the posterior distributions associated with each satellite, we also report ensemble MW mass estimates. The combination of eight low mass satellites with the instantaneous method results in a MW mass of ${M}_{\mathrm{vir},\mathrm{MW}}={1.19}_{-0.21}^{+0.19}\times {10}^{12}$ ${M}_{\odot }$ (${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\,={12.08}_{-0.09}^{+0.06}$). Using the momentum method, the ensemble mass estimate is ${M}_{\mathrm{vir},\mathrm{MW}}={0.96}_{-0.28}^{+0.29}\times {10}^{12}$ (${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })={11.98}_{-0.15}^{+0.11}$). If Sagittarius dSph is included from the latter estimate, it decreases to ${M}_{\mathrm{vir},\mathrm{MW}}={0.85}_{-0.26}^{+0.22}\times {10}^{12}$ ${M}_{\odot }$ (${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\,={11.93}_{-0.16}^{+0.11}$). These ensemble estimates are more precise than the masses inferred by any individual satellite, and they narrow the mass range to less than a factor of two.
• 3.  
  Satellites with high specific orbital angular momentum like Leo II, as opposed to high speed satellites like Leo I or the LMC, now constrain the upper end of the MW's mass range. Satellites in the midst of disruption like the Sagittarius dSph have the lowest specific orbital angular momentum and therefore push the MW's mass to low values. As any one satellite alone can result in a biased MW mass estimate, satellite galaxies that are not undergoing disruption (evidenced by tidal tails or streams) are the most ideal candidates for this method.
• 4.  
  The 68% credible interval of the ensemble MW mass resulting from the momentum method when the Sagittarius dSph is excluded ranges from 0.67 to 1.25 × 10¹² ${M}_{\odot }$. This suggests Leo I could be bound to the MW if its true mass is in the upper half of this range (since ${\text{}}{M}_{\mathrm{vir}}$ = 10¹² ${M}_{\odot }$ corresponds to R${}_{\mathrm{vir}}$ ≈ 260 kpc). Our analysis only includes a fraction of the MW's satellite galaxy population, and while our MW mass estimates are not expected to change drastically as more satellites are included, these results are still preliminary.
• 5.  
  We find there is a linear relationship between host halo mass and the distribution of specific orbital angular momentum for subhalos at a given halo mass. This trend is independent of the presence of massive satellite analogs and their orbital histories (i.e., first infall versus multiple orbits). The presence of a massive satellite analog does, however, change the halo mass range in which low mass subhalos are likely to reside. For example, when we require all halos to harbor a massive satellite analog, subhalos are most abundant in halos with masses ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\approx 11.8\mbox{--}12.3$. When massive satellites are not present, this range broadens and shifts down to ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\approx 10.6\mbox{--}12.1$. The median specific orbital angular momentum of satellites in lower mass host halos is higher if massive satellite analogs are not present, suggesting a massive satellite that is a high fraction of the host halo mass can affect the kinematics of the low mass satellite population. While this trend brackets the orbital angular momenta distribution expected at a given host halo mass, our tests on the conditional independence assumption shows that the individual angular momentum vectors of satellites around a shared host are independent of one another.
• 6.  
  Upon comparing this trend to the distribution of specific orbital angular momenta for MW satellites with measured proper motions, we predict that ultra-faint MW dwarfs will have specific orbital angular momenta values between 5 × 10³ and 5× 10⁴ kpc km s⁻¹. Future proper motion measurements for MW dwarfs (HST-GO-14734, P.I. N. Kallivayalil) and specifically the ultra-faints residing at ∼100 kpc (HST-GO-14236, PI: S. T. Sohn) will allow us to test this hypothesis.

We have shown that combining a majority of the available satellite phase space information has already narrowed the plausible mass range for the MW, and eventually the same can be done for M31 and other nearby galaxies in the era of JWST. As more satellite information becomes available, significant improvements to current cosmological simulations are crucial so a large prior sample of MW analogs, each hosting tens of satellites, can be selected. Together, a larger high precision data set for MW substructures and more advanced simulations will be powerful tools for converging on a precise and accurate MW mass.

E.P. is supported by the National Science Foundation through the Graduate Research Fellowship Program funded by Grant Award No. DGE-1746060. E.P. and G.B. are supported by NSF Grant AST-1714979. K.M. was supported in part by NSF grant AST-1516854. S.T.S. is supported by NASA through grant GO-12966. This research was also funded through a grant for HST program AR-12632. Support for AR-12632 and GO-12966 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The Illustris simulations were run on the Odyssey cluster supported by the FAS Science Division Research Computing Group at Harvard University. We are also grateful to Roeland van der Marel, Laura Watkins, Mark Fardal, Dennis Zaritsky, Vicente Rodriguez-Gomez, Tom McClintock, Shea Garrison-Kimmel, and Hans-Walter Rix for useful discussions that have contributed to this paper.

Software: Example code developed for this project is available from https://github.com/ekta1224/BayesToolsSatDynamics under the MIT open-source software license. This research also utilized IPython (Perez & Granger 2007), numpy (van der Walt et al. 2011), scipy (Jones et al. 2001), and matplotlib (Hunter 2007).

The presence of a galactic disk is known to cause a depletion of subhalos in the inner regions of MW-like halos (e.g., D'Onghia et al. 2010). More specifically, recent work shows that subhalos on radial orbits are most susceptible to this phenomena (Garrison-Kimmel et al. 2017). Here we use the Illustris-1 N-body+hydrodynamical simulation to examine whether there are any significant changes to the intrinsic relationship seen between host halo virial mass and the specific orbital angular momenta of subhalos discussed in Section 5 when baryons are considered. By comparing the results from the two versions of the Illustris simulation, we can also demonstrate whether the resulting MW mass estimates change significantly when both dark matter and baryons are taken into account.

To this end, a new sample of low mass satellite analogs are chosen from the final 20 snapshots of Illustris-1 that satisfy the same exact selection criteria described in Section 3.1. This sample contains 59,156 subhalos residing in 18,858 unique host halos. Notice that one-third fewer subhalos satisfy our selection criteria in Illustris-1 compared to Illustris-Dark, indicating that there is a general depletion of low mass subhalos when disks and other hydrodynamic processes, such as feedback, are included.

Figure 6 shows the distribution of specific orbital angular momenta for subhalos in Illustris-1 binned by host halo virial mass. The gray shaded region is identical to the previous figures and represents the distribution of specific orbital angular momenta for the classical MW satellites. Figure 6 is generally in good agreement with Figures 4 and 5, as indicated by the consistency between the data points and the black solid line (the line of best fit from the dark-matter-only counterpart sample) for host halo masses above ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })=11.75$. The deviation from the line of best fit at host halo masses ${\mathrm{log}}_{10}({M}_{\mathrm{vir}}/{M}_{\odot })\lt 11.75$ is likely due to the depletion of the lowest specific orbital angular momenta satellites and suggests there is a difference in subhalo kinematics between the simulations. The green shaded region indicates the ensemble MW mass estimate when the Illustris-1 data is used. Notice that the latter is lower than that predicted by the dark-matter-only subhalos, indicated by the red dotted line.

Zoom In Zoom Out Reset image size

We have also split the prior sample from Illustris-1 into subsets based on the time of their most recent major merger, where this is defined as the last time there was a 1:4 stellar mass ratio collision such that the Illustris halo-finding algorithm can no longer distinguish between two distinct halos. Hosts with major mergers greater than 6 Gyr ago make up 79% of the sample and the halo mass-specific orbital angular momentum trend for this subset is in good agreement with Figure 6. The scatter associated with various formation histories is therefore implicitly included in our mass estimates, even though our prior selection criteria do no explicitly require any specific formation history. Complementary work by Li et al. (2017) suggest that this scatter can contribute to up to 20% uncertainty to the MW mass estimates.

We conclude that on average, the addition of baryonic physics and other hydrodynamical processes results in lower MW halo mass values. However, the posterior mean MW halo masses using each individual satellite and all satellites simultaneously are still within 2σ of their counterpart MW posterior mean halo masses predicted by the dark-matter-only analysis. It appears that such discrepancy could arise from the depletion of low mass subhalos, which may consequently change the correlation between subhalo properties and their host halo properties. However, a new generation of cosmological simulations with revised feedback formulae and perhaps higher particle mass resolution may be able to address such concerns more confidently.

Since the subhalos considered here all have subhalo masses ≥10⁹ ${M}_{\odot }$ and most have a typical v_max = 20–45 km s⁻¹, their overall depletion is less drastic than expected. At this mass, subhalos are less susceptible to rapid tidal disruption by their hosts. Those subhalos that tend to experience significant depletion due to baryons are generally lower in mass and v_max, such as in the simulation data studied by Garrison-Kimmel et al. (2017). If our method is extended to ultra-faint dwarf galaxies, these issues will need to be reconsidered.

In this analysis, all host halos with a massive satellite analog and one or more low mass satellite analogs were used to infer the halo mass of the MW (see Section 4). The latter constraint, that there should be one or more low mass satellite analogs, rather than 10 per host like the classical satellites around the MW, is necessary due to the current state of cosmological simulations and specifically their resolution limits. Here we calculate the subhalo abundance functions around the host halos used for the MW halo mass calculations presented in Section 4 to demonstrate that on average there are 10–20 subhalos with masses ≥10⁹ in the vicinity of their host's center of mass. However, the general properties of these subhalos, such as their v_max and current distances from their respective hosts, are not all representative of the classical MW satellites and were therefore omitted from Prior 2. In the left panel of Figure 7, we show the cumulative subhalo abundance function for all host halos in Prior 2 that were selected from the final 5 Illustris-Dark snapshots. These halos host exactly one massive satellite analog and at least one or more low mass satellite analogs (see Sections 3.1.1 and 3.1.2 for definitions).

Zoom In Zoom Out Reset image size

In Section 5, we note that the trend between host halo mass and the distribution of satellite specific orbital angular momenta in the presence and absence of a massive satellite analog differs at host halo masses <10^11.5 ${M}_{\odot }$. Below this host halo mass, the median orbital angular momentum is higher when massive satellites are not present. However, inferred MW masses resulting from the ensemble of classical satellites are still in good agreement with each other, independent of this difference. Here, we examine the cumulative subhalo abundance for host halos that do not have a massive satellite analog since 33% of ∼10¹² ${M}_{\odot }$ halos typically host an LMC-mass companion (Patel et al. 2017a).

To create the abundance functions, all subhalos belonging to host halos in Prior 2, regardless of whether they satisfy our selection criteria in Sections 3.1.1 and 3.1.2, with a subhalo mass ≥10⁹ ${M}_{\odot }$ and that reside within their host's virial radius are considered. All included subhalos are not necessarily members of the Prior 2 sample, but this sample encompasses all subhalos above our imposed ∼133 dark matter particle resolution limit. The abundances are plotted as a function of $\mu \equiv {M}_{\mathrm{sub}}/{M}_{\mathrm{vir}}$, or the ratio between subhalo mass and host virial mass. The error bars on the data points indicate the standard deviations in each μ bin. Host halos less massive than 10^11.85 ${M}_{\odot }$ are excluded, as this is beyond the range of MW masses in the literature.

The curves in the left panel of Figure 7 all exhibit an increase near μ = 10⁻² due to the strict one massive satellite criterion (see Section 3.1.1). This is most obvious in the lowest host halo masses (blue curve) where the massive satellites are a significant fraction of their host's mass. The highest mass host halos (>10^12.5 ${M}_{\odot }$, red line) are expected to host the greatest number of subhalos, but the strict one massive satellite criterion causes this trend to drop off more quickly around μ = 10⁻³. Notice that at μ = 10⁻³, all host halos harbor 8–10 subhalos on average. While the results presented in Section 4 use halos that have a massive satellite and one or more low mass satellite analogs, the abundance of subhalos around hosts in Prior 2 are generally in agreement with the observed MW satellite population.

In the right panel of Figure 7, the cumulative abundance functions are shown for all host halos that do not host a massive satellite analog. Notice there is no longer a peak near μ = 10⁻² as the requirement for exactly one massive satellite analog is eliminated. Overall, the shape of the abundance functions changes most significantly at μ > 10⁻³ and the abundance of subhalos is higher when a massive satellite analog is absent. At μ = 10⁻³, hosts tend to have 20–40 subhalos, compared to only 8–10 subhalos when massive satellites are present. The exact cause for the higher subhalo abundance in host halos >10^11.85 ${M}_{\odot }$ is unknown and will be explored in future work. While the overall subhalo abundances differ for these host halo samples, the distribution of satellite specific orbital angular momentum for specifically the low mass satellite analogs in each sample generally remains the same and therefore yields similar MW masses.

Cumulative subhalo abundance functions for the Illustris simulations have been presented by Chua et al. (2017), where they have not imposed any selection criteria on their host halos as we have implemented here. Their samples are chosen from the Illustris-1 simulation, which includes baryons and hydrodynamical processes, so the abundance of subhalos with masses ∼10⁹ ${M}_{\odot }$ is immediately lower (see Appendix A). Using the sample chosen from Illustris-1, they find equivalent matches for each Illustris-1 halo in the Illustris-Dark simulation and formulate their cumulative abundances based on those halo populations. Our Figure 7 data are chosen directly from the Illustris-Dark simulation and do not suffer these depletion effects, leading to higher subhalo abundances overall.