Not so lumpy after all: modeling the depletion of dark matter subhalos by Milky Way-like galaxies

TL;DR

This study demonstrates that the presence of a central galaxy significantly depletes dark matter subhalos near Milky Way-like galaxies, and a simple embedded disk potential in dark matter-only simulations can replicate this effect.

Contribution

The paper shows that adding an embedded disk potential to DMO simulations accurately models baryonic effects on subhalo populations, offering a fast alternative to full baryonic simulations.

Findings

Baryonic physics drastically reduces subhalo numbers near the galaxy center.

Embedded disk potential reproduces subhalo depletion seen in full baryonic simulations.

Subhalo orbits become more tangentially biased due to tidal destruction near the galaxy.

Abstract

Among the most important goals in cosmology is detecting and quantifying small ($M_{\rm halo}\simeq10^{6-9}~\mathrm{M}_\odot$) dark matter (DM) subhalos. Current probes around the Milky Way (MW) are most sensitive to such substructure within $\sim20$ kpc of the halo center, where the galaxy contributes significantly to the potential. We explore the effects of baryons on subhalo populations in $\Lambda$CDM using cosmological zoom-in baryonic simulations of MW-mass halos from the Latte simulation suite, part of the Feedback In Realistic Environments (FIRE) project. Specifically, we compare simulations of the same two halos run using (1) DM-only (DMO), (2) full baryonic physics, and (3) DM with an embedded disk potential grown to match the FIRE simulation. Relative to baryonic simulations, DMO simulations contain $\sim2\times$ as many subhalos within 100 kpc of the halo center; this excess is $\gtrsim5\times$ within 25 kpc. At $z=0$, the baryonic simulations are completely devoid of subhalos down to $3\times10^6~\mathrm{M}_\odot$ within $15$ kpc of the MW-mass galaxy, and fewer than 20 surviving subhalos have orbital pericenters <20 kpc. Despite the complexities of baryonic physics, the simple addition of an embedded central disk potential to DMO simulations reproduces this subhalo depletion, including trends with radius, remarkably well. Thus, the additional tidal field from the central galaxy is the primary cause of subhalo depletion. Subhalos on radial orbits that pass close to the central galaxy are preferentially destroyed, causing the surviving subhalo population to have tangentially biased orbits compared to DMO predictions. Our method of embedding a disk potential in DMO simulations provides a fast and accurate alternative to full baryonic simulations, thus enabling suites of cosmological simulations that can provide accurate and statistical predictions of substructure populations.

Tables

Table 1. Table 1: Properties of the two host halos, m12i   (first presented in Wetzel et al., 2016 ) and m12f   (first presented here and to be presented in more detail in Wetzel et al., in preparation). ‘ -dmo ’ refers to a dark matter-only (DMO) simulation of the same halo, while ‘ -disk ’ indicates a dark matter + embedded disk potential simulation; the entries without an appendix represent the FIRE baryonic simulations. The ‘2 M d subscript 𝑀 d M_{\mathrm{d}} ,’ ‘ M ⋆ + M gas subscript 𝑀 ⋆ subscript 𝑀 gas M_{\star}+M_{\mathrm{gas}} ,’ ‘2 R d subscript 𝑅 d R_{\mathrm{d}} ,’ ‘ z d ​ 0 subscript 𝑧 d 0 z_{\mathrm{d}}0 ,’ and ‘ Hernquist ’ simulations below the break indicate runs discussed in § 3.4 , which include factors of ≲ 2 less-than-or-similar-to absent 2 \lesssim 2 changes to the disk parameters, or, in the final case, a spherically symmetric potential with the same mass. The columns indicate the name of the simulation, the mass of the host halo (using the virial definition of Bryan & Norman, 1998 ) , the maximum circular velocity, V max subscript 𝑉 max V_{\mathrm{max}} , of the host halo (including only DM), the virial velocity, V v subscript 𝑉 v V_{\mathrm{v}} , of the host halo (defined as G ​ M v / R v 𝐺 subscript 𝑀 v subscript 𝑅 v \sqrt{GM_{\mathrm{v}}/R_{\mathrm{v}}} , where R v subscript 𝑅 v R_{\mathrm{v}} is the virial radius), the number of subhalos within 25 25 25 , 50 50 50 , 100 100 100 , and 300 ​ kpc 300 kpc 300~{}\mathrm{kpc} at z = 0 𝑧 0 z=0 with V max > 5 ​ km ​ s − 1 subscript 𝑉 max 5 km superscript s 1 V_{\mathrm{max}}>5~{}{\rm km}\ {\rm s}^{-1} (corresponding roughly to M bound ≃ 5 × 10 6 ​ M ⊙ similar-to-or-equals subscript 𝑀 bound 5 superscript 10 6 subscript M direct-product M_{\rm bound}\simeq 5\times 10^{6}\,\textnormal{M}_{\odot} ), and f bound subscript 𝑓 bound f_{\mathrm{bound}} , the total mass in resolved subhalos divided by the virial mass of the host. M v subscript 𝑀 v M_{\mathrm{v}} and V v subscript 𝑉 v V_{\mathrm{v}} for the baryonic simulations include the contributions of gas and stars. However, we list V max subscript 𝑉 max V_{\mathrm{max}} based only on DM for these runs because the circular velocity curve peaks at small radii ( r ≃ 1.6 ​ kpc similar-to-or-equals 𝑟 1.6 kpc r\simeq 1.6~{}\mathrm{kpc} ) when including baryons in the mass profile; in this case, m12i  reaches V max = 279 ​ km ​ s − 1 subscript 𝑉 max 279 km superscript s 1 V_{\mathrm{max}}=279~{}{\rm km}\ {\rm s}^{-1} while m12f  reaches V max = 283 ​ km ​ s − 1 subscript 𝑉 max 283 km superscript s 1 V_{\mathrm{max}}=283~{}{\rm km}\ {\rm s}^{-1} .

Simulation	$M_v$	$V_{\max }$	$V_v$	$N_{\mathrm{sub}}$	$N_{\mathrm{sub}}$	$N_{\mathrm{sub}}$	$N_{\mathrm{sub}}$	$f_{\mathrm{bound}}$
	(${10}^{12}{\text{M}}_{\odot }$)	($\mathrm{km}{s}^{-1}$)	($\mathrm{km}{s}^{-1}$)
m12f	1.6	183	149	4	62	266	1482	0.037
m12f-dmo	1.6	177	150	42	204	654	2423	0.081
m12f-disk	1.7	195	153	13	97	379	1854	0.049
m12i	1.1	163	134	9	57	370	1432	0.056
m12i-dmo	1.1	162	134	39	200	671	2069	0.078
m12i-disk	1.2	191	136	6	108	473	1712	0.062
m12i-disk-2$M_d$	1.2	220	138	5	96	455	1619	0.051
m12i-disk-$M_{\star }+M_{\mathrm{gas}}$	1.2	188	136	9	45	352	2142	0.062
m12i-disk-2$R_d$	1.2	189	136	16	129	485	1738	0.062
m12i-disk-$z_d$0	1.2	183	136	19	120	529	1733	0.061
m12i-Hernquist	1.2	216	136	10	70	385	2214	0.066

Table 2. Table 2: Properties of the analytic ‘embedded disk’ of the MW-mass central galaxies at several redshifts, along with the properties from the FIRE baryonic simulations that we model the disks on. The first two columns list the total mass in the embedded disk, M d subscript 𝑀 d M_{\mathrm{d}} , and the equivalent stellar mass in the FIRE galaxy; the following two columns list the radial extent of the disk (see text for details). Our analytic disks are typically slightly more compact, particularly at late times, both because the Miyamoto-Nagai potential is an imperfect fit to the galaxies that form in the baryonic simulations, and because we fit the density profile, rather than the mass or the potential. As we show in § 3.4 , the amount of substructure depletion most strongly depends on the disk mass; doubling R d subscript 𝑅 d R_{\mathrm{d}} only weakly affects the surviving subhalo population.

	redshift	$M_d$	$M_{\mathrm{gal}}$	$R_{90,\mathrm{disk}}$	$R_{90,\mathrm{gal}}$
		${10}^{10}{M}_{\odot }$	${10}^{10}{M}_{\odot }$	$\mathrm{kpc}$	$\mathrm{kpc}$
m12i	$z=3$	0.09	0.12	15.75	14.86
	$z=2$	0.32	0.29	15.84	12.57
	$z=1$	1.59	1.44	7.94	10.76
	$z=0.5$	4.28	3.91	4.22	5.93
	$z=0$	8.56	7.70	6.43	9.49
m12f	$z=3$	0.15	0.16	9.33	7.70
	$z=2$	0.71	0.73	6.84	7.43
	$z=1$	2.40	2.18	4.11	7.05
	$z=0.5$	4.86	4.41	4.31	9.01
	$z=0$	10.41	9.39	7.22	10.48

Equations

$$\Phi (R,z)=\frac{G{M}_{\mathrm{d}}}{{\left\{\Delta R^2+{\left[R_{\mathrm{d}}+{\left(\Delta z^2+z_{\mathrm{d}}^2\right)}^{1/2}\right]}^2\right\}}^{1/2}},$$

$$a_x=\frac{-G{M}_{\mathrm{d}}}{{\left\{\Delta R^2+{\left[R_{\mathrm{d}}+{\left(\Delta z^2+z_{\mathrm{d}}^2\right)}^{1/2}\right]}^2\right\}}^{3/2}}\Delta x,$$

$$a_z=\frac{-G{M}_{\mathrm{d}}\left[R_{\mathrm{d}}+{\left(\Delta z^2+z_{\mathrm{d}}^2\right)}^{1/2}\right]}{{\left\{\Delta R^2+{\left[R_{\mathrm{d}}+{\left(\Delta z^2+z_{\mathrm{d}}^2\right)}^{1/2}\right]}^2\right\}}^{3/2}{\left(\Delta z^2+z_{\mathrm{d}}^2\right)}^{1/2}}\Delta z.$$

Full text

††pubyear: 2017††pagerange: Not so lumpy after all: modeling the depletion of dark matter subhalos by Milky Way-like galaxies–References

1 Introduction

