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

E.O. Nadler, R.H. Wechsler, K. Bechtol, Y.-Y. Mao, G. Green, A.

Drlica-Wagner, M. Mcnanna, S. Mau, A.B. Pace, J.D. Simon, et al.

► To cite this version:

E.O. Nadler, R.H. Wechsler, K. Bechtol, Y.-Y. Mao, G. Green, et al., Milky Way Satellite Census – II. Galaxy-Halo Connection Constraints Including the Impact of the Large Magellanic Cloud. The Astrophysical Journal, 2020, 893, pp.48. 10.3847/1538-4357/ab846a hal-02416928

HAL Id: hal-02416928

https://hal.science/hal-02416928v1

Submitted on 9 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

E. O. Nadler<sup>1,2</sup>, R. H. Wechsler<sup>1,2,3</sup>, K. Bechtol<sup>4,5</sup>, Y.-Y. Mao<sup>6,7,8,53</sup>, G. Green<sup>9</sup>, A. Drlica-Wagner<sup>10,11,12</sup>, M. McNanna<sup>4</sup>, S. Mau<sup>11,12</sup>, A. B. Pace<sup>13,14</sup>, J. D. Simon<sup>15</sup>, A. Kravtsov<sup>11</sup>, S. Dodelson<sup>13</sup>, T. S. Li<sup>10,11,15,16,53</sup>, A. H. Riley<sup>14</sup>, M. Y. Wang<sup>13</sup>, T. M. C. Abbott<sup>17</sup>, M. Aguena<sup>18,19</sup>, S. Allam<sup>10</sup>, J. Annis<sup>10</sup>, S. Avila<sup>20</sup>, G. M. Bernstein<sup>21</sup>, E. Bertin<sup>22,23</sup>, D. Brooks<sup>24</sup>, D. L. Burke<sup>2,3</sup>, A. Carnero Rosell<sup>18,25</sup>, M. Carrasco Kind<sup>26,27</sup>, J. Carretero<sup>28</sup>, M. Costanzi<sup>29,30</sup>, L. N. da Costa<sup>18,31</sup>, J. De Vicente<sup>25</sup>, S. Desai<sup>32</sup>, A. E. Evrard<sup>33,34</sup>, B. Flaugher<sup>10</sup>, P. Fosalba<sup>35,36</sup>, J. Frieman<sup>10,11</sup>, J. García-Bellido<sup>20</sup>, E. Gaztanaga<sup>35,36</sup>, D. W. Gerdes<sup>33,34</sup>, D. Gruen<sup>1,2,3</sup>, J. Gschwend<sup>18,31</sup>, G. Gutierrez<sup>10</sup>, W. G. Hartley<sup>24,37</sup>, S. R. Hinton<sup>38</sup>, K. Honscheid<sup>39,40</sup>, E. Krause<sup>41</sup>, K. Kuehn<sup>42,43</sup>, N. Kuropatkin<sup>10</sup>, O. Lahav<sup>24</sup>, M. A. G. Maia<sup>18,31</sup>, J. L. Marshall<sup>14</sup>, F. Menanteau<sup>26,27</sup>, R. Miquel<sup>28,44</sup>, A. Palmese<sup>10,11</sup>, F. Paz-Chinchón<sup>26,27</sup>, A. A. Plazas<sup>16</sup>, A. K. Romer<sup>45</sup>, E. Sanchez<sup>25</sup>, B. Santiago<sup>18,46</sup>, V. Scarpine<sup>10</sup>, S. Serrano<sup>35,36</sup>, M. Smith<sup>47</sup>, M. Soares-Santos<sup>48</sup>, E. Suchyta<sup>49</sup>, G. Tarle<sup>34</sup>, D. Thomas<sup>50</sup>, T. N. Varga<sup>51,52</sup>, and A. R. Walker<sup>17</sup>

(DES Collaboration)

<sup>1</sup>Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA; enadler@stanford.edu

<sup>2</sup>Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA; rwechsler@stanford.edu

<sup>3</sup>SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

<sup>4</sup>Physics Department, 2320 Chamberlin Hall, University of Wisconsin–Madison, 1150 University Avenue, Madison, WI 53706-1390, USA

<sup>5</sup>LSST, 933 North Cherry Avenue, Tucson, AZ 85721, USA

<sup>6</sup>Department of Physics and Astronomy, Rutgers, The State University of New Jersey, Piscataway, NJ 08854, USA

<sup>7</sup>Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA

<sup>8</sup>Pittsburgh Particle Physics, Astrophysics and Cosmology Center (PITT PACC), University of Pittsburgh, Pittsburgh, PA 15260, USA

<sup>9</sup>Max Planck Institute for Astronomy, Königstuhl 17 D-69117, Heidelberg, Germany

<sup>10</sup>Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA

<sup>11</sup>Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

<sup>12</sup>Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

<sup>13</sup>Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15312, USA

<sup>14</sup>George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA

<sup>15</sup>Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA

<sup>16</sup>Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

<sup>17</sup>Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

<sup>18</sup>Laboratório Interinstitucional de e-Astronomia—LineA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

<sup>19</sup>Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil

<sup>20</sup>Instituto de Física Teorica UAM/CSIC, Universidad Autonoma de Madrid, E-28049 Madrid, Spain

<sup>21</sup>Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

<sup>22</sup>CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>23</sup>Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>24</sup>Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

<sup>25</sup>Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

<sup>26</sup>National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, IL 61801, USA

<sup>27</sup>Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 West Green Street, Urbana, IL 61801, USA

<sup>28</sup>Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain

<sup>29</sup>INAF–Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy

<sup>30</sup>Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

<sup>31</sup>Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

<sup>32</sup>Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

<sup>33</sup>Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

<sup>34</sup>Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

<sup>35</sup>Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain

<sup>36</sup>Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

<sup>37</sup>Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

<sup>38</sup>School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

<sup>39</sup>Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA

<sup>40</sup>Department of Physics, The Ohio State University, Columbus, OH 43210, USA

<sup>41</sup>Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

<sup>42</sup>Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

<sup>43</sup>Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

<sup>44</sup>Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

<sup>45</sup>Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

<sup>46</sup>Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS—91501-970, Brazil

<sup>47</sup>School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

<sup>48</sup>Brandeis University, Physics Department, 415 South Street, Waltham, MA 02453, USA

<sup>49</sup>Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

<sup>50</sup>Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

<sup>51</sup>Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany
⁵² Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, D-81679 München, Germany
*Received 2019 December 4; revised 2020 March 13; accepted 2020 March 28; published 2020 April 15*

**Abstract**

The population of Milky Way (MW) satellites contains the faintest known galaxies and thus provides essential insight into galaxy formation and dark matter microphysics. Here we combine a model of the galaxy-halo connection with newly derived observational selection functions based on searches for satellites in photometric surveys over nearly the entire high Galactic latitude sky. In particular, we use cosmological zoom-in simulations of MW-like halos that include realistic Large Magellanic Cloud (LMC) analogs to fit the position-dependent MW satellite luminosity function. We report decisive evidence for the statistical impact of the LMC on the MW satellite population due to an estimated $6 \pm 2$ observed LMC-associated satellites, consistent with the number of LMC satellites inferred from Gaia proper-motion measurements, confirming the predictions of cold dark matter models for the existence of satellites within satellite halos. Moreover, we infer that the LMC fell into the MW within the last 2 Gyr at high confidence. Based on our detailed full-sky modeling, we find that the faintest observed satellites inhabit halos with peak virial masses below $3.2 \times 10^8 M_\odot$ at 95% confidence, and we place the first robust constraints on the fraction of halos that host galaxies in this regime. We predict that the faintest detectable satellites occupy halos with peak virial masses above $10^6 M_\odot$, highlighting the potential for powerful galaxy formation and dark matter constraints from future dwarf galaxy searches.

*Unified Astronomy Thesaurus concepts: Dark matter (353); Milky Way dark matter halo (1049); Galaxy abundances (574)*

**1. Introduction**

The sample of confirmed and candidate Milky Way (MW) satellite galaxies has more than doubled in the last 5 yr. Modern imaging surveys have driven these discoveries; in particular, following the successes of the Sloan Digital Sky Survey (SDSS) in the early 2000s (Willman et al. 2005a, 2005b; Belokurov et al. 2006, 2007, 2008, 2009, 2010; Grillmair 2006, 2009; Sakamoto & Hasegawa 2006; Zucker et al. 2006; Irwin et al. 2007; Walsh et al. 2007), the Dark Energy Survey (DES) and the Panoramic Survey Telescope and Rapid Response System Pan-STARRS1 (PS1) have discovered 17 and three new satellite galaxy candidates, respectively (Bechtol et al. 2015; Drlica-Wagner et al. 2015; Kim & Jerjen 2015; Koposov et al. 2015; Laevens et al. 2015a, 2015b; Luque et al. 2016). These systems are identified as arcminute-scale overdensities of individually resolved stars, and many have already been spectroscopically confirmed. Meanwhile, other surveys with the Dark Energy Camera and VST ATLAS have recently discovered several additional satellites (Martin et al. 2015; Drlica-Wagner et al. 2016; Torrealba et al. 2016a, 2016b, 2018; Koposov et al. 2018; Mau et al. 2019).

Nonetheless, the current census of MW satellites is likely highly incomplete, particularly for faint systems in the outer regions of the MW halo. This is evidenced by the detection of three new ultrafaint satellites in the first ~676 deg² of Hyper Suprime-Cam Strategic Survey Program (HSC-SSP) imaging data (Homma et al. 2016, 2018, 2019) and the discovery of Antlia II, the lowest surface brightness galaxy currently known, using RR Lyrae member stars identified in Gaia DR2 (Torrealba et al. 2019). In the near future, the Legacy Survey of Space and Time (LSST) conducted from the Vera C. Rubin Observatory will be able to detect satellites over the entire southern sky down to a surface brightness of $\mu_V \sim 32$ mag arcsec$^{-2}$ (Ivezić et al. 2008; Tollerud et al. 2008; Hargis et al. 2014; Nadler et al. 2019b).

Interpreting the cosmological and astrophysical implications of these discoveries requires a detailed understanding of the

