A BARYONIC SOLUTION TO THE MISSING SATELLITES PROBLEM

The cosmological paradigm based on cold dark matter (CDM) and dark energy (Λ) has been extremely successful in describing the observed evolution and large-scale structure of our universe. At small scales, however, a number of observations seem to be at odds with the predictions of ΛCDM cosmology. In particular, over a decade ago, it was pointed out by Moore et al. (1999) and Klypin et al. (1999) that the number of high-mass subhalos predicted by high-resolution CDM simulations exceeds the observed number of luminous satellites of the Milky Way (MW) by at least an order of magnitude. This has become known as the "missing satellites problem" (MSP). Although this problem has been mitigated to a degree by the discovery of a number of additional faint satellite galaxies (Willman et al. 2005; Irwin et al. 2007; Liu et al. 2008; Martin et al. 2008; Simon & Geha 2007; Grebel 2000; van den Bergh 2000; Belokurov et al. 2008; Watkins et al. 2009; Belokurov et al. 2010), there remains a considerable discrepancy between the number of observed MW satellites and the number predicted in CDM simulations.

Efforts to resolve this issue have fallen into two broad categories. First, there are proposals in which the star formation rate in satellite galaxies is suppressed, leading to large numbers of low-mass subhalos which are simply unobservable. Possible means for such suppression include photoevaporation resulting from ionizing radiation, e.g., at reionization (Quinn et al. 1996; Thoul & Weinberg 1996; Navarro & Steinmetz 1997; Klypin et al. 1999; Barkana & Loeb 1999; Gnedin 2000; Hoeft et al. 2006; Madau et al. 2008; Alvarez et al. 2009) or lower z due to blazars (Pfrommer et al. 2012), or due to cosmic ray heating (Wadepuhl & Springel 2011). Photoionization is expected to suppress star formation in halos below ∼10⁹ M_☉ (Okamoto et al. 2008). In more massive halos where gas is retained and star formation can begin, further suppression is expected from supernova feedback (Dekel & Silk 1986; Benson et al. 2002; Dekel & Woo 2003; Governato et al. 2007). It is also important to consider the number of faint satellites that remain undetected due to observational incompleteness (Willman et al. 2004; Simon & Geha 2007; Tollerud et al. 2008; Walsh et al. 2009; Koposov et al. 2009; Rashkov et al. 2012). Taken together, it is possible that observational incompleteness combined with suppression of star formation can bring the number of luminous satellites into line with the number of predicted satellites around an MW-mass galaxy.

