THE CONNECTION BETWEEN THE HOST HALO AND THE SATELLITE GALAXIES OF THE MILKY WAY

The Milky Way (MW) galaxy and its satellite galaxies provide an excellent laboratory for constraining cosmology and galaxy formation physics. Accurately modeling the formation of galaxies in the MW system requires stringent constraints on galaxy formation physics and the properties of the dark matter halos in this particular environment. One of the unusual characteristics of the MW galaxy is that it has two massive satellite galaxies, the Large and Small Magellanic Clouds (MCs). Both of them are measured to have maximum circular velocities ${V}_{\max }\geqslant 55\,\mathrm{km}\,{{\rm{s}}}^{-1}$ with magnitudes of ${M}_{{\rm{v}}}=-18.5$ and −17.1, respectively (e.g., van den Bergh 2000; van der Marel et al. 2002; Stanimirović et al. 2004; van der Marel & Kallivayalil 2014). The next brightest satellite is Sagittarius dSph, about 4 mag dimmer, with ${V}_{\max }\sim 20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Strigari et al. 2007). In addition, the MCs have significantly higher stellar masses than the rest of the MW satellite galaxy population. The stellar mass of the LMC is about $1.5\times {10}^{9}\,{M}_{\odot }$, and the stellar mass of the SMC is about $4.6\times {10}^{8}\,{M}_{\odot }$, which are 1.8 and 1.3 dex higher than the stellar mass of the core of Sagittarius, which is $2.1\times {10}^{7}\,{M}_{\odot }$ (McConnachie 2012). Several works have shown that satellites as bright as the LMC are rare around MW-mass galaxies (e.g., Guo et al. 2011; Liu et al. 2011; Tollerud et al. 2011). This rarity suggests that the presence of such high-mass satellites may have important implications for the formation history of the MW halo. In this paper, we will study the connections between the properties of the MW host halo and the existence of MCs, and explore the impact of the host halo properties on the inference of galaxy formation using observations of the satellite galaxies of the MW, including the MCs.

There have been a number of theoretical studies in the literature investigating the implications of high-mass substructures in the MW analogues. Boylan-Kolchin et al. (2011a) used a large sample of MW-mass halos ($4.3\times {10}^{11}\leqslant {M}_{\mathrm{vir}}/{M}_{\odot }\leqslant 4.3\times {10}^{12}$) extracted from the Millennium-II simulation and found that 20% of MW-mass halos host an LMC or SMC, and that only ∼2.5% of such halos host both MCs as a binary. This indicates that the MW system must be rare in a Cold Dark Matter (CDM) universe. A similar conclusion was reached by Busha et al. (2011b) using a different set of simulations.

On the other hand, selecting MW analogues from cosmological simulations by requiring the existence of MCs can be used to constrain the MW's virial mass, the dark matter halo density profile, and its satellite accretion history (e.g., Busha et al. 2011a; Cautun et al. 2014a). In recent studies, Wang et al. (2015) found that halo concentration strongly affects the estimate of halo mass using dynamical tracers of the MW. Mao et al. (2015) demonstrated that halo concentration influences the inferences for the mass and other properties of the MW halo from satellite occupation statistics. It has also been shown that the abundance of subhalos correlates with the mass-assembly history of the halo (e.g., Zentner et al. 2005; Zhu et al. 2006; Mao et al. 2015).

The baryonic processes of galaxy formation also significantly impact the properties of the satellite galaxy population (Bullock et al. 2001; Benson et al. 2002a, 2002b; Font et al. 2011; Gómez et al. 2014; Sawala et al. 2014; Wetzel et al. 2016). A number of authors have adopted various approaches to model the formation of MW satellite galaxies. It has been shown that satellites as bright as the LMC are rarely predicted for MW-mass galaxies, with only ∼10% of the MW-mass galaxies having satellites as bright as the LMC (e.g., Kauffmann & White 1993; Koposov et al. 2009; Liu et al. 2011). Moreover, Koposov et al. (2009) found that an unusually high star formation efficiency was needed in their model to reproduce objects as bright as LMC and SMC, implying again that the MW system is unusual owing to the existence of MCs. Similarly, Okamoto et al. (2010) explored a range of feedback models to add galaxies to some of the high-resolution Aquarius halos (Springel et al. 2008), and again found it difficult to readily reproduce galaxies with luminosities as high as the MC's while simultaneously matching the luminosities of lower mass satellite galaxies. The same conclusion was reached by Starkenburg et al. (2013) using a different model.

While this difficulty could stem from incomplete physics being considered in current galaxy formation models, it is possible that the influence of the underlying dark matter structure is the cause of the issue. As the backbone of galaxy formation, the dark matter halo and its associated subhalos have a strong impact on the properties of the hosted galaxies. For example, Starkenburg et al. (2013 and references therein) showed a strong correlation between the number of satellites and the dark matter mass of the host halo. It has also been shown that the formation of the satellite galaxies of the MW is highly stochastic and sensitively depends on the subhalo population (e.g., Cautun et al. 2014a; Guo et al. 2015). Also, Gómez et al. (2012) found that halo merger histories and galaxy formation physics are degenerate under certain observational data constraints. Thus, to accurately model the formation of a particular galaxy, like the MW, it is important to understand how the prior information about the properties of the dark matter halo, including its concentration, mass-assembly history, and hosted subhalos, affects the model predictions of the MW and its satellite populations. Such an understanding is also important for interpreting observational data on the Andromeda galaxy and MW analogues in the more distant universe (e.g., Guo et al. 2011; Liu et al. 2011; Tollerud et al. 2011; Nierenberg et al. 2013, 2016).

This paper is dedicated to investigating the host halo prior. First, we attempt to gain insight into what halo properties are of interest when modeling the formation of MW satellite galaxies. To achieve this, we study a set of N-body simulations of MW-mass halos, and apply a semi-analytic model to the merger trees extracted from those simulations to study galaxy properties. To understand how the halo prior influences inferences of galaxy formation physics, we employ a Markov-Chain Monte-Carlo (MCMC) machinery that is joint with our SAM to explore the parameter space of a galaxy formation model. Second, using the constrained models, we investigate the effect of the halo prior and gain insight into which aspects of modeling and observational work are needed to further tighten the constraints on the formation of the MW and its satellite galaxies. Specifically, we explore the parameter space in the galaxy formation model to understand the constraining power of observational estimates of the stellar mass–metallicity relation, the kinematics of satellite galaxies, and the number counts of small satellite galaxies brought in by the MCs. In addition, extrapolating the constrained model to central galaxies of MW-mass halos, we demonstrate how the connections between the host halo and the satellite galaxies can also influence the mass of the central galaxies, shedding light on understanding the effects of the halo assembly bias (e.g., Gao et al. 2005; Wechsler et al. 2006; Jing et al. 2007) for galaxy formation.

In this paper, we describe the simulations and the SAM adopted in Section 2. The results on the relation between the host halo properties and the high-mass subhalos from analyzing the dark matter simulations are presented in Section 3. We present the model inferences based on the SAM from the MW satellite stellar mass function in Section 4. We summarize the conclusions of the study and discuss the implications in Section 5.

2.1. The Simulations

2.2. The Semi-analytic Model

In this section, we study the correlation between the existence of a high-mass satellite galaxy like the LMC, and host halo properties. With the simulated dark matter halos, we first study their mass accretion histories (MAHs). The MAH is defined as the virial mass of the main-branch halo as a function of redshift or time, with the virial mass definition following Bryan & Norman (1998). To capture the overall shape of the MAHs, we adopt the two-parameter model proposed by Tasitsiomi et al. (2004). This model combines the exponential form of Wechsler et al. (2002) and the power-law form of van den Bosch (2002),

$$\begin{eqnarray}&&M(z)={M}_{0}{\left(\displaystyle \frac{1+z}{1+{z}_{0}}\right)}^{-\beta }\exp [-\gamma (z-{z}_{0})],\end{eqnarray} \tag{ 1 }$$

where β and γ are model parameters determining the shape of the MAH, and M₀ is the normalization for the halo mass at $z={z}_{0}$. Using this model, McBride et al. (2009) fit the MAHs of halos in the Millennium simulation in a large halo mass range, and found that this two-parameter model captures the halo MAHs remarkably well.⁸ Following the previous work, we also fit the MAHs of the MW halos in our simulations using this model. A merit of this fitting function is that the model parameters characterize the accretion rate of a halo at early and late epochs. In this model, we find that