observational selection effects for each survey under consideration. In a companion paper (Drlica-Wagner et al. 2020, hereafter Paper I), we derive observational selection functions for DES and PS1 based on searches for simulated satellites in each data set. These selection functions encode the probability that satellites in either survey are detectable as a function of their absolute magnitude, heliocentric distance, physical size, and position on the sky. They incorporate realistic photometric error models, selection masks that exclude highly reddened regions near the Galactic disk, and the influence of local stellar density on satellite detectability. Detection sensitivity is linked to sky position because various surveys have imaged different parts of the sky at varying depths, and accurately modeling this effect is crucial in order to disentangle anisotropy in the underlying MW satellite system from selection effects.

In this paper, we combine the observational selection functions derived in Paper I with a detailed model of the galaxy-halo connection and high-resolution cosmological zoom-in simulations of MW-mass host halos to infer the position-dependent MW satellite luminosity function. Although several empirical models have recently been used to study subsets of the MW satellite population (Jethwa et al. 2018; Kim et al. 2018; Newton et al. 2018; Nadler et al. 2019b), this is the first analysis that is directly based on imaging data over more than ~15,000 deg²; indeed, our analysis covers 75% of the high Galactic latitude sky. Moreover, our galaxy-halo connection model allows us to marginalize over astrophysical uncertainties in our fit to the observed DES and PS1 satellite populations. We quantify the impact of the largest MW satellite, the Large Magellanic Cloud (LMC), and its associated satellites on the observed DES and PS1 satellite populations. We find that the satellites accreted with a realistic LMC analog—defined in terms of its mass, heliocentric distance, and infall time—are essential to fit the DES and PS1 luminosity functions simultaneously; this finding constitutes a remarkable confirmation of hierarchical structure formation. We predict that $4.8 \pm 1.7$ (1.1 ± 0.9) of the known satellites observed by DES (PS1) fell into the MW with the LMC, consistent with the number of LMC-associated satellites inferred from Gaia

<footnote n="53">⁵³ NHFP Einstein Fellow.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Visualization of our MW satellite modeling framework. In the first step, we perform high-resolution zoom-in simulations of MW-like halos selected from a larger cosmological volume (Section 3); in the second step, we paint galaxies onto subhalos using a parametric model for the galaxy–halo connection (Section 4); in the third step, we use the observational selection functions derived in Paper I to compute the probability that these satellites would be observed in DES or PS1 imaging data (Section 5); and in the final step, we calculate the likelihood of producing the true DES and PS1 satellite populations given many mock satellite population realizations at fixed galaxy–halo connection model parameters (Section 6). We then iterate this process to constrain our galaxy–halo connection model.</figcaption></figure>

proper-motion measurements (Kallivayalil et al. 2018; Patel et al. 2020).

Our analysis constrains the properties of the lowest-mass halos that host observed satellites, which we infer to have peak virial masses below $3.2 \times 10^8 M_{\odot}$ at 95% confidence. This finding, along with constraints on the faint-end slope of the luminosity function, can be used to inform feedback prescriptions in hydrodynamic simulations (Sawala et al. 2016b; Fitts et al. 2018; Simpson et al. 2018; Munshi et al. 2019; Wheeler et al. 2019). Constraints on the minimum halo mass also hold broad implications for the microphysical properties of dark matter (e.g., Drlica-Wagner et al. 2019; Nadler et al. 2019a). Crucially, our model can be extended to explore the degeneracies between baryonic physics and deviations from the cold dark matter (CDM) paradigm.

This paper is organized as follows. We first provide an overview of our framework in Section 2. We then describe the simulations (Section 3), galaxy–halo connection model (Section 4), observational selection functions (Section 5), and statistical framework (Section 6) used in our analysis. We present our results in Section 7, focusing on the observed DES and PS1 satellite populations (Section 7.1), the impact of the LMC system (Section 7.2), the total MW satellite population

(Section 7.3), galaxy–halo connection model constraints (Section 7.4), the properties of halos that host the faintest observed satellites (Section 7.5), and the implications of our findings for dark matter microphysics (Section 7.6). We discuss the main theoretical uncertainties in our analysis in Section 8, and we conclude in Section 9. Appendices provide additional details on our galaxy–halo connection model (Appendix A) and statistical framework (Appendix B), the robustness of our results to observational systematics (Appendix C) and resolution effects (Appendix D), and the observed DES and PS1 satellite populations (Appendix E).

Throughout, we use the term “galaxy–halo connection model” to refer to a model that describes how the properties of galaxies, including luminosity and size, are related to the properties of halos, such as peak virial mass. Furthermore, “log” refers to the base-10 logarithm.

## 2. Analysis Overview

Using the observed population of MW satellites to constrain our galaxy–halo connection model requires the following ingredients (see Figure 1 for a visualization of each step).
1. A model that predicts the *underlying* MW satellite population.

2. An observational selection function that, convolved with the prediction in the previous step, yields a prediction for the *observed* satellite population.

3. A model for the likelihood of producing the true MW satellite population given the prediction from the previous step.

The first step above can be performed using either a hydrodynamic simulation, in which galaxy formation is modeled at the simulation level, or an empirical prescription for painting galaxies onto halos. We take the latter approach in this paper, which allows for a more flexible modeling framework, as well as the use of simulations that are closer approximations to the MW system. Indeed, our results may help to constrain feedback prescriptions in hydrodynamic simulations. Note that the assumed dark matter model (e.g., cold versus warm dark matter or collisionless versus interacting dark matter) affects the underlying satellite population and often manifests as a cutoff in the abundance of halos—and thus faint galaxies—below a halo mass threshold determined by the microphysical properties of dark matter.

The steps outlined above each rely on tools that have been developed in previous studies and Paper I. Here we simply provide a brief description of each step, and we defer additional details to the Appendices.

## 3. Simulations

### 3.1. General Description

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}$.<fnref n="54" /> 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 \text{ 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_{\text{max}} \approx 9 \text{ km s}^{-1}$ (Mao et al. 2015). To be conservative, we only use subhalos with both $V_{\text{max}} > 9 \text{ km s}^{-1}$ and peak maximum circular velocity $V_{\text{peak}} > 10 \text{ km s}^{-1}$. In Appendix D, we show that these resolution thresholds are sufficient for modeling the satellite populations of interest here.

### 3.2. Host Halo and LMC Analog Selection

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_{\text{max}} \ge 55 \text{ km s}^{-1}$,

2. present-day heliocentric distance $40 \text{ kpc} < D < 60 \text{ 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).<fnref n="55" /> 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

<footnote n="54">⁵⁴ We define virial quantities according to the Bryan–Norman virial overdensity (Bryan & Norman 1998), with $\Delta_{\text{vir}} \simeq 99.2$ as appropriate for the cosmological parameters adopted in our zoom-in simulations: $h = 0.7$, $\Omega_m = 0.286$, $\Omega_b = 0.047$, and $\Omega_{\Lambda} = 0.714$ (Mao et al. 2015).</footnote>

<footnote n="55">⁵⁵ Although detailed exploration of Magellanic Cloud binary systems is beyond the scope of this work, we note that Shao et al. (2018) and Cautun et al. (2019a) found that LMC analogs with SMC-like companions are typically more massive than isolated LMC analogs.</footnote>
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.

## 4. Galaxy–Halo Connection Model

To associate satellite galaxies with subhalos in the simulations described above, we use a modified version of the model developed in Nadler et al. (2019b). This model parameterizes the relationship between satellite and subhalo properties and the effects of baryonic physics on subhalo populations in flexible ways, which allows us to marginalize over the relevant theoretical uncertainties. Additional model details and tests are presented in Appendix A.

### 4.1. Satellite Luminosities

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}$.<fnref n="56" /> This relation is constrained by the GAMA survey (Loveday et al. 2015; Geha et al. 2017) for bright systems ($M_V < -13$ mag) and is extrapolated into the regime of dim satellites by treating the faint-end slope of the satellite luminosity function, $\alpha$, and the lognormal scatter in luminosity at fixed $V_{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$ 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.

### 4.2. Satellite Sizes

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

$$r_{1/2} \equiv A \left( \frac{R_{\text{vir}}}{R_0} \right)^n, \qquad (1)$$

where $A$ and $n$ are free parameters, $R_{vir}$ denotes the subhalo virial radius measured at accretion, and $R_0 = 10$ 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_{\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} > 10$ pc in order to exclude likely star clusters from the analysis. We explore a more conservative cut of $r_{1/2} > 20$ 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.

### 4.3. Subhalo Disruption Due to Baryonic Effects

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 $B$, which is defined such that $B = 1$ corresponds to fiducial hydrodynamic predictions (Nadler et al. 2018) and larger (smaller) values of $B$ correspond to more effective (less effective) subhalo disruption. For each subhalo, we set

$$p_{\text{disrupt}} \equiv (p_{\text{disrupt},0})^{1/B}, \qquad (2)$$

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

### 4.4. Galaxy Formation Efficiency

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),

$$f_{\text{gal}}(M_{\text{peak}}) \equiv \frac{1}{2} \left[ 1 + \text{erf} \left( \frac{M_{\text{peak}} - M_{\text{50}}}{\sqrt{2} \sigma_{\text{gal}}} \right) \right], \qquad (3)$$

where $M_{peak}$ is the largest virial mass a subhalo ever attains, which typically occurs before infall into the MW; $M_{50}$ is the peak halo mass at which 50% of halos host galaxies of any mass; and $\sigma_{gal}$ is the width of the galaxy occupation fraction. In our fiducial model, $M_{50}$ and $\sigma_{gal}$ are free parameters. Note that

<footnote n="56">⁵⁶ We perform abundance matching using $V_{peak}$ to incorporate the effects of halo assembly bias and mitigate the impact of subhalo tidal stripping (Reddick et al. 2013; Lehmann et al. 2017).</footnote>
in the limit $\sigma_{gal} \to 0$, this reduces to a model in which all halos with $M_{peak} > M_{50}$ host a galaxy.

