Milky Way Satellite Census. II. Galaxy–Halo Connection Constraints Including the Impact of the Large Magellanic Cloud

Our model of the underlying MW satellite population is built on high-resolution dark matter–only zoom-in simulations of MW-mass host halos selected from the suite of 45 hosts presented in Mao et al. (2015), which have virial masses between 1.2 and $1.6\times {10}^{12}\ {M}_{\odot }$.⁵⁴ The highest-resolution particles in these simulations have a mass of $3\times {10}^{5}\ {M}_{\odot }\ {h}^{-1}$, and the softening length in the highest-resolution regions is $170\ \mathrm{pc}\ {h}^{-1}$.

Halo catalogs and merger trees were generated using the Rockstar halo finder and the consistent-trees merger code (Behroozi et al. 2013b, 2013c). Subhalos in these simulations are well resolved down to a present-day maximum circular velocity of ${V}_{\max }\approx 9\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Mao et al. 2015). To be conservative, we only use subhalos with both ${V}_{\max }\gt 9\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and peak maximum circular velocity ${V}_{\mathrm{peak}}\gt 10\,\mathrm{km}\,{{\rm{s}}}^{-1}$. In Appendix D, we show that these resolution thresholds are sufficient for modeling the satellite populations of interest here.

The MW might be atypical compared to the average halos of a similar mass (e.g., Boylan-Kolchin et al. 2010; Busha et al. 2011; Rodríguez-Puebla et al. 2013; Fielder et al. 2019); in particular, its satellite population is likely affected by the existence of the LMC system and the "satellites of satellites" that accreted with the LMC into the virial radius of the MW (Lynden-Bell 1976; D'Onghia & Lake 2008; Lu et al. 2016; Dooley et al. 2017). In addition, the detailed merger history of the MW—such as the early accretion of an LMC-mass galaxy inferred from Gaia data—might affect its faint satellite population (Bose et al. 2019).

Thus, we select MW-like host halos that each have an LMC analog with realistic internal and orbital properties; both of these hosts experience an early major merger similar to the Gaia-Enceladus accretion event (see Appendix A.3 for details). We define realistic LMC analogs as subhalos with

• 1.  
  present-day maximum circular velocity ${V}_{\max }\geqslant 55\,\mathrm{km}\,{{\rm{s}}}^{-1}$,
• 2.  
  present-day heliocentric distance $40\ \mathrm{kpc}\lt D\,\lt \,60\ \mathrm{kpc}$, and
• 3.  
  time of accretion onto the MW less than 2 Gyr ago.

These criteria yield two MW-like host halos with virial masses of 1.57 and $1.26\times {10}^{12}\ {M}_{\odot }$, respectively. Both of these hosts were used in the less restrictive host halo set defined in Nadler et al. (2019b), and both have a Navarro–Frenk–White (NFW) concentration parameter that is consistent with constraints set using the combination of satellite and globular cluster dynamics measured by Gaia (Callingham et al. 2019; Watkins et al. 2019). The LMC analogs in these two simulations have present-day virial masses of 1.6 and $2.5\times {10}^{11}\ {M}_{\odot }$, respectively, and both have peak virial masses of $3.0\times {10}^{11}\ {M}_{\odot }$. These LMC analogs accreted onto their host halos 1 and 1.5 Gyr ago, respectively, and their orbital dynamics are consistent with LMC proper-motion measurements (e.g., Kallivayalil et al. 2013).

Our fiducial LMC analogs have masses that are consistent with LMC mass estimates based on stellar stream dynamics, satellite dynamics, and the orbital histories of both Magellanic Clouds (Besla 2015; Peñarrubia et al. 2016; Erkal & Belokurov 2019; Erkal et al. 2019).⁵⁵ However, different studies have adopted various definitions of "LMC mass," and precision in the LMC mass definition (and particularly in the distinction between peak and present-day halo mass) is crucial going forward. We expect that our inference is most sensitive to the peak mass rather than the present-day mass of the LMC because peak mass correlates more directly with the expected abundance of LMC satellites, particularly for recent infall scenarios. Other probes of LMC mass are likely sensitive to these quantities in different ways, and some—including timing arguments (e.g., Peñarrubia et al. 2016) and orbit-rewinding to infer LMC satellite abundances (e.g., Patel et al. 2020)—might be most sensitive to the ratio of the LMC and MW halo masses.

Although the masses of our host halos are consistent with observational constraints for the MW (e.g., Busha et al. 2011; Bland-Hawthorn & Gerhard 2016; Patel et al. 2017), our simulations span a narrower range of host mass relative to the uncertainty on this quantity inferred from Gaia measurements. For example, Callingham et al. (2019) found that the MW host virial mass lies between 1.0 and $1.8\times {10}^{12}\ {M}_{\odot }$ at the 95% confidence level (also see Cautun et al. 2019b; Li et al. 2019a, 2019b). Since subhalo abundance is proportional to host halo mass, predicted satellite abundances scale linearly with MW mass, modulo second-order changes in subhalo disruption due to variations in the mass accretion history of the central galaxy (Kelley et al. 2019; Samuel et al. 2020). Ideally, our analysis would be performed using MW-like host halos—all of which include realistic LMC analogs—that bracket the current range of allowed MW host mass; however, the availability of such MW-like systems is limited by the statistics of our simulations. Thus, we do not marginalize over the full range of allowed MW host masses in this work. We estimate the potential impact of this uncertainty in Appendix A.1.

To associate satellite luminosities with subhalos, we follow Nadler et al. (2019b) by employing an abundance-matching procedure that monotonically relates the absolute V-band magnitude of satellites, M_V, to the peak circular velocity of subhalos, V_peak.⁵⁶ This relation is constrained by the GAMA survey (Loveday et al. 2015; Geha et al. 2017) for bright systems (${M}_{V}\lt -13\ \mathrm{mag}$) and is extrapolated into the regime of dim satellites by treating the faint-end slope of the satellite luminosity function, α, and the lognormal scatter in luminosity at fixed ${V}_{\mathrm{peak}},{\sigma }_{M}$ as free parameters. We assume that this scatter is lognormal and constant as a function of halo properties in our fiducial analysis; we explore a mass-dependent scatter model in Appendix A.2.

Our abundance-matching model is a simple, empirical prescription for assigning satellite luminosities that is not designed to capture the complexities of star formation in ultrafaint dwarf galaxies. For example, Bose et al. (2018) argued that star formation in systems dimmer than ${M}_{V}\sim -5\ \mathrm{mag}$ is effectively shut down by reionization, resulting in two distinct galaxy populations today. While our abundance-matching model is consistent with the current data, which are fit fairly well by a single power-law luminosity function (see Paper I), it will be valuable to investigate more detailed models of stellar mass growth and compare against a wider range of observables, including the inferred star formation histories of MW satellites, in future work.

We assign physical sizes to satellites by extrapolating a modified version of the size–virial radius relation from Kravtsov (2013), which links a galaxy's stellar 3D half-mass radius to its halo's virial radius, into the faint satellite regime. In particular, we set the mean predicted size of each satellite at accretion according to

$$\begin{eqnarray}&&{r}_{1/2}\equiv { \mathcal A }{\left(\displaystyle \frac{{R}_{\mathrm{vir}}}{{R}_{0}}\right)}^{n},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal A }$ and n are free parameters, R_vir denotes the subhalo virial radius measured at accretion, and ${R}_{0}=10\ \mathrm{kpc}$ is a normalization constant. Following Nadler et al. (2019b), we equate the 3D half-mass radii predicted by Equation (1) to azimuthally averaged projected half-light radii; this conversion neglects mass-to-light weighting and projection effects. Nonetheless, this size relation yields reasonable mean sizes when compared to the observed population of classical and SDSS-discovered satellites (Nadler et al. 2019b).

We draw the size of each satellite from a lognormal distribution with a mean given by Equation (1) and a standard deviation of ${\sigma }_{\mathrm{log}R}$, which is a free parameter in our model. When fitting the observed satellite populations, we only compare predicted and mock satellites with ${r}_{1/2}\gt 10\ \mathrm{pc}$ in order to exclude likely star clusters from the analysis. We explore a more conservative cut of ${r}_{1/2}\gt 20\ \mathrm{pc}$ in Appendix C.2.

The size prescription described above assumes that satellite sizes are fixed after accretion onto the MW. However, post-infall effects such as tidal stripping and tidal heating can shrink or enlarge satellites depending on their orbital histories (Peñarrubia et al. 2009; Errani et al. 2015; Fattahi et al. 2018). In Appendix A.4, we show that our key results are not sensitive to these effects using a simple model for satellite size evolution due to tidal stripping.