One of the strongest predictions of $\Lambda$CDM (cosmological constant with cold dark matter) is that dark matter clusters hierarchically: large halos that host Milky Way (MW)-size galaxies are filled with smaller, self-bound clumps known as subhalos. The highest-resolution cosmological simulations of MW-size halos in the $\Lambda$CDM paradigm have demonstrated that dark matter (DM) clumps exist at all resolved masses (e.g. Springel et al., 2008; Kuhlen et al., 2009; Stadel et al., 2009; Garrison-Kimmel et al., 2014a; Griffen et al., 2016).

While at least some of these subhalos are presumed to host faint satellite galaxies, the ‘missing satellites’ problem (Klypin et al., 1999; Moore et al., 1999) points out a sharp discrepancy between the flat luminosity function of observed satellites and the steep, ever-rising mass function of subhalos predicted from numerical simulations. Though the discrepancy can be largely eliminated by invoking gas heating from reionization suppressing star formation in the early Universe (Bullock et al., 2000; Somerville, 2002) and observational incompleteness (Tollerud et al., 2008), this solution demands that subhalos with bound masses smaller than $\sim 10^{8}\,\textnormal{M}_{\odot}$ should be dark and that thousands of $\sim 10^{6}\,\textnormal{M}_{\odot}$ subhalos should be entirely devoid of stars.

Confirming the existence of these tiny, dark subhalos, and further determining their mass function, would simultaneously provide an astounding confirmation of the $\Lambda$CDM theory and rule out large classes of warm dark matter models and inflationary models that predict a cut-off in the power spectrum at low masses (Kamionkowski & Liddle, 2000; Bode et al., 2001; Zentner & Bullock, 2003; Horiuchi et al., 2016; Bozek et al., 2016; Bose et al., 2016). Because these subhalos are dark, however, they must be identified indirectly. Around the MW, the best possibilities for detecting dark substructure are through perturbations to the stellar disk (e.g. Quinn et al., 1993; Feldmann & Spolyar, 2015, and references therein) and via gaps or kinematic distortions in dynamically cold stellar streams formed from disrupting globular clusters, such as Palomar-5 and GD-1 (e.g. Johnston et al., 2002; Koposov et al., 2010; Carlberg et al., 2012; Ngan et al., 2015, and references therein). In fact, Bovy et al. (2017) recently claimed a measurement of $10^{+11}_{-6}$ subhalos of mass $3\times 10^{6}$–$10^{9}\,\textnormal{M}_{\odot}$ within $20~{}\mathrm{kpc}$ of the MW via the Palomar-5 stream. Amorisco et al. (2016), however, argued that giant molecular clouds could also be responsible for some of the fluctuations, and Ibata et al. (2016) reported a null detection using the same stream. Currently undetected stellar streams may provide more information about dark substructure around the MW (Ngan et al., 2016). Around larger, more distant galaxies, dark substructures may be revealed by gravitational lensing anomalies from background sources (Mao & Schneider, 1998; Dalal & Kochanek, 2002; Vegetti et al., 2010; MacLeod et al., 2013; Nierenberg et al., 2014; Hezaveh et al., 2016), particularly with the upcoming instruments on the James Webb Space Telescope (JWST; MacLeod et al., 2013). These lensing studies are sensitive to substructure within a projected Einstein radius that is typically $\sim 5$–$10~{}\mathrm{kpc}$ in size (Dalal & Kochanek, 2002; Fiacconi et al., 2016).

Making predictions for these observations, and thus using them to constrain the properties of DM on small scales, requires a statistical sample of halos simulated in $\Lambda$CDM with sufficient resolution (particle mass $\lesssim 10^{5}\,\textnormal{M}_{\odot}$) to identify the tiny subhalos of interest. While several such simulations exist in the literature (e.g. Springel et al., 2008; Kuhlen et al., 2009; Stadel et al., 2009; Mao et al., 2015; Griffen et al., 2016), the vast majority are purely collisionless (DM-only; DMO), given their low computational cost compared with fully baryonic simulations. However, DMO simulations are problematic because they miss critical baryonic physics, including the simple presence of a central galaxy in the halo, near the very regions where observational probes of substructure are most sensitive. While the smallest subhalos are likely devoid of baryons, they nonetheless should be dynamically influenced by the central galaxy (Taylor & Babul, 2001; Hayashi et al., 2003; Berezinsky et al., 2006; Read et al., 2006a, b; Peñarrubia et al., 2010; D’Onghia et al., 2010). In fact, the mere existence of a stellar halo around the MW implies a significant population of destroyed dwarf galaxies.

Furthermore, observational estimates for the masses of larger subhalos, which host luminous dwarf galaxies and therefore can be identified through direct observations, are well below expectations from DMO simulations of MW-like halos. The magnitude of this discrepancy, which is known as the ‘too-big-to-fail’ problem (Boylan-Kolchin et al., 2011, 2012; see also Tollerud et al., 2014 for similar results for M31 and Garrison-Kimmel et al., 2014b; Papastergis et al., 2015 for field dwarf galaxies), depends on the adopted mass of the MW (di Cintio et al., 2011; Wang et al., 2012; Vera-Ciro et al., 2013), but internal baryonic processes (such as rapid gaseous outflows driven by bursty stellar feedback that can turn the cuspy DM profiles predicted by collisionless simulations into cored profiles with central masses at $r\lesssim 500~{}\mathrm{pc}$ that are consistent with the observations) appear important for resolving TBTF (e.g. Navarro et al., 1996; Read & Gilmore, 2005; Mashchenko et al., 2008; Pontzen & Governato, 2012; Governato et al., 2012; Amorisco et al., 2014; Gritschneder & Lin, 2013; Oñorbe et al., 2015; Chan et al., 2015). As Brooks & Zolotov (2014) showed, the increased tidal forces from the central galaxy also can reduce the central masses of such dwarf galaxies even further, potentially eliminating the disagreement between theory and observations entirely (also see Read et al., 2006a, b). Additionally, the central galaxy potential may completely destroy some of the large, dense subhalos, effectively lowering theoretical predictions for the central masses of the halos expected to host the luminous dwarfs by placing them in correspondingly lower mass, and thus more abundant, subhalos (also see Garrison-Kimmel et al., 2017). However, the amount of destruction that can be unambiguously attributed to the central galaxy remains uncertain.

Recently, the ‘Latte’ simulation (Wetzel et al., 2016), part of the Feedback In Realistic Environments (FIRE) project, achieved baryonic mass resolution of $\sim 7000\,\textnormal{M}_{\odot}$ for a MW-mass halo run to $z=0$. The Latte simulation resolves subhalos down to mass $\sim 10^{6}\,\textnormal{M}_{\odot}$, and as Wetzel et al. (2016) showed, the initial Latte simulation produces a population of satellite dwarf galaxies that agrees with a wide variety of observations around the MW and M31: the distributions of stellar masses, velocity dispersions, and star formation histories, and the relationship between stellar mass and metallicity all agree well those of the MW satellites. Thus, it does not suffer from either the missing satellites or too-big-to-fail problems, at least for resolved satellite dwarf galaxies ($M_{\star}\gtrsim 10^{5}\,\textnormal{M}_{\odot}$). Using the APOSTLE simulations of Local Group-like MW-M31 pairs, Sawala et al. (2017) also recently found good observational agreement for satellite stellar masses and subhalo circular velocities using different treatments of baryonic physics. While these results from baryonic simulations are quite promising, such simulations are sufficiently expensive to prohibit large parameter surveys and statistical samples.

In this work, we extend the initial analysis of Wetzel et al. (2016) to study substructure populations down to the smallest mass scales of relevance for current dark substructure searches ($\sim 10^{6}\,\textnormal{M}_{\odot}$). We also present a second simulation of a MW-mass galaxy in the Latte suite. We will show that properly accounting for the effects of baryons is essential for accurately predicting subhalo populations, even for completely dark subhalos that have no stars. Motivated by previous work that suggested that adding an analytic potential or other inexpensive modification(s) to DMO simulations could capture the key baryonic effects on dark subhalos (e.g. D’Onghia et al., 2010; Zhu et al., 2016; Sawala et al., 2017; Errani et al., 2017; Jethwa et al., 2016), we apply a method for inserting an approximate analytic description of the gravitational potential of the galaxy that forms at the center of each MW-mass halo into a cosmological DMO simulation.

Our approach is particularly interesting for testing the underlying physical drivers, because we both calibrate our input central disk model and benchmark our results against our fully baryonic simulations, which reproduce many observable properties of the satellite populations around the MW and M31. As we will show, many key differences in subhalos populations between DMO and baryonic simulations can be unambiguously attributed to the presence of the central galaxy potential. This success of embedding a central galaxy potential also demonstrates that the substructure populations predicted by cosmological DMO simulations can be significantly improved (i.e. brought into better agreement with fully baryonic simulations) at minimal CPU cost.

In § 2, we discuss the simulations and detail our method of inserting an embedded potential into the center of the host; §3 explores subhalo population statistics with and without a forming galaxy and presents trends with radius in subhalo depletion. We discuss further implications of our results in §4 and conclude in §5.

Throughout this work, we use $h=0.702$, $\Omega_{\mathrm{m}}=0.272$, $\Omega_{\mathrm{b}}=0.0455$, and $\Omega_{\Lambda}=0.728$.

2 Simulations

All of our simulations are cosmological and employ the ‘zoom-in’ technique (Katz & White, 1993; Oñorbe et al., 2014). We run all of our simulations in the same cosmological volume as the AGORA project (Kim et al., 2014), with box length of $60h^{-1}\,\mathrm{Mpc}=85.5~{}\mathrm{Mpc}$. We choose each high-resolution region to contain a single MW-mass ($M_{\mathrm{halo}}\sim 10^{12}\,\textnormal{M}_{\odot}$) halo at $z=0$ that has no neighboring halos of similar or greater mass within $3\,\textnormal{Mpc}$. We focus on two such halos, designated as m12i and m12f, which are part of the Latte sample from the FIRE project. m12i was presented in Wetzel et al. (2016); m12f was simulated with identical parameters using the same pipeline and will be described in detail in Wetzel et al. (in preparation). We chose both halos based only on their virial mass (see Table 1) and not based on their formation/merger history or subhalo population.

We ran all simulations using GIZMO (Hopkins, 2015)111http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html, which uses an updated version of the TREE+PM gravity solver included in GADGET-3 (Springel, 2005). We identify halo centers and create halo catalogs with AHF (Knollmann & Knebe, 2009) and build merger trees usingconsistent-trees (Behroozi et al., 2013a). We generated initial conditions for the DMO and baryonic simulations at $z=100$ using MUSIC (Hahn & Abel, 2011) with second-order Lagrangian perturbation theory. For the embedded disk simulations, we generated initial conditions using the snapshot at $z=3$ from the DMO simulations.