Although our analysis does not constrain $\sigma_{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 $M_{peak}$ is more easily interpretable and connects directly to constraints on alternative dark matter models (e.g., Nadler et al. 2019a).

### 4.5. Orphan Satellites

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

$$p_{disrupt} \equiv (1 - a_{acc})^{\mathcal{O}}, \quad (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.

### 5. Observational Selection Functions

We employ the DES and PS1 survey selection functions derived in Paper I, which have been publicly distributed as machine-learning classifiers that predict satellite detection probability given absolute magnitude, $M_V$; heliocentric distance, $D$; azimuthally averaged projected half-light radius, $r_{1/2}$; and sky position.<fnref n="57" /> The predicted detection probabilities are derived from searches for simulated satellites in catalog-level DES and PS1 data, and they employ geometric cuts that restrict observable satellites to lie within the respective survey footprint and mask regions where satellite detection is challenging due to interstellar extinction, bright nearby stars, and bright extragalactic objects.

We self-consistently apply these position-dependent detection criteria to our predicted satellite populations by matching the on-sky position of our LMC analogs to the true on-sky position of the LMC. In particular, we choose random observer locations 8 kpc from the halo center, and we perform

appropriate rotations to our subhalo populations for each observer location to match the true LMC position. We apply the DES selection function for satellites within the overlap region of the two surveys, and we only count satellites within a fiducial heliocentric distance of 300 kpc.

## 6. Statistical Framework

To fit our galaxy-halo model to the DES and PS1 luminosity functions derived in Paper I, we generate predicted satellite populations given a set of galaxy-halo connection model parameters, $\theta$, by performing mock DES-plus-PS1 surveys using the selection functions described above. For each host halo and each realization of our satellite model, we bin mock-observed satellites according to their absolute magnitude. We further split satellites in each absolute magnitude bin into high ($\mu_V < 28 \text{ mag arcsec}^{-2}$) and low ($\mu_V \ge 28 \text{ mag arcsec}^{-2}$) surface brightness samples to incorporate satellite size information in our fit.<fnref n="58" /> We list the DES and PS1 satellites used in this analysis in Table C1.

Next, we calculate the number of predicted satellites in each bin $i$ via

$$n_i = \sum_{s_i} p_{\text{detect},s_i} \times (1 - p_{\text{disrupt},s_i}) \times f_{\text{gal},s_i}, \quad (5)$$

where $s_i$ indexes the satellites in bin $i$, $p_{\text{detect}}$ is the detection probability returned by the appropriate observational selection function, $p_{\text{disrupt}}$ is the disruption probability due to baryonic effects (Equation (2)), and $f_{\text{gal}}$ is the galaxy occupation fraction (Equation (3)). For objects that lie in the overlap region of the DES and PS1 footprints, we calculate $p_{\text{detect}}$ using the DES selection function.

We note that detection probability mainly depends on surface brightness and present-day heliocentric distance (see Paper I), disruption probability mainly depends on orbital properties (Nadler et al. 2018), and galaxy occupation depends on $M_{peak}$ according to Equation (3). Thus, our model for satellite detectability is coupled to our galaxy occupation fraction model, since surface brightness is directly linked to $M_{peak}$ due to our abundance-matching assumption. Nonetheless, our results are largely unaffected if we exclude the galaxy occupation fraction from our model, confirming that $f_{\text{gal}}$ can be interpreted as the probability that a halo hosts a satellite brighter than $M_V = 0$ mag, corresponding to the faintest satellite in our observational sample.

We assume that satellites populate each bin in absolute magnitude-surface brightness parameter space according to an independent Poisson point process with a rate parameter $\lambda$ that depends on absolute magnitude, surface brightness, and our galaxy-halo connection model parameters. Because our model yields noisy estimates of $\lambda$, we marginalize over its range of possible values in each bin, following Nadler et al. (2019b). The likelihood of observing the set of DES and PS1 satellites, $s_{\text{DES}}$ and $s_{\text{PS1}}$ (specified by their absolute magnitudes and surface brightnesses), given a set of model parameters $\theta$ is then

$$P(s_{\text{DES}}, s_{\text{PS1}}|\theta) = \prod_{\text{bins } i} P(n_{\text{DES},i}|\hat{n}_{\text{DES},i}) \times P(n_{\text{PS1},i}|\hat{n}_{\text{PS1},i}), \quad (6)$$

<footnote n="57">⁵⁷ The DES Y3A2 and PS1 DR1 selection functions are publicly available at https://github.com/des-science/mw-sats.</footnote>

<footnote n="58">⁵⁸ In particular, we calculate the effective surface brightness averaged within the half-light radius as $\mu_V = M_V + 36.57 + 2.5 \log(2\pi r_{1/2}^2)$, where $r_{1/2}$ is measured in units of kpc.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Predicted DES and PS1 satellite luminosity functions resulting from a joint fit to these satellite populations. Dark (light) blue bands correspond to 68% (95%) confidence intervals from our fiducial eight-parameter galaxy-halo connection model, dashed red lines show the 68% confidence interval for a model using host halos without LMC analogs ("No LMC"), and black lines show the observed luminosity functions within each survey footprint. Our fiducial model, which includes realistic LMC analogs, is decisively favored over the No LMC scenario, with a Bayes factor of ~$10^4$.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

where $n_{DES,i}$ ($n_{PS1,i}$) is the observed number of DES (PS1) satellites in bin $i$, and $\hat{n}_{DES,i}$ ($\hat{n}_{PS1,i}$) is a vector of the number of mock DES (PS1) satellites in bin $i$ from several realizations of our model at fixed $\theta$. These realizations include draws over host halos, observer locations, and our galaxy-halo connection model, which is stochastic at fixed $\theta$. Note that steps 1–3 in Figure 1 generate mock satellite populations $\hat{n}_{DES}$ and $\hat{n}_{PS1}$, and step 4 compares these to the observed populations $n_{DES}$ and $n_{PS1}$. The explicit forms of $P(n_{DES,i}|\hat{n}_{DES,i})$ and $P(n_{PS1,i}|\hat{n}_{PS1,i})$ are given in Equation (B1).

Finally, we use Bayes's theorem to compute the posterior distribution over galaxy-halo connection model parameters,

$$P(\theta|\mathbf{s}_{\text{DES}}, \mathbf{s}_{\text{PS1}}) = \frac{P(\mathbf{s}_{\text{DES}}, \mathbf{s}_{\text{PS1}}|\theta)P(\theta)}{P(\mathbf{s}_{\text{DES}}, \mathbf{s}_{\text{PS1}})}, \quad (7)$$

where $P(\theta)$ is our prior on the galaxy-halo connection model parameters (given in Appendix B.2), $P(\mathbf{s}_{\text{DES}}, \mathbf{s}_{\text{PS1}})$ is the Bayesian evidence, and $P(\mathbf{s}_{\text{DES}}, \mathbf{s}_{\text{PS1}}|\theta)$ is given by Equation (6). To sample from this posterior, we run $10^5$ iterations of the Markov Chain Monte Carlo (MCMC) sampler `emcee` (Foreman-Mackey et al. 2013) to sample the eight free parameters $\theta = (\alpha, \sigma_M, M_{50}, B, \sigma_{\text{gal}}, A, \sigma_{\text{log}R}, n)$ using 32 walkers. We discard a burn-in period of 20 autocorrelation lengths, corresponding to ~$10^4$ steps, which leaves more than 100 independent samples.

## 7. Results

We now present our results, focusing on the observed DES and PS1 satellite populations (Section 7.1), the impact of the LMC system (Section 7.2), the total MW satellite population (Section 7.3), the galaxy-halo connection model constraints (Section 7.4), the properties of the halos that host faint satellites (Section 7.5), and the implications for dark matter microphysics (Section 7.6).

### 7.1. Observed Satellite Populations

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.

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_{\text{peak}} > 10 \text{ km s}^{-1}$ and $V_{\text{max}} > 9 \text{ km s}^{-1}$ does not significantly affect our inference.

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 \text{ 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.

### 7.2. The Impact of the LMC

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.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Size distributions derived by fitting our galaxy-halo connection model to the DES and PS1 satellite populations. Dark (light) blue bands correspond to 68% (95%) confidence intervals from our fiducial eight-parameter model, dashed red lines show the 68% confidence interval for a model using host halos without LMC analogs ("No LMC"), and black lines show the observed size distributions.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

(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^4$, $10^4$, and $10^3$, respectively. In addition, both host halos in the Early LMC Infall case are individually disfavored with Bayes factors of ~$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_{peak} > 10$ km s⁻¹ and $V_{max} > 9$ km s⁻¹) 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 < 0$ mag and $r_{1/2} > 10$ 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).<fnref n="59" /> In addition, the orbital

<footnote n="59">⁵⁹ A recent analysis based on Gaia proper-motion measurements and hydrodynamic simulations suggests that two bright satellites in or near DES, Fornax and Carina, may also be LMC-associated (Pardy et al. 2020; however, see Patel et al. 2020).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Left panel: total MW satellite luminosity function inferred from our joint fit to the DES and PS1 satellite populations (blue) compared to the current census of confirmed and candidate MW satellites (black) and the empirical estimate derived in Paper I (gray), which assumes an isotropic satellite distribution and a cored NFW radial satellite distribution. The 68% confidence intervals from hydrodynamic simulations of the Local Group using the FIRE feedback prescription are shown in red (Garrison-Kimmel et al. 2019). Luminosity function slopes predicted from hydrodynamic simulations with (solid green line) and without (dashed green line) H₂-based star formation are shown for comparison (Munshi et al. 2019); these predictions do not account for subhalo disruption due to the Galactic disk. Note that the Paper I prediction (gray) differs from the “All Known Satellites” curve (black) at the bright end because it does not include the LMC, SMC, or Sagittarius. Right panel: surface brightness distribution of MW satellites with $M_V < 0$ mag and $r_{1/2} > 10$ pc as a function of the limiting observable surface brightness of an all-sky survey. Arrows indicate approximate detection limits for current surveys. Note that LSST Y1 is expected to have similar detection sensitivity to HSC (Ivezić et al. 2008; Tollerud et al. 2008; Hargis et al. 2014; Nadler et al. 2019b).</figcaption></figure>

<figure><img src="image_9.png" /></figure>

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.

### 7.3. The Total MW Satellite Population

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).

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$ mag arcsec⁻² 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
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Posterior distribution from our fit to the DES and PS1 satellite populations. Dark (light) shaded contours represent 68% (95%) confidence intervals. Shaded areas in the marginal distributions and parameter summaries correspond to 68% confidence intervals. Note that $\sigma_M$, $\sigma_{gal}$, and $\sigma_{log R}$ are reported in dex, $M_{50}$ is reported as $\log(M_{50}/M_\odot)$, and $A$ is reported in pc.</figcaption></figure>

luminosity function at $M_V < -16$ 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$ mag (corresponding to stellar mass $M_* \lesssim 10^5 M_\odot$), can be used to inform subgrid star formation prescriptions.

### 7.4. Galaxy–Halo Connection Model Constraints

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\sigma$ relative to the
<table><caption>**Table 1**
Galaxy–Halo Connection Model Constraints Derived from Our Fit to the DES and PS1 Satellite Populations</caption><thead><tr><th>Parameter</th><th>Physical Interpretation</th><th>95% Confidence Interval</th></tr></thead><tbody><tr><td>Faint-end slope</td><td>Power-law slope of satellite luminosity function</td><td>-1.46 &lt; α &lt; -1.39</td></tr><tr><td>Luminosity scatter</td><td>Scatter in luminosity at fixed V<sub>peak</sub></td><td>0 dex* &lt; σ<sub>M</sub> &lt; 0.19 dex</td></tr><tr><td>50% occupation mass</td><td>Mass at which 50% of halos host galaxies</td><td>7.5* &lt; log(M<sub>50</sub>/M<sub>☉</sub>) &lt; 7.93</td></tr><tr><td>Baryonic effects</td><td>Strength of subhalo disruption due to baryons</td><td>0.3 &lt; B &lt; 2.1</td></tr><tr><td>Occupation scatter</td><td>Scatter in galaxy occupation fraction</td><td>0 dex* &lt; σ<sub>gal</sub> &lt; 0.67 dex</td></tr><tr><td>Size amplitude</td><td>Amplitude of galaxy–halo size relation</td><td>0* pc &lt; A &lt; 110 pc</td></tr><tr><td>Size scatter</td><td>Scatter in half-light radius at fixed halo size</td><td>0 dex* &lt; σ<sub>log R</sub> &lt; 1.2 dex</td></tr><tr><td>Size power-law index</td><td>Power-law index of galaxy–halo size relation</td><td>0* &lt; n &lt; 1.8</td></tr></tbody></table>

**Note.** Asterisks mark prior-driven constraints.

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 < 0.19$ 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 $\sigma_M$ is consistent with previous studies of faint galaxy samples (e.g., Lehmann et al. 2017). Meanwhile, large values of $\sigma_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 $\sigma_M$, finding that $\sigma_M = 0.15$ dex ($\sigma_M = 0.2$ dex) is disfavored relative to $\sigma_M = 0$ dex with a Bayes factor of 30 (100).<fnref n="60" /> These upper limits are comparable to the scatter typically inferred from abundance-matching analyses of brighter systems ($\sigma_M \sim 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_{☉}$ at 95% confidence. Note that this summary statistic depends on the lower

limit of our prior on $M_{50}$, since the $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 $M_{50}$ to confirm that this summary statistic is robust. We find that $M_{50} = 8.5 \times 10^7 M_{☉}$ ($M_{50} = 1.5 \times 10^8 M_{☉}$) is disfavored relative to arbitrarily low values of $M_{50}$ with a Bayes factor of 50 (100). The current data are not able to place a lower limit on $M_{50}$, which would correspond to the detection of a cutoff in galaxy formation.

4. Our posterior is consistent with $B = 1$, corresponding to our fiducial baryonic disruption model. Although a large spread in disruption strength is allowed by the data, extremely efficient ($B > 2.1$) and inefficient ($B < 0.3$) subhalo disruption relative to hydrodynamic simulations is strongly disfavored. These constraints widen when our lognormal prior on $B$ is relaxed; however, zero subhalo disruption (corresponding to $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 $M_{50}$. Models with large scatter in the occupation fraction ($\sigma_{gal} < 0.67$ 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 $\sigma_{gal}$ posterior is driven by the lower limit of our $M_{50}$ prior; as this limit decreases, the $\sigma_{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_{vir} = 10$ kpc at accretion, is constrained to lie between 12 and 73 pc at 68% confidence. For larger values of $A$, satellites are too large to be detected with high probability, and the DES and PS1 luminosity functions are underpredicted; for smaller values of $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_{log R}$, faint satellites upscatter to large sizes too frequently, which results in underpredicted luminosity functions. Our 68% confidence lower limit on $\sigma_{log R}$ of 0.33 dex is consistent with the value estimated in Kravtsov (2013). Lower values of $\sigma_{log R}$ lead to slightly

<footnote n="60">⁶⁰ We provide details on our Bayes factor calculations in Appendix B.3.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Left panel: fraction of halos that host galaxies, inferred from our fit to the DES and PS1 satellite populations. The solid line shows the median inferred galaxy occupation fraction, and dark (light) shaded contours represent 68% (95%) confidence intervals. The resolution limit of our simulations is indicated by the dashed vertical line. Right panel: SMHM relation inferred from our fit to the DES and PS1 satellite populations. An extrapolation of the mean SMHM relation derived from more luminous field galaxies is shown in gray (Behroozi et al. 2013a). Stars illustrate the mean of the predicted $M_{\text{peak}}$ range corresponding each observed DES and PS1 satellite, and top ticks indicate the corresponding present-day virial masses of the halos that host these systems.</figcaption></figure>

too many predicted DES and PS1 satellites; however, our results are consistent with $\sigma_{\log R} = 0$ 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.

### 7.5. Properties of Halos that Host the Faintest Satellites

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 \text{ km 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 $M_{\text{min}} < 3.2 \times 10^8 M_{\odot}$ and $V_{\text{peak,min}} < 21 \text{ km s}^{-1}$ at 95% confidence. Furthermore, we predict that the faintest observed satellite inhabits a halo with $M_{\text{peak}} = 1.5 \times 10^8 M_{\odot}$, on average.<fnref n="61" />

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_{\text{peak}} \approx 20 \text{ km 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 $M_{\text{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_{\text{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 $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 $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 $M_{50} \sim 10^8 M_{\odot}$ (Katz et al. 2019).

Our galaxy occupation fraction constraint implies that models with $M_{50} > 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

<footnote n="61">⁶¹ The faintest observed satellite in our analysis, Cetus II, is detected by DES with $M_V = 0.02$ mag (Drlica-Wagner et al. 2015; Table C1).</footnote>
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 $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 > -6$ 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.<fnref n="62" /> 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 $M_{\text{peak}} > 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.

### 7.6. Implications for Dark Matter Microphysics

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_{\text{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 $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_{\text{hm}}$ equal to our upper limit on the minimum halo mass, i.e., $M_{\text{hm}} < 3.2 \times 10^8 M_{\odot}$, which corresponds to a lower limit on the half-mode scale of $k_{\text{hm}} > 36 h \text{ Mpc}^{-1}$. Using the relations in Nadler et al. (2019a) with the cosmological parameters $h = 0.7$ and $\Omega_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 \times 10^{-30} \text{ cm}^2$ 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.

## 8. Theoretical Uncertainties

We aim to present a thorough galaxy-halo connection model that allows us to marginalize over the most important theoretical uncertainties when modeling the MW satellite population. Nonetheless, our modeling choices necessarily affect the predicted number of detected low-mass halos and thus our upper limit on the minimum halo mass, $M_{\text{min}}$. In this section, we briefly discuss the main uncertainties in this analysis and their impact on our $M_{\text{min}}$ constraint.

To do so, we consider the upper limit on $M_{\text{min}}$ as a function of modeling assumptions, starting with the most conservative model possible and adding one assumption at a time. We illustrate the results of this exercise in Figure 7 for upper limits calculated as follows.

(i) Minimal CDM. Assuming a maximally massive MW halo given Gaia constraints (i.e., a virial mass of $1.8 \times 10^{12} M_{\odot}$; Callingham et al. 2019; Cautun et al. 2019b; Li et al. 2019a, 2019b), count the subhalos within

<footnote n="62">⁶² For many MW satellites, this will likely require spectroscopy with giant segmented mirror telescopes (Simon et al. 2019).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Impact of modeling assumptions on the minimum subhalo mass inferred from the observed DES and PS1 satellite populations. The first three models match the number of subhalos to the number of confirmed DES and PS1 satellites, and the last two models populate subhalos with galaxies to fit the position-dependent MW satellite luminosity function and size distribution.</figcaption></figure>

the virial radius of the MW in order of decreasing $V_{\text{peak}}$ until the number of kinematically confirmed DES and PS1 satellites is matched, and set the lowest corresponding value of $M_{\text{peak}}$ equal to the upper limit on $M_{\text{min}}$.

(ii) MW host mass. Repeat the previous step with the MW host mass fixed to its average value in our two fiducial simulations (i.e., an average virial mass of $1.4 \times 10^{12} M_{\odot}$).

(iii) Subhalo disruption. Repeat the previous step many times with the subhalo number weighted by disruption probability, sampling $\mathcal{B}$ from our fiducial posterior, to calculate an upper limit on $M_{\text{min}}$ at 95% confidence.

(iv) Satellites (confirmed). Repeat the previous step including the observational detection probabilities for mock satellites in the DES and PS1 footprints by drawing satellite properties from our fiducial posterior.

(v) Satellites (unconfirmed). Repeat the previous step including the unconfirmed candidate satellites detected by DES and PS1 in the observed tally.

This yields $M_{\text{min}} < (17, 14, 12, 6.5, 2.5) \times 10^8 M_{\odot}$ for models (i)–(v), respectively. Note that models (i)–(iii) are extremely conservative, since subhalos are counted in order of decreasing $V_{\text{peak}}$; however, these models do not reproduce the observed position-dependent MW satellite luminosity function or radial distribution. Model (iv) yields the conservative limit presented in Appendix C.1, and model (v) yields our fiducial constraint, uncorrected for MW host halo mass. Although we have not explicitly considered artificial subhalo disruption in this list of theoretical uncertainties (e.g., van den Bosch & Ogiya 2018; van den Bosch et al. 2018), our fiducial orphan satellite model effectively assumes that subhalo disruption in dark matter-only simulations is entirely artificial, which is a conservative choice.

Figure 7 shows that both fitting the satellite luminosity function and including the population of faint, kinematically unconfirmed satellite galaxies in our fit yield significant increases in constraining power. We emphasize that our galaxy-halo connection model is conservative from the perspective of upper limits on the minimum halo mass, since we assume that high-mass halos host the brightest observed satellite galaxies. Moreover, we marginalize over many uncertainties in the connection between low-mass halos and faint galaxies. Thus, the largest gain in constraining power likely results from our detailed use of observational selection functions, i.e., from the fact that some satellites are not detected

in DES or PS1 data. Given our extensive validation of the DES and PS1 selection functions in Paper I, we are therefore confident in our minimum halo mass constraints.

## 9. Conclusions

We have presented the results of a forward-modeling framework for MW satellites applied to recent searches for satellites in photometric surveys over nearly the entire high Galactic latitude sky. Our analysis includes position-dependent observational selection effects that faithfully represent satellite searches in DES and PS1 imaging data, and our galaxy-halo connection model allows us to marginalize over theoretical uncertainties in the relationship between galaxy and halo properties, the effects of baryonic physics on subhalo populations, and the stochastic nature of galaxy formation in low-mass halos. By performing a Bayesian analysis of the observed DES and PS1 satellite populations, we find decisive statistical evidence for the following.

1. The LMC impacts the observed MW satellite population, contributing $4.8 \pm 1.7$ (1.1 $\pm$ 0.9) LMC-associated satellites to the DES (PS1) satellite populations.

2. The LMC fell into the MW within the last 2 Gyr.

3. The faintest satellites currently known occupy halos with peak virial masses less than $3.2 \times 10^8 M_{\odot}$.

4. The faintest detectable satellites (i.e., dark matter-dominated systems with $M_{\star} > 100 M_{\odot}$) occupy halos with peak virial masses greater than $10^6 M_{\odot}$.

These results have broad implications for galaxy formation and dark matter physics. For example, comparing our inferred luminosity function and galaxy occupation fraction to predictions from hydrodynamic simulations will help break degeneracies subgrid star formation and feedback models. Meanwhile, extending our model to study the evolution of the luminosity function will shed light on high-redshift faint galaxy populations (e.g., Boylan-Kolchin et al. 2015; Weisz & Boylan-Kolchin 2017) and the MW’s local reionization history (e.g., Busha et al. 2010; Lunnan et al. 2012; Katz et al. 2019).

Finally, our statistical detection of low-mass halos translates directly into constraints on a suite of dark matter properties, including warmth in thermal production scenarios, initial velocity distribution in nonthermal production scenarios, self-interaction cross section, interaction strength with the Standard Model, formation redshift, stability, and quantum mechanical behavior on astrophysical scales. Exploring the interplay between galaxy formation physics and alterations to the standard CDM paradigm will be crucial in order to extract these signals from upcoming observations of ultrafaint galaxies, and forward-modeling approaches like the one developed here will drive these studies forward.

This paper has gone through internal review by the DES collaboration. Our code and subhalo catalogs are available at github.com/eonadler/subhalo Satellite connection; please contact the authors with additional data requests. We thank Ralf Kaehler for providing the visualizations of our simulations used in Figure 1. We thank Shea Garrison-Kimmel and Andrew Wetzel for sharing data from the ELVIS on FIRE simulations. We thank Susmita Adhikari, Arka Banerjee, Ekta Patel, and Andrew Pontzen for useful conversations. Finally, we are grateful to the anonymous referee for constructive feedback.
This research received support from the National Science Foundation (NSF) under grant No. NSF AST-1517422, grant No. NSF PHY17-48958 through the Kavli Institute for Theoretical Physics program “The Small-Scale Structure of Cold(?) Dark Matter,” and grant No. NSF DGE-1656518 through the NSF Graduate Research Fellowship received by E.O.N. Support for Y.-Y.M. was provided by the Pittsburgh Particle Physics, Astrophysics and Cosmology Center through the Samuel P. Langley PITT PACC Postdoctoral Fellowship and NASA through NASA Hubble Fellowship grant No. HST-HF2-51441.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.

This research made use of computational resources at SLAC National Accelerator Laboratory, a U.S. Department of Energy Office; the authors are thankful for the support of the SLAC computational team. This research made use of the Sherlock cluster at the Stanford Research Computing Center (SRCC); the authors are thankful for the support of the SRCC team. This research made use of https://arXiv.org and NASA's Astrophysics Data System for bibliographic information.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the NSF's National Optical-Infrared Astronomy Laboratory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

This work is based in part on observations at Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Laboratory, which is operated by the Association of

Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under grant Nos. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. The IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013), including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by the Fermi Research Alliance, LLC, under contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.

*Software:* Astropy (Astropy Collaboration et al. 2013), ChainConsumer (Hinton 2016), emcee (Foreman-Mackey et al. 2013), Healpy (healpy.readthedocs.io), Pandas (McKinney 2010), Scikit-Learn (Pedregosa et al. 2011), IPython (Pérez & Granger 2007), Jupyter (jupyter.org), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001), Seaborn (seaborn.pydata.org).

## Appendix A
### Galaxy–Halo Connection Model Details

In this appendix, we examine several components of our galaxy–halo connection model to determine whether any of our assumptions significantly impact the results presented above.

#### A.1. MW Host Halo Mass

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_{\text{MW}}$, to impact our constraints. Consider a toy model in which $N_{\text{sat}}$ satellites brighter than a limiting magnitude $M_{V,\text{min}}$ must be predicted in order to match the observed luminosity function. In this toy model, the predicted number of satellites is given by

$$N_{\text{sat}}(<M_{V,\text{lim}}) = \int_{M_{V,\text{lim}}}^{-\infty} \frac{dN_{\text{sat}}}{dM_{V}} dM_{V} \\ = \int_{M_{\text{min}}(M_{V,\text{lim}})}^{\infty} f\left(\frac{dN_{\text{sub}}}{dM_{\text{peak}}}, \theta\right) dM_{\text{peak}}, \quad (\text{A1})$$

where $dN_{\text{sub}}/dM_{\text{peak}}$ is the subhalo mass function, $M_{\text{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 $\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_{\text{sub}}/dM_{\text{peak}} \propto M_{\text{peak}}^{-2}$ (e.g.,
Mao et al. 2015), we have

$$N_{\text{sat}} \propto M_{\text{MW}} \int_{M_{\text{min}}}^{\infty} M_{\text{peak}}^{-2} dM_{\text{peak}} \propto \frac{M_{\text{MW}}}{M_{\text{min}}}. \quad (\text{A2})$$

Thus, for a fixed observed satellite count $N_{\text{sat}}$, we expect our 95% confidence level upper limit on $M_{\text{min}}$ to scale linearly with host mass. In addition, because the error on MW mass is independent of the error on $M_{\text{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 $M_{\text{min}} < 3.2 \times 10^8 M_{\odot}$ ($M_{\text{min}} < 2 \times 10^8 M_{\odot}$) for a maximally high-mass (maximally low-mass) host halo given the current $2\sigma$ observational uncertainty on the MW virial mass of $1.0 \times 10^{12} < M_{\text{MW}}/M_{\odot} < 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_{\text{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 < -13$ 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_{\text{MW}}$. We find that the total number of satellites inferred within the virial radius scales almost exactly linearly with $M_{\text{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 $M_{\text{min}}$ is roughly independent of $M_{\text{MW}}$ for our fiducial simulations.

## A.2. Mass-dependent Scatter

Here we test a model where the abundance-matching scatter in luminosity at fixed $V_{\text{peak}}$, $\sigma_M$, depends on peak halo mass. Motivated by the model in Garrison-Kimmel et al. (2017a), we set

$$\sigma_M \equiv \sigma_{M,0} - \gamma_M (\log M_{\text{peak}} - \log M_1), \quad (\text{A3})$$

where $\sigma_{M,0}$ is a free parameter that captures the amplitude of the luminosity scatter, $\gamma_M$ is a free parameter that captures its mass dependence, and $M_1 = 10^{11} M_{\odot}$ is fixed. By rerunning our fit with $\gamma_M$ as an additional ninth free parameter, we find that large values of $\gamma_M$ are ruled out by the DES and PS1 satellite populations at high statistical significance, with $\gamma_M < 0.07$ at 95% confidence. Large values of $\gamma_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 $\sigma_M$ apply to $\gamma_M$, so this

upper limit should be interpreted with caution. Introducing mass-dependent scatter does not significantly affect our inferred upper bound on $M_{50}$, implying that our fiducial minimum halo mass constraint does not depend on the details of our luminosity scatter model.

## A.3. Radial Scaling

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

$$r_{\text{sat}} \equiv \chi r_{\text{sub}}, \quad (\text{A4})$$

where $r_{\text{sat}}$ is a satellite's distance from the center of its host halo, which we equate to its galactocentric distance, and $r_{\text{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 $\chi = 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 $\chi$, we also show the 68% confidence interval for our fiducial posterior evaluated with $\chi = 0.5$.

To test the impact of radial scaling, we refit the DES and PS1 satellite populations with $\chi$ as an additional ninth free parameter. As expected, decreasing $\chi$ 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 $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 $\chi = 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 (~$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.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A1.** Radial distributions derived from our fit to the DES and PS1 satellite populations. Our fiducial eight-parameter galaxy occupation fraction model is shown in blue. Dark (light) blue bands correspond to 68% (95%) confidence intervals, dashed red lines show the 68% confidence interval for a model using host halos without LMC analogs (No LMC), and black lines show the observed radial distributions. Dotted-dashed blue lines show the 68% confidence interval for a model with a radial scaling parameter of $\chi = 0.5$.</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A2.** Size distributions derived by fitting to the DES and PS1 satellite populations. Our fiducial eight-parameter galaxy occupation fraction model is shown in blue. Dark (light) blue bands correspond to 68% (95%) confidence intervals, dashed red lines show the 68% confidence interval for a model with a concentration-dependent galaxy-halo size relation, and dotted-dashed blue lines show the 68% confidence interval for a model with an extreme dependence of satellite size on tidal stripping.</figcaption></figure>

### A.4. Tidal Stripping

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

$$r'_{1/2} \equiv r_{1/2} \left( \frac{V_{\max}}{V_{\text{acc}}} \right)^{\beta}, \quad (A5)$$

where $r'_{1/2}$ 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_{\text{acc}}$) is the maximum circular velocity of a subhalo today (at accretion), and $\beta > 0$ is a parameter that controls the strength of size reduction due to tidal stripping. We set $\beta = 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 < \beta < 2$ describes the results of high-resolution simulations well.

In Figure A2, we illustrate predicted size distributions for our fiducial posterior evaluated with $\beta = 3$; a large value of $\beta$ 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.

### A.5. Concentration-dependent Satellite Sizes

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

$$r_{1/2} \equiv A \left(\frac{c}{10}\right)^{\gamma} \left(\frac{R_{\text{vir}}}{R_0}\right)^n, \qquad (\text{A6})$$

with $A = 0.02$, $n = 1$, $\gamma = -0.7$, $R_0 = 1$ 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 $\gamma$ 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_{\text{logR}}$ (0.88 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 $\gamma$, and our analysis does not place an upper limit on $A$ in this case. Here $\gamma$ 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 $\gamma = -0.7$ in Figure A2).

### A.6. Orphan Satellite Contribution

To test the importance of orphan satellites, we refit the DES and PS1 satellite populations with $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 $O$. In particular, the 95% confidence upper limit on $M_{50}$ increases by less than $1\sigma$ 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 \pm 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.

## Appendix B
Statistical Framework Details

Next, we provide additional details on our statistical framework, and we discuss several caveats.

### B.1. Poisson Process Likelihood

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, \dots, \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

$$P(n_i|\hat{n}_{i,1}, \dots, \hat{n}_{i,\hat{N}}) = \left(\frac{\hat{N} + 1}{\hat{N}}\right)^{-(\hat{n}_{i,1} + \dots + \hat{n}_{i,\hat{N}} + 1)} \times (\hat{N} + 1)^{-n_i} \frac{\Gamma(\hat{n}_{i,1} + \dots + \hat{n}_{i,\hat{N}} + n_i + 1)}{\Gamma(n_i + 1)\Gamma(\hat{n}_{i,1} + \dots + \hat{n}_{i,\hat{N}} + 1)}, \quad (\text{B1})$$

where the dependence on galaxy-halo connection model parameters $\theta$ is implicit, and we assumed (i) a flat prior on $\lambda_i$ for $\lambda_i \ge 0$ and (ii) that $n_i$ and all $\hat{n}_{i,\nu}$ are drawn from the same Poisson distribution with rate parameter $\lambda_i$. Note that our method yields noninteger numbers of mock satellites by counting each system as $p_{\text{detect}} \times (1 - p_{\text{disrupt}}) \times f_{\text{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.

### B.2. Priors

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 ~0.25 dex. For $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 $M_{50}$ posterior is flat below ~5 × $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.<fnref n="63" /> Priors for $\mathcal{B}$ and $n$ are set based on studies that identify the preferred values of these parameters, and priors for $\sigma_{\text{gal}}$, $A$, and $\sigma_{\text{logR}}$ are chosen to be uniform with conservative upper bounds.

### B.3. Bayes Factor Calculation

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 $\theta$ within a fixed Mahalanobis distance of a point $\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.

### B.4. Caveats and Future Work

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

<footnote n="63">⁶³ However, as noted in Section 7.4, our reported Bayes factors are independent of this choice.</footnote>
<table><caption>**Table B1**
Prior Distributions for the Parameters Varied in Our Fiducial Eight-parameter Fit to the DES and PS1 Satellite Populations</caption><thead><tr><th>Free Parameter</th><th>Prior Distribution</th><th>Motivation</th></tr></thead><tbody><tr><td>Faint-end slope</td><td>arctan α ~ unif(-1.1, -0.9)</td><td>Jeffreys prior for -2 &lt; α &lt; -1.2</td></tr><tr><td>Luminosity scatter</td><td>σ<sub>M</sub> ~ unif(0, 2) dex</td><td>Conservative upper limit (Garrison-Kimmel et al. 2017a; Lehmann et al. 2017)</td></tr><tr><td>50% occupation mass</td><td>log(M<sub>50</sub>/M<sub>☉</sub>) ~ unif(7.5, 11)</td><td>Lower limit corresponds to simulation resolution limit (Mao et al. 2015)</td></tr><tr><td>Baryonic effects</td><td>ln(B) ~ N(μ = 1, σ = 0.5)</td><td>Hydrodynamic simulations (Nadler et al. 2018, 2019b)</td></tr><tr><td>Occupation scatter</td><td>σ<sub>gal</sub> ~ unif(0, 1) dex</td><td>Hydrodynamic simulations (Fitts et al. 2018; Graus et al. 2019)</td></tr><tr><td>Size amplitude</td><td>A ~ unif(0, 0.5) kpc</td><td>Empirical galaxy-halo size relation (Kravtsov 2013)</td></tr><tr><td>Size scatter</td><td>σ<sub>log R</sub> ~ unif(0, 2) dex</td><td>Empirical galaxy-halo size relation (Kravtsov 2013)</td></tr><tr><td>Size power-law index</td><td>n ~ N(μ = 1, σ = 0.5)</td><td>Empirical galaxy-halo size relation (Kravtsov 2013)</td></tr></tbody></table>

Note. Here $N(μ, σ)$ denotes a normal distribution with mean $μ$ and standard deviation $σ$.

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 $T$ of observable space $x$ is given by a Poisson distribution with rate $λ_T = \int_T λ(x)dx$, where $λ(x)$ is referred to as the “rate function.” Given a rate function $λ(x)$, the likelihood of observing $N$ events at a set of points $\{x_i\}_{i=1}^N$ is

$$p(\{x_i\}_{i=1}^N | \lambda) = \exp \left[ - \int_T \lambda(x) dx \right] \prod_{i=1}^N \lambda(x_i), \quad (B2)$$

where we suppressed the dependence of the rate on our galaxy-halo connection model parameters $θ$. In our case, the “events” $\{x_i\}_{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 $λ$,

$$p(\{x_i\} | \{\hat{x}_j\}) = \frac{\int \mathcal{D}\lambda \, p(\{x_i\} | \lambda) p(\{\hat{x}_j\} | \lambda) p(\lambda)}{\int \mathcal{D}\lambda \, p(\{\hat{x}_j\} | \lambda) p(\lambda)}. \quad (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(λ)$, 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

$$p(\{n_j\}_{j=1}^N | \{\lambda_j\}) = \exp \left( - \sum_{\text{bins } j} \lambda_j V_j \right) \prod_{\text{bins } k} \frac{\lambda_k^{n_k}}{n_k!}, \quad (B4)$$

where $λ_j$ is the rate in bin $j$, $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 $λ(x)$ 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.

## Appendix C
### Robustness to Observational Systematics

We now present a set of tests in order to verify the robustness of our key results to various observational systematics.

#### C.1. Kinematically Unconfirmed Satellites

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
<figure><img src="image_17.png" /><figcaption>**Figure C1.** Posterior distribution from our fit to the kinematically confirmed DES and PS1 satellite populations. Dark (light) shaded contours represent 68% (95%) confidence intervals. Shaded areas in the marginal distributions and parameter summaries correspond to 68% confidence intervals. Note that $\sigma_M$, $\sigma_{gal}$, and $\sigma_{logR}$ are reported in dex, $M_{50}$ is reported as $\log(M_{50}/M_{\odot})$, and $A$ is reported in pc. Note that $\sigma_{gal}$ is not constrained at 68% confidence in this fit.</figcaption></figure>

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 $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 ~1σ to $150 \pm 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.
<table><caption>**Table C1**
MW Satellites Used in Our Analysis</caption><thead><tr><th>Name</th><th><i>M</i><sub>V</sub><br>(mag)</th><th><i>D</i><br>(kpc)</th><th><i>r</i><sub>1/2</sub><br>(pc)</th></tr></thead><tbody><tr><td></td><td>DES</td><td></td><td></td></tr><tr><td>Fornax</td><td>-13.46</td><td>147</td><td>707</td></tr><tr><td>Sculptor</td><td>-10.82</td><td>84</td><td>223</td></tr><tr><td>Reticulum II</td><td>-3.88</td><td>30</td><td>31</td></tr><tr><td>Eridanus II<sup>a</sup></td><td>-7.21</td><td>380</td><td>158</td></tr><tr><td>Tucana II</td><td>-3.8</td><td>58</td><td>165</td></tr><tr><td>Grus II<sup>*</sup></td><td>-3.9</td><td>53</td><td>92</td></tr><tr><td>Horologium I</td><td>-3.55</td><td>79</td><td>31</td></tr><tr><td>Tucana III<sup>*</sup></td><td>-2.4</td><td>25</td><td>44</td></tr><tr><td>Tucana IV</td><td>-3.5</td><td>48</td><td>128</td></tr><tr><td>Phoenix II</td><td>-3.30</td><td>83</td><td>21</td></tr><tr><td>Horologium II<sup>*</sup></td><td>-2.6</td><td>78</td><td>33</td></tr><tr><td>Tucana V<sup>*</sup></td><td>-1.6</td><td>55</td><td>16</td></tr><tr><td>Pictor I<sup>*</sup></td><td>-3.45</td><td>114</td><td>18</td></tr><tr><td>Columba I<sup>*</sup></td><td>-4.2</td><td>183</td><td>98</td></tr><tr><td>Cetus II<sup>*</sup></td><td>0.02</td><td>30</td><td>17</td></tr><tr><td>Grus I<sup>*</sup></td><td>-3.47</td><td>120</td><td>21</td></tr><tr><td>Reticulum III<sup>*</sup></td><td>-3.31</td><td>92</td><td>64</td></tr><tr><td></td><td>PS1</td><td></td><td></td></tr><tr><td>Leo I</td><td>-11.78</td><td>254</td><td>226</td></tr><tr><td>Leo II</td><td>-9.74</td><td>233</td><td>165</td></tr><tr><td>Draco</td><td>-8.71</td><td>76</td><td>180</td></tr><tr><td>Ursa Minor</td><td>-9.03</td><td>76</td><td>272</td></tr><tr><td>Sextans</td><td>-8.72</td><td>86</td><td>345</td></tr><tr><td>Canes Venatici I</td><td>-8.8</td><td>218</td><td>338</td></tr><tr><td>Boötes I</td><td>-6.02</td><td>66</td><td>160</td></tr><tr><td>Ursa Major II</td><td>-4.25</td><td>32</td><td>85</td></tr><tr><td>Coma Berenices</td><td>-4.38</td><td>44</td><td>57</td></tr><tr><td>Sagittarius II</td><td>-5.2</td><td>69</td><td>32</td></tr><tr><td>Willman 1</td><td>-2.53</td><td>38</td><td>20</td></tr><tr><td>Canes Venatici II</td><td>-5.17</td><td>160</td><td>55</td></tr><tr><td>Segue 1</td><td>-1.30</td><td>23</td><td>20</td></tr><tr><td>Segue 2<sup>*</sup></td><td>-1.86</td><td>35</td><td>34</td></tr><tr><td>Crater II</td><td>-8.2</td><td>117</td><td>1066</td></tr><tr><td>Draco II<sup>*</sup></td><td>-0.8</td><td>22</td><td>17</td></tr><tr><td>Triangulum II<sup>*</sup></td><td>-1.60</td><td>30</td><td>13</td></tr><tr><td>Hercules</td><td>-5.83</td><td>132</td><td>120</td></tr><tr><td>Cetus II<sup>b</sup></td><td>0.02</td><td>30</td><td>17</td></tr></tbody></table>

**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.

<sup>a</sup> Eridanus II is not included because it lies outside our fiducial 300 kpc heliocentric distance cut.

<sup>b</sup> Cetus II is detected in both PS1 and DES; in our analysis, we only count this system in the observed DES population.

## C.2. Satellite Size Criterion

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 $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.

## C.3. Biases in Measured Satellite Properties

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 $\Delta M_V = +1$ 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 $\Delta M_V = -1$ 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 $\Delta M_V = -1$ mag ($\Delta M_V = +1$ mag); however, the total predicted number of MW satellites with $M_V < 0$ mag and $r_{1/2} > 10$ pc and our 95% confidence upper limit on $M_{50}$ are not significantly affected in either case.

## Appendix D
Resolution and Sample Variance

To assess the impact of resolution effects on our fiducial simulations and results, we compare the subhalo maximum circular velocity function, radial distribution, and size distribution from one of our fiducial host halos (excluding LMC satellites) to those from a higher-resolution resimulation of the same host. In particular, we resimulate this halo with a $4 \times 10^4 M_{\odot} h^{-1}$ high-resolution particle mass and an $85$ pc $h^{-1}$ minimum softening length. We find that the distributions of all relevant subhalo properties are not significantly affected above the resolution limit of our fiducial simulations. Moreover, by rerunning our analysis, we find that none of our galaxy-halo connection model constraints are significantly affected when using a higher-resolution simulation.

We also assess the impact of sample variance on our fiducial subhalo and satellite populations, since the final positions of LMC satellites might be sensitive to the realizations of small-scale density fluctuations in our fiducial simulations. In particular, we resimulate both of our fiducial host halos at standard resolution with different random seeds for small-scale phases in the matter power spectrum below $60$ kpc $h^{-1}$. We find that the properties of the MW host halo and LMC halo are not significantly affected in these resimulations, and that the resulting subhalo populations are nearly identical in terms of their distributions of $M_{\text{peak}}$, $V_{\text{peak}}$, halo size at accretion, and present-day heliocentric distance, implying that our results are robust to sample variance in the phases of the matter power spectrum on small scales.

## Appendix E
Observed Satellite Data Vectors

The confirmed and candidate DES and PS1 satellites that pass the detection criteria defined in Paper I are listed in Table C1. Note that, although Kim 2 (DES) and Laevens 1 (PS1) formally pass these detection criteria, we do not include these systems in our analysis or Table C1 because they are suspected to be star clusters (Paper I).

## ORCID iDs

E. O. Nadler https://orcid.org/0000-0002-1182-3825

R. H. Wechsler https://orcid.org/0000-0003-2229-011X

K. Bechtol https://orcid.org/0000-0001-8156-0429
Y.-Y. Mao https://orcid.org/0000-0002-1200-0820

G. Green https://orcid.org/0000-0001-5417-2260

A. Drlica-Wagner https://orcid.org/0000-0001-8251-933X

M. McNanna https://orcid.org/0000-0001-5435-7820

S. Mau https://orcid.org/0000-0003-3519-4004

A. B. Pace https://orcid.org/0000-0002-6021-8760

A. Kravtsov https://orcid.org/0000-0003-4307-634X

S. Dodelson https://orcid.org/0000-0002-8446-3859

A. H. Riley https://orcid.org/0000-0001-5805-5766

M. Y. Wang https://orcid.org/0000-0002-8226-6237

S. Allam https://orcid.org/0000-0002-7069-7857

J. Annis https://orcid.org/0000-0002-0609-3987

A. Carnero Rosell https://orcid.org/0000-0003-3044-5150

M. Carrasco Kind https://orcid.org/0000-0002-4802-3194

J. Carretero https://orcid.org/0000-0002-3130-0204

J. De Vicente https://orcid.org/0000-0001-8318-6813

S. Desai https://orcid.org/0000-0002-0466-3288

A. E. Evrard https://orcid.org/0000-0002-4876-956X

J. Frieman https://orcid.org/0000-0003-4079-3263

J. García-Bellido https://orcid.org/0000-0002-9370-8360

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

G. Gutierrez https://orcid.org/0000-0003-0825-0517

K. Kuehn https://orcid.org/0000-0003-0120-0808

N. Kuropatkin https://orcid.org/0000-0003-2511-0946

J. L. Marshall https://orcid.org/0000-0003-0710-9474

F. Menanteau https://orcid.org/0000-0002-1372-2534

R. Miquel https://orcid.org/0000-0002-6610-4836

A. A. Plazas https://orcid.org/0000-0002-2598-0514

A. K. Romer https://orcid.org/0000-0002-9328-879X

E. Sanchez https://orcid.org/0000-0002-9646-8198

M. Smith https://orcid.org/0000-0002-3321-1432

M. Soares-Santos https://orcid.org/0000-0001-6082-8529

E. Suchyta https://orcid.org/0000-0002-7047-9358

G. Tarle https://orcid.org/0000-0003-1704-0781

A. R. Walker https://orcid.org/0000-0002-7123-8943

## References

Adams, R. P., Murray, I., & MacKay, D. J. C. 2009, in Proc. 26th Annual Int. Conf. on Machine Learning, ICML '09 (New York: ACM), 9

Agertz, O., Pontzen, A., Read, J. I., et al. 2020, MNRAS, 491, 1656

Alvarez, M. A., Busha, M., Abel, T., & Wechsler, R. H. 2009, ApJL, 703, L167

Amorisco, N. C. 2019, MNRAS, 489, L22

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bechtol, K., Drlica-Wagner, A., Balbinot, E., et al. 2015, ApJ, 807, 50

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013a, ApJ, 770, 57

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013b, ApJ, 762, 109

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013c, ApJ, 763, 18

Belokurov, V., Walker, M. G., Evans, N. W., et al. 2008, ApJ, 686, L83

Belokurov, V., Walker, M. G., Evans, N. W., et al. 2009, MNRAS, 397, 1748

Belokurov, V., Walker, M. G., Evans, N. W., et al. 2010, ApJL, 712, L103

Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJL, 647, L111

Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2007, ApJ, 654, 897

Besla, G. 2015, arXiv:1511.03346

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bose, S., Deason, A. J., Belokurov, V., & Frenk, C. S. 2019, arXiv:1909.04039

Bose, S., Deason, A. J., & Frenk, C. S. 2018, ApJ, 863, 123

Boylan-Kolchin, M., Springel, V., White, S. D. M., & Jenkins, A. 2010, MNRAS, 406, 896

Boylan-Kolchin, M., Weisz, D. R., Johnson, B. D., et al. 2015, MNRAS, 453, 1503

Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, ApJ, 765, 22

Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80

Busha, M. T., Alvarez, M. A., Wechsler, R. H., Abel, T., & Strigari, L. E. 2010, ApJ, 710, 408

Busha, M. T., Marshall, P. J., Wechsler, R. H., Klypin, A., & Primack, J. 2011, ApJ, 743, 40

Callingham, T. M., Cautun, M., Deason, A. J., et al. 2019, MNRAS, 484, 5453

Cautun, M., Benitez-Llambay, A., Deason, A. J., et al. 2019b, arXiv:1911.04557

Cautun, M., Deason, A. J., Frenk, C. S., & McAlpine, S. 2019a, MNRAS, 483, 2185

Côté, B., Silvia, D. W., O'Shea, B. W., Smith, B., & Wise, J. H. 2018, ApJ, 859, 67

D'Onghia, E., & Lake, G. 2008, ApJL, 686, L61

Deason, A. J., Wetzel, A. R., Garrison-Kimmel, S., & Belokurov, V. 2015, MNRAS, 453, 3568

Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415

Dooley, G. A., Peter, A. H. G., Carlin, J. L., et al. 2017, MNRAS, 472, 1060

Drlica-Wagner, A., Bechtol, K., Allam, S., et al. 2016, ApJL, 833, L5

Drlica-Wagner, A., Bechtol, K., Mau, S., et al. 2020, ApJ, 893, 47

Drlica-Wagner, A., Bechtol, K., Rykoff, E. S., et al. 2015, ApJ, 813, 109

Drlica-Wagner, A., Mao, Y.-Y., Adhikari, S., et al. 2019, arXiv:1902.01055

Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS, 487, 2685

Erkal, D., & Belokurov, V. A. 2019, MNRAS, 487, 2685

Errani, R., Penarrubia, J., & Tormen, G. 2015, MNRAS, 449, L46

Fattahi, A., Navarro, J. F., Frenk, C. S., et al. 2018, MNRAS, 476, 3816

Fielder, C. E., Mao, Y.-Y., Newman, J. A., Zentner, A. R., & Licquia, T. C. 2019, MNRAS, 486, 4545

Fitts, A., Boylan-Kolchin, M., Bullock, J. S., et al. 2018, MNRAS, 479, 319

Flaxman, S., Wilson, A., Neill, D., Nickisch, H., & Smola, A. 2015, in Proc. Machine Learning Research, Vol. 37, Proc. 32nd Int. Conf. on Machine Learning, ed. F. Bach & D. Blei (Lille: PMLR), 607

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Garrison-Kimmel, S., Bullock, J. S., Boylan-Kolchin, M., & Bardwell, E. 2017a, MNRAS, 464, 3108

Garrison-Kimmel, S., Hopkins, P. F., Wetzel, A., et al. 2019, MNRAS, 487, 1380

Garrison-Kimmel, S., Wetzel, A., Bullock, J. S., et al. 2017b, MNRAS, 471, 1709

Geha, M., Wechsler, R. H., Mao, Y.-Y., et al. 2017, ApJ, 847, 4

Graus, A. S., Bullock, J. S., Kelley, T., et al. 2019, MNRAS, 488, 4585

Grillmair, C. J. 2006, ApJL, 645, L37

Grillmair, C. J. 2009, ApJ, 693, 1118

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101

Hargis, J. R., Willman, B., & Peter, A. H. G. 2014, ApJL, 795, L13

Hinton, S. R. 2016, JOSS, 1, 00045

Homma, D., Chiba, M., Komiyama, Y., et al. 2019, PASJ, 71, 94

Homma, D., Chiba, M., Okamoto, S., et al. 2016, ApJ, 832, 21

Homma, D., Chiba, M., Okamoto, S., et al. 2018, PASJ, 70, S18

Hunter, J. D. 2007, CSE, 9, 90

Irwin, M. J., Belokurov, V., Evans, N. W., et al. 2007, ApJL, 656, L13

Jahn, E. D., Sales, L. V., Wetzel, A., et al. 2019, MNRAS, 489, 5348

Jethwa, P., Erkal, D., & Belokurov, V. 2016, MNRAS, 461, 2212

Jethwa, P., Erkal, D., & Belokurov, V. 2018, MNRAS, 473, 2060

Jiang, F., Dekel, A., Kneller, O., et al. 2019, MNRAS, 488, 4801

John, S. T., & Hensman, J. 2018, arXiv:1804.01016

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, v. 1.2.3 (USA: Trelgol Publishing), https://scipy.org/

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2008, ApJ, 873, 111

Kallivayalil, N., Sales, L. V., Zivick, P., et al. 2018, ApJ, 867, 19

Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J., & Alcock, C. 2013, ApJ, 764, 161

Katz, H., Ramsoy, M., Rosdahl, J., et al. 2019, arXiv:1905.11414

Kelley, T., Bullock, J. S., Garrison-Kimmel, S., et al. 2019, MNRAS, 487, 4409

Kennedy, R., Frenk, C., Cole, S., & Benson, A. 2014, MNRAS, 442, 2487

Kim, D., & Jerjen, H. 2015, ApJL, 808, L39

Kim, S. Y., Peter, A. H. G., & Hargis, J. R. 2018, PhRvL, 121, 211302

Koposov, S. E., Belokurov, V., Torrealba, G., & Evans, N. W. 2015, ApJ, 805, 130

Koposov, S. E., Walker, M. G., Belokurov, V., et al. 2018, MNRAS, 479, 5343

Kravtsov, A. V. 2013, ApJL, 764, L31

Laevens, B. P. M., Martin, N. F., Bernard, E. J., et al. 2015a, ApJ, 813, 44

Laevens, B. P. M., Martin, N. F., Ibata, R. A., et al. 2015b, ApJL, 802, L18

Lehmann, B. V., Mao, Y.-Y., Becker, M. R., Skillman, S. W., & Wechsler, R. H. 2017, ApJ, 834, 37

Li, Z.-Z., Qian, Y.-Z., Han, J., Wang, W., & Jing, Y. P. 2019a, ApJ, 886, 69

Li, Z.-Z., Qian, Y.-Z., Han, J., et al. 2019b, arXiv:1912.02086

Loveday, J., Norberg, P., Baldry, I. K., et al. 2015, MNRAS, 451, 1540
Lovell, M. R., Frenk, C. S., Eke, V. R., et al. 2014, MNRAS, 439, 300

Lovell, M. R., Gonzalez-Perez, V., Bose, S., et al. 2017, MNRAS, 468, 2836

Lu, Y., Benson, A., Mao, Y.-Y., et al. 2016, ApJ, 830, 59

Lunnan, R., Vogelsberger, M., Frebel, A., et al. 2012, ApJ, 746, 109

Luque, E., Queiroz, A., Santiago, B., et al. 2016, MNRAS, 458, 603

Lynden-Bell, D. 1976, MNRAS, 174, 695

Macciò, A. V., & Fontanot, F. 2010, MNRAS, 404, L16

Mao, Y.-Y., Williamson, M., & Wechsler, R. H. 2015, ApJ, 810, 21

Martin, N. F., Nidever, D. L., Besla, G., et al. 2015, ApJL, 804, L5

Mau, S., Cerny, W., Pace, A. B., et al. 2019, arXiv:1912.03301

McKinney, W. 2010, in Proc. 9th Python in Science Conf., ed. S. van der Walt & J. Millman, 56

Munshi, F., Brooks, A. M., Christensen, C., et al. 2019, ApJ, 874, 40

Murray, I. 2007, PhD thesis, Univ. College London

Nadler, E. O., Gluscevic, V., Boddy, K. K., & Wechsler, R. H. 2019a, ApJL, 878, L32

Nadler, E. O., Mao, Y.-Y., Green, G. M., & Wechsler, R. H. 2019b, ApJ, 873, 34

Nadler, E. O., Mao, Y.-Y., Wechsler, R. H., Garrison-Kimmel, S., & Wetzel, A. 2018, ApJ, 859, 129

Newton, O., Cautun, M., Jenkins, A., Frenk, C. S., & Helly, J. C. 2018, MNRAS, 479, 2853

Pardy, S. A., D'Onghia, E., Navarro, J. F., et al. 2020, MNRAS, 492, 1543

Patel, E., Besla, G., & Mandel, K. 2017, MNRAS, 468, 3428

Patel, E., Kallivayalil, N., Garavito-Camargo, N., et al. 2020, ApJ, in press

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825

Peñarrubia, J., Gómez, F. A., Besla, G., Erkal, D., & Ma, Y.-Z. 2016, MNRAS, 456, L54

Peñarrubia, J., Navarro, J. F., McConnachie, A. W., & Martin, N. F. 2009, ApJ, 698, 222

Pérez, F., & Granger, B. E. 2007, CSE, 9, 21

Reddick, R. M., Wechsler, R. H., Tinker, J. L., & Behroozi, P. S. 2013, ApJ, 771, 30

Rey, M. P., Pontzen, A., Agertz, O., et al. 2019, ApJL, 886, L3

Rodríguez-Puebla, A., Avila-Reese, V., & Drory, N. 2013, ApJ, 773, 172

Sakamoto, T., & Hasegawa, T. 2006, ApJL, 653, L29

Sales, L. V., Navarro, J. F., Kallivayalil, N., & Frenk, C. S. 2017, MNRAS, 465, 1879

Samuel, J., Wetzel, A., Tollerud, E., et al. 2020, MNRAS, 491, 1471

Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016a, MNRAS, 457, 1931

Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016b, MNRAS, 456, 85

Shao, S., Cautun, M., Deason, A. J., Frenk, C. S., & Theuns, T. 2018, MNRAS, 479, 284

Simon, J., Bechtol, K., Drlica-Wagner, A., et al. 2019, BAAS, 51, 409

Simpson, C. M., Grand, R. J. J., Gómez, F. A., et al. 2018, MNRAS, 478, 548

Tollerud, E. J., Bullock, J. S., Strigari, L. E., & Willman, B. 2008, ApJ, 688, 277

Tollerud, E. J., & Peek, J. E. G. 2018, ApJ, 857, 45

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2018, MNRAS, 475, 5085

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2019, MNRAS, 488, 2743

Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016a, MNRAS, 459, 2370

Torrealba, G., Koposov, S. E., Belokurov, V., et al. 2016b, MNRAS, 463, 712 van den Bosch, F. C., & Ogiya, G. 2018, MNRAS, 475, 4066

van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Walsh, S. M., Jerjen, H., & Willman, B. 2007, ApJL, 662, L83

Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2019, ApJ, 873, 118

Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435

Weisz, D. R., & Boylan-Kolchin, M. 2017, MNRAS, 469, L83

Wetzel, A. R., Hopkins, P. F., Kim, J.-H., et al. 2016, ApJL, 827, L23

Wheeler, C., Hopkins, P. F., Pace, A. B., et al. 2019, MNRAS, 490, 4447

Willman, B., Blanton, M. R., West, A. A., et al. 2005a, AJ, 129, 2692

Willman, B., Dalcanton, J. J., Martinez-Delgado, D., et al. 2005b, ApJL, 626, L85

Wright, A. H., Robotham, A. S. G., Driver, S. P., et al. 2017, MNRAS, 470, 283

Zhang, D., Luo, Y., & Kang, X. 2019, MNRAS, 486, 2440

Zucker, D. B., Belokurov, V., Evans, N. W., et al. 2006, ApJL, 650, L41