To incorporate the effects of baryonic physics—and particularly the tidal influence of the Galactic disk—on our simulated subhalo populations, we apply a random forest algorithm trained on hydrodynamic simulations to predict the probability that each subhalo will be disrupted in a hydrodynamic resimulation based on its orbital and internal properties (Garrison-Kimmel et al. 2017b; Nadler et al. 2018). We model the strength of this disruption effect using the free parameter ${ \mathcal B }$, which is defined such that ${ \mathcal B }=1$ corresponds to fiducial hydrodynamic predictions (Nadler et al. 2018) and larger (smaller) values of ${ \mathcal B }$ correspond to more effective (less effective) subhalo disruption. For each subhalo, we set

$$\begin{eqnarray}&&{p}_{\mathrm{disrupt}}\equiv {\left({p}_{\mathrm{disrupt},0}\right)}^{1/{ \mathcal B }},\end{eqnarray} \tag{ 2 }$$

where ${p}_{\mathrm{disrupt},0}$ is the fiducial disruption probability returned by the machine-learning algorithm in Nadler et al. (2018).

The stochastic, nonlinear nature of galaxy formation in low-mass halos likely leads to a smoothly varying fraction of occupied halos, rather than a sharp cutoff in the efficiency of galaxy formation (Sawala et al. 2016b; Fitts et al. 2018; Munshi et al. 2019; Wheeler et al. 2019). Thus, in our fiducial model, we parameterize the fraction of halos that host galaxies of any mass, referred to as the galaxy occupation fraction, following Graus et al. (2019),

$$\begin{eqnarray}&&{f}_{\mathrm{gal}}({{ \mathcal M }}_{\mathrm{peak}})\equiv \displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(\displaystyle \frac{{{ \mathcal M }}_{\mathrm{peak}}-{{ \mathcal M }}_{50}}{\sqrt{2}{\sigma }_{\mathrm{gal}}}\right)\right],\end{eqnarray} \tag{ 3 }$$

where ${{ \mathcal M }}_{\mathrm{peak}}$ is the largest virial mass a subhalo ever attains, which typically occurs before infall into the MW; ${{ \mathcal M }}_{50}$ is the peak halo mass at which 50% of halos host galaxies of any mass; and σ_gal is the width of the galaxy occupation fraction. In our fiducial model, ${{ \mathcal M }}_{50}$ and ${\sigma }_{\mathrm{gal}}$ are free parameters. Note that in the limit ${\sigma }_{\mathrm{gal}}\to 0$, this reduces to a model in which all halos with ${{ \mathcal M }}_{\mathrm{peak}}\gt {{ \mathcal M }}_{50}$ host a galaxy.

Although our analysis does not constrain σ_gal, Equation (3) is a simple, physically motivated form of the occupation fraction that will be interesting to explore in future work. Note that we parameterize the occupation fraction in terms of peak halo mass (rather than, e.g., V_peak) because ${{ \mathcal M }}_{\mathrm{peak}}$ is more easily interpretable and connects directly to constraints on alternative dark matter models (e.g., Nadler et al. 2019a).

Because we model faint galaxies that can potentially inhabit subhalos near the resolution limit of our simulations, it is important to account for artificially disrupted subhalos that might host "orphan" satellites (Guo et al. 2011; see Bose et al. 2019 for a recent example of the importance of orphans in satellite modeling). To model orphans, we follow the prescription in Nadler et al. (2019b), which identifies disrupted subhalos in each simulation, interpolates their orbits to z = 0 using a softened gravitational force law and a dynamical friction model, and accounts for tidal stripping with a mass-loss model calibrated on high-resolution test simulations. We parameterize the effective abundance of orphans by setting their disruption probabilities equal to

$$\begin{eqnarray}&&{p}_{\mathrm{disrupt}}\equiv {\left(1-{a}_{\mathrm{acc}}\right)}^{{ \mathcal O }},\end{eqnarray} \tag{ 4 }$$

where a_acc is the final scale factor at which each subhalo enters the virial radius of the MW, and ${ \mathcal O }$ is a parameter that captures deviations from disruption probabilities in hydrodynamic simulations, which are well fit by ${ \mathcal O }=1$ (Nadler et al. 2019b). Thus, larger (smaller) values of ${ \mathcal O }$ correspond to a greater (smaller) contribution from orphan satellites.

Following Nadler et al. (2019b), we include orphan satellites by fixing ${ \mathcal O }=1$ in our fiducial model. Thus, we effectively assume that subhalo disruption in dark matter–only simulations is a numerical artifact (van den Bosch & Ogiya 2018; van den Bosch et al. 2018) but that subhalo disruption in hydrodynamic simulations is a physical effect. We show that our results are largely insensitive to the value of ${ \mathcal O }$ in Appendix A.6.

Figure 2 shows the 68% and 95% confidence intervals for the observed DES and PS1 satellite luminosity functions given by draws from the posterior of our fiducial model, which is consistent with both data sets. We note that the DES and PS1 likelihoods individually yield consistent results.

Zoom In Zoom Out Reset image size

Figure 3 shows the corresponding satellite size distributions drawn from our fiducial posterior. Our model is consistent with the sizes of both the observed DES and PS1 satellites. It very slightly overpredicts the sizes of observed DES systems; however, we reiterate that our size model does not allow for size reduction due to tidal stripping or size enlargement due to tidal heating, which affect satellites with close pericentric passages to the Galactic disk (e.g., Amorisco 2019). Our findings in Appendix A.4 suggest that the post-infall size evolution of satellites in subhalos with ${V}_{\mathrm{peak}}\gt 10\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${V}_{\max }\gt 9\,\mathrm{km}\,{{\rm{s}}}^{-1}$ does not significantly affect our inference.

Zoom In Zoom Out Reset image size

Our fiducial model is consistent with the outer radial distributions of both DES and PS1 satellites, but it slightly underpredicts the number of satellites near the center of the MW ($D\lesssim 100\ \mathrm{kpc}$), particularly in PS1. We explore this minor discrepancy in Appendix A.3, where we show that our galaxy–halo connection model constraints and inferred total MW satellite population are largely unaffected if the radial distribution is forced to match the data.

To assess the impact of the LMC and its satellites on the MW satellite population, we test the following models in addition to our fiducial model, which includes a realistic, recently accreted LMC system by construction.

• (i)  
  No LMC: a model with four host halos that have the same mean concentration as our fiducial hosts but no LMC analog.
• (ii)  
  Misplaced LMC: a model with our fiducial host halos in which subhalo positions are reflected, effectively placing the DES footprint in the northern hemisphere.
• (iii)  
  Early LMC Infall: a model with two host halos that have the same mean concentration as our fiducial hosts with LMC analogs that pass our LMC V_max and heliocentric distance cuts but fall into the MW 2 and 6 Gyr ago, respectively.

For each alternative LMC scenario listed above, we refit the observed DES and PS1 satellite populations, sampling over the same eight parameters used in our fiducial analysis.

Our fiducial model is favored over the No LMC, Misplaced LMC, and Early LMC Infall scenarios with Bayes factors of ∼10⁴, 10⁴, and 10³, respectively. In addition, both host halos in the Early LMC Infall case are individually disfavored with Bayes factors of $\sim {10}^{3}$. Thus, we find decisive statistical evidence for the impact of the LMC on the MW satellite population, particularly within and near the DES footprint. Moreover, we infer that the LMC system fell into the MW within the last 2 Gyr at high confidence. We also note that, in our fiducial host with more massive MW and LMC halos, the LMC reaches pericenter near the second-to-last simulation snapshot (i.e., ∼150 Myr ago). Performing our analysis using the final snapshot for this host noticeably degrades the fit due to the dispersal and disruption of LMC satellites during the LMC's pericentric passage. Thus, we use the second-to-last snapshot for this host in our fiducial analysis, and we remark that satellite abundances can potentially constrain the number of allowed pericentric passages for the LMC.

The alternative LMC scenarios defined above are strongly disfavored relative to our fiducial model because they cannot produce a sufficient number of dim satellites in the DES footprint without overpredicting the number of observed PS1 satellites. This is a direct consequence of the spatial overdensity of subhalos near the LMC analogs in our fiducial simulations; in particular, the projected density of resolved subhalos within 50^◦ of the LMC on the sky is enhanced by ∼50% relative to the density on a random patch of sky.

To quantify the number of satellites in our fiducial model that are associated with the LMC, we explore the following definitions of LMC-associated subhalos.

• (i)  
  Fiducial definition. A subhalo is associated with the LMC if it is within the virial radius of the LMC halo at the time of LMC infall into the MW.