The properties of the two host halos, along with the number of resolved subhalos identified by AHF within several radial cuts of that host, are listed in Table 1. The first column lists the names of the simulations: those identified with ‘dmo’ are purely collisionless, while the simulations labeled as ‘disk’ are DMO with an embedded galactic potential; the other rows show the FIRE baryonic simulations. The remaining columns list the Bryan & Norman (1998) virial masses, maximum circular velocities, $V_{\mathrm{max}}$, and virial velocities ($\sqrt{GM_{\mathrm{v}}/R_{\mathrm{v}}}$) of the hosts, here calculated from the full particle distribution (including gas and stars in the baryonic simulations), the total number of subhalos with $V_{\mathrm{max}}>5{\rm km}\ {\rm s}^{-1}$ that survive to $z=0$ within $25$, $50$, $100$, and $300~{}\mathrm{kpc}$ of the halo center, and the fraction of the host mass that resides in self-bound subhalos.

2.1 Baryonic simulations

The baryonic simulation of m12i analyzed here is the same run presented in Wetzel et al. (2016); m12f was simulated with identical code and parameters. The physics and numerical prescriptions are therefore given in Wetzel et al. (2016). Briefly, however, the baryonic simulations are part of the Feedback In Realistic Environments (FIRE; Hopkins et al., 2014) project222http://fire.northwestern.edu. Specifically, they use the updated FIRE-2 code, which features identical physics as FIRE-1 but incorporates several numerical improvements. In particular, FIRE-2 adopts the new mesh-free finite mass (‘MFM’) method for more accurate hydrodynamics (Hopkins, 2015). We model radiative heating and cooling from $10-10^{10}\,$K (following CLOUDY tabulations; Ferland et al., 1998), and accounting for self-shielding and photo-heating both by a UV background (from Faucher-Giguère et al., 2009) and by local sources. Star formation occurs in self-gravitating gas (according to the criterion in Hopkins et al., 2013) that is also molecular and self-shielding (following Krumholz & Gnedin, 2011), Jeans unstable, and exceeds a minimum density threshold $n_{\rm sf}>1000\,{\rm cm^{-3}}$. A star particle is then spawned probabilistically from a gas particle, inheriting the same mass and metallicity. The simulations follow several stellar feedback mechanisms, including (1) local and long-range momentum flux from radiation pressure (in the initial UV/optical single-scattering, and from re-radiated light in the IR), (2) energy, momentum, mass, and metal injection from core-collapse and Ia supernovae, as well as stellar mass loss from OB and AGB stars, and (3) photo-ionization and photo-electric heating. Every star particle is treated as a single stellar population with a mass, age, and metallicity. We tabulate all feedback event rates, luminosities and energies, mass-loss rates, and other quantities directly from stellar evolution models (STARBURST99 v7.0; Leitherer et al., 1999) assuming a Kroupa (2001) initial mass function (IMF).

Full details of FIRE-2 are provided in Hopkins et al. (2017). The source code and numerical parameters of our baryonic simulations are exactly identical to those in all FIRE-2 simulations (Wetzel et al., 2016; Su et al., 2016; Fitts et al., 2016).

The FIRE simulations have been shown to reproduce a wide variety of observables, including the relationships between stellar mass and halo mass, the Kennicutt-Schmidt law, bursty star formation histories, the star forming main sequence (Hopkins et al., 2014), galactic winds (Muratov et al., 2015, 2017), the gas and stellar phase $M_{\star}$-metallicity relations (Ma et al., 2016), the $M_{\star}$-size relation (El-Badry et al., 2016), the HI content of galaxy halos at both low and high redshift (Faucher-Giguère et al., 2015, 2016; Hafen et al., 2017), and the structure and star formation histories of isolated dwarf galaxies (Oñorbe et al., 2015; Chan et al., 2015; Fitts et al., 2016). Moreover, in simulations of MW-mass halos, in addition to forming a realistic MW-like galaxy in terms of stellar mass and disk morphology (Wetzel et al., 2016; Ma et al., 2017), the FIRE model yields reasonable populations of dwarf galaxies around those galaxies, in terms of the distributions of stellar masses and velocity dispersions, as well as a wide range of star formation histories that agree well with those of the actual MW satellites.

Both m12i and m12f form thin, radially extended stellar disks with $M_{\star}(R<R_{90},z<z_{90})=6.2\times 10^{10}\,\textnormal{M}_{\odot}$ and $7.5\times 10^{10}\,\textnormal{M}_{\odot}$, respectively, where $R_{90}$ and $z_{90}$ are the radius and height that contain 90% of the mass. Thus, these galaxies are comparable to, if slightly more massive than, the MW in stars (Bland-Hawthorn & Gerhard, 2016). At $z=0$, the total gas fraction, $M_{\rm gas}/(M_{\star}+M_{\rm gas})$, within $R_{90}$ and $z_{90}$ is 13% for m12i and 15% for m12f.

The gravitational force softenings and kernel smoothing lengths for gas particles are fully adaptive and conservative (following Price & Monaghan, 2007). Hydrodynamic smoothings and gravitational force softenings are always self-consistently matched. The minimum gas smoothing/Plummer equivalent softening achieved in both simulations is $\epsilon_{\mathrm{gas,\,min}}=1~{}\mathrm{pc}$ (corresponding to a density of $n_{\rm gas}\approx\!10^{7}\,\mathrm{cm}^{-3}$), thus ensuring that dense, star-forming regions are well resolved. We choose softenings for the DM particles to be comparable to the typical gas softening in the host galaxy’s disk: $\epsilon_{\rm DM}=20~{}\mathrm{pc}$. The softenings for the stars are $\epsilon_{\rm stars}=8~{}\mathrm{pc}$, chosen to match the gas softening at the density threshold for star-forming regions, $n_{\rm sf}>1000~{}\mathrm{cm}^{-3}$. All (minimum) softening lengths quoted here are fixed in physical units after $z=9$, and evolve comoving with the scale factor prior to that redshift. Each simulations is initialized at an ‘effective’ resolution of $2\times 8192^{3}$ particles within the box, resulting in a dark matter particle mass of $m_{\rm DM}=3.5\times 10^{4}\,\textnormal{M}_{\odot}$ and a (initial) gas or star particle mass of $m_{\rm gas}=7.1\times 10^{3}\,\textnormal{M}_{\odot}$.

2.2 Dark matter-only simulations

The dark matter-only (DMO) simulations are identical to the baryonic simulations, except that they include only dark matter particles, and the baryonic mass is included in the dark matter particles. Consequently, the particle masses are larger by a factor of $1/(1-f_{\mathrm{b}})$, where $f_{\mathrm{b}}=\Omega_{b}/\Omega_{\mathrm{m}}\simeq 0.17$ is the cosmic baryon fraction: $m_{\rm DM,\,DMO}=4.2\times 10^{4}\,\textnormal{M}_{\odot}$. This increased particle mass has non-trivial effects on comparisons of DMO and baryonic simulations, both directly, as in the case of subhalo mass functions, and indirectly, through implied resolution cuts. This difference is particularly relevant for low-mass subhalos, which have lost essentially all of their baryonic mass by $z=0$ from cosmic reionization, feedback-driven gas heating, and ram-pressure stripping within the host halo. Thus, a given subhalo may contain an identical number of DM particles in the collisionless and baryonic simulations, but be more massive in the former. Therefore, for low-mass subhalos, which are expected to have lost (nearly) all of their baryonic mass, the most physical way to quote their masses in DMO simulations is to correct for this (presumed) baryonic mass loss. We thus multiply all particle and subhalo masses by $1-f_{\mathrm{b}}$ in all post-processing (after halo finding), and we similarly suppress the maximum circular velocities of all subhalos $V_{\mathrm{max}}$ by $\sqrt{1-f_{\mathrm{b}}}$ (similar to Zolotov et al., 2012). Finally, we perform halo finding only on the DM particles in the baryonic simulations, to achieve the fairest comparison. Therefore, the differences between baryonic and DMO simulations that we quote here are somewhat smaller than quoted in Wetzel et al. (2016), who did not include this correction.

2.3 Embedded disk potentials

To include the effects of the disk of the central galaxy – which grows naturally within baryonic simulations – in our DMO simulations, we add an additional gravitational acceleration to every particle active during each timestep, as given by a Miyamoto & Nagai (1975) potential:

[TABLE]

where $\Delta R$ and $\Delta z$ indicate the relative position from the center of the potential ($x_{\rm disk},y_{\rm disk},z_{\rm disk}$) in cylindrical coordinates: $\Delta R^{2}=(x-x_{\mathrm{disk}})^{2}+(y-y_{\mathrm{disk}})^{2}$; $\Delta z^{2}=(z-z_{\mathrm{disk}})^{2}$. This potential has three parameters: $R_{\mathrm{d}}(z)$, the disk scale length, $z_{\mathrm{d}}(z)$, the disk scale height, and $M_{\mathrm{d}}(z)$, the total mass in the potential. Importantly, the acceleration from the Miyamoto-Nagai disk is analytic: in the plane of the disk, taken here to be the $x-y$ plane, it is

[TABLE]

where $\Delta x=(x-x_{\mathrm{disk}})$, and similarly in the $y$ direction. The acceleration along the minor axis, taken to be the $z$ direction, is

[TABLE]

Within the DMO simulation, we track the center of the MW-mass host halo across cosmic time using a single massive ($m_{\mathrm{p}}=10^{8}\,\textnormal{M}_{\odot}/h$) particle inserted at the center of the main branch of the host at $z=3$, using a large softening length ($\epsilon=3.7$ physical $\mathrm{kpc}$). This effectively acts as a small bulge in the center of the galaxy. Because this particle is significantly more massive than the high-resolution particles that comprise the remainder of the halo, dynamical friction acts to keep it near the center of the host.333At $z=0$, the particle is within $12~{}\mathrm{pc}$ of the center of m12i found by AHF and within $235~{}\mathrm{pc}$ of m12f, even without including the particle in the halo finding. The position of this particle at each timestep determines the center of the disk potential. To minimize computational cost, $(x,y,z)$ offsets from this particle are computed for all other particles while traversing the gravity tree.