$$\begin{eqnarray}&&\displaystyle \frac{d\,\mathrm{log}\,M}{{dz}}=-\left(\displaystyle \frac{\beta }{1+z}+\gamma \right),\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{d\,\mathrm{log}\,M}{d\,\mathrm{log}\,a}=\beta +\gamma (1+z).\end{eqnarray} \tag{ 3 }$$

The derivations show that at the present time, when z = 0 and a = 1, the accretion rate $d\,\mathrm{log}\,M/d\,\mathrm{log}\,a$ is characterized by $\beta +\gamma$, and at early times, the accretion rate is characterized by γ. Motivated by these indications, we will use $\beta +\gamma$ and γ extracted from the best-fit model to characterize the accretion rate at late and early times for each halo, and to study the correlations between the MAH and the satellite population.

In Panel (A) of Figure 1, we show the distribution of all the simulated MW host halos in the parameter space defined by the MAH fitting parameters, $\beta +\gamma$ and γ. The halos from the cosmological box (c125-2048) are shown as circles, and the halos from the zoom-in simulation are shown as stars with a black circle at the center. All the halos populate a particular region of the parameter space with a mode located around $\beta +\gamma =0.3$ and $\gamma =0.7$. This distribution is in agreement with what is found in other simulations (McBride et al. 2009; Taylor 2011). In the diagram, halos in the upper left branch have rapid late-time accretion, while those in the lower-right branch have rapid early-time accretion.

Zoom In Zoom Out Reset image size

The size of the symbols in Figure 1 is linearly proportional to the scale radius, ${r}_{{\rm{s}}}$, of the NFW profile, and inversely proportional to the concentration parameter $c={R}_{\mathrm{vir}}/{r}_{{\rm{s}}}$, where ${R}_{\mathrm{vir}}$ is the virial radius of the halo. From the figure, it can be seen that the concentrations are correlated with the MAH fitting parameters. In general, halos with faster late-time accretion tend to have lower concentration (or larger ${r}_{{\rm{s}}}$ for given ${R}_{\mathrm{vir}}$). This is more clearly seen in panel (B) of Figure 1, where we plot the distribution of the simulated host halos in a diagram defined by halo concentration c and $\beta +\gamma$. In the diagram, one can find that rapid accreting halos also tend to have a lower concentration and hence tend to be located in the upper left part of the diagram.

In the figure, we also color code each halo according to the mass of its most massive subhalo, ${M}_{\mathrm{sub}1}$, which is defined by the total mass confined in the phase space of the subhalo with the Rockstar halo finder (Behroozi et al. 2013b). The scale of $\mathrm{log}{M}_{\mathrm{sub}1}$ and the corresponding color are shown in the color bar on the right side of Figure 1. An interesting trend shown in the figure is that the mass of the most massive subhalo is also correlated with the MAH parameters and the concentration of the host. Hosts that have more rapid late-time accretion and lower concentration are more likely to host a high mass subhalo. We note that when we choose to use other subhalo properties as a proxy for high-mass satellites, for instance the peak subhalo mass (${M}_{\mathrm{peak}}$), the subhalo mass at accretion (${M}_{\mathrm{acc}}$), or the maximum subhalo circular velocity (${V}_{\max }$), the trends we describe here still hold, but become weaker (see also Mao et al. 2015). This is because these quantities more closely reflect the accretion of a high-mass subhalo at a much earlier time. The massive stellar streams or debris of tidally disrupted satellite galaxies may also contain information about the properties of the MW host halo, but the instantaneous mass of high-mass subhalos is more closely related to the present-day concentration and the recent accretion history of the host.

Adjacent to panels (A) and (B), we show the marginalized distribution of γ, $\beta +\gamma$, and c in panels (A1), (A2), and (B1), respectively. The distribution of $\beta +\gamma$ is repeated in panel (B2). The solid color bars represent the distribution of the cosmological simulation halos. The colors of each bar indicates the median mass of the most massive subhalo for all the hosts in the bin. Again, we can see a clear sequence that halos with higher $\beta +\gamma$ and a lower concentration tend to host a higher mass subhalo. The gray hashed bars represent the zoom-in simulation halos. The small sample size of the zoom-in simulation covers the distribution of the population reasonably well with some discrepancy in the halo concentration distribution.

To quantify how the concentration parameter and the MAH fitting parameter $\beta +\gamma$ affect the probability for an MW-mass halo to host MCs, we count the number of subhalos that can possibly host the MCs. Observational estimates suggest that the maximum circular velocity, ${V}_{\max },$ is $\sim 50\mbox{--}60\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for the SMC (Stanimirović et al. 2004; Harris & Zaritsky 2006) and $\gt 80\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for the LMC (Olsen et al. 2011; Kallivayalil et al. 2013; van der Marel & Kallivayalil 2014). We choose halos that host at least two subhalos with ${V}_{\max }\geqslant 55\,\mathrm{km}\,{{\rm{s}}}^{-1}$ from the cosmological simulation. Overall, 28% of the MW-mass halos in the simulation satisfy this criterion, broadly consistent with predictions of other simulations (Boylan-Kolchin et al. 2010; Busha et al. 2011b; Rodríguez-Puebla et al. 2013; Cautun et al. 2014b). In the upper panels of Figure 2, we show the cumulative distribution functions of the concentration parameter c and the MAH fitting parameter $\beta +\gamma$ of those halos selected from the cosmological simulation. In comparison, dashed lines show the same cumulative distribution functions for all the MW-mass hosts. The comparison clearly shows that, at a fixed mass, halos that host high-mass subhalos tend to have a lower concentration and more rapid late-time accretion. When selecting an MW host with both host mass and subhalo ${V}_{\max }$, the median concentration decreases to c = 8.2 as opposed to the median concentration c = 10.4 for the entire simulated sample selected at a fixed mass. A Kolmogorov-Smirnov test on the distribution functions safely rejects the hypothesis that the distributions are the same with a p-value less than $6.9\times {10}^{-18}$, confirming that the difference in the cumulative distribution function is statistically significant.

Zoom In Zoom Out Reset image size

In the lower panels of Figure 2, we show the fraction of the halos that host at least two subhalos with ${V}_{\max }\geqslant 55\,\mathrm{km}\,{{\rm{s}}}^{-1}$ as a function of host concentration (lower left), and as a function of $\beta +\gamma$ (lower right). We split the entire halo sample into 10 bins according to their concentration or $\beta +\gamma$. The circles denote the mean fraction of such hosts for a given bin. The error bar for each bin represents the standard deviation estimated assuming the Poisson distribution. As one can see, the fraction of halos that host at least two ${V}_{\max }\geqslant 55\,\mathrm{km}\,{{\rm{s}}}^{-1}$ subhalos decreases with concentration and increases with $\beta +\gamma$. In this sample of a fixed halo's mass, for host halos with concentrations lower than five, more than 67% of them host two or more MC analogs. This probability drops below 20% when the concentration is higher than 10. For the MAH fitting parameter $\beta +\gamma$, when $\beta +\gamma$ is higher than one, more than half of those halos can host at least two MC-mass subhalos. The result shown in the lower panels of Figure 2 demonstrates that the probability for a fixed mass halo to host MCs strongly depends on the concentration parameter and its late-time accretion. An MW-mass halo can have a high probability to host MCs if its concentration is low and it has rapid late accretion.

To explore how varying accretion history and concentration affects the satellite galaxy populations, we split the halo sample into two groups according to the location of a halo in the c–$(\beta +\gamma )$ diagram of Figure 1 (panel (B)). We use a division line, which is shown by a dashed line in the right panel of Figure 1, to split the halo samples into two groups. The upper left part of the diagram is named "Group 1" (GP1), and the lower-right part is named "Group 2" (GP2). Halos in Group 1 are less concentrated, have more rapid recent accretion, and tend to host high-mass subhalos, while halos in Group 2 are more concentrated, have slower recent accretion, and do not tend to host high-mass subhalos. The division defined here is not intended to be a physical classification, but merely to roughly split the host halo samples into two subsamples with equal size so that we can contrast the differences between them in the following studies. With a larger sample, one could split the host population into finer bins to better elucidate the trends we explore in this paper. For the current sample from zoom-in simulations, to which we will apply our SAM, we have 20 halos in Group 1, and 18 halos in Group 2.