• (ii)  
  Gravitationally influenced definition. A subhalo is associated with the LMC if it has ever passed within the virial radius of the LMC halo.

Here LMC infall is defined as the snapshot at which the center of the LMC halo crosses within the MW virial radius. Note that nearly all systems that satisfy our strict definition are bound to the LMC at the time of LMC infall.

Under the fiducial (gravitationally influenced) definitions above, we predict that 52 ± 8 (181 ± 25) total LMC-associated subhalos (above our cuts of ${V}_{\mathrm{peak}}\gt 10\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${V}_{\max }\gt 9\,\mathrm{km}\,{{\rm{s}}}^{-1}$) exist within the virial radius of the MW today, where the 95% confidence interval is estimated by drawing from our fiducial posterior. We predict that 48 ± 8 (164 ± 25) of these subhalos form galaxies with ${M}_{V}\lt 0\ \mathrm{mag}$ and ${r}_{1/2}\gt 10\ \mathrm{pc}$ (in agreement with an earlier estimate by Jethwa et al. 2018), and that 41 ± 7 (118 ± 21) of these satellites survive tidal disruption due to the Galactic disk. Of these surviving LMC-associated satellites, we predict that 4.8 ± 1.7 (11 ± 3.6) are currently observed by DES and 1.1 ± 0.9 (6.1 ± 2.1) are currently observed by PS1.

Our statistical probe of LMC satellite association is remarkably consistent with the number of observed LMC satellites inferred from Gaia proper-motion measurements, which indicate that four galaxies in or near the DES footprint—excluding the Small Magellanic Cloud (SMC)—are associated with the LMC, and that two satellites in or near the PS1 footprint are potentially associated with the LMC (Kallivayalil et al. 2018; Patel et al. 2020).⁵⁹ In addition, the orbital dynamics of our predicted LMC satellites are consistent with Gaia proper-motion measurements for these likely LMC-associated systems. These predictions are also consistent with other empirical models (Deason et al. 2015; Jethwa et al. 2016; Dooley et al. 2017; Sales et al. 2017; Kallivayalil et al. 2018; Erkal & Belokurov 2019; Zhang et al. 2019) and with hydrodynamic simulations of isolated LMC analogs (Jahn et al. 2019).

In Appendix D, we show that the properties of our LMC-like systems are not significantly affected by the realizations of small-scale power in our fiducial simulations. However, we caution that the number of predicted LMC satellites observed by DES and PS1 depends on the particular properties of our two LMC analogs. Thus, exploring the robustness of these results using a suite of zoom-in simulations selected to contain realistic LMC systems with a range of internal and orbital properties is an important avenue for future work.

Figure 4 shows the total MW satellite luminosity function and surface brightness distribution resulting from our fit to the DES and PS1 satellite populations. We infer that a total of 220 ± 50 satellites with M_V < 0 mag and r_1/2 > 10 pc exist within the virial radius of the MW, where uncertainties correspond to 68% confidence intervals calculated by sampling from our fiducial posterior. Thus, we predict that ∼150 satellites remain undiscovered in a standard CDM scenario, roughly one-fourth of which are associated with the LMC. This is larger than the fraction of satellites that have ever fallen into the MW that are associated with the LMC because our fiducial LMC analogs accreted recently, making their satellites less likely to be disrupted. Our prediction for the total number of MW satellites is consistent with several recent studies (Jethwa et al. 2018; Kim et al. 2018; Newton et al. 2018; Nadler et al. 2019b), and it is lower than the empirical estimate in Paper I, which was recognized to be inflated due to the assumption of an isotropic satellite distribution. This prediction will be tested by upcoming deep imaging surveys; indeed, HSC-SSP has already started to probe this population of distant, low surface brightness MW satellites by discovering three new ultrafaint satellite candidates in ∼676 deg² of imaging data (Homma et al. 2016, 2018, 2019).

Zoom In Zoom Out Reset image size

To estimate whether our predictions are consistent with HSC-SSP observations, we draw realizations of the MW satellite population from our fiducial posterior and calculate the number of systems within the DES or PS1 footprints that would not be observed by the appropriate survey. We then estimate the number of these systems currently observed by an HSC-like survey covering 676 deg² that detects all satellites (i.e., systems with M_V < 0 mag and r_1/2 > 10 pc) down to a surface brightness of ${\mu }_{V}=32\ \mathrm{mag}\ {\mathrm{arcsec}}^{-2}$ and out to a heliocentric distance of 300 kpc, assuming an isotropic satellite distribution at high Galactic latitudes and accounting for subhalo disruption. There are six known satellites in the HSC footprint, but two of the six (Sextans and Leo IV) are detected at high significance in PS1 by at least one of the search algorithms in Paper I. We find that our mock HSC survey detects 1.75 ± 0.6 satellites, which is in slight tension with the number of systems detected by HSC (four, after discounting Sextans and Leo IV).

Figure 4 illustrates several predictions from hydrodynamic simulations of isolated and satellite dwarf galaxies. Our results are largely consistent with the luminosity function of bright MW satellites in hydrodynamic simulations of the Local Group using the Feedback In Realistic Environments (FIRE) feedback prescription, down to the FIRE resolution limit of ∼−6 mag (Garrison-Kimmel et al. 2019). Note that these FIRE simulations do not include LMC or SMC analogs, which accounts for the discrepancy with both our predictions and the observed luminosity function at ${M}_{V}\lt -16\ \mathrm{mag}$. Interestingly, other recent hydrodynamic simulations indicate that different star formation prescriptions significantly impact the amplitude and faint-end slope of the luminosity function for satellites of isolated LMC-like halos (Munshi et al. 2019). Thus, our constraints on the faint-end slope, which are driven by satellites with ${M}_{V}\gtrsim -6\ \mathrm{mag}$ (corresponding to stellar mass ${M}_{* }\,\lesssim {10}^{5}\ {M}_{\odot }$), can be used to inform subgrid star formation prescriptions.

The posterior distribution for our fiducial model is shown in Figure 5, and the corresponding galaxy–halo connection model constraints are listed in Table 1. Note that we obtain statistically consistent results when fitting the DES and PS1 satellite populations with either of our two fiducial simulations individually in terms of both the Bayesian evidence and the galaxy–halo connection model constraints. In particular, no parameter constraints shift by more than 1σ relative to the fiducial values reported below when the fit is performed using either simulation individually. We now discuss each constraint in detail.

• 1.  
  The inferred faint-end slope of the satellite luminosity function is steeper than that reported in a previous study based on classical and SDSS satellites (Nadler et al. 2019b). Our constraint is consistent with the faint-end slope derived from higher-luminosity field galaxies in the GAMA survey (Loveday et al. 2015; Wright et al. 2017), even though it is based on a sample that extends nearly 10 mag fainter than that used in GAMA. We note that α is the most sensitive parameter in our analysis to modeling assumptions and details of the observed satellite population.
• 2.  
  The scatter in luminosity at fixed V_peak is constrained to ${\sigma }_{M}\lt 0.19\ \mathrm{dex}$ at 95% confidence, which may inform hydrodynamic feedback prescriptions that predict a steep increase in luminosity scatter at low masses (e.g., see Wechsler & Tinker 2018). Our lack of a lower limit on σ_M is consistent with previous studies of faint galaxy samples (e.g., Lehmann et al. 2017). Meanwhile, large values of σ_M are not allowed because too many low-V_peak satellites upscatter to observable luminosities, resulting in overpredicted luminosity functions. To confirm that this upper limit is robust, we calculate Bayes factors by drawing samples from the posterior in bins of σ_M, finding that σ_M = 0.15 dex (σ_M = 0.2 dex) is disfavored relative to σ_M = 0 dex with a Bayes factor of 30 (100).⁶⁰ These upper limits are comparable to the scatter typically inferred from abundance-matching analyses of brighter systems (σ_M ∼ 0.2 dex) and smaller than that from hydrodynamic simulations of dwarf galaxies (e.g., Rey et al. 2019); however, we caution that our constraint might be impacted by the use of only two independent realizations of the MW satellite population. In addition, it is potentially misleading to compare global constraints on scatter to those derived from the MW alone. Both of these caveats are important to explore in future work.