In our embedded disk simulations, we allow this disk potential to evolve over time to match the stellar disk that forms in the corresponding baryonic simulation of the same system. We first embed the disk potential at $z=3$ ($\sim 11.5~{}\mathrm{Gyr}$ ago), initializing with parameters ($M_{\mathrm{d}}$, $R_{\mathrm{d}}$, and $z_{\mathrm{d}}$) obtained by jointly fitting the average density profile of each baryonic simulation’s stellar disk along the major and minor axes, defined by solving for the eigenvectors of the moment of inertia tensor of all stars within 20 kpc. We then linearly interpolate each parameter in scale factor between fits performed at twelve additional snapshots, equally spaced in $z$, to $z=0$. In all fits, we bound the total mass in the potential to be less than $0.1$ dex greater than the stellar mass of the disk in the baryonic simulation. This method yields an excellent match between the stellar mass in the baryonic simulations and the mass in the embedded disks within fixed physical radii.444At the thirteen steps that anchor the time evolution, the disk potential matches the stellar mass within 20 kpc to within $10\%$ in m12f and within $20\%$ in m12i at worst; the average agreement over those thirteen steps is $\sim 7\%$ across both simulations.

Table 2 lists the physical parameters of the analytic disk potential, and equivalent quantities for the galaxies that form in the FIRE simulations. The first two columns quote the mass, quantified as $M_{\mathrm{d}}$ for the disk potential (the total mass when integrated to infinity) and by $M(<r_{90})/0.9$ for the FIRE simulations, where $r_{90}$ is the 3D radius that contains 90% of the mass. The last two columns give the radial extent of the disk, quantified here by the 2D radius $R_{90}$, where $M(R<R_{90})=0.9\times M(z<z_{90})$ and $z_{90}$ similarly contains 90% of the total stellar mass of the galaxy.

While we carefully design the analytic potential to match the parameters of the stellar disk in the baryonic simulations, we note several limitations. For simplicity, we fix the orientation of the disk to be along the $x-y$ axis plane (rather than allowing it to rotate through arbitrary angles). Additionally, in the baryonic simulations, m12i and m12f do not form permanent well-ordered disks until $z\approx 0.5$ and $z\approx 0.6$, respectively, so our assumed potential overestimates the thinness of the disk at early times. However, as we show below, the ‘thinness’ of the disk is not important for capturing its tidal effects on subhalos – in fact, replacing the disk with a spherical Hernquist (1990) profile of the same mass produces nearly the same effect (see §3.4). Furthermore, we add this disk potential to the DMO simulation without adjusting the mass of dark matter particles, thus slightly increasing the overall mass within the host halo. However, this error is likely small to the overall system because the disks comprise only $\sim 6\%$ of the total mass within $R_{\mathrm{v}}$ and $10\%$ within 80 kpc.

We also fit our fiducial disks to only the stellar mass from the baryonic simulation; we do not try to fit the gaseous component. This is a reasonable approximation at $z=0$, where the gas fractions in the disks of our simulated host galaxies are $\approx 15\%$, comparable to the MW, M31, and similar galaxies (Yin et al., 2009; Catinella et al., 2010, 2012, 2013). At higher redshifts, however, when the disks were more gas-rich, our method underestimates the total baryonic mass within the disk. At these redshifts, stellar feedback drives significant gas flows in and out of galaxy on short timescales ($\lesssim 100$ Myr), compromising the accuracy of any simple, constant analytic description of the resultant potential. In principle, this lack of incorporating gas implies that our stellar-disk-only model represents a lower-limit to the level of substructure depletion from the central galaxy, though, as we show in §3.4 by increasing the mass of the embedded disk to include the gas in the galaxy, the purely gravitational contribution from gas is small.

However, as we demonstrate below, even this simple model is remarkably successful at reproducing the statistical properties of surviving subhalos as compared with the fully baryonic simulations. Thus, while there is room for further progress, our method represents a substantial improvement over DMO simulations at essentially the same CPU cost.

3 Results

Figure 1 shows images of dark matter in our three types of simulations, highlighting the qualitative effect of the central galaxy on the dark matter distribution relative to the DMO simulation. The top panels visualize, from left to right, the dark matter density within a cube 500 kpc on a side in the DMO, FIRE baryonic, and embedded disk simulations of m12i; the lower panels show a zoomed-in view of a cube 100 kpc across. The log-scaled colormap changes from the upper to lower panels, but is identical across the simulations. The dark matter near the central galaxies (lower panels) clearly responds to the extra mass present in both the baryonic and the embedded disk simulations. However, the dark matter takes on a slightly more disk-like shape around the analytic embedded potential; we posit that the relatively spherical cusp around the central galaxy in the baryonic simulation arises because of the continually evolving orientation of the central stellar disk, along with the fluctuating gas distribution, which is more spherical in a time-average sense because of time-dependent gas inflows and outflows.

Even more striking than the enhancement in the central dark matter density, however, is the severe reduction in the number of subhalos within the central $r\lesssim 50~{}\mathrm{kpc}$. The central galaxy has destroyed an enormous fraction of the satellites that the DMO simulation predicts in the central regions, where observational probes are most sensitive.

3.1 Dependence on mass and radius at $z=0$

Figure 2 quantifies the differences in the subhalo populations. The top row shows the cumulative counts of subhalos within $100~{}\mathrm{kpc}$ of m12i (left) and m12f (right) as a function of $V_{\mathrm{max}}$. The lower sub-panels show the ratio between the cumulative counts in the collisionless simulations (with and without the embedded potential) and the FIRE baryonic simulations. The DMO simulations, plotted in thin dashed lines, overpredict the subhalo $V_{\mathrm{max}}$ function relative to the baryonic simulations, plotted as solid lines, at all $V_{\mathrm{max}}$. The light dashed lines, which plot counts in the DMO simulations without correcting for the difference in particle masses between DMO and baryonic simulations, show that DMO predictions improve relative to the FIRE simulations after correcting for this effect, but only slightly – the disk reduces counts at fixed $V_{\mathrm{max}}$ by much more, particularly in m12f. In fact, the embedded disk simulations provide a significantly better match to the FIRE baryonic simulations at all $V_{\mathrm{max}}$. As demonstrated in the bottom panels, the disk simulations agree with the baryonic simulations to within $\sim 20\%$ at nearly all $V_{\mathrm{max}}$.

The bottom row in Figure 2 illustrates the radial distribution of the subhalos included in the top row by showing the cumulative number of subhalos as a function of 3D distance from the halo center, $r$. As expected, the depletion is greatest at the center of the host: both the m12i and m12f baryonic simulations have no subhalos within $15~{}\mathrm{kpc}$ and only $1$–$2$ within $20~{}\mathrm{kpc}$ at $z=0$, while the DMO simulations predict $\sim 7$ subhalos within $15~{}\mathrm{kpc}$ and $\gtrsim 10$ within $20~{}\mathrm{kpc}$ (also see Table 1). The embedded potential captures this effect well, particularly around m12f: only two subhalos remain within $20~{}\mathrm{kpc}$. The match is slightly worse around m12i, where three subhalos remain within $20~{}\mathrm{kpc}$ and two within $15~{}\mathrm{kpc}$, but the embedded disk simulation still improves upon the DMO simulation by a factor of $\sim 3$. Though we do not explicitly plot it here, we also note that, even after accounting for the increased destruction via the disk, the radial distributions of the surviving subhalos are more extended than predicted from DMO simulations.

This dearth of substructure within $\sim 20~{}\mathrm{kpc}$ has strong implications for attempts to discover dark subhalos via stellar streams around the MW. Thus far, the best studied streams, Palomar 5 at $19~{}\mathrm{kpc}$ (indicated by the grey vertical band; Carlberg et al., 2012) and GD-1 at $\sim 15~{}\mathrm{kpc}$ (plotted as the gold band; Koposov et al., 2010), are within this distance, suggesting that baryonic simulations predict far fewer gaps in stellar streams than DMO simulations. The null result from Ibata et al. (2016) therefore may be in line with predictions from $\Lambda$CDM including baryonic effects. However, those gaps could be created by subhalos at an earlier point in their orbits than their positions at $z=0$; we therefore examine the distribution of pericentric distances $d_{\mathrm{peri}}$ in §3.3.

While $V_{\mathrm{max}}$ can be measured robustly in simulations, the degree to which a subhalo disrupts a stellar stream is more closely related to its bound mass (along with the interaction distance and velocity; e.g. Yoon et al., 2011; Carlberg, 2013; Sanders et al., 2016; Sanderson et al., 2016). Therefore, the left panel of Figure 3 plots the cumulative subhalos count, again within $100~{}\mathrm{kpc}$, as a function of the bound mass assigned to them by AHF, $M$. Because we run AHF only on dark matter particles, $M$ therefore represents the bound dark matter mass. The results are similar to those using $V_{\mathrm{max}}$ as in Figure 2, though the discrepancies between the DMO and baryonic simulations are less severe: the DMO simulations overpredict subhalo counts above fixed $M$ by a factor of $\gtrsim 1.7$.

The right panel of Figure 3 plots the cumulative radial distribution of subhalos with $M>3\times 10^{6}\,\textnormal{M}_{\odot}$, which are reliably resolved in our simulations with 85 particles (Appendix B presents an explicit resolution test). The distributions are again similar to those in the lower plots in Figure 2: simulations with a central galaxy predict only $1$–$3$ subhalos in the central $20~{}\mathrm{kpc}$, while the DMO simulations overpredict that count by a factor of $>5$. Therefore, the lack of substructure near the galaxy is independent of whether subhalos are selected by $V_{\mathrm{max}}$ or $M$.

Because correcting for the reduction in the particle mass from the lack of baryons is more straightforward for $M$ than for $V_{\mathrm{max}}$, which depends on the mass profile of the subhalo, we select subhalos based on $M$ for the remainder of the paper, using the resolution cut of $M>3\times 10^{6}\,\textnormal{M}_{\odot}$.