For each group, we show the general trend of the MAH and the concentration of the host halos as a function of redshift in panels (A) and (B) of Figure 3, respectively. We show the median and the 20th, 50th, and 80th percentiles of the distributions of these functions predicted by the cosmological simulation. In the comparisons, the Group 1 hosts show relatively more rapid mass accretion at late times ($z\lt 1$), while the Group 2 hosts have flat MAHs at late times. The concentration also shows a clear difference between the two groups. The median values of the concentration for both groups stay around five from high redshift until $z\sim 1$ for both halo groups. Group 1 halos slowly increase their concentration, with the $80\mathrm{th}$ percentile being around 5–10 at z = 0. In contrast, Group 2 halos increase their concentration rapidly at $z\lt 1$. At z = 0, the 80th percentile of the concentration distribution for Group 2 hosts is between 10 and 15.

Zoom In Zoom Out Reset image size

In our simulated halo samples, the median concentration is 8.2 for Group 1 hosts and 12.5 for Group 2 hosts. The typical values of the concentration parameter of the simulated halos, especially those in Group 1, appear to be much lower than some values derived from fitting an NFW profile to observational kinematic data of the MW (e.g., Battaglia et al. 2005; Deason et al. 2012). One important difference to note is that the concentration values we report here are for halos in dark-matter-only simulations, while the MW in reality has baryons. To make a fair comparison with existing kinematic constraints, we predict the circular velocity as a function of radius for the simulated halos by employing a toy model, which accounts for the gravitational potential of the disk and the bulge of the MW and the effect of halo contraction as a response to the formation of the MW galaxy. In this toy model, we assume that the mass of the MW disk is ${M}_{\mathrm{disk}}=5.2\times {10}^{10}\,{M}_{\odot }$, and the mass of the MW bulge is ${M}_{\mathrm{bulge}}=9\times {10}^{9}\,{M}_{\odot }$ (Licquia & Newman 2015). We also assume that the disk follows an exponential profile with ${r}_{d}=3\,\mathrm{kpc}$ and that the bulge follows a Hernquist profile (Hernquist 1990) with a scale radius ${r}_{{\rm{b}}}=0.6\,\mathrm{kpc}$ (Smith et al. 2007). We adopt the adiabatic contraction model proposed by Blumenthal et al. (1986) to numerically compute the final dark matter density profile, ignoring other detailed dynamical effects, which can affect the circular velocity at about the 5% level (Gnedin et al. 2004; Choi et al. 2006).

Panel (C) of Figure 3 shows the predictions for each halo group. The color bands and the central solid lines show the 20th, 50th, and 80th percentiles of the distribution and the median of the predictions for each halo group. To show the effect of the baryonic components, we also plot in dotted lines the rotation curves that resulted from only the dark matter halo with an NFW profile with the median value of concentration for each group. In addition, the dashed line shows the rotation curves for the case where the MW disk and bulge are added, but the halo contraction effect is ignored. By comparing these predictions, we find that the gravitational effect of the baryonic components of the MW substantially boosts the circular velocity in the inner halo $r\lt 60\,\mathrm{kpc}$. We overplot observational constraints from McMillan & Binney (2010), McMillan (2011), Wilkinson & Evans (1999), Xue et al. (2008), Gnedin et al. (2010), and Watkins et al. (2010) on the same panel to compare with the predictions. This comparison shows that, when the baryonic component and the contraction effect are taken into account, both halo groups, even with $c\sim 8$ as predicted by dark-matter-only simulations, are consistent with most of the current kinematic data constraints, though only Group 2 is consistent with the McMillan (2011) constraint.

We also show the distribution of the cumulative subhalo mass function of each host group in panel (D) of Figure 3. The mass function for Group 1 hosts is systematically shifted to higher masses relative to the mass function of Group 2 hosts. This is because the first two most massive subhalos in Group 1 hosts typically have higher mass than those of Group 2 hosts. The median mass of the most massive subhalo in Group 1 hosts is about ${10}^{10.75}\,{M}_{\odot }$, and ${10}^{10.25}\,{M}_{\odot }$ for Group 2 hosts. For this reason, the cumulative mass function of Group 2 hosts appears to be steeper than that of the Group 1 hosts in the regime of $N(\gt {M}_{\mathrm{sub}})\lt 3$. We shift down the median of the cumulative mass function of Group 1 hosts by two and plot it as a red dashed line in the same figure, demonstrating that the two mass functions (with one of them being shifted) agree very well for $N\geqslant 3$. This comparison indicates that the main difference in the subhalo population between hosts with different MAH and concentration is the masses of the two most massive subhalos. When those rare high-mass subhalos are excluded, the rest of the subhalo population is basically the same from halo to halo in terms of the mass distribution. We verify this statement for other properties of the satellite population when we apply our SAM on the host halos in the following sections.

4.1. Satellite Galaxy Stellar Mass Function

We have shown that, at a fixed halo mass, the formation history and concentration vary and span a wide range. Correlated with those variations, the subhalo population also varies. We have split the simulated MW halos into two subsamples with roughly equal sizes based on their concentration and formation history to represent two different halo priors for a given final mass. In this section, we test the effects of halo priors defined by the two halo groups on the inference of galaxy formation model parameters by performing a Bayesian inference from the MW satellite stellar mass function.

The data we use in this paper is the stellar mass function of the MW satellite galaxies, down to the $11\mathrm{th}$ most massive satellite galaxy (${M}_{* }=2.9\times {10}^{5}\,{M}_{\odot }$). We adopt the stellar masses and memberships of MW satellite galaxies compiled in McConnachie (2012). Tollerud et al. (2008) have shown that the incompleteness of the MW satellite galaxy count becomes important only for fainter dwarfs ${M}_{{\rm{v}}}\gt -7$ or $L\lt {10}^{5}{L}_{\odot }$ unless there is a significant low-surface-brightness population of satellites like Crater 2 (Torrealba et al. 2016). We restrict this study to higher masses to avoid uncertainties in incompleteness corrections.

The theoretical prediction for the mass function is straightforward, but a likelihood function of the satellite mass function for the data given model is needed to perform a Bayesian inference. Supported by the tests presented in the Appendix, we adopt the Negative Binomial Distribution (NBD) as a model for the likelihood function of the satellite stellar mass function,

$$\begin{eqnarray}&&L{(D| \theta )={{\rm{\Pi }}}_{i}P({N}_{i}| r,{p}_{i})={{\rm{\Pi }}}_{i}\displaystyle \frac{{\rm{\Gamma }}({N}_{i}+r)}{{\rm{\Gamma }}(r){\rm{\Gamma }}({N}_{i}+1)}{p}_{i}^{r}(1-{p}_{i})}^{{N}_{i}},\end{eqnarray} \tag{ 4 }$$

where N_i is the observed number of satellite galaxies for a given stellar mass bin i per MW halo; the two parameters r and p_i are determined by the model as $r=1/{s}_{{\rm{I}}}^{2}$, ${p}_{i}=1/(1+{s}_{{\rm{I}}}^{2}{\mu }_{i});$ ${\mu }_{i}$ is the expected number for the $i\mathrm{th}$ mass bin predicted by the model. ${s}_{{\rm{I}}}$ is the fractional scatter from the intrinsic scatter, ${\sigma }_{{\rm{I}},i}$, with respect to the Poisson scatter, ${\mu }_{i}$, defined as ${s}_{{\rm{I}}}\equiv {\sigma }_{{\rm{I}},i}/{\mu }_{i}$ (see the Appendix). We note that the value of ${s}_{{\rm{I}}}$ may vary as a function of mass bin and can be simulated for any given model if a large number of merger trees are utilized. In this paper, however, we assume that it is a constant for the mass bins and treat it as a nuisance parameter to be sampled with MCMC because the number of the simulated hosts is not large enough to yield an accurate estimate for this parameter. We marginalize the posterior distribution of this parameter in the analysis.