• 3.  
  The peak mass at which 50% of halos host galaxies is inferred to be less than $8.5\times {10}^{7}\ {M}_{\odot }$ at 95% confidence. Note that this summary statistic depends on the lower limit of our prior on ${{ \mathcal M }}_{50}$, since the ${{ \mathcal M }}_{50}$ posterior flattens near its lower limit, which is chosen based on the resolution of our simulations. Thus, we also calculate Bayes factors by drawing from the posterior in bins of ${{ \mathcal M }}_{50}$ to confirm that this summary statistic is robust. We find that ${{ \mathcal M }}_{50}=8.5\times {10}^{7}\ {M}_{\odot }$ (${{ \mathcal M }}_{50}=1.5\times {10}^{8}\ {M}_{\odot }$) is disfavored relative to arbitrarily low values of ${{ \mathcal M }}_{50}$ with a Bayes factor of 50 (100). The current data are not able to place a lower limit on ${{ \mathcal M }}_{50}$, which would correspond to the detection of a cutoff in galaxy formation.
• 4.  
  Our posterior is consistent with ${ \mathcal B }=1$, corresponding to our fiducial baryonic disruption model. Although a large spread in disruption strength is allowed by the data, extremely efficient (${ \mathcal B }\gt 2.1$) and inefficient (${ \mathcal B }\lt 0.3$) subhalo disruption relative to hydrodynamic simulations is strongly disfavored. These constraints widen when our lognormal prior on ${ \mathcal B }$ is relaxed; however, zero subhalo disruption (corresponding to ${ \mathcal B }=0$) is robustly ruled out.
• 5.  
  The scatter in the galaxy occupation fraction is consistent with zero, which makes sense given our lack of a lower limit on ${{ \mathcal M }}_{50}$. Models with large scatter in the occupation fraction (${\sigma }_{\mathrm{gal}}\lt 0.67\ \mathrm{dex}$), corresponding to extremely stochastic galaxy formation, are disfavored relative to a step function occupation fraction at 95% confidence. Note that the slope of the σ_gal posterior is driven by the lower limit of our ${{ \mathcal M }}_{50}$ prior; as this limit decreases, the σ_gal posterior flattens.
• 6.  
  The amplitude of the galaxy–halo size relation, defined as the typical size of a satellite in a halo with ${R}_{\mathrm{vir}}=10\ \mathrm{kpc}$ at accretion, is constrained to lie between 12 and 73 pc at 68% confidence. For larger values of ${ \mathcal A }$, satellites are too large to be detected with high probability, and the DES and PS1 luminosity functions are underpredicted; for smaller values of ${ \mathcal A }$, many predicted satellites do not pass our r_1/2 > 10 pc cut, and the luminosity functions are underpredicted.
• 7.  
  The scatter in the galaxy–halo size relation is constrained to lie between 0.33 and 0.91 dex at 68% confidence. For larger values of ${\sigma }_{\mathrm{log}R}$, faint satellites upscatter to large sizes too frequently, which results in underpredicted luminosity functions. Our 68% confidence lower limit on ${\sigma }_{\mathrm{log}R}$ of 0.33 dex is consistent with the value estimated in Kravtsov (2013). Lower values of ${\sigma }_{\mathrm{log}R}$ lead to slightly too many predicted DES and PS1 satellites; however, our results are consistent with ${\sigma }_{\mathrm{log}R}=0\ \mathrm{dex}$ at 95% confidence.
• 8.  
  The power-law index of the galaxy–halo size relation is constrained to lie between 0.5 and 1.45 at 68% confidence. For shallower power-law slopes, satellite sizes do not change appreciably as a function of halo size, which results in a worse joint fit to the observed absolute magnitude and surface brightness distribution. We note that the posterior widens when our Gaussian prior on n is relaxed.

Zoom In Zoom Out Reset image size

Table 1.  Galaxy–Halo Connection Model Constraints Derived from Our Fit to the DES and PS1 Satellite Populations

Parameter	Physical Interpretation	95% Confidence Interval
Faint-end slope	Power-law slope of satellite luminosity function	−1.46 < α < −1.39
Luminosity scatter	Scatter in luminosity at fixed Vpeak	$0\ {\mathrm{dex}}^{* }\lt {\sigma }_{M}\lt 0.19\ \mathrm{dex}$
50% occupation mass	Mass at which 50% of halos host galaxies	${7.5}^{* }\lt \mathrm{log}({{ \mathcal M }}_{50}/{M}_{\odot })\lt 7.93$
Baryonic effects	Strength of subhalo disruption due to baryons	$0.3\lt { \mathcal B }\lt 2.1$
Occupation scatter	Scatter in galaxy occupation fraction	$0\ {\mathrm{dex}}^{* }\lt {\sigma }_{\mathrm{gal}}\lt 0.67\ \mathrm{dex}$
Size amplitude	Amplitude of galaxy–halo size relation	${0}^{* }\ \mathrm{pc}\lt { \mathcal A }\lt 110\ \mathrm{pc}$
Size scatter	Scatter in half-light radius at fixed halo size	$0\ {\mathrm{dex}}^{* }\lt {\sigma }_{\mathrm{log}R}\lt 1.2\ \mathrm{dex}$
Size power-law index	Power-law index of galaxy–halo size relation	${0}^{* }\lt n\lt 1.8$

Note. Asterisks mark prior-driven constraints.

Download table as:  ASCIITypeset image

We now explore the properties of the lowest-mass halos inferred to host MW satellites. The left panel of Figure 6 shows the galaxy occupation fraction derived from our statistical inference, where uncertainties are estimated by drawing from our fiducial posterior. By sampling from our fiducial posterior, we infer that halos with a peak virial mass below $2.5\times {10}^{8}\ {M}_{\odot }$ and peak circular velocity below $19\,\mathrm{km}\,{{\rm{s}}}^{-1}$ host at least one of the faintest observed satellites. To convert these into maximally conservative upper limits, we account for the uncertainty in MW host halo mass using the procedure described in Appendix A.1, which yields limits on the minimum halo mass and peak circular velocity of ${{ \mathcal M }}_{\min }\lt 3.2\times {10}^{8}\ {M}_{\odot }$ and ${V}_{\mathrm{peak},\min }\lt 21\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at 95% confidence. Furthermore, we predict that the faintest observed satellite inhabits a halo with ${{ \mathcal M }}_{\mathrm{peak}}=1.5\times {10}^{8}\ {M}_{\odot }$, on average.⁶¹

Zoom In Zoom Out Reset image size

These results improve upon the minimum halo mass constraint derived from classical and SDSS satellites (Nadler et al. 2019b) by a factor of 2, and they are consistent with the constraints reported in Jethwa et al. (2018). Moreover, these upper limits are not in significant tension with the expected atomic cooling limit of ${V}_{\mathrm{peak}}\approx 20\,\mathrm{km}\,{{\rm{s}}}^{-1}$, contrary to recent studies based on the radial MW satellite distribution (e.g., Graus et al. 2019) and consistent with the findings in Bose et al. (2019).

We caution that the median galaxy occupation fraction shown in Figure 6 is driven by the assumed functional form in Equation (3) and is therefore arbitrary. Although the functional form in Equation (3) is consistent with results from hydrodynamic simulations for ${{ \mathcal M }}_{\mathrm{peak}}\gtrsim {10}^{9}\ {M}_{\odot }$, this particular functional form is not required to fit the DES and PS1 luminosity functions. Rather, we have evidence that f_gal > 50% above a peak virial mass of $\sim {10}^{8}\ {M}_{\odot }$. To verify that the assumed form of the galaxy occupation fraction does not impact our constraints, we also test a binned model in which we fit for ${{ \mathcal M }}_{50}$ and a corresponding 90% occupation mass. We find that the resulting 50% and 90% occupation constraints are consistent with those inferred from our fiducial analysis.

A wide range of galaxy occupation fractions have been reported in hydrodynamic simulations, with some placing ${{ \mathcal M }}_{50}$ as high as $\sim {10}^{9}\ {M}_{\odot }$ (Sawala et al. 2016b; Fitts et al. 2018). However, recent hydrodynamic simulations run at higher resolution result in efficient galaxy formation in significantly lower-mass halos, and some claim that all halos down to the resolution limit consistently host star particles (Wheeler et al. 2019). In addition, simulations of galaxy formation at early pre-reionization epochs show that stellar systems form in halos with masses as low as $\sim {10}^{7}\ {M}_{\odot }$ (e.g., see Figure 13 in Côté et al. 2018 for a compilation of recent simulation results). Most recently, high-resolution simulations of high-redshift galaxy formation that include the effects of spatially and temporally inhomogeneous reionization find ${{ \mathcal M }}_{50}\sim {10}^{8}\ {M}_{\odot }$ (Katz et al. 2019).