In general, the differences between the DMO and FIRE baryonic simulations are the largest among the most massive subhalos ($V_{\mathrm{max}}\gtrsim 25~{}{\rm km}\ {\rm s}^{-1}$), which form a non-negligible amount of stars and are therefore influenced by both internal and external baryonic effects. These are the subhalos that are important for the too-big-to-fail problem, and the relative lack of these subhalos in the baryonic simulations is key for solving the too-big-to-fail problem in the satellite populations of those hosts. Surprisingly, even without including the internal changes driven by feedback and bursty outflows, the embedded disk simulations capture these quantitative trends remarkably well: the accuracy of the embedded disk simulations is largely independent of subhalo mass at $r\lesssim 100~{}\mathrm{kpc}$. This agreement suggests that subhalo survival near the center of the host halo is dominated by gravitational interactions, such that these subhalos are likely to suffer the same fate of disruption independent of their internal structure (on $\sim 500~{}\mathrm{pc}$ scales). However, when including all subhalos out to $300~{}\mathrm{kpc}$ (see Appendix A), the disk simulations do overpredict the number of subhalos with $V_{\mathrm{max}}\gtrsim 25~{}{\rm km}\ {\rm s}^{-1}$ by a factor of $3-5$ relative to the baryonic simulations, indicating that internal baryonic processes remain important in regulating the population of subhalos at larger distances. These results generally agree with those of Brooks & Zolotov (2014), who found that tidal effects dominate for subhalos that pass near the galactic disk.

Therefore, the tidal field from the central galaxy is partially responsible for eliminating some of the subhalos that a too-big-to-fail analyses would identify as problematic ‘massive failures’ in a DMO simulation. As shown in Wetzel et al. (2016), m12i does not suffer from a too-big-to-fail problem when simulated with FIRE baryonic physics, indicating that a combination of internal feedback and tidal interactions resolve TBTF around that galaxy. The embedded disk potential alone reduces the number of subhalos with $V_{\mathrm{max}}>25~{}{\rm km}\ {\rm s}^{-1}$ within $300~{}\mathrm{kpc}$ from $12$ to $7$, suggesting that the central galaxy is responsible for roughly half of the necessary changes for this system, with internal feedback driving lower inner densities in the surviving subhalos and accounting for the remaining discrepancy. We draw similar conclusions from m12f, which hosts $9$ subhalos with $V_{\mathrm{max}}>25~{}{\rm km}\ {\rm s}^{-1}$ in the DMO run, but only $5$ when simulated with an embedded disk, and identically zero when simulated with full physics. Of course, the exact relative contribution from internal feedback versus external tidal forces depends strongly on the individual histories of those large subhalos, their spatial distribution within the host and, as we will show below, their orbital characteristics.

3.2 Destruction versus mass stripping?

In general, the tidal force from the central disk potential can either strip off a portion of the outer mass of a subhalo, shifting it to lower $M$ and $V_{\mathrm{max}}$, or it can completely destroy the subhalo, either through tidal shocking (Gnedin et al., 1999) or repeated stripping events, removing it from the population entirely. In order to distinguish between subhalo mass loss from total destruction, we measure subhalos using $M_{\mathrm{peak}}$, defined as the largest virial mass that the main branch of each subhalo’s merger tree reached. Thus, this quantity is preserved for a subhalo after infall and is not susceptible to mass loss, provided the subhalo survives to $z=0$.

Figure 4 shows the cumulative counts of resolved subhalos within $100~{}\mathrm{kpc}$ as a function of $M_{\mathrm{peak}}$. If the central galaxy typically strips subhalos, but does not entirely destroy them, then the $M_{\mathrm{peak}}$ functions from the FIRE and disk simulations will agree with those of the DMO runs. Instead, at fixed $M_{\mathrm{peak}}$, the number of subhalos is at least $2-3\times$ larger in the DMO simulation, meaning that the central galaxy has destroyed at least $50-70\%$ of the substructure that currently resides within $100~{}\mathrm{kpc}$ of the halo center. As before, the embedded disk captures all but $\sim 10$–$20\%$ of this destruction, independent of $M_{\mathrm{peak}}$ except for the few largest subhalos in m12f.

3.3 Which subhalos are destroyed?

Thus far, we have examined subhalos only as a function of their position at $z=0$. However, we also can ask: out to what orbital pericenter distance are subhalos significantly affected by the presence of the central disk? The left panel of Figure 5 shows the cumulative distribution of the pericentric distance, $d_{\mathrm{peri}}$, that subhalos within $100~{}\mathrm{kpc}$ that survive to $z=0$ experienced, normalized to the total number of subhalos within $100~{}\mathrm{kpc}$. We define $d_{\mathrm{peri}}$ as the smallest physical distance reached between the main branch of the host and that of a given subhalo.555The positions of the two halos are interpolated in scale factor with a third-order spline to improve the time resolution, as in Fillingham et al. (2015). The curves, which use the same colors and styles as previous Figures, are thus required to equal $1$ at $100~{}\mathrm{kpc}$. The lower panel plots the ratio of the cumulative, non-normalized distributions, again relative to the FIRE baryonic simulations. The right panel of Figure 5, meanwhile, plots the (non-normalized) differential distribution of these same subhalos.

Clearly, when a central galaxy is present, the orbit of a subhalo is important to its survival, and the differences in the distributions are stark: the median pericentric distances of the DMO samples are $20~{}\mathrm{kpc}$ smaller than when the disk is included. Nearly every subhalo that reaches the central $10~{}\mathrm{kpc}$ is destroyed by the central galaxy: only 2 (3) subhalos that have passed within $10~{}\mathrm{kpc}$ of the center of m12i (m12f) survive to $z=0$ in the baryonic simulations, whereas $\sim 100$ such objects exist in the DMO runs. Similarly, $50\%$ of the surviving subhalos have pericentric distances $\lesssim 20~{}\mathrm{kpc}$ in the DMO simulation of m12f; that fraction drops to $\sim 5\%$ in both the baryonic and embedded disk simulations. The impact of the central galaxy is much weaker for subhalos that have never passed within $\sim 30~{}\mathrm{kpc}$. This is in qualitative agreement with previous simulations by D’Onghia et al. (2010), who reached similar conclusions by inserting a disk potential into a halo extracted from a cosmological simulation.

The embedded disk simulations match the baryonic simulations well: the median pericentric distances are only $\sim 5~{}\mathrm{kpc}$ smaller in the former, and the lower panels demonstrate that the distributions match better than a factor of $2$ down to $d_{\mathrm{peri}}\approx 20~{}\mathrm{kpc}$, though the embedded disk does overestimate the number of subhalos with $d_{\mathrm{peri}}\sim 15~{}\mathrm{kpc}$. The disk simulations remain an improvement over purely DMO simulations, however: the latter over-predict the number of surviving subhalos that have passed within $10~{}\mathrm{kpc}$ of a MW-like galaxy by a factor of $\sim 50$. The remaining excess in the disk simulations is at least partially caused by the reduced central mass at early times due to our exclusion of gas from the disk model: doubling the mass of the disk at all times (§3.4) completely eliminates this discrepancy within $10~{}\mathrm{kpc}$, though such a change does overcorrect for the total mass in the FIRE simulations after $z\sim 1.5$.

The distributions plotted in Figure 5 are largely independent of bound mass $M$: increasing the mass cut to only include subhalos with $M>3\times 10^{7}\textnormal{M}_{\odot}$ yields nearly identical trends. Therefore, the sharp contrast between the DMO and embedded disk simulations further suggests that the primary role of the disk is to completely destroy subhalos and remove them from the population, particularly for subhalos that come near to the galaxy. In fact, the differences in the count of resolved subhalos within $100~{}\mathrm{kpc}$ of both hosts can be completely accounted for by the difference in the number of subhalos with $d_{\mathrm{peri}}<20~{}\mathrm{kpc}$.

The distributions in Figure 5 also inform our understanding of the differences between m12i and m12f. Specifically, the relative depletion in m12f is larger than in m12i, likely because the pericenter distribution in the former is skewed to smaller radii, making its subhalo population more sensitive to the presence of a central galaxy.666While m12i and m12f have a similar total mass within a 2 Mpc sphere, m12f is within a more filamentary environment, which may explain why its satellites have more radial orbits. This difference motivates the need to run many such simulations to explore the full range of subhalo orbital characteristics between hosts.

Figure 2 showed that m12i and m12f baryonic simulations have no resolved subhalos within $17-20\,\mathrm{kpc}$ of the central galaxy at $z=0$. However, Figure 5 shows that $\approx 20$ subhalos have orbited within this distance at some point, they simply are not located within this region by $z=0$. Furthermore, Figure 5 plots only pericentric distances of surviving subhalos, and those that have been destroyed also could have induced dynamical effects in the inner halo, such as creating holes in stellar streams, prior to their destruction. That said, our results show that the total timescales over which such subhalos could act is significantly shorter. We leave the broader question of how streams form and evolve in fully cosmological simulations that include a central galaxy for later work.

Figure 5 clearly shows that $\gtrsim 90\%$ of subhalos that pass within $10$–$20~{}\mathrm{kpc}$ of the central galaxy are destroyed. This dependence on pericenter implies that the velocity distribution of surviving subhalos should be significantly biased relative to DMO simulations, such that systems on radial, plunging orbits with low specific angular momentum should be preferentially destroyed.

The left panel of Figure 6 shows the cumulative distribution of the tangential velocities, $V_{\mathrm{tan}}$, of subhalos within $100~{}\mathrm{kpc}$, normalized to the total number within $100~{}\mathrm{kpc}$, demonstrating the strength of this effect. While subhalos in the DMO simulation have a mean $V_{\mathrm{tan}}$ of $\sim 150~{}{\rm km}\ {\rm s}^{-1}$, the presence of a central galaxy increases it to $\sim~{}200-250~{}{\rm km}\ {\rm s}^{-1}$. The lower panel, which again plots the ratio in the absolute distributions, demonstrates that the DMO simulations over predict the number of subhalos within $100~{}\mathrm{kpc}$ with $V_{\mathrm{tan}}<100~{}{\rm km}\ {\rm s}^{-1}$ by a factor of ten. The embedded disk simulations, however, agree to within a factor of $2$–$3$ at all $V_{\mathrm{tan}}$. While the suppression of subhalos with low $V_{\mathrm{tan}}$ in runs with a central galaxy is primarily caused by the destruction of those systems, the deeper potential well in the central regions also leads to an increase in the velocities of surviving systems near the galaxy.