However, while this combination of effects may explain the faint or low-mass regime of the MSP, Boylan-Kolchin et al. (2011) recently pointed out that it still fails at the massive end. Dubbed the "too big to fail" (TBTF) problem, this aspect of the MSP arises because the most massive subhalos in ultra-high resolution dark matter (DM)-only simulations of MW-analog galaxies (Diemand et al. 2008; Springel et al. 2008) are too dense to host the observed satellites of the MW (see also Wolf & Bullock 2012; Hayashi & Chiba 2012). The simulations always contain a population of subhalos (6–22, varying within the errors of the MW's measured mass) that are more massive than any of the dwarf spheroidals observed in the MW (Boylan-Kolchin et al. 2012). In other words, while the abundance of lower luminosity satellites may be made consistent with CDM predictions, the MSP remains a puzzle because simulations predict too many massive satellites.

A second class of possible solutions to the MSP considers departures from the standard assumptions of cold and collisionless DM. Whereas typical particle DM models predict that weakly interacting massive particles (WIMPs) will form halos with masses as small as 10⁻⁶ M_☉ or so (depending on the temperature at which the WIMPs undergo kinetic decoupling with the cosmic neutrino background), the formation of small-scale structure can be strongly suppressed if the DM particles are not entirely cold (Colín et al. 2000; Macciò & Fontanot 2010; Lovell et al. 2012). DM in the form of sterile neutrinos with masses of ∼1–10 keV have received attention within this context (Dodelson & Widrow 1994), although many other warm dark matter (WDM) scenarios could potentially accommodate a similar suppression of small-scale power (Pagels & Primack 1982; Hooper et al. 2007). WDM primarily addresses the low-mass end of the MSP by suppressing the formation of small-scale structure. It can also help at the high-mass end, however, because the delayed halo collapse times in a WDM cosmology result in lower concentrations and hence reduced central densities (Lovell et al. 2012).

Self-interacting dark matter (SIDM) models (Carlson et al. 1992; Spergel & Steinhardt 2000; Loeb & Weiner 2011) provide another mechanism to achieve the same effect. In this case, the interactions prevent the formation of the steep central density cusps that are the hallmark of CDM halos. Density profiles in SIDM instead exhibit a central core (Vogelsberger et al. 2012; Rocha et al. 2012; Zavala et al. 2012), which helps to address the MSP in two ways: cored halos are more susceptible to tidal disruption, potentially removing many of the excess low-mass halos; second, in the surviving halos, the core reduces the central densities, alleviating the TBTF problem (Vogelsberger et al. 2012). Furthermore, there is considerable observational evidence for the existence of DM cores in the centers of low surface brightness galaxies (Kuzio de Naray et al. 2008; de Blok 2010; Oh et al. 2011) and in at least two MW dwarf spheroidal satellites (Walker & Peñarrubia 2011).

In this work, we consider a new type of solution to the MSP in which the shape of satellite DM profiles are altered, but as a result of baryonic physics, rather than through modifications in the particle physics sector. This model incorporates suppression of star formation from photoionization and supernova feedback, but additionally considers the tidal effects due to the presence of a baryonic disk, which is not found in DM-only simulations. The loss of gas in satellite halos through tidal stripping has been proposed as a means to limit star formation (Strigari et al. 2007). However, tidal stripping in a cuspy halo by itself does not reduce the central densities of the most massive subhalos enough to bring them into agreement with the observed kinematics of the MW satellites (Read et al. 2006a; Boylan-Kolchin et al. 2012). An additional modification to the central densities of the most massive satellites appears to be required.

It has become broadly accepted in recent years that baryonic processes can alter the distributions of DM within halos, potentially steepening (Blumenthal et al. 1986; Gnedin & Zhao 2002; Gnedin et al. 2011) or flattening (Navarro et al. 1996; Read & Gilmore 2005; Mashchenko et al. 2006, 2008; Governato et al. 2010, 2012; Pasetto et al. 2010; de Souza et al. 2011; Cloet-Osselaer et al. 2012; Macciò et al. 2012; Pontzen & Governato 2012; Teyssier et al. 2012) the profile depending on the mass of the halo in question, and on the strength of various feedback mechanisms. Recently, the effects of baryons on the DM halo profiles of luminous satellite galaxies was studied by Zolotov et al. (2012). These authors found that satellites with M_star ≳ 10⁷ M_☉ tend to develop "cored" density distributions through supernova feedback prior to infall (here, cored refers to any inner density slope, γ, shallower than −1, because γ becomes flatter with increasing stellar mass, see Governato et al. 2012). The flatter density profiles make them more vulnerable to tidal effects of the baryonic disk (see also Taylor & Babul 2001; Read et al. 2006a; Choi et al. 2009; Wetzel & White 2010; Romano-Díaz et al. 2010; D'Onghia et al. 2010a; Peñarrubia et al. 2010). Brooks & Zolotov (2012) showed that the presence of the disk increases tidal stripping for galaxies of all masses compared to the DM-only case, but most strongly reduces the central densities of cored satellites (see also Stoehr et al. 2002; Hayashi et al. 2003; Peñarrubia et al. 2010), which are the most luminous. In other words, baryonic physics significantly reduces the observable mass of satellites. Brooks & Zolotov (2012) demonstrated that these combined baryonic effects produce a z = 0 satellite distribution with luminosities and kinematics comparable to the MW and M31.

In this work, we apply the results of Zolotov et al. (2012, hereafter Z12) to the subhalo populations found in the DM-only Via Lactea II (VL2) simulation. We find that after taking into account baryonic effects, the number of surviving massive subhalos is strongly reduced. Adopting a stellar mass to halo mass relation, we show that the number of luminous satellites is also compatible with observations. This suggests that a proper accounting of baryonic effects reduces the predicted number of massive, luminous subhalos that should survive in CDM, bringing the observations and theory into agreement without the need of invoking WDM or DM with exotic properties or interactions.

We make use of the publicly available subhalo catalogs⁵ from the VL2 DM-only simulation (Diemand et al. 2008), supplemented by additional information extracted from the simulation outputs, consisting of (1) mass enclosed within 1 kpc at z = 0, (2) infall times and the maximum value in their rotation curves, v_max, at infall, and (3) their full orbital information (Kuhlen et al. 2012), including number and distance of apo- and pericenter passages. The VL2 halo was run using a Wilkinson Microwave Anisotropy Probe (WMAP) year 3 cosmology (Spergel et al. 2007). The halo was selected to have no low z major mergers, thought to be similar to the merging history of the MW. The z = 0 halo has roughly 450 million particles with masses of 4100 M_☉ within its virial⁶ radius of 402 kpc at z = 0, yielding a halo mass of 1.9 × 10¹² M_☉.

Starting from all subhalos within R_vir at z = 0, we restrict our study to only those subhalos with v_peak (the largest value of v_max over the subhalo's entire history) >13.3 km s⁻¹. This v_peak value was chosen to correspond to 10⁵ M_☉ in stellar mass, or a V-band absolute magnitude of −7 (see Section 3.3). This allows our sample to cover the full range of luminosities of the classical dwarf satellites. Selecting for v_peak > 13.3 km s⁻¹ yields a sample of 410 subhalos. It is commonly assumed that all halos with v_peak < 20 km s⁻¹ are strongly affected by UV heating, unbinding their gas (e.g., Okamoto et al. 2008). This process is expected to either leave them entirely devoid of stars ("dark"), or very faint if they were able to form a few stars prior to reionization (e.g., Bovill & Ricotti 2009; Salvadori & Ferrara 2009; Li et al. 2010; Wolf et al. 2010), and prone to be missed by current surveys. In Section 3.2, we consider the impact of UV heating on the VL2 subhalos, and show that, indeed, most of these low-mass subhalos remain dark. However, a few of the satellites in this mass range should be capable of forming stars.

The v_max-function of our final subhalo sample (i.e., those halos that survive to z = 0 and had v_peak > 13.3 km s⁻¹) is shown in Figure 1, for various times. The solid line shows the cumulative v_peak-function for the halos in our sample, while the short-dashed curve shows v_max at infall. We note that v_peak and v_max at infall, v_infall, can be slightly different for these subhalos. The majority of the VL2 subhalos in our sample reach v_peak at high z and then grow very little, with an overall slight decrease (<10%) in v_max. A few subhalos undergo a substantial reduction in v_max between high z and infall, presumably due to encounters with other halos that strip their mass. The resulting velocity function at z = 0 is shown in Figure 1 as the dotted line. The circular velocity at 1 kpc, v_1 kpc, is shown for comparison (long-dashed line). In the following, we adopt v_1 kpc to examine the effect of baryons on the central mass distribution of satellite galaxies.

Zoom In Zoom Out Reset image size

Z12 proposed a correction to the v_1 kpc values of z = 0 DM-only subhalos to account for the missing baryonic physics that lowers the central masses of luminous subhalos,

$$\begin{equation} \Delta (v_{\rm 1\,kpc}) = 0.2v_{\rm infall}-0.26 \, {\rm \,{\rm km\, s}^{-1} }. \end{equation} \tag{ 1 }$$

For subhalos with v_infall < 30 km s⁻¹, the Z12 correction is designed to account for a reduction in central mass due to (1) loss of gas, either due to UV heating before infall or stripping after infall and (2) the tidal effect of the baryonic disk (which does not exist in DM-only runs). These processes should act on a subhalo even if it is too low mass to retain gas at reionization, or massive enough to be luminous. These processes should also occur irrespective of whether the subhalo has a cuspy or a cored DM density profile. However, neither are typically accounted for in DM-only simulations. In subhalos with v_infall > 30 km s⁻¹, Z12 found that enough star formation takes place to significantly flatten the DM density profiles prior to infall. Hence, at v_infall > 30 km s⁻¹, the Z12 correction accounts for (1) and (2) as in the lower mass case, but also for an additional reduction in the central masses of the satellites due to supernova feedback and the enhanced tidal stripping that occurs as the central density profiles of massive satellites become more shallow.

In Z12, Equation (1) was derived using halos with 20 km s⁻¹ < v_infall < 50 km s⁻¹. Hence, it has not been tested down to the lower v_infall of our lowest mass VL2 subhalos. However, we will apply it to all subhalos with v_infall < 50 km s⁻¹, and in Section 3 we show that in all cases where the correction yields unphysical (negative) results, the unphysical halos lose more than 97% of their mass at infall, and should be considered completely destroyed in the tidal presence of a baryonic disk. At v_infall > 50 km s⁻¹, subhalos are Magellanic-like and gas-rich at accretion, possibly including an additional effect of adiabatic contraction that is not accounted for in the correction. There are five massive VL2 subhalos with v_infall > 50 km s⁻¹, for which the Z12 correction is not applied.

Most of the MW and M31 dwarf spheroidal galaxies (dSphs) have central v_c values <20 km s⁻¹ (measured at the half-light radii, which are all ≲1 kpc; McConnachie 2012). At z = 0, the VL2 host halo contains 28 satellites with v_1 kpc > 20 km s⁻¹, grossly inconsistent with the observational results. After applying the correction of Z12, there are only five satellites with v_1 kpc > 20 km s⁻¹. This is a first indication that baryonic physics, at least as implemented in the simulations of Z12, appears to be quite capable of reducing the mismatch between the central densities of the most massive subhalos in DM-only simulations and the observed kinematics in the classical MW dwarf satellite galaxies.

It is clear that the baryonic correction proposed by Z12 dramatically reduces the number of massive halos expected in a DM-only MW-mass run (e.g., from 28 to 5 for satellites with v_1 kpc > 20 km s⁻¹ in VL2). However, it is less clear if it reduces the number of luminous subhalos. That is, does the correction take 20 Fornax-like satellites that initially have v_1 kpc > 20 km s⁻¹ at z = 0 and simply shift them to v_1 kpc < 20 km s⁻¹? If so, then there may be fewer massive halos, but the overall number of luminous satellites would still be much larger than observed in the MW or M31. In this section, we explore whether both the mass and luminosity of MW-mass galaxies can be reproduced by accounting for baryonic effects.

3.1. Destruction Rates

We must first consider whether the VL2 halos should have all survived to z = 0, or whether it is likely that baryonic physics would have led to their destruction. Peñarrubia et al. (2010) examined the effect of a baryonic disk on both cored and cuspy subhalos. They showed that tides from the disk are dominant in the center of cored subhalos, leading to much more mass loss than DM-only runs predict. Even cuspy subhalos will undergo more mass loss when a baryonic disk is present if they are on highly elliptical orbits. Importantly, cuspy halos tend to survive even after substantial (99.99%) mass loss, while cored halos can be completely disrupted. This implies that there could be a large number of satellites that survive in the DM-only VL2 run that would not survive in a baryonic run.

In what follows, we assume that a VL2 subhalo that has lost a certain fraction of its mass since infall should have been fully disrupted if baryons had been included, and should no longer appear as a bound luminous satellite at z = 0. To set the limiting fraction, we refer to Figure 2 in Peñarrubia et al. (2010), where tidal effects of a disk on cored dwarf galaxies are compared to those on cuspy halos with no disk present. Adopting their "mixed" model, in which a cored (γ = 0) subhalo evolves within a parent halo with a cuspy (γ = −1) density profile, we find that cored subhalos that have pericenters ≲20 kpc undergo 99.9% mass loss, which we consider to be completely disrupted. Without the disk presence, the results in Peñarrubia et al. (2010) show that a cuspy subhalo with a pericenter of 20 kpc would experience 90% mass loss. Hence, we assume that any VL2 subhalos with pericenters <20 kpc should experience strong disk tides that need to be accounted for. We conclude that if a VL2 DM-only satellite has lost more than 90% of its mass after infall, and has had pericentric passages that take it within 20 kpc of the parent halo's center, then the same satellite in a baryonic run is extremely likely to be fully disrupted by tides.

To be conservative, we apply this limit only to subhalos that have v_infall > 30 km s⁻¹. Z12 and Governato et al. (2012) demonstrated that cored DM density profiles exist in satellites with a stellar mass of more than 10⁷ M_☉, corresponding to halos with v_infall > 30 km s⁻¹, and consistent with the energy arguments explored in both Boylan-Kolchin et al. (2012) and Peñarrubia et al. (2012). While it is possible that halos at lower masses have cores (but at smaller radii that are below the resolution limit of current cosmological simulations), we assume they may remain cuspy, and 90% mass loss may then not be enough to fully destroy them. Instead, for the remainder of the halos with v_infall < 30 km s⁻¹, we adopt the results of Wetzel & White (2010) and assume that a halo must have lost 97% of its mass since infall to be fully stripped, and no longer appear as a bound luminous satellite at z = 0. We verified that the majority of the halos that have lost 97% of their mass have tidal radii less than 1 kpc, suggesting that their inner luminous regions should indeed be stripped. Adopting 97% is conservative, as both Z12 and Peñarrubia et al. (2008) found that halos begin to have their stars stripped after losing 90% of their halo mass after infall.

We use the difference in v_max at infall and v_max at z = 0 for a given VL2 subhalo to estimate the amount of mass that it has lost since infall. Following the results of Peñarrubia et al. (2010),

$$\begin{equation} \frac{v_{{\rm max}}(z=0)}{v_{{\rm infall}}} = \frac{2^{0.4}x^{0.3}}{(1+x)^{0.3}} \end{equation} \tag{ 2 }$$

for a subhalo with a DM density slope, γ = −1.0, where x ≡ mass(z = 0)/mass(z = infall). A slope of −1 is roughly the slope in the central regions for all of the halos formed in DM-only simulations (Navarro et al. 1997; Springel et al. 2008), and we verified that the inner density slopes in the VL2 subhalos are also consistent with −1. We use the above equation to find (1) all of the subhalos that lose more than 97% of their infall mass (x = 0.03) and (2) all of the subhalos with v_infall > 30 km s⁻¹ that lose more than 90% of their infall mass (x = 0.1) and have pericentric passages under 20 kpc. We consider these two populations of subhalos to be "destroyed."

3.2. Identifying Dark Subhalos

As we will show below, even after considering destruction due to tidal effects on the VL2 subhalos, 307 satellites remain at z = 0 from our original sample of 410 with v_peak >13.3 km s⁻¹. Roughly 100 of the surviving subhalos should be as bright as the Draco dSph (see next section), and thus luminous enough to have been detected, but 100 such luminous satellites is much larger than the number observed in the MW or M31 to date. Yet many of these satellites are in the halo mass range that is expected to be strongly affected by UV heating, and some of these subhalos should therefore be inefficient at forming stars and will remain dark. In this section, we apply the results of Okamoto et al. (2008) to identify the subhalos that may remain dark.

Okamoto et al. (2008) use hydrodynamical simulations to identify the characteristic halo mass, M_char, that retains 50% of the cosmic baryon fraction, f_bar, as a function of z. Their simulations adopt a uniform ionizing background that accounts for H and He i reionization at z = 9, and He ii reionization at z = 3.5. We have converted their M_char(z) results into v_char(z), adopting a WMAP3 cosmology and using an overdensity of 200ρ_crit. We assume that if a VL2 subhalo has a v_peak value above v_char(z_peak), then it retains enough baryons to be luminous. Note that v_char(z) is a virial quantity, but we do not have the virial masses of the VL2 subhalos. Rather, we have v_max values. Most satellites have 1.2 < v_max/v_vir < 1.8, corresponding to a range of concentrations⁷10 < c < 40, with a mean v_max/v_char = 1.4 for a halo with c = 20 (Bullock et al. 2001; Prada et al. 2012).

Figure 2 illustrates the affect of applying this model to the VL2 subhalo population. Diamonds show the v_peak(z_peak) values of the VL2 satellites that survive after our tidal stripping considerations. To put the v_char(z) of Okamoto et al. (2008) into v_max space for comparison with the VL2 subhalos, we have multiplied v_char(z) by 1.4 to derive the solid line in Figure 2 (i.e., we assume v_max/v_char = 1.4, a typical value for a subhalo). The shaded region surrounding the solid line shows the full range of 1.2 < v_max/v_char < 1.8.

Zoom In Zoom Out Reset image size

Clearly, the number of subhalos above v_char depends sensitively on the concentrations (and hence v_max/v_char) of individual subhalos. Adopting v_max/v_char = 1.4 (1.2) leads to 40 (70) surviving, luminous satellites at z = 0. Certainly, scatter in concentration is expected among the individual VL2 halos, which we have not accounted for. Furthermore, concentrations at these halo masses are expected to decrease with increasing z (Bullock et al. 2001; Eke et al. 2001; Zhao et al. 2003; Macciò et al. 2008; Prada et al. 2012). Again, we do not account for the change in concentration with z, but note that the trend would allow for lower v_peak halos to lie above v_char with increasing z. The overall determination of whether a particular halo is likely to be luminous will depend on the halo's formation and evolution. For the purposes of this work, we assume that all of the halos above the solid line in Figure 2 are likely to be luminous.

3.3. Assigning Luminosities

Finally, we assess whether the surviving satellites are luminous by using the v_infall − M_star relation from Z12 ⁸ to assign stellar masses to the VL2 subhalos,

$$\begin{equation} \frac{{M}_{\rm star}}{{M}_{\rm \odot }} = 0.018 \left(\frac{v_{\rm infall}}{\,{\rm km\, s}^{-1} } \right)^6. \end{equation} \tag{ 3 }$$

We use the tight log (M_star) − M_V for the subhalos,

$$\begin{equation} \log _{10} \left(\frac{{M}_{\rm star}}{{M}_{\rm \odot }} \right) = 2.37\hbox{--}0.38M_V, \end{equation} \tag{ 4 }$$

in Z12 to further assign V-band magnitudes, M_V. Z12 showed that this relation produced simulated satellite luminosity functions that were in good agreement with the classical dSph populations of the MW and M31. Additionally, Munshi et al. (2012) have shown that the simulations that Equations (3) and (4) are drawn from yield an excellent match to the z = 0 stellar-to-halo mass relation from Moster et al. (2013). Most importantly, Governato et al. (2012) demonstrated that these mass-to-luminosity relations produced excellent agreement with field galaxies in the same luminosity range as the satellites we examine here. In other words, despite the fact that past simulations have overproduced stars (Zolotov et al. 2009; Guo et al. 2010; Sawala et al. 2011; Brooks et al. 2011; Leitner 2012; Moster et al. 2013), the simulations used to derive Equations (3) and (4) are in excellent agreement with observed stellar-to-halo mass relations. There is no indication that the luminosities predicted by these relations are too bright, though we stress that it has not been tested fainter than the luminosity range of the classical dSphs.

Zoom In Zoom Out Reset image size

Figure 3 shows the resulting M_V and v_1 kpc for the VL2 satellites. The top panel shows the results directly from the VL2 catalog at z = 0, while the bottom panel shows the results after the destruction, heating, and velocity corrections considered in this paper. Filled red data points are those subhalos that are likely to survive tidal effects and retain enough baryons to be luminous, circles with "x" through them are those that are likely to have been destroyed by tidal effects in the presence of baryons, and empty circles are those that fall below the characteristic mass to retain baryons and form stars (see Figure 2). The filled black circles in Figure 3 identify a population of satellites that have lost more than 90% of their mass since infall, but do not meet the destruction criteria that we outlined above. Peñarrubia et al. (2008) found that halos that lose more than 90% of their mass begin to strip stars (see also Z12). Hence, the luminosities of the black points in Figure 3 should be considered upper limits, as these halos will likely have had stars stripped and be fainter than at infall.

Most of the brightest (M_V < −12) subhalos meet our criteria for destruction. We verified that with the exception of one massive satellite accreted at z = 0.5, all of these luminous, destroyed satellites were accreted z > 1.5. The majority have v_peak > 40 km s⁻¹, and by definition have orbits that take them within 20 kpc of the center of the parent halo. This is consistent with the idea that early, massive satellites contribute to the growth of the inner stellar halo in MW-mass galaxies (Bullock & Johnston 2005; Zolotov et al. 2009). In DM-only runs that neglect the presence of the disk, some of the cuspy inner remnants of these early, massive satellites are capable of surviving to z = 0.

We note that M_V = −7 corresponds to v_infall = 13.3 km s⁻¹. The subhalos fainter than M_V = −7 in Figure 3 are those that had v_peak > 13.3 km s⁻¹ but a v_infall < 13.3 km s⁻¹. Had we adopted v_peak to assign stellar masses, all subhalos in Figure 3 would be brighter than M_V = −7. Thus, the luminosities below M_V ∼ −9 should be taken with caution, but we verified that using v_peak instead of v_infall had almost no impact on the luminosities of the satellites brighter than M_V = −9. Hence, whether we adopt v_peak or v_infall has no effect on the number, masses, or luminosities of surviving, luminous subhalos brighter than M_V ∼ −9. This luminosity corresponds to the faintest of the classical MW dwarf satellites. Hence, use of v_infall does not change our conclusions regarding whether baryonic physics can address the massive MSP.

Before applying the Z12 correction, the VL2 run contains more than 20 luminous satellites that have v_1 kpc > 20 km s⁻¹, completely inconsistent with the satellite population of either the MW or M31. After applying the Z12 correction and considering satellites that are likely to be destroyed by baryonic physics or remain dark, only three satellites with v_1 kpc > 20 km s⁻¹ remain. These three satellites are all more luminous than Fornax, the MW's brightest dSph, which we discuss further in the next section.

It can be seen from Figure 3 that applying the Z12 correction overcorrects a small number of subhalos, resulting in negative velocities. These are a population of halos that have lost so much mass after infall that they have very low v_1 kpc values at z = 0, and result in an overcorrection. While Figure 3 already identifies these as halos that have lost at least 97% of their mass since infall, we verified that they actually lost more than 99.9% of their mass after infall, and that the majority have tidal radii less than 1 kpc, and thus can safely be associated with destroyed subhalos.

Two conclusions can be drawn from Figure 3. The first conclusion is that a correction such as that suggested by Z12 is required to bring the masses (and hence velocities) of the predicted subhalo population in line with observational results. Neglecting the baryonic effects that reduce the central masses of subhalos will inevitably lead to a population of subhalos substantially more massive (v_1 kpc > 20 km s⁻¹) than any of the dSphs observed around the MW. The second conclusion is that both tidal destruction in the presence of a baryonic disk and UV heating must be considered in order to bring the total number of luminous satellites in line with observations. All subhalos with M_V brighter than −10 should be bright enough to be detected around our MW. The uncorrected VL2 catalog would suggest that more than 100 detectable subhalos should exist, grossly inconsistent with the classical dSph population (roughly a dozen in the MW and two dozen in M31). Considering destruction mechanisms alone does not bring the observable number into line with observations. An additional correction that assumes suppression of star formation is necessary to further reduce the number of luminous satellites.

Considering the effects of tidal destruction and heating can substantially reduce the number of luminous satellites that are predicted to survive at z = 0 around an MW-mass galaxy. The exact number depends on the assumptions adopted, particularly for the number of dark satellites. The model adopted in this work is not intended to be conclusive, but rather to motivate more rigorous work on this topic. Ideally, future work will adopt a semi-analytic model that follows the growth and merger history of individual halos to determine if they are massive enough to retain baryons and form stars. More work is also needed to understand the influences of the disk on the survival of substructure. Previous work has examined the influence of tidal stripping, though not always with the added presence of a galaxy disk (Taylor & Babul 2001; Stoehr et al. 2002; Hayashi et al. 2003; Kravtsov et al. 2004; Kazantzidis et al. 2004; Read et al. 2006b; Peñarrubia et al. 2008, 2010; Romano-Díaz et al. 2010; Nickerson et al. 2011). In this work, we have emphasized the effect of tidal stripping on satellites in the presence of a host with a baryonic disk (or equivalently, any concentration of centralized baryons, e.g., at high z the central baryons may not be a fully stable disk, but the concentration will still tidally impact the substructure; see Chang et al. 2012). The Z12 correction, however, neglects disk shocking that occurs when a subhalo passes directly through the baryonic disk (Taylor & Babul 2001; D'Onghia et al. 2010a), and will lead to even faster disruption of a subhalo.⁹ Tidal heating of subhalos should also occur (Gnedin et al. 1999; Mayer et al. 2001; D'Onghia et al. 2010b; Kazantzidis et al. 2011), but requires very high resolution to capture (Choi et al. 2009) and is unlikely to be accounted for properly in the Z12 correction. Hence, all of these processes require more thorough study to fully understand the influence of baryons on the evolution of satellites. Finally, one galaxy simulation alone is not sufficient to understand whether baryonic processes can solve the MSP. A statistical sample of MW-mass realizations must be used to quantify the theoretical predictions.

By the same token, one observed galaxy alone is not sufficient to understand whether baryonic processes can solve the MSP. Looking at two galaxies, the MW and M31, it is already clear that the subhalo population of MW-mass galaxies can vary significantly. The surviving luminous VL2 satellites (shown as red circles in Figure 3) include three satellites with luminosities and velocities larger than the dSph population of the MW. However, M31 does have such galaxies, and the surviving VL2 subhalos would more closely resemble the luminous satellite population of M31. M31 has at least four luminous satellites with magnitudes brighter than M_V < −14. Three of these are dwarf ellipticals, a population that is missing in the MW. All three of the VL2 subhalos were accreted at z > 1, have orbital pericenters under 50 kpc, and are unlikely to retain gas, making them look probably very similar to the dwarf ellipticals of M31. We also note that many of the VL2 satellites identified as luminous at z = 0 reach v_peak at relatively low z (see Figure 2), and may have extended star formation to quite low z. Z12 also found that satellites with M_V < −8 had extended star formation histories, consistent with observational estimates (Grebel & Gallagher 2004; Dellenbusch et al. 2008; Weisz et al. 2011). We anticipate that the more recently discovered ultra-faint satellites should have had their star formation truncated at high z (Brown et al. 2012).

In comparing our results for the luminous satellites of VL2 to those of the MW and M31, we should consider the impact of halo host mass on satellite populations. It has been shown that the number of subhalos at a given mass scales with host mass, since the mass assembly of DM halos is self-similar (Stewart et al. 2008; Fakhouri et al. 2010). Several works have argued that the halo mass of M31 is nearly twice as massive as the MW's halo mass (e.g., Kallivayalil et al. 2009; Guo et al. 2010; Watkins et al. 2010), suggesting that M31 should have both a brighter satellite population, and more massive satellites, than the MW. While the exact halo mass of the MW is still unknown, recent estimates (e.g., Smith et al. 2007; Xue et al. 2008; Guo et al. 2010; Gnedin et al. 2010) find a halo mass range of 0.7–2.0 × 10¹² M_☉. The virial mass of VL2 (1.7 × 10¹² M_☉) is near the upper range of the MW's mass estimates, but perhaps much closer to the virial mass of M31 (Watkins et al. 2010), particularly if M31 is twice as massive as the MW. Therefore, the resemblance of the VL2 bright satellite population, i.e., its three satellites with M_V < −14, to the bright satellites of M31 is perhaps indicative that the halo mass of VL2 more closely matches the halo mass of M31 than the MW (see also Wang et al. 2012; Vera-Ciro et al. 2013).

Because this simulated population is a better match to M31 than the MW, we should also consider how common the MW satellite distribution is. A few authors have recently begun exploring this question. On the observational side, searches for satellites in the Sloan Digital Sky Survey around galaxies as luminous as the MW find ∼10 or fewer satellites more luminous than Fornax (Busha et al. 2011; Guo et al. 2011; Lares et al. 2011; Strigari & Wechsler 2012). Adopting a semi-analytic model to assign stellar masses, Vera-Ciro et al. (2013) find that 2–5 subhalos brighter than M_V < −14 consistently exist in the six MW-mass halos in the high-resolution Aquarius DM-only runs. Using a larger statistical sample, both Wang et al. (2012) and Purcell & Zentner (2012) find a 10%–20% probability of finding a subhalo population similar to the MW's at a halo mass of 10¹² M_☉. This suggests that the MW is somewhat rare, though not exceedingly so. The rarity of the MW's subhalo luminosity function seems to be linked to a gap in luminosities between Fornax (at M_V = −13.4) and the Small Magellanic Cloud (at M_V = −16.8), though more work needs to be done to quantify how common the gap is, and if having a pair of galaxies as bright as the Magellanic Clouds has any influence on the existence of such a gap.

Despite the caveats listed above, it is clear that baryons offer a promising solution to solving the MSP, without an additional form of WDM or SIDM. Certainly, heating and destruction have long been considered potential solutions to bring the number of predicted satellites in line with the number of observed satellites (Bullock et al. 2000; Somerville 2002; Kravtsov et al. 2004; Simon & Geha 2007; Koposov et al. 2009; Nickerson et al. 2011; Peñarrubia et al. 2012). The key point, however, is that even if the number of luminous satellites can be brought into line, the overall distribution in masses (and observed velocities) of those satellites is usually too high (Boylan-Kolchin et al. 2011, 2012). The Z12 model builds on this previous work by considering the additional impact of baryons on the DM density slopes of satellites, finding shallow profiles in satellites brighter than M_V < −12, which further reduces their mass. When combined with more effective tides after infall (Taylor & Babul 2001; Stoehr et al. 2002; Read et al. 2006b; Peñarrubia et al. 2010), both the number and masses of the subhalo population can be brought into agreement with observations for the first time.

In fact, the model utilized by Z12 goes beyond solving the MSP or TBTF problems. The model in Z12 has also been shown to create bulgeless dwarf disk galaxies (Governato et al. 2010) and smaller bulges in higher mass galaxies (Brook et al. 2012). This is because winds driven by supernovae preferentially remove low angular momentum material from galaxies (Christensen et al.  2012; Brook et al. 2011). Importantly, these winds are naturally driven (rather than artificially inserted) when supernovae deposit thermal energy into a high density medium, and create highly overpressurized bubbles. Hence, a model must be able to reproduce the high densities found in star-forming molecular clouds (>10–100 atoms/cc) to naturally drive such winds. The ability to model such high density regions is what allowed for DM core creation in the dSph satellites simulated in Z12. Previous studies of dSphs in MW-mass satellites allowed for star formation at lower densities (<1 amu cc⁻¹), and hence found that baryons either steepened the inner densities or had little affect (Sales et al. 2007; Okamoto et al. 2010; Nickerson et al. 2011; Wadepuhl & Springel 2011; Parry et al. 2012; Di Cintio et al. 2012; Sawala et al. 2012). Critically, the expanding overpressurized bubbles are able to drive rapid fluctuations in the potential well of a galaxy. The cumulative effect is to expand the central orbits of DM, transforming an initially steep, cuspy DM density profile into a shallow, cored profile (e.g., Read & Gilmore 2005; Mashchenko et al. 2006; Pontzen & Governato 2012; Macciò et al. 2012; Pontzen & Governato 2012; Teyssier et al. 2012). Governato et al. (2012) and Oh et al. (2011) showed that the density profiles of simulated galaxies within this model are in excellent agreement with the high-resolution observations of field galaxies.

In summary, the model adopted by Z12 not only yields a population of satellite galaxies in agreement with MW or M31 observations, potentially solving both the "missing satellites" and "TBTF" problems, but it also simultaneously addresses other small-scale problems within CDM. Hence, a unified baryonic solution remains viable to solve the small-scale crisis of CDM. Future work must quantify the impact of baryons on the DM structure of galaxies, and make testable predictions across all galaxy masses.

It has long been recognized that the observed population of MW satellites is at odds with the distribution predicted by DM-only simulations. In this article, we have argued that effects associated with baryonic physics can reconcile the results of such simulations with observations, without the need for DM that is warm, self-interacting, or with properties that are otherwise different from those of the standard cold and collisionsless paradigm. In particular, supernova feedback in luminous satellites can reduce the density of DM in the inner volumes of these systems, while the presence of a baryonic disk can enhance the degree of tidal stripping that takes place. The combination of these effects leads to a reduction of the masses of the predicted subhalo population, bringing the overall number of satellites into concordance with that observed. In particular, in the VL2 sample considered, we found that these effects reduced the number of subhalos with masses larger than seen in the MW satellites from more than 20 to only a few.

After determining the distribution of massive DM subhalos predicted to be present in an MW-mass galaxy, we turned our attention on the question of how many of these objects are likely to host luminous satellites. For reasonable assumptions regarding the stellar-to-halo mass relationship and for the criteria for destruction via tidal stripping, we predict a luminous satellites population that is in adequate agreement with both the MW and M31.

While the work presented here is not intended to represent the final word on this topic, we have shown that baryonic effects can lead to a population of satellites around MW-mass galaxies that is in good agreement with observations. This strongly reduces the motivation for WDM or SIDM scenarios and represents yet another success for CDM.

A.B. acknowledges support from The Grainger Foundation. A.Z. acknowledges support from the Lady Davis Foundation. A.Z.'s work was partially supported by the ISF grant 6/08, by GIF grant G-1052-104.7/2009, and by the DFG grant STE1869/1-1.GE625/15-1. This work was supported in part by the U.S. National Science Foundation, grants OIA-1124453 (PI: P. Madau) and OIA-1124403 (PI: A. Szalay).