Our galaxy occupation fraction constraint implies that models with ${{ \mathcal M }}_{50}\gt {10}^{8}\ {M}_{\odot }$ are in significant tension with the observed MW satellite population, as long as MW satellite formation is representative of galaxy formation at this halo mass scale, on average. This assumption may not be true if the reionization history of the MW's Lagrangian volume differs from the average reionization history of an MW-mass halo hosting dwarf galaxies of the masses considered here (however, see Alvarez et al. 2009; Busha et al. 2010). Note that analyses based on H i surveys of Local Group dwarfs indicate a suppression mass scale similar to our ${{ \mathcal M }}_{50}$ constraint (Tollerud & Peek 2018).

Due to our abundance-matching assumption, the lowest-mass halos in our model host the faintest galaxies, on average. Thus, our constraints on the masses of these halos are conservative, since the most massive halos in our simulations are forced to host more easily observable satellites at fixed distance and size, modulo baryonic disruption effects and abundance-matching scatter. In other words, our abundance-matching model yields a testable prediction: the faintest galaxies should inhabit the halos with the lowest pre-infall virial masses. We expect this correlation to be weakened by post-infall effects, including tidal stripping, but we can nevertheless infer the present-day joint distribution of halo mass and satellite luminosity or stellar mass. We illustrate this stellar mass–halo mass (SMHM) relation in the right panel of Figure 6. Our inferred SMHM relation is generally consistent with recent results (e.g., Jethwa et al. 2018). Like the faint-end slope of the luminosity function, the SMHM relation can be used to discriminate between different subgrid models of star formation and stellar feedback (Munshi et al. 2019). As in previous studies, we find that the SMHM relation in the ultrafaint regime falls off more steeply than extrapolations of abundance-matching relations derived using higher-mass field galaxies (Behroozi et al. 2013a). Interestingly, Agertz et al. (2020) found that full on-the-fly radiative transfer is necessary to match the steepness and normalization of our inferred SMHM relation for a fixed hydrodynamic feedback prescription.

Ultimately, our predictions must be confronted with the dynamical mass function of observed satellites, measurements of which will improve significantly in the era of upcoming spectroscopic facilities and giant segmented mirror telescopes (Simon et al. 2019). A preliminary comparison of our joint predicted distribution of luminosity and V_max with the measured stellar velocity dispersions of DES and PS1 satellites suggests that our model is consistent with the inferred central densities of low-luminosity satellites (${M}_{V}\gt -6\ \mathrm{mag}$). Although there is a systematic discrepancy between observed and predicted values of V_max for brighter systems (the "too big to fail" problem), our simple comparison does not account for the conversion from line-of-sight velocity dispersion measured within observed half-light radii to V_max or the tidal effects of the Galactic disk on the density profiles of surviving subhalos. Moreover, the systems for which predicted densities are higher than those inferred observationally are susceptible to baryonic feedback processes that core the inner regions of halos (Di Cintio et al. 2014), and this effect has been shown to alleviate the too big to fail problem (Brooks et al. 2013; Sawala et al. 2016a; Wetzel et al. 2016; Lovell et al. 2017; Garrison-Kimmel et al. 2019).

Finally, we explore the properties of the halos inferred to host the faintest potentially detectable galaxies. In particular, we calculate the minimum peak halo mass necessary for halos to contain a stellar population of at least $100\ {M}_{\odot }$, chosen to represent the approximate threshold for which it would be possible to observationally confirm a stellar system as a dark matter–dominated dwarf galaxy.⁶² By populating a higher-resolution version of one of our fiducial simulations and sampling from the posterior of our abundance-matching relation, we find that systems at the observational threshold occupy halos with ${{ \mathcal M }}_{\mathrm{peak}}\gt {10}^{6}\ {M}_{\odot }$ at 95% confidence. To detect even lower-mass halos, gravitational probes of dark matter that are independent of baryonic content, e.g., gravitational lensing or stellar streams, must be employed.

Many deviations from CDM lead to a cutoff in the abundance of low-mass halos. Several authors have used MW satellite abundances to constrain a free-streaming cutoff induced by warm dark matter (e.g., Macciò & Fontanot 2010; Kennedy et al. 2014; Lovell et al. 2014; Jethwa et al. 2018). Nadler et al. (2019a) showed that similar constraints apply to other dark matter models, resulting in limits on the velocity-independent scattering cross section between dark matter and baryons. Our statistical detection of halos with peak virial masses below $3.2\times {10}^{8}\ {M}_{\odot }$ therefore translates directly into constraints on various microphysical properties of dark matter.

Nadler et al. (2019a) found that the minimum halo mass inferred in this fashion is comparable to the limit on the half-mode mass M_hm, which corresponds to the scale at which the matter power spectrum is suppressed by a critical amount relative to CDM due to dark matter free streaming or interactions. Performing a statistical inference in which the half-mode mass is varied will constrain it to lie below our upper limit on ${{ \mathcal M }}_{50}$, since the abundance of halos at and above the half-mode mass is reduced relative to CDM.

Thus, for a simple and conservative estimate of the dark matter constraints resulting from our analysis, we set the upper limit on M_hm equal to our upper limit on the minimum halo mass, i.e., ${M}_{\mathrm{hm}}\lt 3.2\times {10}^{8}\ {M}_{\odot }$, which corresponds to a lower limit on the half-mode scale of ${k}_{\mathrm{hm}}\gt 36\ h\ {\mathrm{Mpc}}^{-1}$. Using the relations in Nadler et al. (2019a) with the cosmological parameters h = 0.7 and Ω_m = 0.286 corresponding to the simulations used in our analysis (see Section 3.1), this yields a lower limit of 3.4 keV on the mass of thermal relic warm dark matter and an upper limit of 6 × 10⁻³⁰ cm² on the velocity-independent dark matter–baryon scattering cross section for a 10 keV dark matter particle mass, both at 95% confidence. We leave a detailed investigation of dark matter constraints to future work.

The analysis in this paper is restricted to two fiducial MW-like hosts with virial masses of 1.57 and $1.26\times {10}^{12}\ {M}_{\odot }$. However, we expect the uncertainty in the MW host halo mass, M_MW, to impact our constraints. Consider a toy model in which N_sat satellites brighter than a limiting magnitude ${M}_{V,\min }$ must be predicted in order to match the observed luminosity function. In this toy model, the predicted number of satellites is given by

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{sat}}(\lt {M}_{V,\mathrm{lim}}) & = & {\int }_{{M}_{V,\mathrm{lim}}}^{-\infty }\displaystyle \frac{{{dN}}_{\mathrm{sat}}}{{{dM}}_{V}}{{dM}}_{V}\\ & = & {\int }_{{{ \mathcal M }}_{\min }({M}_{V,\mathrm{lim}})}^{\infty }f\left(\displaystyle \frac{{{dN}}_{\mathrm{sub}}}{d{{ \mathcal M }}_{\mathrm{peak}}},{\boldsymbol{\theta }}\right)d{{ \mathcal M }}_{\mathrm{peak}},\end{array}\end{eqnarray} \tag{ A1 }$$

where ${{dN}}_{\mathrm{sub}}/d{{ \mathcal M }}_{\mathrm{peak}}$ is the subhalo mass function, ${{ \mathcal M }}_{\min }$ is the lowest halo mass populated by an observed satellite, and f encapsulates the observational selection, subhalo disruption, and galaxy occupation effects that determine whether each halo hosts an observable satellite, all of which depend on galaxy–halo connection model parameters ${\boldsymbol{\theta }}$. Neglecting the dependence of the latter effects on host mass (which we expect to be subdominant compared to the overall rescaling of subhalo abundances), using the standard linear relationship between subhalo abundance and host mass, and assuming a standard subhalo mass function slope of ${{dN}}_{\mathrm{sub}}/d{{ \mathcal M }}_{\mathrm{peak}}\propto {{ \mathcal M }}_{\mathrm{peak}}^{-2}$ (e.g., Mao et al. 2015), we have

$$\begin{eqnarray}&&{N}_{\mathrm{sat}}\propto {M}_{\mathrm{MW}}{\int }_{{{ \mathcal M }}_{\min }}^{\infty }{{ \mathcal M }}_{\mathrm{peak}}^{-2}\ d{{ \mathcal M }}_{\mathrm{peak}}\propto \displaystyle \frac{{M}_{\mathrm{MW}}}{{{ \mathcal M }}_{\min }}.\end{eqnarray} \tag{ A2 }$$

Thus, for a fixed observed satellite count N_sat, we expect our 95% confidence level upper limit on ${{ \mathcal M }}_{\min }$ to scale linearly with host mass. In addition, because the error on MW mass is independent of the error on ${{ \mathcal M }}_{\min }$, we expect these uncertainties to add in quadrature.