We similarly examine the distributions of subhalo radial velocities in the right panel of Figure 6. However, in contrast to the tangential distributions, we find that they agree across all simulations to within $\sim 50\%$ at most $V_{\mathrm{rad}}$. We do find that the relative fraction of subhalos with large, positive radial velocities ($V_{\mathrm{rad}}\gtrsim 75~{}{\rm km}\ {\rm s}^{-1}$) is slightly suppressed in the simulations with a central galaxy relative to the DMO simulations, as expected if subhalos that have passed close to the halo center and would be moving away from the host are preferentially destroyed. However, overall these results indicate that a subhalo’s likelihood of being destroyed by the central galaxy does not depend significantly on its radial velocity. As a result, the primary effect of the central galaxy is on the resultant distribution of tangential, and not radial, velocities of the surviving subhalo population. There is a slight trend wherein the pure DMO simulations underpredict the fraction of high $\left|V_{\mathrm{rad}}\right|$ subhalos (i.e., overpredict the fraction of systems with $V_{\mathrm{rad}}\sim 0~{}{\rm km}\ {\rm s}^{-1}$), consistent with an overall increase in subhalo velocities in the presence of a galaxy due to the added mass at the center of the host.

This tangential velocity bias may have direct consequences on the expected perturbations induced in the cold stellar streams that are suggested as best able to constrain dark subhalo counts in the MW, which stretch along the plane of the sky (e.g. Newberg & Carlin, 2016). Specifically, if both the observed streams and the subhalo population are tangentially biased, one would expect more ‘glancing’ blows, wherein subhalos are either co-rotating or counter-rotating relative to the orbit of the stream, and fewer perpendicular interactions, where the subhalo punches directly through a cold stream, relative to what is expected for an isotropic population of subhalos (the standard assumption in most work on stream-subhalo interactions, e.g. Yoon et al., 2011; Carlberg, 2012; Erkal et al., 2016). Counter-rotating interactions will have a small impact on streams due to their relatively short interaction times, but a larger fraction of the encounters may by co-rotating than otherwise expected. These occur over long enough timescales that the impulse approximation typically used to estimate changes in the stellar velocities (e.g. Sanders et al., 2016) may break down. We defer a full analysis of the degree to which this effect is observable to future work.

Similarly, the preferential destruction of subhalos on highly radial orbits has important implications for the structure of the stellar halo, particularly close to the MW. Specifically, the majority of the destroyed satellites that make up the stellar halo (Bullock et al., 2001; Bullock & Johnston, 2005; Bell et al., 2008) were likely on plunging, highly radial orbits, such that their resultant stellar streams are likely to be radially extended and reach to greater radii. Therefore, streams from destroyed satellites may prove useful for detecting clear, coherent gaps from the surviving dark subhalos, which are largely on tangential orbits that are perpendicular to the expected orientations of the streams from destroyed satellites.

3.4 Varying the disk parameters

In order to explore the dependence on the specific disk potential that we add to the dark matter simulations, we additionally simulate m12i with modified disk parameters. Specifically, we hold two of the three disk parameters, $R_{\mathrm{d}}(z)$, $z_{\mathrm{d}}(z)$, and $M_{\mathrm{d}}(z)$ (see Equation 1) at their fiducial values, while varying the third. We investigate: $R_{\mathrm{d}}(z)\rightarrow 2R_{\mathrm{d}}(z)$, spreading the same amount of mass over a larger surface area; $z_{\mathrm{d}}\rightarrow 1~{}\mathrm{pc}$ at all times, forcing a nearly infinitesimally thin disk; and $M_{\mathrm{d}}(z)\rightarrow 2M_{\mathrm{d}}(z)$, making the disk more massive. We also test the dynamical importance of the gas in the galaxy by holding $R_{\mathrm{d}}$ and $z_{\mathrm{d}}$ fixed at their fiducial values, but setting $M_{\mathrm{d}}$ equal to the fiducial (stellar) mass plus the mass in gas within $r_{90}$ (i.e., $M_{\mathrm{d}}\rightarrow M_{\star}+M_{\rm gas})$. In order to further test the importance of the shape and orientation of the galactic potential, we also simulate m12i with a Hernquist (1990) sphere, rather than the Miyamoto-Nagai disk potential.

Figure 7 show the bound mass functions (left) and cumulative radial distributions (right) of these simulations, along with the three versions of m12i presented previously. The solid, dashed, and dotted black lines are identical to those in the left panels of Figure 2, while the green, blue, and orange lines indicate the results of the simulations with $2R_{\mathrm{d}}$, $z_{\mathrm{d}}=1~{}\mathrm{pc}$, and $2M_{\mathrm{d}}$, respectively. The lavender lines plot counts obtained by additionally including the gas in the mass in the disk, and the light red dashed line plots subhalo counts with a Hernquist potential of identical mass and radial extent, quantified by $r_{90}$.

In summary, subhalo depletion most directly correlates with the mass of the central disk. The simulation with the $2M_{\mathrm{d}}$ disk results in fewer subhalos, especially at $r\lesssim 70\,\mathrm{kpc}$. The $2R_{\mathrm{d}}$ disk yields slightly more subhalos at small $r$, while the dependence on $z_{\mathrm{d}}$ is weak, with the $z_{\mathrm{d}}=1~{}\mathrm{pc}$ disk predicting slightly fewer subhalos than the fiducial disk. The shape, and consequently the orientation, of the potential appears to be of secondary importance, with the $z_{\mathrm{d}}=1~{}\mathrm{pc}$ disk and the Hernquist sphere both yielding very similar results to the fiducial disk. However, the slight excess in the simulation with a Hernquist potential, relative to the fiducial disk, suggests that some subhalos are destroyed by disk shocking that would otherwise survive the enhanced tidal forces (also see D’Onghia et al., 2010). Therefore, small gains may be achieved by matching the orientation of the real galaxy, but we emphasize again that even a disk with an imperfect, fixed orientation yields subhalo populations with distributions that typically agree to within $\sim 25\%$, without the additional complexity of determining the correct orientation.

Interestingly, at most masses and radii, the $2M_{\mathrm{d}}$ disk in fact yields a slightly better match to the FIRE baryonic simulation than the fiducial disk simulation.777 The $2M_{\mathrm{d}}$ simulation also provides a better match to the $d_{\mathrm{peri}}$, $V_{\mathrm{tan}}$, and $V_{\mathrm{rad}}$ distributions of the FIRE baryonic simulations than the fiducial disk. As noted in Section 2.1, we match the mass of the fiducial disk to the stellar disk that forms in the baryonic simulation, but this neglects the (complex) contribution from (fluctuating) gas in/near the disk. Therefore, we may posit that this improved match is a consequence of more accurately matching the total baryonic mass at the center of the halo at early times, when the disk is dominated by gas ($z\gtrsim 1$).

However, the simulation with $M_{\mathrm{d}}\rightarrow M_{\star}+M_{\mathrm{gas}}$ differs only slightly from the results of the fiducial disk: adding the gas to the disk mass yields marginally fewer subhalos at fixed mass and radius. Because the differences in the disk mass between this run and the fiducial run occur primarily at early times,888While adding the gas mass drastically increases the mass of the disk at early times (by a factor of $\sim$10 at $z=3$), the contribution after $z=1$ is less than a factor of two, and the total mass is within $\sim 15\%$ of the fiducial (stellar) mass by $z=0.5$. the close agreement between the two runs suggests that the mass at early times is relatively unimportant, likely because the majority of subhalos fall into the MW-size host halo after $z=1$ (e.g. Wetzel et al., 2015a). The lingering discrepancy between the $M_{\star}+M_{\mathrm{gas}}$ run and the fully hydrodynamical simulation further indicates that the gas does have important, if second-order, dynamical effects on subhalos due to a combination of internal feedback, ram pressure and viscous stripping in the hot halo, and interactions with inflows and outflows.

Overall, we emphasize that factors of $\sim 2$ changes to the parameters of the embedded disks have a relatively modest effect on the subhalo population at $z=0$, as compared with the much stronger difference from not using an embedded disk. Moreover, our experiments with the shape of the potential show that the detailed geometry and orientation of the disk (difficult to predict in non-fully baryonic simulations) are of small, secondary importance: what matters (to leading order) is a central potential of the correct baryonic mass and radial size. Thus, even if the exact galaxy that would form in a given halo is not completely constrained, one still can improve upon the predictions of DMO simulation by inserting a central galaxy with parameters drawn from empirical scaling relations.

3.5 Comparison with Sawala et al. (2017)

Recently, Sawala et al. (2017, hereafter S17) demonstrated that the subhalo populations of four of the hosts in the APOSTLE project (Sawala et al., 2016; Fattahi et al., 2016), which simulates MW/M31-like pairs of halos using the EAGLE prescriptions for star formation and feedback (Schaye et al., 2015), are suppressed relative to their DMO counterparts. Here we explicitly compare their results to the differences that arises between the FIRE simulations and their DMO analogues.

Figure 1 of S17 indicates that at the lowest masses resolved by the APOSTLE simulations, the DMO counterparts overpredict subhalo counts within 300 kpc by a factor of $\sim$1.3. While the DMO analogues of the Latte simulations overpredict subhalo counts within 100 kpc by a factor of $\sim$2-3 (Figure 2), the ratio of $N(>M)_{\mathrm{DMO}}/N(>M)_{FIRE}$ within 300 kpc is roughly 1.2 - 1.5 (Figure 9), in close agreement with the results of S17. At smaller radii, however, Figure 3 of S17 suggests that the DMO simulations typically overpredict counts within 50 kpc by a factor of $\sim$1.75 - 2, while our results (Figure 8) suggest that DMO simulations overpredict counts within 50 kpc by a factor of 3 - 4. Accordingly, we predict a larger discrepancy (between DMO and hydrodynamical simulations) in the radial distributions at small distances: subhalo densities at $r\sim 30$ kpc are a factor of $\sim$4 times higher in the DMO analogues than in the FIRE simulations, while Figure 4 of S17 indicates that their DMO analogues only predict a factor of $\sim$1.6 discrepancy at the same radius.

Therefore, while our general conclusions regarding depleted radial distributions agree well, results from the FIRE simulations suggest that the discrepancy grows larger at small radius ($r<100$ kpc), but the results of S17 indicate a roughly constant offset within $\sim$50 kpc. While this disagreement is second-order, it may be due to differences in the central masses of the galaxies that form in the two sets of simulations. Per Figure 7 of Fattahi et al. (2016), the APOSTLE disks tend to lie slightly below expectations from abundance matching analyses. Both m12i and m12f, however, fall slightly above those same relationships (Hopkins et al., 2017). As evidenced by the parameter survey presented in §3.4, increasing the mass density at the center of host destroys more satellites at small radii without strongly altering results at $r\gtrsim 100~{}\mathrm{kpc}$.