We use MCMC to sample the posterior probability density distribution. Two separate runs are performed with the Group 1 hosts and the Group 2 hosts, respectively. In both runs, we allow six model parameters and the nuisance parameter ${s}_{{\rm{I}}}$ to vary within their prior ranges. A brief description of the model parameters and priors are listed in Table 1. For a detailed explanation of these parameters, readers are referred to Lu et al. (2014). For each halo prior, we run the MCMC for 20,000 iterations with 144 parallel chains using the differential evolution algorithm (Ter Braak 2006). The convergence test is done with the Gelman–Rubin test (Gelman & Rubin 1992), requiring the potential scale reduction factor $\hat{R}\lt 1.2$. After removing outliers and pre-burn-in states, we obtain ∼2,000,000 posterior samples from the MCMC for each run.

Table 1.  Summary of Semi-analytic Model Parameters

Notation	Meaning of Parameter	Prior	Best Fit for GP1	Best Fit for GP2
${{\rm{\Sigma }}}_{\mathrm{SF}}$	the threshold gas surface density for star formation in units of ${M}_{\odot }\,{\mathrm{pc}}^{-2}$	$\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{SF}}\in [-1.0,2.0]$	no constraint	no constraint
${\alpha }_{\mathrm{LD}}$	the normalization of the mass-loading factor for feedback outflow	$\mathrm{log}{\alpha }_{\mathrm{LD}}\in [-9.0,1.0]$	−0.98	−5.85
${\beta }_{\mathrm{LD}}$	the power index for the circular velocity dependence of the mass-loading factor	${\beta }_{\mathrm{LD}}\in [0.0,15.0]$	5.32	11.7
${V}_{\mathrm{out}}$	characteristic halo circular velocity in $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, below which all outflow mass leaves the host halo	$\mathrm{log}{V}_{\mathrm{out}}\in [1.5,3.0]$	no constraint	no constraint
${\beta }_{\mathrm{out}}$	the steepness of the transition from total outflow for halos with ${V}_{c}\lt {V}_{\mathrm{out}}$ to no outflow for halos with ${V}_{c}\gt {V}_{\mathrm{out}}$	${\beta }_{\mathrm{out}}\in [0.0,8.0]$	no constraint	no constraint
${\gamma }_{\mathrm{RI}}$	the fraction of outflow mass reincorporated back into the halo	$\mathrm{log}{\gamma }_{\mathrm{RI}}\in [-2.0,0.0]$	no constraint	no constraint
${s}_{{\rm{I}}}$	the fractional intrinsic scatter of the satellite mass function	${s}_{{\rm{I}}}\in [0.05,0.5]$	no constraint	no constraint

Download table as:  ASCIITypeset image

We show the 2D marginalized posterior distribution for the Group 1 halo prior and the Group 2 halo prior in Figure 4. We also list the "best-fit" model (maximum likelihood) parameters for each group in Table 1 for comparison. As one can see, under the same data constraints, using a different halo prior results in different posteriors of the model parameters. An obvious change is the normalization (${\alpha }_{\mathrm{LD}}$) and the power index (${\beta }_{\mathrm{LD}}$) for the parameterization of the mass-loading factor of outflow. In the model, the outflow mass-loading factor is parameterized as