Given our fiducial minimum halo mass of $2.5\times {10}^{8}\ {M}_{\odot }$ derived for an average host mass of $1.4\times {10}^{12}\ {M}_{\odot }$, we therefore expect ${{ \mathcal M }}_{\min }\lt 3.2\times {10}^{8}\ {M}_{\odot }$ (${{ \mathcal M }}_{\min }\lt 2\times {10}^{8}\ {M}_{\odot }$) for a maximally high-mass (maximally low-mass) host halo given the current 2σ observational uncertainty on the MW virial mass of $1.0\times {10}^{12}\lt {M}_{\mathrm{MW}}/{M}_{\odot }\lt 1.8\times {10}^{12}$ (Callingham et al. 2019; Cautun et al. 2019b; Li et al. 2019a, 2019b). We expect the remaining galaxy–halo connection model parameters and associated errors to remain largely unchanged, although rerunning our analysis using additional simulations is required to confirm this hypothesis. We expect the inferred total satellite count to scale linearly with MW mass; thus, given our fiducial prediction of 220 total MW satellites with M_V < 0 mag and r_1/2 > 10 pc, we expect 280 (170) such satellites for a maximally high-mass (maximally low-mass) host halo.

Above, we implicitly assumed that our galaxy–halo connection model is capable of adjusting the number of faint satellites, which make the largest contribution to ${N}_{\mathrm{sat}}$, while simultaneously matching the bright end of the observed satellite luminosity functions. This assumption holds because we have fixed the abundance-matching prescription to the relation derived from GAMA data for ${M}_{V}\lt -13\ \mathrm{mag}$ while allowing the faint-end slope to vary. We note that the results of Newton et al. (2018) suggest that the inferred number of satellites within a fixed physical radius is independent of M_MW. We find that the total number of satellites inferred within the virial radius scales almost exactly linearly with M_MW, as expected from the linear scaling of subhalo abundance with host halo mass, and further confirming the consistency of the results among our two fiducial simulations. In addition, we find that ${{ \mathcal M }}_{\min }$ is roughly independent of M_MW for our fiducial simulations.

Here we test a model where the abundance-matching scatter in luminosity at fixed V_peak, σ_M, depends on peak halo mass. Motivated by the model in Garrison-Kimmel et al. (2017a), we set

$$\begin{eqnarray}&&{\sigma }_{M}\equiv {\sigma }_{M,0}-{\gamma }_{M}(\mathrm{log}{{ \mathcal M }}_{\mathrm{peak}}-\mathrm{log}{{ \mathcal M }}_{1}),\end{eqnarray} \tag{ A3 }$$

where ${\sigma }_{M,0}$ is a free parameter that captures the amplitude of the luminosity scatter, γ_M is a free parameter that captures its mass dependence, and ${{ \mathcal M }}_{1}={10}^{11}\ {M}_{\odot }$ is fixed. By rerunning our fit with γ_M as an additional ninth free parameter, we find that large values of γ_M are ruled out by the DES and PS1 satellite populations at high statistical significance, with γ_M < 0.07 at 95% confidence. Large values of γ_M are disfavored because abundant, low-mass halos host satellites that upscatter to observable luminosities too often to match the observed DES and PS1 luminosity functions; however, the same caveats noted in Section 7.4 for our constraint on σ_M apply to γ_M, so this upper limit should be interpreted with caution. Introducing mass-dependent scatter does not significantly affect our inferred upper bound on ${{ \mathcal M }}_{50}$, implying that our fiducial minimum halo mass constraint does not depend on the details of our luminosity scatter model.

To account for potential biases in our radial subhalo distributions due to artificial disruption and halo finder incompleteness, we define the parameter χ by

$$\begin{eqnarray}&&{r}_{\mathrm{sat}}\equiv \chi {r}_{\mathrm{sub}},\end{eqnarray} \tag{ A4 }$$

where ${r}_{\mathrm{sat}}$ is a satellite's distance from the center of its host halo, which we equate to its galactocentric distance, and r_sub is the galactocentric distance of the corresponding subhalo.

In our main analysis, we take subhalo positions directly from the simulation data and therefore assume χ = 1. However, as noted above, our fiducial model slightly underpredicts the observed radial distribution of satellites close to the center of the MW in the PS1 footprint. We plot the predicted DES and PS1 radial distributions for our fiducial model in Figure A1; to illustrate the effect of varying χ, we also show the 68% confidence interval for our fiducial posterior evaluated with χ = 0.5.

Zoom In Zoom Out Reset image size

To test the impact of radial scaling, we refit the DES and PS1 satellite populations with χ as an additional ninth free parameter. As expected, decreasing χ reduces the tension between the predicted and observed inner radial distribution of PS1 satellites; however, doing so does not significantly affect the goodness of fit for the observed luminosity functions and size distributions. Moreover, our key constraints, including the upper limit on ${{ \mathcal M }}_{50}$, and our conclusions regarding the impact of the LMC system are not affected. In particular, the Bayes factors in favor of our fiducial LMC model relative to the alternative LMC scenarios defined in Section 7.2 are unchanged. We note that, since we have only fit to observed absolute magnitudes and surface brightnesses, the discrepancy with the observed radial distribution for χ = 1 might not persist for a fit that includes galactocentric distance; we comment on the technical difficulties associated with such a fit in Appendix B.4.

Bose et al. (2019) suggested that the Gaia-Enceladus accretion event, in which an LMC-mass galaxy merged with the MW 8–11 Gyr ago, might lead to a relative overabundance of ultrafaint satellites in the MW. Because of dynamical friction, this overabundance would be particularly evident in the innermost regions of the MW and might affect the observed radial satellite distribution. Interestingly, both host halos used in this work experience a Gaia-Enceladus-like accretion event, following the definition in Bose et al. (2019) of a massive ($\sim {10}^{11}\ {M}_{\odot }$) halo merging with the MW halo between z = 1 and 2. Given that we still predict an underabundance of observed satellites in the inner regions of both fiducial host halos, the Gaia-Enceladus-like events they experience do not seem to be sufficient to ease the tension with the observed radial distribution. Nonetheless, exploring the relationship between the mass accretion history of the MW and the present-day radial distribution of observed ultrafaint satellites in detail is an interesting avenue for future work.

Following Nadler et al. (2019b), we test a model for the evolution of satellite sizes by changing the mean sizes predicted by Equation (1) to

$$\begin{eqnarray}&&{r}_{1/2}^{{\prime} }\equiv {r}_{1/2}\ {\left(\displaystyle \frac{{V}_{\max }}{{V}_{\mathrm{acc}}}\right)}^{\beta },\end{eqnarray} \tag{ A5 }$$

where ${r}_{1/2}^{{\prime} }$ denotes the satellite half-light radius at z = 0, r_1/2 is the half-light radius at accretion predicted by Equation (1), V_max (V_acc) is the maximum circular velocity of a subhalo today (at accretion), and β > 0 is a parameter that controls the strength of size reduction due to tidal stripping. We set β = 0 in our fiducial analysis, meaning that satellite sizes are fixed based on halo sizes at accretion. However, tidal stripping after infall can shrink satellite sizes; for example, Peñarrubia et al. (2009) found that 1 < β < 2 describes the results of high-resolution simulations well.

In Figure A2, we illustrate predicted size distributions for our fiducial posterior evaluated with β = 3; a large value of β was chosen to test an extreme dependence of satellite sizes on tidal stripping. We find that even this extreme model does not impact the observed satellite size distributions, indicating that our results are robust to assumptions about tidal stripping. Our simulations lack the spatial resolution to test whether the Peñarrubia et al. (2009) prescription holds in detail and alters observed satellite size distributions, but this—along with an exploration of size enlargement due to tidal heating—is an interesting avenue for future work.

Zoom In Zoom Out Reset image size

Jiang et al. (2019) found that galaxy sizes in two hydrodynamic simulations follow a size relation similar to that in Kravtsov (2013), with an additional dependence on halo concentration. In particular, the size relation

$$\begin{eqnarray}&&{r}_{1/2}\equiv { \mathcal A }{\left(\displaystyle \frac{c}{10}\right)}^{\gamma }{\left(\displaystyle \frac{{R}_{\mathrm{vir}}}{{R}_{0}}\right)}^{n},\end{eqnarray} \tag{ A6 }$$

with ${ \mathcal A }=0.02,n=1,\gamma =-0.7,{R}_{0}=1\ \mathrm{kpc}$, and halo concentration c measured as a function of redshift, fits the hydrodynamic simulation results in Jiang et al. (2019) with a residual scatter of ∼0.15 dex. This relation implies that more concentrated halos host less extended stellar systems at a fixed virial radius in these simulations.