We caution, however, that the results in S17 are based on the subhalo populations of only four hosts, and that our results include only two hosts. Because some degree of halo-to-halo scatter is expected, it’s not clear that our results should agree exactly with those of S17. It is encouraging, therefore, that all six hosts demonstrate a similar trend: a clear reduction in subhalo populations, with the largest differences occuring within $\sim 50$ - 100 kpc and with subhalos on radial orbits (low tangential velocities) preferentially depleted (Figure 8 of S17).

4 Implications

Dark matter-only simulations have been used extensively to interpret data and to make predictions for subhalo observables. Here we briefly discuss several such investigations that are likely affected by our results and speculate how the effects of a central galaxy potential in more massive halos could affect the broader use of dark matter simulations to interpret data.

As discussed in §1 and §3, the significant depletion that we see in our baryonic and embedded disk simulations at small radii has important implications for ongoing searches for substructure in the MW halo. Perhaps the most promising method involves looking for gaps and other anomalies in the stellar streams of Palomar-5 and GD-1. Both of these streams exist within $20$ kpc of the Galaxy, therefore they sit within the region that is most severely depleted (see Figure 2). Our results imply that there may be no dark substructures in such a region today to affect these streams, and that even if they orbited in this region in the past, their lifetime is much shorter than predicted from DMO simulations. However, further work is needed to sample full subhalo evolutionary histories, and as we discuss in the next section, a more focused effort on simulating tuned MW analogues (set to the mass of the MW galaxy, rather than being matched to a comparison hydrodynamical simulation) with a statistical sample of initial conditions will be required to determine the range of subhalo counts that we expect around our Galaxy within this region.

More generally, any analyses that depends strongly on the subhalo mass function, the radial distribution, or the velocity distribution, particularly near the center of the host, will be severely impacted by the presence of a central galaxy. One example is the expected completeness corrections to the faint-end of the MW stellar mass function (Tollerud et al., 2008; Hargis et al., 2014). This depends non-trivially on the radial distribution of subhalos, because the correction is based on how many more subhalos we expect to exist within some large volume (for example, $r\leq 400$ kpc) given an observed number within a smaller (for example, $r\leq 100$ kpc), observationally complete volume. The enhanced central depletion in our simulations suggests that there are fractionally even more satellites waiting to be discovered at large radii, given the number that we currently observe within the completeness limits of SDSS, DES, or Pan-STARRS. As pointed out by Ahmed et al. (2017), who also found that subhalo populations are more radially extended in baryonic simulations than in DMO, depleting the central region of substructure has the additional effect of increasing the statistical significance of any potential “planes of satellites.”

The missing satellites problem is clearly reduced in severity by the destruction of subhalos with $M\sim 10^{9}$–$10^{10}\,\textnormal{M}_{\odot}$. Similarly, the tension in comparing the number of dwarf galaxies containing old stars around the MW with expectations from the observed steep luminosity functions during the reionization era, as discussed by Boylan-Kolchin et al. (2014), will be reduced; many of those ancient stars instead would be dispersed into the stellar halo via the enhanced disruption brought on by the central disk. Similarly, the constraints on stellar mass functions at $z\sim 5$ based on MW galaxy counts presented in Graus et al. (2016) will be modified.

Limits on warm dark matter models (e.g. Horiuchi et al., 2014, 2016) associated with requiring enough subhalos to host the known satellites will become even tighter. Predictions for the substructure boost for dark matter annihilation signals in the Galactic center (Kuhlen et al., 2008; Bovy, 2009) will decrease. However, subhalos would still be expected to contribute to a boost at larger radii, which could be important for indirect detection searches around the M31 halo or in the all-sky background from the MW halo itself.

The reduction in subhalo counts also modifies the results of using counts of satellite galaxies in the Local Group to constrain the $M_{\star}-M_{\mathrm{halo}}$ relation at low masses (e.g. Garrison-Kimmel et al., 2014a; Brook et al., 2014; Garrison-Kimmel et al., 2017). Specifically, matching to fewer dark matter halos, as implied by our results, requires a flatter log-slope, which also shifts the relation closer towards one that alleviates the ‘too-big-to-fail’ problem (Boylan-Kolchin et al., 2011, 2012; Garrison-Kimmel et al., 2014b; Papastergis et al., 2015, and also see Jethwa et al., 2016). The effect of enhanced disruption from a central disk in alleviating the too-big-to-fail problem was first emphasized by Brooks & Zolotov (2014).

The preferential destruction of subhalos on plunging orbits also has consequences for analyses that utilize the velocity distributions of subhalos in DMO simulations. For example, the number of ‘backsplash galaxies’ (galaxies beyond the virial radius of a more massive host that orbited within the virial radius in the past; e.g. Teyssier et al., 2012; Wetzel et al., 2014; Garrison-Kimmel et al., 2014a) that have passed close to the central galaxy should be drastically reduced, though the number that have had only glancing encounters with the host halo will be largely unaltered. Similarly, the works of Fillingham et al. (2015) and Wetzel et al. (2015b) explored the timescales over which the satellite galaxies of the MW and M31 had their star formation quenched after infall using subhalo catalogs taken from the ELVIS simulations. While they found that the environmental quenching timescales of satellite dwarf galaxies are short ($\lesssim 2$ Gyr; much shorter than for more massive satellites), our results imply even shorter quenching timescales, because many of the subhalos that fell in at early times have preferentially smaller pericenters (Wetzel et al., 2015a) and therefore should be destroyed by the central galaxy. Because satellites with low tangential velocities also are destroyed with high efficiency, however, the constraints on ram-pressure stripping from Fillingham et al. (2016) may become less stringent. Similarly, the recent discrepancy pointed out by Cautun & Frenk (2017), wherein the MW satellites are on more tangentially biased orbits than expected from $\Lambda$CDM, is alleviated by the destruction of these low $V_{\mathrm{tan}}$ subhalos.

While our analysis focused on MW-size halos, substructure depletion caused by a central galaxy is likely important around larger ($M_{\mathrm{halo}}\sim 10^{13}\,\textnormal{M}_{\odot}$) galaxies as well. This is particularly relevant for searches for dark subhalos via lensing anomalies (e.g. Vegetti et al., 2010; MacLeod et al., 2013; Nierenberg et al., 2014; Hezaveh et al., 2016). However, the strength of the impact of these more massive galaxies on their subhalos is unclear. The MW represents the mass at which the ratio of central galaxy stellar mass to host halo mass is highest (Leauthaud et al., 2012; Behroozi et al., 2013b; Moster et al., 2013), so this is the mass scale at which the central galaxy most strongly affects the halo potential. Thus, the effect in larger hosts may be weaker. Furthermore, more massive halos generally assemble at later times (with satellites falling in later), suggesting that subhalo destruction may be further suppressed relative to our findings here.

In fact, Fiacconi et al. (2016) recently presented two baryonic simulations of galaxies in this mass regime, and found that baryonic contraction actually increases the number of massive subhalos near the center of the host relative to DMO simulations. However, Graus et al. (in preparation) perform a similar analysis using the Illustris and Illustris Dark simulations (Vogelsberger et al., 2014b, a) and find the opposite effect: the baryonic simulations have fewer subhalos at $\sim 10^{8}\,\textnormal{M}_{\odot}$ around lens hosts, with the largest deficit near the central regions. Despali & Vegetti (2017) similarly compared the DMO and baryonic versions of the Illustris and EAGLE (Schaye et al., 2015) simulations, and demonstrated that subhalo counts around hosts of mass $10^{12.5}$–$10^{14}\,\textnormal{M}_{\odot}$ are suppressed by $\sim 30\%$ at subhalo masses $10^{9}\,\textnormal{M}_{\odot}$ in the baryonic simulations, with the largest suppression near the center of the hosts (also see Chua et al., 2016, who reached similar conclusions using the Illustris simulations). The impact on the lowest-mass subhalos, however, remains largely unexplored due to the difficulty of simulating galaxies of this mass at high resolution with baryonic physics.

Finally, if a similar degree of substructure depletion occurs for roughly LMC to MW-size subhalos ($M_{\mathrm{halo}}\sim 10^{11}$–$10^{12}\,\textnormal{M}_{\odot}$) around massive hosts ($M_{\mathrm{halo}}\gtrsim 10^{14}\,\textnormal{M}_{\odot}$), it could have important implications for the use of DMO simulations to interpret small-scale (‘one-halo’) clustering statistics through subhalo abundance matching (e.g. Conroy et al., 2006; Reddick et al., 2013), as compared with studies that examined these trends for subhalo disruption without modeling a central galaxy (e.g. Wetzel & White, 2010). Specifically, if the subhalos that are massive enough to host bright galaxies ($M_{\star}\gtrsim 10^{10}\,\textnormal{M}_{\odot}$) are preferentially depleted in galaxy groups owing to the potential of the brightest cluster galaxy, then abundance matching analyses would tend to assign too little stellar mass to subhalos of a given $M$ to avoid overproducing those galaxies. If this effect is important in groups and clusters, then it may also influence the use of similar approaches to understand trends between galaxy color and age with subhalo accretion times (e.g. Mahajan et al., 2011; Wetzel et al., 2013; Hirschmann et al., 2014; Oman & Hudson, 2016). However, the increased concentrations of MW-mass subhalos in baryonic simulations may also make these systems more resistant to stripping (e.g. Chua et al., 2016), such that even if there are fewer satellites overall, a greater proportion of those remaining would have $M_{\mathrm{halo}}\sim 10^{12}\textnormal{M}_{\odot}$ than DMO simulations would predict.

5 Conclusions

The interplay between dark matter and baryons, and the subsequent importance of baryons in correctly predicting the properties of DM halos with cosmological simulations has a rich history in this field (e.g. di Cintio et al., 2011; Governato et al., 2012; Pontzen & Governato, 2012; Zolotov et al., 2012; Brooks & Zolotov, 2014; Arraki et al., 2014; Sawala et al., 2015; Di Cintio et al., 2014; Oñorbe et al., 2015; Chan et al., 2015; Sawala et al., 2016; Cui et al., 2016; Wetzel et al., 2016). In this paper, we explored the role of baryons in affecting the low-mass subhalos of interest for observational searches for dark subhalos: dark matter mass $M>3\times 10^{6}\,\textnormal{M}_{\odot};\>V_{\rm max}>5~{}{\rm km}\ {\rm s}^{-1}$.