$$\begin{eqnarray}&&\eta ={\alpha }_{\mathrm{LD}}{\left(\displaystyle \frac{{V}_{{\rm{c}}}}{220\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-{\beta }_{\mathrm{LD}}},\end{eqnarray} \tag{ 5 }$$

where ${V}_{{\rm{c}}}$ is the circular velocity of a halo, ${\alpha }_{\mathrm{LD}}$ and ${\beta }_{\mathrm{LD}}$ are model parameters to be constrained. The normalization parameter ${\alpha }_{\mathrm{LD}}$ for Group 1 is higher than that for Group 2. Correspondingly, the power index for Group 1 is lower than that for Group 2. The normalization for the mass-loading factor is defined as the mass-loading factor for halos with a circular velocity ${V}_{{\rm{c}}}=220\,\mathrm{km}\,{{\rm{s}}}^{-1}$. A higher ${\alpha }_{\mathrm{LD}}$ requires a higher fraction of supernova energy to power outflows to keep the baryon fraction of a halo low. The Group 1 halos are those having high probability to host massive subhalos. The high-mass subhalos typically have higher baryon mass to start with before being accreted into the host. To keep the baryon mass fraction low in those high-mass subhalos, the model is required to have stronger feedback to suppress star formation. The Group 2 halos, which do not tend to have massive subhalos, typically require lower feedback to allow a relatively higher stellar mass fraction to fit the satellite mass function.

Zoom In Zoom Out Reset image size

As shown in Figure 4, it is clear that the two parameters for the outflow mass-loading factor are strongly degenerate, regardless of which halo prior is adopted. The degeneracy can be approximately described by a linear function as ${\beta }_{\mathrm{LD}}=A\mathrm{log}{\alpha }_{\mathrm{LD}}+B$, where A and B are the slope and the intercept of the linear function. The straight orange line in the insert panel of Figure 4 shows this linear function with $A=-1.2$ and B = 4. Recall that the outflow mass-loading factor is parameterized as Equation (5). Hence, if ${\alpha }_{\mathrm{LD}}$ and ${\beta }_{\mathrm{LD}}$ follow the aforementioned linear relationship, then we have $\mathrm{log}{V}_{220}=1/A$ and $\mathrm{log}\eta =-B/A$. Using the values of A and B, we derived from the degeneracy, we find that the outflow mass-loading factor $\eta \approx 2100$, when ${V}_{{\rm{c}}}=32\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This is the generic feature for all models that are on the ridge of the degeneracy. This velocity scale is interesting because it is between the circular velocity of SMC and less massive satellites, such as Canis Major and Sagittarius dSph. The degenerate models collectively require extremely high outflow mass-loading factor for subhalos with circular velocity $\leqslant 32\,\mathrm{km}\,{{\rm{s}}}^{-1}$. It would be interesting to investigate the significance of this velocity scale using other data sets of galaxy populations.

We marginalize over the posterior to show how the models, together with their predetermined halo priors, reproduce the MW satellite stellar mass function. The reproduced stellar mass function using the two underlying dark matter halo priors are shown in the left panel of Figure 5. The red bands denote the posterior predictive distribution using the Group 1 halo prior. The solid line in the middle of the bands denotes the median of the predictive distribution, and the bands from darker to lighter color enclose 20%, 50%, and 80% of the predictive probability distribution. Similarly, the blue bands show the predictive distribution using the Group 2 halo prior with the same three levels of confidence range. In addition, we also show the observed satellite mass function of the MW in the same figure. We find that when the model is applied to Group 1 hosts, it can achieve a very good match to the observed stellar mass function for the most massive 11 satellite galaxies. When the model is applied to Group 2 hosts, while the model matches the mass function equally well as using Group 1 hosts for satellites with masses lower than the SMC, it still tends to predict lower stellar masses for the two most massive satellites. It is because the hosts do not host high-mass subhalos, precluding solutions that can make sufficiently high-mass galaxies to match the stellar mass of MCs. We compute the Bayes Factor, which is defined as the ratio of the marginalized likelihoods of the models based on the Group 1 and Group 2 halo priors. We find that the Bayes Factor ${ \mathcal B }={{ \mathcal M }}_{1}/{{ \mathcal M }}_{2}=1.42$, which indicates that the halo prior of a Group 1 type MW host is only weakly preferred by the data over the halo prior of Group 2 type host. Due to large uncertainties, the model is not able to rule out either of the two halo priors using the observed MW satellite galaxy stellar mass function.

Zoom In Zoom Out Reset image size

Using the constrained models, we show the biases in the model predicted MW satellite mass function if a mismatched halo prior is adopted. To illustrate the extreme case, we apply the model constrained using Group 1 halos to the merger trees of Group 2 halos, and the model constrained using Group 2 halos to the Group 1 halo merger trees to predict the satellite galaxy stellar mass function. The posterior predictive distributions are shown in the right panel of Figure 5. The model constrained for Group 1 halos predicts a significantly lower stellar mass function when Group 2 halo merger trees are employed. On the other hand, the model constrained for Group 2 halos predicts a too high mass function when Group 1 halo merger trees are applied. The observed stellar mass function is excluded from the 50% confidence range of the predictive distributions of both models. The MCs are even outside the 80% confidence range. This result again demonstrates the strong influence of the underlying dark matter halo formation history on the hosted galaxies.

4.2. Satellite-subhalo Matching

Using the posteriors, we can study the connections between MW satellite galaxies and their host subhalos based on the two halo priors and the corresponding models specifically tuned for each halo prior. Instead of choosing one single model (the best or any arbitrary one), we use the posterior samples obtained with MCMC to produce the posterior predictive distribution (Lu et al. 2012). The predictive distribution marginalizes over the model uncertainties that are allowed by the uncertainties of the constraining data and quantifies the confidence levels of the predictions for the given model and prior. Figure 6 shows the maximum circular velocity, ${V}_{\max },$ of subhalos as a function of their stellar mass, as predicted by the models using the two halo priors. The color bands with decreasing intensity denote the 20%, 50%, and 80% confidence ranges of the predictive distribution for each model, and the solid lines show the median predictions. The error bars on each predictions represent the $1\sigma$ galaxy-to-galaxy scatter averaged over all the posterior sample models. At high mass, $\mathrm{log}{M}_{* }/{M}_{\odot }\gt 6.5$, the model based on the Group 1 halo prior predicts higher ${V}_{\max }$ than the model based on Group 2 halos for given stellar masses. This is because Group 1 hosts typically have higher mass subhalos than Group 2 hosts in this regime. To fit the same satellite stellar masses, the model for Group 1 hosts has to populate satellite galaxies into the few relatively higher mass subhalos. At lower mass, $\mathrm{log}{M}_{* }/{M}_{\odot }\lt 6.5$, the ${V}_{\max }\mbox{--}{M}_{* }$ relation predicted by the two models becomes similar. This stellar mass scale corresponds to ${V}_{\max }\sim 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and the halo mass is about $3\times {10}^{9}\,{M}_{\odot }$. The ${V}_{\max }\mbox{--}{M}_{* }$ relation converges at this mass because this is where the subhalo mass functions of the two host groups become similar, as Figure 3 showed. This is also reflected in the degenerate behavior of the feedback parameters as we discussed in Section 4. The models based on different halo priors require different strengths of feedback for higher mass halos, but similar outflow mass-loading factors for halos with circular velocities lower than $\sim 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$. It is also worth pointing out that while the constrained models vary substantially across the parameter space, the monotonic trend between the predicted stellar mass and subhalo ${V}_{\max }$ is preserved. This means that while the detailed subhalo-galaxy matching may vary depending on the host halo prior and specific model adopted, the rank order of the satellites based on their stellar mass still generally follows the depth of the gravitational potential of the subhalos.

Zoom In Zoom Out Reset image size

We also compare the model prediction of the ${V}_{\max }\mbox{--}{M}_{* }$ relation with that derived from observations by Peñarrubia et al. (2008), Boylan-Kolchin et al. (2012), and Sawala et al. (2016) for a number of MW classical dwarf spheroidals in Figure 6. These widely cited results have quite different ${V}_{\max }$ estimates for the same galaxy. Comparing them with our model predictions, we stress the importance of accurate estimates of ${V}_{\max }$ for model inference. Among the three different estimates, the Peñarrubia et al. (2008) results have significantly higher ${V}_{\max }$ for given stellar masses and thus agree more with the model based on the Group 2 halo prior. As argued in Boylan-Kolchin et al. (2012), Peñarrubia et al. (2008) assumed the NFW density profile and the size–mass relation of field halos for subhalos, which result in the overestimation of ${V}_{\max }$. We also find that the assigned NFW profiles in Peñarrubia et al. (2008) have significantly larger sizes than what we find in the SAM predictions. The radii where the rotation curve peaks, ${R}_{\mathrm{vmax}}$, for the classic dwarfs derived by Peñarrubia et al. (2008) are typically as large as 3–5 kpc, which are rare and are only relevant for the subhalos of the MCs in our simulations. Boylan-Kolchin et al. (2012) matched the MW dwarfs with simulated subhalos based on their abundances and found much lower ${V}_{\max }$ for given stellar mass. Using both dark matter only simulations and hydrodynamical simulations, Sawala et al. (2016) also found lower ${V}_{\max }$ values. These results are overplotted in Figure 6, with symbols noted as "DMO" for dark-matter-only simulations and "Hydro" for hydrodynamical simulations. Compared to these results, our models, regardless of which halo prior is adopted, predict much higher subhalo ${V}_{\max }$ values at a given stellar mass. All the data points from Boylan-Kolchin et al. (2012) and Sawala et al. (2016) are below the lower bound of the 80% confidence range of the predictions. The higher ${V}_{\max }$ values in our model indicate that the simulated subhalos are too dense compared to the observed classic dwarfs, as Boylan-Kolchin et al. (2011b) showed (also known as the "too-big-to-fail" problem). By varying the halo prior, we show that the variation in halo concentration and MAH for fixed halo mass results in large scatter in the subhalo ${V}_{\max }\mbox{--}{M}_{* }$ relation. However, the effect of the increased scatter is only at high mass and not at stellar masses relevant to the too-big-to-fail problem. Our results demonstrate that, at this halo mass, if no baryonic processes are included, the probability for CDM subhalos to be consistent with MW dwarf kinematics is low. At a lower mass scale, this discrepancy will have even stronger statistical significance (Jiang & van den Bosch 2015).

We stress, however, the importance of understanding the baryonic processes, because recent high-resolution baryonic simulations of MW-mass halos have shown promise in addressing the "too-big-to-fail" problem (Brooks & Zolotov 2014; Sawala et al. 2016; Wetzel et al. 2016), largely because (1) tidal effects with the MW-mass stellar disk reduce the number of surviving subhalos (e.g., Zolotov et al. 2012; Wetzel et al. 2016; Zhu et al. 2016) and (2) stellar feedback can reduce the inner densities of subhalos at these masses (e.g., Read & Gilmore 2005; Pontzen & Governato 2014; Brook & Di Cintio 2015; Chan et al. 2015).

4.3. Metallicity Relation of MW Satellites

We now make predictions for the stellar-phase metallicity as a function of stellar mass for the MW satellite galaxies. The posterior predictive distributions for both host halo priors are shown in Figure 7. In the figure, the red bands are predicted using the posterior constrained for the Group 1 halo prior, and the blue bands are predictions using the posterior constrained with the Group 2 halo prior. The bands with different intensity of each color show the 20%, 50%, and 80% predictive distribution. We find that while the model can be tuned to reproduce the satellite stellar mass function using either halo prior, they make rather different predictions for the stellar-phase metallicity-stellar mass relation. For the Group 2 prior, because the hosts typically do not have high-mass subhalos, the model needs weaker feedback and outflow to yield a relatively higher stellar mass for given subhalo mass. The consequence of the weaker outflow is to leave more metals in the galaxies, resulting in a higher metallicity for a given stellar mass. In contrast, the Group 1 prior imposes host halos that typically host high-mass subhalos. The high-mass subhalos require stronger outflows, which eject more metals from the galaxy. The result suggests that for a given galaxy formation scenario, when it is calibrated to the stellar mass function, the properties of the underlying dark matter halo, including the formation history, concentration and subhalo masses, can leave imprints in the satellite galaxy metallicity relation. Because the high-mass end of the subhalo mass function is correlated with the formation history and concentration of the host halo, the metallicity may also reflect the formation history and other properties of the MW host halo and its analogues.

Zoom In Zoom Out Reset image size

The stars with error bars are observational data compiled by McConnachie (2012), and the mean observational relation is shown as solid and dashed magenta lines (Kirby et al. 2013). We note that our mass–metallicity relation using either halo prior has a steeper slope than the observed relation. Using this relation to constrain our galaxy formation models is beyond the scope of this paper, but the difference in the relations using the different priors may provide useful comparisons with data in future work.

4.4. Satellites Accreted with LMC

In a CDM cosmology, the higher amplitude of the power spectrum at small scales gives rise to a large number of substructures in an MW-mass halo (e.g., Klypin et al. 1999; Moore et al. 1999; Kravtsov et al. 2004). Wetzel et al. (2015) studied the fraction of satellites of MW-mass halos accreted as part of a group in the ELVIS simulation (Garrison-Kimmel et al. 2014)⁹ and found that higher mass satellite galaxies have higher probability of hosting satellites when they are finally accreted into MW. Using a similar method for tracking the accretion history of subhalos, we explore the number of satellite galaxies that are accreted into the MW host halo together with the most massive subhalo. These smaller halos are subhalos of the subhalo of the MCs and are accreted into the MW host with the MCs. These sub-subhalos may remain gravitationally bound with the MC subhalo, or may have been tidally stripped from the MC subhalo potential and live in the vicinity of the MCs. Using the models that are calibrated to match the satellite mass function, we test if the satellite galaxies that are accreted with the MCs can provide information to discriminate the halo priors.

In the SAM prediction, we identify the most massive satellite galaxy in each MW-mass halo as its "LMC" and track the merger history of its subhalo using the merger tree. In the merger tree, we find the time when it is just accreted into the MW host and all subhalos hosted by the "LMC" host subhalo. We then track the history of all galaxies hosted by the subhalo since this epoch and identify the ones that are still present in the MW halo at z = 0. We repeat this calculation for every model in the posterior samples to produce a posterior predictive distribution. In Figure 8, we show the stellar mass function of the satellite galaxies that are accreted into the MW host together with the "LMC" as its satellites predicted by the constrained models at z = 0. We show the predictions using the two different host halo priors. We find that Group 1 hosts have systematically higher numbers of satellites that are accreted with the most massive satellite than Group 2 hosts. Down to ${10}^{3}\,{M}_{\odot }$, a Group 1 host typically has more than eight satellite galaxies accreted with the most massive satellite, but a Group 2 host has no more than five satellite galaxies with stellar mass higher than ${10}^{3}\,{M}_{\odot }$ accreted with its "LMC." The reason is simply that the subhalos of the "LMC" in Group 1 hosts are typically more massive than those in the Group 2 hosts, and thus they bring in more low-mass satellites. While many of these small satellites have been stripped out of the "LMC" subhalo, they are still located near it in the simulation. Transforming their positions into a coordinate system centered on the host halo, we find that on average of 57% of the small satellites are located within 40° of the most massive satellite, and about 80% of them are located within 100° (see also Deason et al. 2015; Yozin & Bekki 2015; Jethwa et al. 2016; Sales et al. 2016, for detailed analyses). In the year of 2015, 21 new dwarf galaxy candidates were discovered around the MW, and a significant fraction of these are within 50° of either the LMC or the SMC (e.g., Drlica-Wagner et al. 2015). Although it is not yet clear how many of these candidates are true satellites of the MCs, our result suggests that identifying low-mass satellite galaxies that are associated with MCs may provide useful constraints on the subhalo mass of the MC and the MW host prior.

Zoom In Zoom Out Reset image size

4.5. Central Galaxy Stellar Mass

The halo prior discussed in this paper is characterized by the concentration, MAH, and the high-mass subhalo content of the host halo. These host halo properties are expected to affect the mass of the central galaxy (e.g., Lehmann et al. 2015; Zentner et al. 2016). In this subsection, we extrapolate the models constrained to the MW satellite stellar mass function for different halo priors to study the effect of the halo prior on the stellar mass of the MW galaxy. As we discussed in Section 2.2, the modeling scheme adopted for previous sections does not follow the formation of the central galaxy to focus computation on the satellite galaxy populations. In this section, we employ the full model (Lu et al. 2011), which has the additional ability to make predictions for the central galaxy. We run the full model with the same posterior parameter samples as in previous sections to predict the probability distribution of the central galaxy stellar mass that is implied by the constrained models. In Figure 9, the red solid line shows the prediction of the constrained model using the Group 1 halo prior, and the blue solid line shows the prediction of the constrained model using the Group 2 halo prior. When both models are best constrained to the observed MW satellite stellar mass function, the two models predict different stellar masses for the central galaxy. As we show in the figure, the Group 2 halo prior systematically predicts higher central galaxy stellar masses than the Group 1 halo prior.

Zoom In Zoom Out Reset image size

The difference in the predicted central galaxy stellar mass stems from two effects. First, the halos of these two host halo groups have different formation histories and subhalo content. Second, the models that are best constrained to match the MW satellite galaxy stellar mass function for each halo prior are different.

By construction, Group 1 halos typically form later and host more high-mass subhalos than the Group 2 counterparts. To show the effect of different halos on the central galaxy mass, we apply a fixed model to different merger trees to predict the central galaxy stellar masses. In Figure 9, the dashed red line shows the central galaxy stellar mass distribution predicted by the model constrained using Group 1 halo prior but applied to Group 2 merger trees. Similarly, the blue dashed line shows the same prediction made by the model constrained using Group 2 halos but applied to Group 1 merger trees. We find that, for a given model, Group 1 halos systematically form lower mass central galaxies than Group 2 halos. The effect could stem from the fact that Group 1 halos host more high-mass subhalos, which can lock up a considerable fraction of baryons that would have been available to fuel star formation in the central galaxy. However, we expect that this effect is relatively small because even the most massive subhalo in Group 1 halos typically makes up a rather small fraction (typically $\lt 15 \%$) of the total mass. To test if the difference in the central galaxy mass between the two halo groups still exists when the effect of locking baryons in subhalos is excluded, we apply the models to the main-branch halo MAHs of the merger trees, in which galaxy formation in subhalos is completely ignored as the baryons associated with the halo in this treatment can cool only into the central galaxy. The red dotted line in Figure 9 shows the predicted central galaxy stellar mass distribution of Group 1 halo main-branch MAHs, and the blue dotted line shows the prediction of Group 2 halo MAHs. One can find that while all central galaxies are predicted to have higher stellar masses than when we use the full merger trees (in this case, baryonic mass in halos is unable to cool and to form galaxies in subhalos), Group 2 halos still form systematically more massive central galaxies than Group 1 halos, even after removing the effect of locking up baryons in subhalos. In summary, Group 2 halos, which have higher concentrations and slower recent accretion, form higher stellar mass for the central galaxy than Group 1 halos mainly because Group 2 halos typically form earlier and are denser in the center. Star formation, at least in the model (e.g., Lu et al. 2014), is more efficient at higher redshift and denser halos. Consequently, the central galaxies of Group 2 halos form more stars over the history than those of Group 1 halos.

We also test the effect of different models. As we demonstrated in Section 4.1, different levels of feedback are required by the two different halo priors to fit the MW satellite stellar mass function for each halo prior. Group 1 halos require a higher level of feedback than Group 2 halos. When the model is extrapolated to work on the MW host halo, the stronger feedback in Group 1 halos prevents stars from forming more efficiently than in Group 2 halos. We can see this effect by applying the two different models on the same merger tree set (comparing the blue dashed with red solid lines or the blue solid with red dashed lines). The blue dashed line in Figure 9 shows the prediction of the model constrained for the Group 2 halo prior but applied to Group 1 halos. The difference between this line and the red solid line shows the effect of changing galaxy formation models that are constrained for different halo priors. This is the same effect that can be seen by comparing the blue solid line and the red dashed line.

In Figure 9, we also include the observational estimate of the MW galaxy stellar mass of Licquia & Newman (2015), who constrained the stellar mass of the MW galaxy to ${M}_{* }=(6.08\pm 1.14)\times {10}^{10}\,{M}_{\odot }$. Both of the model predictions (the solid red and blue lines) agree with the observational result quite well. The Group 2 halo prior seems to agree with the observational estimate better, as its predicted stellar mass distribution peaks right at the observational constraint, but the Group 1 halo prior is also consistent with the observational result. We stress that the predictions for the central galaxies are made with models that are specifically constrained to the MW satellite galaxies. No information from higher mass scales is used to constrain the models. The prediction could change if the physics that acts on the MW-sized galaxies is different from what we infer from the MW satellite galaxies. We also note that we have fixed the halo mass in this paper. The Group 1 halo prior would achieve a better match with the data if the host halo of the MW has a higher mass than we have assumed here.

In this paper, we have analyzed a suite of N-body cosmological simulations of MW-mass halos with various mass-assembly histories and concentrations. The analysis shows that these properties of a host halo are correlated with the probability for a halo to host high-mass subhalos. For fixed halo mass, halos that have rapid recent accretion, characterized by the fitting parameter $\beta +\gamma$ of the mass-assembly history, tend to host high-mass subhalos that still exist in the present day. These halos also tend to have a lower concentration. The result suggests that the existence of high-mass satellite galaxies may provide a clue to the formation history and the structure of the underlying host dark matter halo. By applying a semi-analytic model to halos with different mass-assembly histories, we show herein that the best-fit galaxy formation models depend sensitively on the host mass-assembly history. Therefore, it may be possible to infer the mass-assembly history of hosts, and the properties of hosts that correlate with the assembly history, based on properties of their satellite galaxies.

In this work, we focused specifically on the MW and its satellite system, especially the unusually high-mass satellites, the Magellanic Clouds. As suggested by recent observational estimates using different probes, the dynamical mass of the LMC is at least $1.7\times {10}^{10}\,{M}_{\odot }$ (van der Marel & Kallivayalil 2014) or even higher (e.g., Vera-Ciro et al. 2013; Peñarrubia et al. 2016). The high masses of the MCs have important implications for the structure of the MW's dark halo because this would suggest that the host halo of the MW is more likely to have experienced rapid accretion since $z=0.5\sim 1$ (e.g., Besla et al. 2010) and that it may have relatively lower concentration than typical halos at the same mass (i.e., Group 1).

This argument of a low value for the MW halo concentration seems to be, at face value, in tension with some results in the literature. Observations of kinematic tracers of the MW can provide direct constraints of the dynamical mass distribution of the MW halo. Some results suggest that the MW halo appears to have a concentration value as high as or even higher than 18 (e.g., Battaglia et al. 2005; Deason et al. 2012). One possible solution to the tension is that the MW halo is a significant outlier or an early formed halo with the existence of the MCs as a "transient" coincidence  (e.g., Deason et al. 2016). We stress, however, that the formation of a central galaxy in the halo center can cause a contraction of the dark matter in the inner halo. As we have shown in Section 3, it is still plausible that the MW has formed in a low-concentration halo, but the effects of the baryonic components of the MW give rise to a centrally peaked rotation curve, which can be mimicked by a halo with a much higher concentration value (though this may be in tension with Binney & Piffl (2015), who argue against halo contraction). Moreover, if the mass of the MW host halo is higher than that assumed in our simulations, the low-concentration rapid-accretion hosts could achieve better fits to the MW kinematic constraints (e.g., Wang et al. 2012; Cautun et al. 2014a).

These disparate conclusions about the MW halo's concentration can be solved with more data and better modeling of the dynamical effects of the Magellanic Clouds on the MW. Our analysis suggests that more accurate measurements of the MW halo kinematics within $r\lt 60\,\mathrm{kpc}$ are still needed to tighten the constraints on the radial mass distribution of the MW halo. We also note that, with dark-matter-only simulations, Busha et al. (2011a) found that MW analogs, when selected based on not only the virial mass but also the phase-space position of the MCs, tend to have slightly higher concentration values ($c\sim 11$). It is worth noting, however, given the close distances of the MCs to the Galactic center, the observed speeds of the MCs are likely affected by the baryonic component of the MW Galaxy. Also, the observed speeds of the MCs can be biased if the MCs are not yet in equilibrium orbits. As such, more detailed modeling of the kinematics of the MCs is needed to understand what constraints on the host halo concentration can be derived.

The connection between the host halo properties and the subhalo populations also has interesting implications for the properties of the satellites of the MW and the stellar mass of the MW galaxy itself. To understand these implications, we constrained a flexible galaxy formation model with the observed MW stellar mass function of classical dwarf galaxies with different priors on the formation history and concentration of the host halo. We studied how different halo priors affect the inference of galaxy formation physics and, more interestingly, inferred what observations can help distinguish between different host halo priors and constrain the feedback processes.

First, we found that when assuming different host halo priors for the model, the observed stellar mass function requires a different strength of and circular velocity dependence for feedback (Figure 4 and Table 1). After the uncertainties in relevant processes of the model are marginalized, we find that host halos with rapid recent accretion require stronger feedback than those without rapid recent accretion to fit the satellite galaxy mass function. The reason for this is because halos with rapid recent accretion tend to host high-mass subhalos, which then requires a lower stellar-mass-to-halo-mass ratio to match the observed mass function. Moreover, our model inference from the MW stellar mass function weakly prefers the host halo prior that the MW host halo had rapid recent accretion and has lower concentration than typical halos with a similar virial mass. However, due to large uncertainties in the model, especially the uncertainties in the feedback processes, our model inference cannot completely rule out the alternative prior that the MW host halo has slow recent accretion and higher concentration. This result suggests that, in order to determine the properties of the host halos using galaxy formation models, one has to better understand the feedback processes in galaxy formation.

Second, we found that to match the observed stellar mass function, different host halos require a significantly different mapping between subhalo mass or circular velocity and stellar mass, or the ${V}_{\max }\mbox{--}{M}_{* }$ relation (Figure 6). The difference in the subhalo maximum circular velocity, ${V}_{\max },$ for a given stellar mass can be as large as $\sim 25\,\mathrm{km}\,{{\rm{s}}}^{-1}$ between different simulated host halos. The scatter from the host halo prior, however, is not large enough to account for the observationally inferred low ${V}_{\max }$ values of the MW classical dwarfs (Boylan-Kolchin et al. 2011b), at least for the fixed halo mass employed in this paper. Moreover, this difference in the ${V}_{\max }\mbox{--}{M}_{* }$ relation between different simulated host halos mainly exists in satellites more massive than ${M}_{* }\sim {10}^{7}\,{M}_{\odot }$, and vanishes for satellites with masses lower than $\sim {10}^{6.5}\,{M}_{\odot }$. This is because the difference in the subhalo mass function of different host halos of a given mass only exists in the very high-mass end, where statistical fluctuation matters. When looking at the number counts of low-mass satellite galaxies, the MW is not a significant statistical outlier because it is for high-mass satellites (e.g., Strigari & Wechsler 2012). Ideally, a stronger test of the CDM cosmology with local group satellite galaxy kinematics should use very-low-mass galaxies, in which the subhalo population suffers less from the statistical fluctuations inherent from the halo-to-halo scatter (Purcell & Zentner 2012; Jiang & van den Bosch 2015).

Third, we find that the feedback processes assumed in the model strongly affect the mass–metallicity relation of dwarf galaxies (Figure 7). This is because feedback not only suppresses star formation but also governs the flow of baryons, including metals, produced during galaxy formation (e.g., Lu et al. 2015a, 2015b). When a model requires stronger outflow to overcome the deeper gravitational potential of larger subhalos, it inevitably results in a lower metallicity in dwarf galaxies unless the metallicity of the outflow is different from that of the ISM (Dalcanton 2007). Therefore, one can find clues about the properties of the dark matter halo potential by modeling the strength and metallicity of galactic outflow.

Fourth, we find that halos with different MAHs, concentrations, and subhalo populations predict different numbers of small satellite galaxies that are accreted together with the MCs. When constrained to reproduce the observed satellite galaxy stellar mass function, the host halos that have recent rapid accretion and lower concentration tend to host more satellites associated with the most massive satellite. The result suggests that finding and identifying low-mass satellite galaxies down to ${M}_{* }\sim {10}^{3}\,{M}_{\odot }$ that are physically or historically associated with high-mass satellites may provide useful constraints on the property of the host halo (see also Deason et al. 2015). The ongoing surveys for low-mass satellite galaxies within the MW and especially around the vicinity of the MCs (e.g., Bechtol et al. 2015; Drlica-Wagner et al. 2015) may provide useful data for further theoretical investigations.

Finally, assuming the model constrained to observations of MW satellite galaxies can be extrapolated to the MW host halo, we find that the host halo properties that are connected to the satellite galaxy populations also affect the formation of central galaxies. Since halos with lower concentration and more rapid recent accretion typically form later and require stronger feedback to match the satellite stellar mass function, these halos form a central galaxy later and with lower stellar mass than halos with higher concentration and slower recent accretion. These effects explain the phenomenology of halo occupation and abundance matching models in which at fixed mass low-concentration halos are required to host a central galaxy with lower stellar masses (e.g., Reddick et al. 2013; Lehmann et al. 2015; Hearin et al. 2016; Zentner et al. 2016), which is attributed to the effect of the assembly bias (Gao et al. 2005; Wechsler et al. 2006; Jing et al. 2007). While the effects of the assembly bias on galaxy populations have been studied with SAMs on large-scale survey data in previous works (e.g., Croton et al. 2007; Wang et al. 2013), we have demonstrated in this paper that, at a fixed halo mass, the formation history of the host halo can affect both the central galaxy and the satellite populations simultaneously.

The interplay between MAH, host halo properties, and galaxy evolution has important consequences beyond the MW. Interestingly, the Andromeda galaxy also has massive satellites, and thus may have rapid recent accretion too. It would be interesting to investigate how probable it is that in CDM cosmology two nearby L* galaxies both host similar high-mass satellite galaxies, and how the accretion of high-mass subhalos correlates with the large-scale environment. These studies will shed further light on the position of our own galaxy, the MW, in the cosmic web.

We thank Matthew Becker for providing the c125-2048 simulation, which was run using computing resources at SLAC. The resimulations were performed using computational resources at SLAC and at NERSC. We thank the SLAC computational team for their consistent support. We also acknowledge the Ahmanson Foundation for providing computational resources used in this work. Support for programs HST-AR-13270.003-A, HST-AR-13896.005-A, and HST-AR-13896.009-A 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. Y.L. thanks Evan Kirby, Josh Simon, Robyn Sanderson, Shude Mao, and Simon White for useful discussions. S.T. was supported by the Alvin E. Nashman Fellowship in Theoretical Astrophysics. A.R.W. was supported by a Moore Prize Fellowship through the Moore Center for Theoretical Cosmology and Physics at Caltech and by a Carnegie Fellowship in Theoretical Astrophysics at Carnegie Observatories.

We study the probability distribution of the MW satellite mass function predicted by the N-body simulations and the SAM adopted in this paper. Two alternative models for the probability distribution function are tested here. One is the Poisson probability distribution, and the other is the NBD. Boylan-Kolchin et al. (2010) found that the NBD provides a more accurate description for the distribution of the subhalo mass function predicted by N-body simulation. Similar findings have also been reported in independent studies (e.g., Busha et al. 2011a; Cautun et al. 2014b). To choose an accurate model for the distribution function of the satellite galaxies predicted by a galaxy formation model, we apply a fiducial SAM to a set of halo merger trees of MW-mass halos extracted from the cosmological simulation c125-2048.

We first extract the subhalo mass function from the simulation. Using a fiducial SAM, which was tuned to match the field galaxy stellar mass function of the local universe (Lu et al. 2014), we then make a prediction for the satellite galaxy stellar mass function. In Figure 10, we show the subhalo mass function on the left panel and the SAM predicted satellite stellar mass function on the right panel. Each blue horizontal box covers the lower to upper quartile values of the mass of the $i\mathrm{th}$ subhalo (satellite), with a middle line indicating the median of the mass distribution. The whiskers extending from the boxes show the range of the masses in the sample. On the other hand, the gray shaded region covers the 1σ range of the cumulative distribution function $N(\gt M)$ as a function of $\mathrm{log}\,M$. In comparison, we also show the expected 1σ scatter of the mass functions with red vertical error bars assuming the mass functions are Poisson distributed. We will show later that the real distribution of the mass functions quantitatively deviates from Poisson. Moreover, we also overplot the observed MW satellite stellar mass function in the right panel. Without further tuning the model against the data, the fiducial model prediction is in good agreement with the data, with most of the data points encompassed by the $1\sigma$ region of the predicted distribution, leaving only two outliers with a small deviation. Although the fiducial model does not perfectly match the data, it is sufficiently close. We can use it as a representative model to explore the statistical behavior of predicted satellite populations.

Zoom In Zoom Out Reset image size

The probability distribution function of the number count of subhalos is found to be better described by the NBD (Boylan-Kolchin et al. 2010),

$$\begin{eqnarray}&&P{(N| r,p)=\displaystyle \frac{{\rm{\Gamma }}(N+r)}{{\rm{\Gamma }}(r){\rm{\Gamma }}(N+1)}{p}^{r}(1-p)}^{N},\end{eqnarray} \tag{ 6 }$$

where N is the number of subhalos per host in a given mass range, ${\rm{\Gamma }}(x)\equiv (x-1)!$ is the Gamma function, and r and p are two parameters. This distribution function has also been used to describe the probability of the number of satellite galaxies in HOD models (e.g., Berlind & Weinberg 2002). Boylan-Kolchin et al. (2010) argued that the reason for this model to better describe the distribution of the subhalo mass function than the Poisson distribution is because it captures the intrinsic non-Poisson scatter of the mass function. The authors defined a parameter ${s}_{{\rm{I}}}\equiv {\sigma }_{{\rm{I}}}/\mu$, the fractional scatter from the intrinsic scatter, ${\sigma }_{{\rm{I}}}$, with respect to the Poisson scatter, μ. The two parameters in the NBD are then determined as

$$\begin{eqnarray}&&p=\displaystyle \frac{1}{1+{s}_{{\rm{I}}}^{2}\mu },r=\displaystyle \frac{1}{{s}_{{\rm{I}}}^{2}}.\end{eqnarray} \tag{ 7 }$$

In Figure 11, we show the distribution of the subhalo number count of the MW-mass hosts in three different subhalo mass bins, which are noted in each panel. In agreement with previous studies, we also find that the NBD describes the subhalo number count distribution much better than the Poisson distribution. Boylan-Kolchin et al. (2010) found that ${s}_{{\rm{I}}}\approx 0.18$ yields a good fit to the subhalo mass function distribution predicted by the Millennium-II Simulation. We find, however, that a larger value for ${s}_{{\rm{I}}}$ is needed to fit the distribution predicted by our simulations. In the figure, we have used ${s}_{{\rm{I}}}=0.26$. When the expected number counts become smaller ($\bar{N}\lt 4$), the NBD approaches the Poisson distribution, making it hard to distinguish between the two models. We also show the satellite galaxy number count distribution at three different stellar masses. Same as the subhalo mass function distribution, the number count distribution for given stellar mass can also be accurately described by the NBD. Busha et al. (2011a) found that adding an exponential tail to the NBD can better capture the distribution at very high N. Because the tail only covers a small fraction of the probability distribution, we ignore this part in order to keep the model simple without losing accuracy above the level our inference can capture. Using the fiducial SAM, we predict stellar masses for the satellite galaxies and compute the stellar mass function for each MW host. We show the distribution of the number of satellite galaxies in three stellar mass bins in the lower panel of Figure 11. Again, we use ${s}_{{\rm{I}}}=0.26$ for the NBD plotted in the figure to compare with the model predictions and find that the NBD matches the simulated distribution of the stellar mass function remarkably well.

Zoom In Zoom Out Reset image size

We now quantitatively test which model better matches the moments of the simulation predicted distribution of the subhalo and satellite mass functions. Using the large sample provided by the cosmological simulations, we compute the second, third, and fourth order central moments of the subhalo mass function and the stellar mass function for given mass bins. The $n\mathrm{th}$ order central moment is defined as

$$\begin{eqnarray}&&{\mu }_{n}=E[{(x-E[x])}^{n}],\end{eqnarray} \tag{ 8 }$$

where function $E[x]$ is the expectation of x. 1133 MW halos from the simulation are used. The simulation predicted moments for the subhalo mass function and the stellar mass function as functions of masses are shown as green lines in the upper panels and the lower panels of Figure 12, respectively. The moments monotonically increase with decreasing mass. Based on the hypothesis that the distribution function is Poisson or NBD, we generate the same number of Monte Carlo mass functions, assuming each mass function is a random realization of the assumed distribution function. We adopt the simulation predicted means to assign the expectation value for the Monte Carlo mass functions. Using the Monte Carlo mass functions, we can compute the same central moments for each given mass bin. We then replicate the Monte Carlo simulation 10,000 times to obtain 10,000 samples of the moments for each mass bin. We use these 10,000 values to construct a reference distribution for each moment for any given mass bin. In Figure 12, we show the distribution of the Monte Carlo subhalo mass function and the satellite stellar mass function assuming the samples are Poisson distributed (blue) or follow the NBD (red). As one can see, the simulation predicted moments are largely deviated from the moments based on the Poisson distribution, especially in the low-mass bins. These Monte Carlo simulations suggest that while the simulation predicted distribution is still consistent with Poisson at high-mass bins, it is clearly inconsistent with Poisson for low-mass bins. In contrast, the NBD moments match the moments of the simulated mass functions remarkably well.

Zoom In Zoom Out Reset image size