To test whether a concentration-dependent size model is favored by the DES and PS1 data, we refit these satellite populations with γ as an additional ninth free parameter. Because the concentration of subhalos after infall into the MW is difficult to measure accurately in our simulations, we measure the concentration at the time of accretion when implementing Equation (A6). We find that our galaxy–halo connection model constraints are largely unchanged in this case, although the upper limits on ${\sigma }_{\mathrm{log}R}$ ($0.88\ \mathrm{dex}$) and n (1.7) are more stringent than in our fiducial model. We find that the amplitude of the size relation is degenerate with γ, and our analysis does not place an upper limit on ${ \mathcal A }$ in this case. Here γ itself is constrained to lie between −1.5 and −0.2 at 95% confidence. The predicted luminosity functions and size distributions are nearly identical to those from our fiducial analysis (we illustrate the size distribution for our fiducial posterior evaluated with γ = −0.7 in Figure A2).

To test the importance of orphan satellites, we refit the DES and PS1 satellite populations with ${ \mathcal O }=0$, which adds zero orphans to our fiducial subhalo populations and effectively assumes that there is no artificial subhalo disruption in our simulations. Our constraints are virtually unaffected by this extreme variation in ${ \mathcal O }$. In particular, the 95% confidence upper limit on ${{ \mathcal M }}_{50}$ increases by less than 1σ to ${10}^{8}\ {M}_{\odot }$, and the rest of our galaxy–halo connection model constraints are also not significantly affected. The total number of MW satellites with M_V < 0 mag and r_1/2 > 10 pc decreases to 190 ± 50, as expected from the absence of an orphan satellite population. Thus, ∼15% of the systems in our best-fit model are orphan satellites; these satellites might be associated with heavily stripped or disrupting subhalos.

The Poisson point process likelihood in our statistical comparison to observed satellites is implemented as follows. Suppose we observe n_i real satellites and ${\hat{n}}_{i,\nu }$ mock satellites in an absolute magnitude bin i, where $\nu =1,...,\hat{N}$ runs over all model realizations, including different host halos, and draws from our stochastic galaxy–halo connection model. The likelihood of observing these satellites given our model realizations, which enters Equation (6), is then

$$\begin{eqnarray}&&\begin{array}{l}P({n}_{i}| {\hat{n}}_{i,1},\ldots ,{\hat{n}}_{i,\hat{N}})={\left(\displaystyle \frac{\hat{N}+1}{\hat{N}}\right)}^{-({\hat{n}}_{i,1}+...+{\hat{n}}_{i,\hat{N}}+1)}\\ \quad \times \,{\left(\hat{N}+1\right)}^{-{n}_{i}}\displaystyle \frac{{\rm{\Gamma }}({\hat{n}}_{i,1}+...+{\hat{n}}_{i,\hat{N}}+{n}_{i}+1)}{{\rm{\Gamma }}({n}_{i}+1){\rm{\Gamma }}({\hat{n}}_{i,1}+...+{\hat{n}}_{i,\hat{N}}+1)},\end{array}\end{eqnarray} \tag{ B1 }$$

where the dependence on galaxy–halo connection model parameters ${\boldsymbol{\theta }}$ is implicit, and we assumed (i) a flat prior on λ_i for λ_i ≥ 0 and (ii) that n_i and all ${\hat{n}}_{i,\nu }$ are drawn from the same Poisson distribution with rate parameter λ_i. Note that our method yields noninteger numbers of mock satellites by counting each system as ${p}_{\mathrm{detect}}\times (1-{p}_{\mathrm{disrupt}})\times {f}_{\mathrm{gal}}$ objects according to Equation (5), so we have replaced factorials in the Poisson likelihood with appropriate gamma functions. Our results are unaffected if we enforce integer satellite counts by performing a binary mock observation of each predicted satellite according to its detection probability.

We list the prior distributions used in our fiducial analysis in Table B1, several of which are informed by previous work. The prior on the faint-end slope is a noninformative Jeffreys prior (Jethwa et al. 2018). The upper limit on the luminosity scatter is chosen to be very conservative; for example, Lehmann et al. (2017) found that abundance-matching scatter at the luminosity scale of the brightest systems used in our analysis is less than $\sim 0.25\ \mathrm{dex}$. For ${{ \mathcal M }}_{50}$, we set the lower limit of the prior based on the resolution limit of our simulations, which is a maximally conservative choice from the perspective of the inferred upper limit on this quantity. In particular, while we can decrease the lower limit of this prior because the ${{ \mathcal M }}_{50}$ posterior is flat below $\sim 5\times {10}^{7}\ {M}_{\odot }$ due to the limited sensitivity of the DES and PS1 satellite searches, doing so would artificially decrease the inferred 95% confidence upper limit.⁶³ Priors for ${ \mathcal B }$ and n are set based on studies that identify the preferred values of these parameters, and priors for ${\sigma }_{\mathrm{gal}},{ \mathcal A }$, and ${\sigma }_{\mathrm{log}R}$ are chosen to be uniform with conservative upper bounds.

Table B1.  Prior Distributions for the Parameters Varied in Our Fiducial Eight-parameter Fit to the DES and PS1 Satellite Populations

Free Parameter	Prior Distribution	Motivation
Faint-end slope	$\arctan \alpha \sim \mathrm{unif}(-1.1,-0.9)$	Jeffreys prior for $-2\lt \alpha \lt -1.2$
Luminosity scatter	${\sigma }_{M}\sim \mathrm{unif}(0,2)\ \mathrm{dex}$	Conservative upper limit (Garrison-Kimmel et al. 2017a; Lehmann et al. 2017)
50% occupation mass	$\mathrm{log}({{ \mathcal M }}_{50}/{M}_{\odot })\sim \mathrm{unif}(7.5,11)$	Lower limit corresponds to simulation resolution limit (Mao et al. 2015)
Baryonic effects	$\mathrm{ln}({ \mathcal B })\sim { \mathcal N }(\mu =1,\sigma =0.5)$	Hydrodynamic simulations (Nadler et al. 2018, 2019b)
Occupation scatter	${\sigma }_{\mathrm{gal}}\sim \mathrm{unif}(0,1)\ \mathrm{dex}$	Hydrodynamic simulations (Fitts et al. 2018; Graus et al. 2019)
Size amplitude	${ \mathcal A }\sim \mathrm{unif}(0,0.5)\ \mathrm{kpc}$	Empirical galaxy–halo size relation (Kravtsov 2013)
Size scatter	${\sigma }_{\mathrm{log}R}\sim \mathrm{unif}(0,2)\ \mathrm{dex}$	Empirical galaxy–halo size relation (Kravtsov 2013)
Size power-law index	$n\sim { \mathcal N }(\mu =1,\sigma =0.5)$	Empirical galaxy–halo size relation (Kravtsov 2013)

Note. Here ${ \mathcal N }(\mu ,\sigma )$ denotes a normal distribution with mean μ and standard deviation σ.

Download table as:  ASCIITypeset image

To calculate Bayes factors, we estimate the Bayesian evidence using the bounded harmonic mean method described in Nadler et al. (2019b). In particular, for a given posterior, we select samples of galaxy–halo connection model parameters ${\boldsymbol{\theta }}$ within a fixed Mahalanobis distance of a point ${{\boldsymbol{\theta }}}_{0}$ in a high-density region of the posterior. We then average the inverse of the posterior probabilities for these samples, and we normalize by the volume of the sampled region. We repeat this procedure for high-density regions that contain 10%–25% of the total number of MCMC samples, and we average over these percentiles to obtain the mean Bayesian evidence.

In this work, we fit to observed MW satellites in an observable parameter space x that consists of absolute magnitude and two large surface brightness bins. However, it would be more constraining to perform our inference in a higher-dimensional space that includes galactocentric distance. There are two main difficulties inherent in our statistical modeling.

• 1.  
  We have binned observed and modeled satellites assuming that the unknown Poisson process rate in each bin is independent from the rate in other bins. This assumption is unphysical, as the rate should vary smoothly in observable parameter space.
• 2.  
  As the number of bins increases, the number of satellites per bin decreases, which causes the uncertainty in the rate parameter to increase and our model to become increasingly unconstrained. This is a particularly challenging problem as we move to higher-dimensional parameter spaces, since the number of bins increases rapidly with dimensionality.

To address these issues, it is possible to connect rates in nearby regions of parameter space in an unbinned fashion using a correlated prior. This is equivalent to imposing that our galaxy–halo model should produce satellite abundances that vary smoothly as a function of observable quantities. We now lay out the mathematical formalism necessary for introducing this prior.