Our exploration relied on two MW-mass dark matter halos from the Latte simulation suite (Wetzel et al., 2016), each simulated (1) with full baryonic physics from the FIRE project, (2) with dark matter only, and (3) with dark matter plus an embedded disk potential for the central galaxy that evolves to match the stellar disk in the corresponding baryonic simulation. As shown in Figure 2, relative to the DMO simulations, subhalo counts in the baryonic simulations are lower by a factor of $\sim 5$ within $25$ kpc of the halo center. Both of the baryonic simulations are completely devoid of substructure within $15$ kpc at $z=0$, though $\sim 20$ subhalos with pericentric distances $<20~{}\mathrm{kpc}$ survive to the present day. Subhalo depletion becomes less important with increasing radius: within $\sim 300$ kpc (roughly the halo virial radius), subhalo counts in the baryonic simulations are lower by only $\sim 15$–$30\%$ compared to the DMO simulations. This depletion is driven by the preferential destruction of subhalos on radial orbits that get to $d_{\rm peri}\lesssim 30\,\mathrm{kpc}$ from the halo center. This in turn biases the orbital velocities of surviving subhalos at $z=0$ to be more tangential in simulations that have central galaxies than they are in DMO simulations.

Importantly, the simulations with embedded disk potentials reproduce at least $\sim 75\%$ of the depletion evident in the mass functions and radial distributions. This good agreement is important for two reasons. First, because we carefully matched the evolving disk potential to the stellar disk in the baryonic simulations, this comparison provides a clear physical explanation for the origin of subhalo depletion in the baryonic simulations: the bulk of the subhalo depletion arises simply from the tidal field of the central galaxy itself. These results also imply that the majority of subhalo depletion is independent of the exact model of feedback or star formation, because our embedded disk simulations include none of these. However, we emphasize that we matched the evolving disk potential to the stellar mass in the baryonic simulations; this therefore provides a conservative under-estimate of the total baryonic mass of the central galaxy, especially at early times when the simulated galaxies were gas rich. As Figure 7 showed, though, adding the approximate gas mass to the disk yields only marginally better agreement between the embedded disk and baryonic simulations. Therefore, a small fraction of the substructure depletion is not due to interactions with the central galaxy, and cannot be captured with this technique.

Second, given its significantly improved accuracy as compared with DMO simulations, our method of embedding a central disk potential provides an inexpensive way to significantly improve the accuracy of DMO simulations. The population statistics of surviving subhalos in the embedded disk simulations display much better agreement with those in the FIRE baryonic simulations in every statistic that we have checked, including infall times, pericentric distances, radial velocities, total orbital velocities, radial profiles, and counts as a function of mass at $z=0$ or at infall. We also find that the disk simulations yield more accurate counts around the hosts at higher redshift, though the differences are less dramatic than at $z=0$. Thus, one can vastly improve upon predictions from purely DMO simulations, at nearly the same CPU cost, by simply including a central galactic potential.

Moreover, as Figure 7 showed, subhalo depletion is largely insensitive to the central disk thickness, and it depends only mildly on disk radius and the detailed shape of the potential. The most important parameter is simply the central galaxy mass. Thus, one can (to reasonable approximation) adopt a simple, spherically-symmetric analytic potential with parameters taken from observations (e.g. abundance-matching) or large-volume simulations.

The differences that we see between simulations with and without central galaxies are particularly important for dark substructure searches that rely on cold stellar streams within 20 kpc of the Galaxy. However, it is difficult to estimate the expected halo-to-halo variance based on just two halos. As shown in Figure 5, the intrinsic pericenter distribution for a given subhalo population can vary considerably, which likely affects subhalo disruption considerably. A larger suite of hosts simulated with embedded disks that match the MW is required to make more concrete statistical statements (Kelley et al., in preparation).

We demonstrated that baryonic effects are crucially important for interpreting ongoing substructure searches. Galaxies exist within the centers of halos and are dynamically important within the vicinity of disks like the MW. Our method of embedding galactic potentials in cosmological zoom-in simulations provides an avenue to producing more accurate substructure predictions, relative to baryonic simulations, without the millions of CPU hours required by those simulations. These embedded galaxies therefore hold the promise of improved predictions for the statistical properties of subhalo populations, which are necessary to fully interpret the results of many upcoming observations, including those aimed at detecting tiny, dark subhalos.

Acknowledgements

The authors thank Jo Bovy, Mia Bovill, and Brandon Bozek for valuable discussions and the anonymous referee for several suggestions that have improved the manuscript. The authors also thank Erik Tollerud for providing software used to create visualizations and Alexander Knebe, Peter Behroozi, and Oliver Hahn, respectively, for making AHF, consistent-trees, and MUSIC publicly available.

Support for SGK was provided by NASA through Einstein Postdoctoral Fellowship grant number PF5-160136 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. AW was supported by a Caltech-Carnegie Fellowship, in part through the Moore Center for Theoretical Cosmology and Physics at Caltech, and by NASA through grant HST-GO-14734 from STScI. JSB, and TK were supported by NSF grant AST-1518291 and by NASA through HST theory grants (programs AR-13921, AR-13888, and AR-14282.001) awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under NASA contract NAS5-26555. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. MBK acknowledges support from the National Science Foundation (grant AST-1517226) and from NASA through HST theory grants (programs AR-12836, AR-13888, AR-13896, and AR-14282) awarded by STScI. CAFG was supported by NSF through grants AST-1412836 and AST-1517491, and by NASA through grant NNX15AB22G. DK acknowledges support from NSF grant AST-1412153 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. EQ was supported by NASA ATP grant 12-APT12-0183, a Simons Investigator award from the Simons Foundation, and the David and Lucile Packard Foundation. Support for ASG was provided by NSF grant AST-1009973.

Numerical calculations were run on the Caltech compute cluster ‘Zwicky’ (NSF MRI award #PHY-0960291) and allocation TG-AST130039 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF. Resources supporting this work were also provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This work also made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy (Jones et al., 01), ipython (Perez & Granger, 2007), Mayavi (Ramachandran & Varoquaux, 2011), and NASA’s Astrophysics Data System.

Appendix A: Distributions within $50$ and $300~{}\mathrm{kpc}$

The Figures presented in the main body of the paper plot subhalo counts, as a function of $V_{\mathrm{max}}$, $d_{\mathrm{peri}}$, and $V_{\mathrm{tan}}$, within $100~{}\mathrm{kpc}$. For completeness, we plot here similar distributions, but including subhalos within either $50~{}\mathrm{kpc}$ (Figure 8) or $300~{}\mathrm{kpc}$ (Figure 9) of the host centers. The fractional difference in the purely DMO simulations are lower on $\sim 300~{}\mathrm{kpc}$ scales than within either $50$ or $100~{}\mathrm{kpc}$, but the impact of the central galaxy remains present in the lack of subhalos with low $V_{\mathrm{tan}}$ or small $d_{\mathrm{peri}}$. Within $50~{}\mathrm{kpc}$, the central galaxy is extremely destructive: only $\sim 25\%$ of the subhalos present in the DMO simulations survive in the presence of a disk. Similarly, no subhalos with $V_{\mathrm{tan}}\lesssim 100~{}{\rm km}\ {\rm s}^{-1}$ remain within $50~{}\mathrm{kpc}$ of the galaxy.

Appendix B: Resolution

Because of the cost of simulating the systems at even higher resolution than that presented in the main body, we instead establish that we reliably identify substructures at the same particle count using lower-resolution simulations. Specifically, we compare counts as a function of mass and radius to simulations of m12i with particle masses $8$ times larger than our fiducial simulations ($m_{\rm p}=3.4\times 10^{5}\,\textnormal{M}_{\odot}$) and with softening lengths a factor of two larger. A detailed discussion of resolution in the FIRE simulations can be found in Hopkins et al. (2017).

The left panel of Figure 10 shows cumulative counts as a function of bound mass assigned by AHF, $M$, within $100~{}\mathrm{kpc}$ of m12i, corrected for $f_{\mathrm{b}}$. The black lines present counts in simulations with the fiducial disk parameters, while the brown lines plot the DMO counterparts; fiducial resolution simulations are plotted in solid, and those at lower resolution are dashed. The lower sub-panel shows the ratio between the low and high resolution versions of the run with and without the disk. The dotted vertical line represents the resolution cut of $3\times 10^{6}\textnormal{M}_{\odot}\sim 85$ particles adopted in the main text based on inspection of the differential mass functions. The dashed vertical line indicates a factor of $8$ larger mass, corresponding to an identical number of particles in the lower resolution simulation ($M=2.4\times 10^{7}\,\textnormal{M}_{\odot}$, after multiplying by $1-f_{\mathrm{b}}$). Counts in the disk simulations agree well at that mass (though there are fluctuations at higher mass), and counts in the DMO simulations only differ by $\sim 15\%$. Convergence is generally even better on larger scales: counts agree to within $10\%$ at nearly all $M>2.4\times 10^{7}\textnormal{M}_{\odot}$. The low-resolution DMO simulation does, however, underpredict the subhalo count within $300~{}\mathrm{kpc}$ by $\sim 10\%$, and only agrees for $\sim 300$ particles; at that particle count, however, the low resolution disk simulation overpredicts the subhalo count by $\sim 10\%$.

The right panel of Figure 10 shows the radial distribution of subhalos with $M>2.4\times 10^{7}\textnormal{M}_{\odot}$ in the same four simulations. Counts in the disk runs agree to within a few percent at nearly all radii, and the DMO simulations are within $\sim 10$–$15\%$ of one another. The offset in the latter is roughly constant with radius, indicating that the deficiency discussed above is relatively independent of distance. The most significant deviations occur at $r\lesssim 40\mathrm{kpc}$ or $M>3\times 10^{8}\,\textnormal{M}_{\odot}$, where the number of subhalos is small ($<10$) and therefore subject to significant scatter, depending the subhalos’ orbital phases. In general, the differences between the resolution levels are small compared to the differences between the DMO and disk simulations.