Our model of the distribution of satellites in observable space is an inhomogeneous Poisson process, where the number of "events" in any region ${ \mathcal T }$ of observable space x is given by a Poisson distribution with rate ${\lambda }_{{ \mathcal T }}={\int }_{{ \mathcal T }}\lambda \left(x\right){dx}$, where $\lambda \left(x\right)$ is referred to as the "rate function." Given a rate function $\lambda \left(x\right)$, the likelihood of observing N events at a set of points ${\left\{{x}_{i}\right\}}_{i=1}^{N}$ is

$$\begin{eqnarray}&&p\left({\left\{{x}_{i}\right\}}_{i=1}^{N}| \lambda \right)=\exp \left[-\int \lambda \left(x\right){dx}\right]\prod _{i=1}^{N}\lambda \left({x}_{i}\right),\end{eqnarray} \tag{ B2 }$$

where we suppressed the dependence of the rate on our galaxy–halo connection model parameters ${\boldsymbol{\theta }}$. In our case, the "events" ${\left\{{x}_{i}\right\}}_{i=1}^{N}$ are the locations of detected satellites in an observable parameter space. Note that in this formulation, there is no binning in x.

Calculating this likelihood exactly is challenging because, in order to compare observed and modeled satellite populations, we must integrate over the unknown rate function λ,

$$\begin{eqnarray}&&p\left(\left\{{x}_{i}\right\}| \left\{{\hat{x}}_{j}\right\}\right)=\displaystyle \frac{\int { \mathcal D }\lambda \,p\left(\left\{{x}_{i}\right\}| \lambda \right)p\left(\left\{{\hat{x}}_{j}\right\}| \lambda \right)p\left(\lambda \right)}{\int { \mathcal D }\lambda \,p\left(\left\{{\hat{x}}_{j}\right\}| \lambda \right)p\left(\lambda \right)}.\end{eqnarray} \tag{ B3 }$$

Here, both the numerator and denominator contain functional integrals over the rate; these integrals are performed over an infinite-dimensional space consisting of the rate at each point in observable parameter space. Further, this rate is a stochastic function in our galaxy–halo connection model due to satellite luminosity and size scatter. This makes our model an inhomogeneous Poisson process with a stochastic rate function, which is known as a "Cox process." The prior on the rate function, $p\left(\lambda \right)$, must admit only positive rates; one possible choice is to treat the logarithm of the rate as a Gaussian process. Models involving Cox processes are often termed "doubly intractable" due to the presence of intractable integrals over the rate function (Murray 2007).

There are, however, several approaches to make Cox processes tractable. As noted above, we bin satellites in absolute magnitude and split the sample into two large surface brightness bins, so that our likelihood is over the number of counts in each bin, rather than the locations of the points. This is equivalent to assuming that the rate function is constant in each bin and leads to the likelihood

$$\begin{eqnarray}&&p\left({\left\{{n}_{j}\right\}}_{j=1}^{{N}_{\mathrm{bins}}}| \left\{{\lambda }_{j}\right\}\right)=\exp \left(-\sum _{\mathrm{bins}\ j}{\lambda }_{j}{{ \mathcal V }}_{j}\right)\prod _{\mathrm{bins}\ k}\displaystyle \frac{{\lambda }_{k}^{{n}_{k}}}{{n}_{k}!},\end{eqnarray} \tag{ B4 }$$

where λ_j is the rate in bin j, ${{ \mathcal V }}_{j}$ is the volume of bin j, and n_j is the number of events in bin j. Binning turns the functional integral over $\lambda \left(x\right)$ in Equation (B3) into a finite-dimensional integral over the value of λ in each bin. Choosing Cartesian bins in observable parameter space then renders the problem tractable (Flaxman et al. 2015). There also exist approaches that avoid binning the observable space (Adams et al. 2009; John & Hensman 2018), which we intend to explore in future work.

To assess possible systematic uncertainties associated with the observed set of DES and PS1 satellites presented in Paper I, we rerun the entire analysis using only satellites that are confirmed to exhibit dark matter–dominated internal kinematics. The candidate satellites excluded from this reanalysis are indicated in Table C1. As shown in Figure C1, our galaxy–halo connection model constraints are largely unaffected by refitting the DES and PS1 satellite populations under the conservative assumption that all unconfirmed systems are star clusters. Most importantly, the upper limit on ${{ \mathcal M }}_{50}$ only increases by ∼1σ, to $5\times {10}^{8}\ {M}_{\odot }$ at 95% confidence, and the minimum halo mass increases to $6.5\times {10}^{8}\ {M}_{\odot }$, similar to the minimum halo mass inferred from classical and SDSS satellites in Nadler et al. (2019b). In addition, the total predicted number of MW satellites decreases by $\sim 1\sigma$ to 150 ± 60. These shifts are expected, since unconfirmed satellite candidates constitute many of the faintest systems in our fiducial sample. Thus, we conclude that our key constraints and predictions are not highly sensitive to the nature of kinematically unconfirmed satellite candidates.

Zoom In Zoom Out Reset image size

Table C1.  MW Satellites Used in Our Analysis

Name	MV	D	r1/2
	(mag)	(kpc)	(pc)
	DES
Fornax	−13.46	147	707
Sculptor	−10.82	84	223
Reticulum II	−3.88	30	31
Eridanus IIa	−7.21	380	158
Tucana II	−3.8	58	165
Grus II*	−3.9	53	92
Horologium I	−3.55	79	31
Tucana III*	−2.4	25	44
Tucana IV	−3.5	48	128
Phoenix II	−3.30	83	21
Horologium II*	−2.6	78	33
Tucana V*	−1.6	55	16
Pictor I*	−3.45	114	18
Columba I*	−4.2	183	98
Cetus II*	0.02	30	17
Grus I*	−3.47	120	21
Reticulum III*	−3.31	92	64
	PS1
Leo I	−11.78	254	226
Leo II	−9.74	233	165
Draco	−8.71	76	180
Ursa Minor	−9.03	76	272
Sextans	−8.72	86	345
Canes Venatici I	−8.8	218	338
Boötes I	−6.02	66	160
Ursa Major II	−4.25	32	85
Coma Berenices	−4.38	44	57
Sagittarius II	−5.2	69	32
Willman 1	−2.53	38	20
Canes Venatici II	−5.17	160	55
Segue 1	−1.30	23	20
Segue ${2}^{* }$	−1.86	35	34
Crater II	−8.2	117	1066
Draco II*	−0.8	22	17
Triangulum II*	−1.60	30	13
Hercules	−5.83	132	120
Cetus II*b	0.02	30	17

Notes. Properties of confirmed and candidate DES and PS1 satellites used in our analysis, listed in order of detection significance (Paper I). Asterisks mark kinematically unconfirmed systems.

^aEridanus II is not included because it lies outside our fiducial 300 kpc heliocentric distance cut. ^bCetus II is detected in both PS1 and DES; in our analysis, we only count this system in the observed DES population.

Download table as:  ASCIITypeset image

Next, we test whether a more conservative satellite size criterion impacts our results. For this test, we self-consistently exclude all observed and predicted satellites with r_1/2 > 20 pc from our statistical inference, rather than the r_1/2 > 10 pc cut used in our fiducial analysis. Our key constraints are not significantly affected; for example, the 95% confidence level upper limit on ${{ \mathcal M }}_{50}$ increases slightly, to $1.5\times {10}^{8}\ {M}_{\odot }$. The upper limit on the amplitude of the galaxy–halo size relation, which was 110 pc in our fiducial analysis, increases to 220 pc, as we might expect from excluding small satellites in the fit.

Finally, we test whether systematic offsets in measured satellite properties could affect our conclusions. In particular, we assume that every measured DES and PS1 satellite absolute magnitude is offset from the fiducial value listed in Table C1 by ${\rm{\Delta }}{M}_{V}=+1\ \mathrm{mag}$, which is similar to the width of the absolute magnitude bins used in our fiducial analysis. We rerun the entire analysis with these shifted magnitudes, and we repeat this procedure for ${\rm{\Delta }}{M}_{V}=-1\ \mathrm{mag}$. In both cases, we still obtain a good joint fit to the DES and PS1 luminosity functions. As expected, the inferred faint-end slope is steeper (shallower) than that obtained from our fiducial analysis for ${\rm{\Delta }}{M}_{V}=-1\ \mathrm{mag}$ (${\rm{\Delta }}{M}_{V}=+1\ \mathrm{mag}$); however, the total predicted number of MW satellites with ${M}_{V}\lt 0\ \mathrm{mag}$ and ${r}_{1/2}\gt 10\ \mathrm{pc}$ and our 95% confidence upper limit on ${{ \mathcal M }}_{50}$ are not significantly affected in either case.