# Milky Way Satellite Census. IV. Constraints on Decaying Dark Matter from Observations of Milky Way Satellite Galaxies

S. Mau, E. O. Nadler, R. H. Wechsler, A. Drlica-Wagner, K. Bechtol, G.

Green, D. Huterer, T. S. Li, Y. -Y. Mao, C. E. Martínez-Vázquez, et al.

► To cite this version:

S. Mau, E. O. Nadler, R. H. Wechsler, A. Drlica-Wagner, K. Bechtol, et al. Milky Way Satellite Census. IV. Constraints on Decaying Dark Matter from Observations of Milky Way Satellite Galaxies. The Astrophysical Journal, 2022, 932, 10.3847/1538-4357/ac6e65. insu-03748259

HAL Id: insu-03748259

https://insu.hal.science/insu-03748259

Submitted on 9 Aug 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
Milky Way Satellite Census. IV. Constraints on Decaying Dark Matter from Observations of Milky Way Satellite Galaxies

S. Mau<fnref n="1" />,<fnref n="2" /> ID, E. O. Nadler<fnref n="3" />,<fnref n="4" /> ID, R. H. Wechsler<fnref n="1" />,<fnref n="2" />,<fnref n="5" /> ID, A. Drlica-Wagner<fnref n="6" />,<fnref n="7" />,<fnref n="8" /> ID, K. Bechtol<fnref n="9" /> ID, G. Green<fnref n="10" /> ID, D. Huterer<fnref n="11" /> ID, T. S. Li<fnref n="12" /> ID, Y.-Y. Mao<fnref n="13" />,⁵³ ID, C. E. Martínez-Vázquez<fnref n="14" /> ID, M. McNanna<fnref n="9" /> ID, B. Mutlu-Pakdil<fnref n="8" /> ID, A. B. Pace<fnref n="15" /> ID, A. Peter<fnref n="16" /> ID, A. H. Riley<fnref n="17" /> ID, L. Strigari<fnref n="17" /> ID, M.-Y. Wang<fnref n="15" /> ID, M. Aguena<fnref n="18" /> ID, S. Allam<fnref n="7" /> ID, J. Annis<fnref n="7" /> ID, D. Bacon<fnref n="19" /> ID, E. Bertin<fnref n="20" />,<fnref n="21" /> ID, S. Bocquet<fnref n="22" /> ID, D. Brooks<fnref n="23" /> ID, D. L. Burke<fnref n="2" />,<fnref n="5" /> ID, A. Carnero Rosell<fnref n="18" />,<fnref n="24" />,<fnref n="25" /> ID, M. Carrasco Kind<fnref n="26" />,<fnref n="27" /> ID, J. Carretero<fnref n="28" /> ID, M. Costanzi<fnref n="29" />,<fnref n="30" />,<fnref n="31" /> ID, M. Crocce<fnref n="32" />,<fnref n="33" /> ID, M. E. S. Pereira<fnref n="11" />,<fnref n="34" /> ID, T. M. Davis<fnref n="35" /> ID, J. De Vicente<fnref n="36" /> ID, S. Desai<fnref n="37" /> ID, P. Doel<fnref n="23" />, I. Ferrero<fnref n="38" /> ID, B. Flaugher<fnref n="7" /> ID, J. Frieman<fnref n="7" />,<fnref n="8" /> ID, J. García-Bellido<fnref n="39" /> ID, M. Gatti<fnref n="40" /> ID, G. Giannini<fnref n="28" /> ID, D. Gruen<fnref n="22" />,<fnref n="41" /> ID, R. A. Gruendl<fnref n="26" />,<fnref n="27" /> ID, J. Gschwend<fnref n="18" />,<fnref n="42" /> ID, G. Gutierrez<fnref n="7" /> ID, S. R. Hinton<fnref n="35" /> ID, D. L. Hollowood<fnref n="43" /> ID, K. Honscheid<fnref n="16" />,<fnref n="44" /> ID, D. J. James<fnref n="45" /> ID, K. Kuehn<fnref n="46" />,<fnref n="47" /> ID, O. Lahav<fnref n="23" /> ID, M. A. G. Maia<fnref n="18" />,<fnref n="42" /> ID, J. L. Marshall<fnref n="17" /> ID, R. Miquel<fnref n="28" />,<fnref n="48" /> ID, J. J. Mohr<fnref n="22" />,<fnref n="49" /> ID, R. Morgan<fnref n="9" /> ID, R. L. C. Ogando<fnref n="42" /> ID, F. Paz-Chinchón<fnref n="26" />,<fnref n="50" /> ID, A. Pieres<fnref n="18" />,<fnref n="42" /> ID, M. Rodriguez-Monroy<fnref n="36" /> ID, E. Sanchez<fnref n="36" /> ID, V. Scarpine<fnref n="7" /> ID, S. Serrano<fnref n="32" />,<fnref n="33" /> ID, I. Sevilla-Noarbe<fnref n="36" /> ID, E. Suchyta⁵¹ ID, G. Tarle<fnref n="11" /> ID, C. To<fnref n="16" /> ID, D. L. Tucker<fnref n="7" /> ID, and J. Weller<fnref n="49" />,⁵² ID

(DES Collaboration)

<footnote n="1">¹Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA; smau@stanford.edu</footnote>

<footnote n="2">²Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="3">³Carnegie Observatories, 813 Santa Barbara Street, Pasadena, CA 91101, USA; enadler@carnegiescience.edu</footnote>

<footnote n="4">⁴Department of Physics & Astronomy, University of Southern California, Los Angeles, CA, 90007, USA</footnote>

<footnote n="5">⁵SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="6">⁶Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="7">⁷Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="8">⁸Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="9">⁹Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue Madison, WI 53706-1390, USA</footnote>

<footnote n="10">¹⁰Max Planck Institute for Astronomy, Königstuhl 17, D-69117, Heidelberg, Germany</footnote>

<footnote n="11">¹¹Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="12">¹²Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto ON, M5S 3H4, Canada</footnote>

<footnote n="13">¹³Department of Physics and Astronomy, Rutgers, The State University of New Jersey, Piscataway, NJ 08854, USA</footnote>

<footnote n="14">¹⁴Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile</footnote>

<footnote n="15">¹⁵Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15312, USA</footnote>

<footnote n="16">¹⁶Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="17">¹⁷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</footnote>

<footnote n="18">¹⁸Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil</footnote>

<footnote n="19">¹⁹Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK</footnote>

<footnote n="20">²⁰CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="21">²¹Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="22">²²Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="23">²³Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK</footnote>

<footnote n="24">²⁴Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="25">²⁵Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="26">²⁶Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA</footnote>

<footnote n="27">²⁷Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA</footnote>

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

<footnote n="29">²⁹Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="30">³⁰INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="31">³¹Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

<footnote n="32">³²Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="33">³³Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="34">³⁴Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="35">³⁵School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

<footnote n="36">³⁶Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain</footnote>

<footnote n="37">³⁷Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="38">³⁸Institute of Theoretical Astrophysics, University of Oslo. P.O. Box 1029 Blindern, NO-0315 Oslo, Norway</footnote>

<footnote n="39">³⁹Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="40">⁴⁰Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="41">⁴¹Excellence Cluster Origins, Boltzmannstr. 2, D-85748 Garching, Germany</footnote>

<footnote n="42">⁴²Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil</footnote>

<footnote n="43">⁴³Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="44">⁴⁴Department of Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="45">⁴⁵Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="46">⁴⁶Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia</footnote>

<footnote n="47">⁴⁷Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA</footnote>

<footnote n="48">⁴⁸Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="49">⁴⁹Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany</footnote>

<footnote n="50">⁵⁰Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>
⁵¹ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

⁵² Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, D-81679 München, Germany
*Received 2022 January 27; revised 2022 May 4; accepted 2022 May 9; published 2022 June 27*

**Abstract**

We use a recent census of the Milky Way (MW) satellite galaxy population to constrain the lifetime of particle dark matter (DM). We consider two-body decaying dark matter (DDM) in which a heavy DM particle decays with lifetime $\tau$ comparable to the age of the universe to a lighter DM particle (with mass splitting $\epsilon$) and to a dark radiation species. These decays impart a characteristic “kick velocity,” $V_{\text{kick}} = \epsilon c$, on the DM daughter particles, significantly depleting the DM content of low-mass subhalos and making them more susceptible to tidal disruption. We fit the suppression of the present-day DDM subhalo mass function (SHMF) as a function of $\tau$ and $V_{\text{kick}}$ using a suite of high-resolution zoom-in simulations of MW-mass halos, and we validate this model on new DDM simulations of systems specifically chosen to resemble the MW. We implement our DDM SHMF predictions in a forward model that incorporates inhomogeneities in the spatial distribution and detectability of MW satellites and uncertainties in the mapping between galaxies and DM halos, the properties of the MW system, and the disruption of subhalos by the MW disk using an empirical model for the galaxy-halo connection. By comparing to the observed MW satellite population, we conservatively exclude DDM models with $\tau < 18$ Gyr (29 Gyr) for $V_{\text{kick}} = 20 \text{ km s}^{-1}$ (40 km s$^{-1}$) at 95% confidence. These constraints are among the most stringent and robust small-scale structure limits on the DM particle lifetime and strongly disfavor DDM models that have been proposed to alleviate the Hubble and $S_8$ tensions.

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

**1. Introduction**

The $\Lambda$ cold dark matter (CDM) paradigm—in which dark energy is a cosmological constant and dark matter (DM) is stable, cold, and collisionless—has been established as the concordance model of cosmology, accurately predicting the detailed content and structure of the universe throughout cosmic history (DES Collaboration 2018; Planck Collaboration 2020). However, several recent cosmological tensions— namely the “Hubble tension,” concerning the present-day expansion rate of the universe, $H_0$, and the “$\sigma_8$ tension,” concerning the amplitude of matter clustering on quasi-linear scales—potentially point to new physics beyond $\Lambda$CDM (see Verde et al. 2019 for a review). To address these tensions, it is essential to stress-test every assumption underlying the $\Lambda$CDM model, including the nature of the DM particle.

Any viable DM model must accurately predict both the expansion history of the universe and the formation of structure throughout cosmic history. To date, all cosmological observations are consistent with a DM particle that is stable against decays. However, there is a class of models that allow DM to decay, either to other dark sector particles or to Standard Model particles, which is consistent with current observational constraints. These decays often transfer energy from DM to radiation-like species in the late-time universe, increasing the present-day expansion rate relative to its extrapolated evolution based on CMB measurements (e.g., Vattis et al. 2019; Clark et al. 2021b). Intriguingly, recent studies of models where the DM particle decays with a lifetime comparable to the Hubble time suggest that DM decays can potentially alleviate the $H_0$ and/or $\sigma_8$ tensions (e.g., Enqvist et al. 2015; Vattis et al. 2019; Abellán et al. 2022; Blinov et al. 2020; Pandey et al. 2020;

Haridasu & Viel 2020; Abellán et al. 2021; Clark et al. 2021a, 2021b).

In this paper, we focus specifically on DM models where late-time decays into other dark sector particles alter the velocity distribution of DM, which we refer to generically as decaying dark matter (DDM). While these DDM models have predominantly been studied using expansion history and large-scale structure probes (e.g., see Wang & Zentner 2012; Chen et al. 2021), cosmic structure on smaller scales offers a unique window into the microphysics of DM decays. In particular, because viable DDM lifetimes are comparable to the age of the universe, the effects of DM decays accumulate over time as matter perturbations grow in the post-recombination epoch. Late-time decays thus suppress the clustering of matter on small scales at late times; this effect has been leveraged to constrain DDM models using $\text{Ly}\alpha$ forest data (e.g., McDonald et al. 2006; Wang et al. 2013).

DM decays impart momentum to the daughter particles, injecting kinetic energy into virialized DM halos. These momentum “kicks”—which depend on the details of the decay mechanism—gradually reduce the central density of cuspy DM halos and deplete halos of mass (e.g., Peter et al. 2010b). This process also makes subhalos more susceptible to tidal disruption as they orbit within a host halo, leading to the suppression of the abundance of low-mass halos and subhalos in DDM models at late times. This manifests as a deficit of small-scale structure, including faint satellite galaxies, relative to $\Lambda$CDM predictions (Peter et al. 2010b; Peter & Benson 2010; Wang et al. 2014) and also alleviates potential tensions between the predicted and inferred central densities of bright dwarf galaxies (the “too big to fail” problem; Boylan-Kolchin et al. 2011, 2012; Garrison-Kimmel et al. 2014; Papastergis et al. 2015).

The population of Milky Way (MW) satellite galaxies, which contains the faintest observed galaxies in the universe, therefore offers a unique testing ground for DDM. Faint dwarf galaxies orbiting the MW occupy the smallest DM halos directly associated with galaxies and are the most dark matter-

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

<footnote n="53">⁵³ NHFP Einstein Fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
dominated systems known (see a recent review by Simon 2019). DM microphysics that impacts the formation, late-time abundance, and internal structure of small DM halos and subhalos can therefore be inferred from the abundance and properties of MW satellites.

Here, we derive robust and stringent constraints on the DM particle lifetime using a state-of-the-art census of the MW satellite population. Although DDM has previously been studied using MW satellites (e.g., Peter & Benson 2010; Wang et al. 2014), our work is the first to constrain it by leveraging observations over nearly the full sky, including the population of ultra-faint dwarf galaxies recently discovered by the Dark Energy Survey (DES; DES Collaboration 2005, 2016). Following the procedure developed in Drlica-Wagner et al. (2020) and Nadler et al. (2020b, 2021c), we combine (i) the observed population of MW satellites in DES and Pan-STARRS1 (PS1; Chambers et al. 2016), (ii) observational selection functions derived from MW satellite searches in DES and PS1 data, (iii) a detailed forward model of the connection between MW satellite galaxies and DM subhalos, (iv) high-resolution cosmological simulations of MW-mass halos in DDM cosmologies, and (v) a Bayesian inference framework to constrain the DM particle lifetime. This method was used to constrain microphysical DM properties including its primordial velocity distribution, Standard Model coupling, and particle mass in Nadler et al. (2021c).

Unlike previous DDM studies, our model only relies on the inferred abundance of low-mass subhalos above a minimum mass threshold, rather than their inferred central densities, which makes our constraints complementary to previous studies. This minimum mass threshold corresponds to the peak mass of the smallest halo inferred to host MW satellites observed in the DES and PS1 data for CDM (Nadler et al. 2020b), and therefore represents a conservative upper limit on this quantity for DDM, which further reduces halo masses as explored below. In turn, our results provide a foundation for future work that combines satellite abundances with stellar velocity dispersion measurements to search for signatures of DM decays. Similar approaches will be useful to test other DM properties that simultaneously modify satellite abundances and density profiles at an observable level, including DM self-interactions (e.g., Vogelsberger et al. 2016; Tulin & Yu 2018; Drlica-Wagner et al. 2019; Salucci 2019; Kim & Peter 2021).

This paper is outlined as follows. In Section 2, we describe our fiducial DDM model. We describe the impact of DDM physics and analytically estimate subhalo disruption timescales in Section 3. In Section 4, we derive a fitting function for the suppression of the subhalo mass function using cosmological zoom-in simulations of MW-like halos in DDM cosmologies. We incorporate this suppression in a forward model of the MW satellite population to derive DDM constraints in Section 5. We compare our results to other small-scale and cosmological probes of DDM in Section 6.1, and we conclude in Section 7. With the exception of Section 3, we work in units with $c = 1$. Furthermore, “log” refers to the base-10 logarithm.

## 2. Decaying Dark Matter Overview

In this section, we outline classes of DM models that undergo decays, connect common model assumptions to cosmological observables, and describe the particular model of DDM that we consider in this work. We emphasize that the DDM models we describe are phenomenological, rather than first-principles

particle physics descriptions; models that feature scattering between DM and other particles are broadly categorized under self-interacting and/or interacting DM.

### 2.1. Decaying Dark Matter Phenomenology

DM models featuring decays typically involve massive parent DM particles, a fraction of which undergo particle decay to some number of daughter particles with a given decay lifetime. Accordingly, decaying DM models can be broadly differentiated based on the following characteristics:

1. The lifetime of the decay;

2. The number of particles involved in each DM decay, and their masses; and

3. Whether the decays exclusively involve dark sector particles or involve any Standard Model particles.

DDM is typically regarded as a modification to the CDM paradigm rather than a completely distinct model; for instance, superWIMPs (Feng et al. 2003a, 2003b; Ichiki et al. 2004) are an example of decaying CDM (Cen 2001). However, each of the assumptions listed are associated with specific cosmological signatures that distinguish them from stable CDM. We now consider the cosmological effects of each assumption in turn.

The lifetime of the DM decay sets the onset and rapidity of decays. For extremely short lifetimes compared to the age of the universe, with $\tau \lesssim \mathcal{O}(\text{yr})$ (Kaplinghat 2005), decays before recombination can transition a fraction of the dark matter into a “warm” state by introducing an additional nonthermal velocity component; this can change the initial conditions for structure formation and increase the effective number of relativistic degrees of freedom in the early universe (e.g., Blinov et al. 2020 and references therein). For decays that occur after recombination but on a shorter timescale than the age of the universe, cosmological observables are generally affected by the resulting changes in DM composition at late times. For extremely long decay lifetimes—of the order the age of the universe or longer—the observable impacts relative to a stable DM model are minimal. Furthermore, note that the lifetime of the decay is occasionally parameterized according to the fraction of DM that has undergone decays by $z=0$ (e.g., Cen 2001; Sánchez-Salcedo 2003).

The number of particles involved in the decay has several consequences for cosmological observables. By allowing for the parent particle to decay into $N$ products, it is possible to introduce up to $N$ new and potentially unique dark sector particles. The case of $N=2$ allows for the decay products to take on a range of mass splitting values; our fiducial model (detailed below) is an example of such a decay in which one of the decay products is a dark radiation component, and the other is a massive particle of comparable mass to the parent particle. For $N > 2$, the complexity of the model increases (e.g., as in Blackadder & Koushiappas 2014, 2016; Haridasu & Viel 2020), leading to cosmological signatures that are generally more difficult to disentangle than the $N=2$ case.

Finally, it is interesting to consider whether or not the decay products include Standard Model particles. Particle collider experiments have been used to constrain dark matter decays that operate via Standard Model portals and scattering mechanisms, including the Higgs portal (e.g., ATLAS Collaboration 2021) or a more general coupling in the case of pseudo-Dirac dark matter (e.g., Jordan et al. 2018; Bhattacharya & Slatyer 2019; González & Toro 2021). Similarly, gamma-ray and X-ray telescopes
constrain WIMP-like models that decay into Standard Model particles (e.g., Ackermann et al. 2015; Gaskins 2016; Ando et al. 2021).

## 2.2. Fiducial Decaying Dark Matter Model

In this paper, our fiducial DDM model consists of a cold, massive DM particle $\chi$ that decays to a slightly less massive daughter DM particle $\chi'$ and a massless dark radiation species $\gamma'$,

$$ \chi \rightarrow \chi' + \gamma'. \quad (1) $$

In our fiducial cosmology, $\chi$ composes the entirety of the initial DM. Denoting the mass of $\chi$ by $M$ and the mass of $\chi'$ by $m$, this model can be parameterized by the decay lifetime $\tau$ and the mass splitting (Wang et al. 2014)

$$ \epsilon \equiv \frac{M - m}{M}. \quad (2) $$

These decays impart a recoil velocity

$$ V_{\text{kick}} \approx \epsilon, \quad (3) $$

assuming $V_{\text{kick}} \ll 1$, on the daughter DM particle in the center-of-mass rest frame of the parent particle (recall that we work in natural units with $c = 1$). Following Peter (2010) and Wang et al. (2014), we consider a parameter space in which the decays occur at late times with small recoil velocities, with $\tau \sim O(\text{Gyr})$ and $\epsilon \ll 1$. Such models only mildly affect large-scale structure (e.g., Poulin et al. 2016) while yielding potentially observable effects on smaller scales, and particularly on low-mass DM halos. Note that our fiducial DDM model makes no assumptions about interactions between the dark sector (i.e., any of the particles $\chi$, $\chi'$, or $\gamma'$) and the Standard Model.

In Section 4.3, we demonstrate that the late-time suppression of DDM subhalo abundances can be fit reasonably well using a functional form similar to that used to describe the suppression of the subhalo mass function in warm dark matter (WDM). However, we emphasize that late-time subhalo disruption is a dynamical effect in our DDM models while, for WDM-like models, the effect is imprinted on the linear matter power spectrum at early times. This distinction between DDM and WDM models can also lead to differences in the relative abundance of low-mass isolated halos and subhalos at a given mass scale; in particular, the abundance of isolated halos and subhalos is roughly equally suppressed in WDM-like models (Stücker et al. 2022), while subhalo abundances can be preferentially suppressed in DDM.

It is informative to consider the limiting cases of our fiducial two-body decay model. Longer DDM lifetimes lead to fewer DM particle decays, and the model becomes more similar to CDM at arbitrarily later times. Similarly, a smaller mass splitting between parent and daughter particles leads to smaller kick velocities—and less energy carried away by the dark radiation component—and the model again approaches CDM. We will consider lifetimes comparable to the Hubble time (or longer) and kick velocities comparable to the circular velocities of subhalos that host MW satellite galaxies; thus, our analysis only directly constrains decays that change the distribution of DM structure relative to CDM at late times.

# 3. Effects of Dark Matter Decays on Milky Way Satellites

To develop intuition for the effects of dark matter decays on MW satellite galaxy abundances, we consider a simple toy model wherein we compute the maximum initial halo mass that remains above a $z=0$ mass threshold corresponding to the subhalos that host the faintest galaxies we use in our likelihood analysis. This calculation provides a rough estimate of the halo mass scales that are significantly affected in the DDM models we consider; however, we use the cosmological zoom-in simulations presented in Section 4 to make precise predictions for DDM subhalo abundances.

We note that, although the toy model presented in this section applies equally to subhalos and isolated halos, subhalos experience tidal effects (after infall) that accelerate their mass loss and that are exacerbated for DDM subhalos with reduced central densities (Peter 2010; also see Figure 3). Thus, our estimates of the halo masses that are significantly affected by DM decays represent conservative lower limits for the subhalo masses that are impacted, and the combination of decays and tidal stripping shapes the simulated DDM subhalo populations we study in Section 4. We reinstate units of $c$ in this section for clarity.

For a halo with an initial mass $M_i$ at time $t_i$, it is possible to determine the time needed for decays to reduce the halo mass below a given threshold. By choosing this threshold to correspond to the minimum mass of subhalos that host observed MW satellites, this calculation yields a rough estimate of the impact of DM decays on MW satellite abundances. In addition, it highlights an important difference between DDM and other non-CDM models, including WDM, which are commonly considered on dwarf galaxy scales. In DDM, low-mass halos can form at early times but evaporate due to decays by late times, whereas, in models that only suppress the linear matter power spectrum (e.g., WDM), subhalos below a mass threshold never form.

For a halo with an initial number of particles $N_i$ that decay with a lifetime $\tau$, the number of particles that have not undergone decay within a time interval $\Delta t \equiv t_f - t_i$ is given by

$$ N(t) = \exp\left[-\frac{\Delta t}{\tau}\right]N_i, \quad (4) $$

assuming that all daughter particles are ejected from the halo due to the recoil kick velocity. Because we will analyze DDM models with $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$, this is a reasonable assumption for the subhalos that drive our constraints, which typically have peak virial masses of $\sim 10^8 M_{\odot}$ and peak maximum circular velocities of $V_{\text{peak}} \approx 20 \text{ kms}^{-1}$ (Nadler et al. 2019b, 2020b), corresponding to escape velocities of $V_{\text{esc}} \approx 28 \text{ kms}^{-1}$. In the regime of $V_{\text{kick}} \gtrsim V_{\text{peak}}$ (which guarantees that $V_{\text{kick}} \gtrsim V_{\text{max}}$), the impact of decays on low-mass halos is effectively determined by the DM lifetime alone (Wang et al. 2014).$^{54}$

We also assume that all parent and all daughter DM particles, respectively, have the same mass; hence, the total

<footnote>$^{54}$ This is consistent with a calculation of the energy injection due to DM decays following Abdelgader & Melia (2008), which indicates that decays are unlikely to unbind halos of the masses we consider within a Hubble time.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Final subhalo mass (in units of peak virial mass) versus dark matter decay lifetime for a subhalo that undergoes decays for the age of the universe ($t_H = 13.7$ Gyr). For each peak subhalo mass, dashed horizontal lines indicate a mass threshold $M_{\text{thresh}} = 2 \times 10^7 M_{\odot}$ that we assume corresponds to the minimum mass for subhalos that host satellites galaxies. This mass threshold corresponds to subhalos resolved with roughly 100 particles in our zoom-in simulations (see Section 4.1). For decay lifetimes of $\lesssim 10$ Gyr, the abundance of low-mass ($M_{\text{peak}} \sim 10^8 M_{\odot}$) subhalos above the mass threshold is significantly affected by decays alone (gray region).</figcaption></figure>

mass of the halo after time interval $\Delta t$ is then given by

$$M_{\text{sub}}(\Delta t) = \exp\left[-\frac{\Delta t}{\tau}\right]M_i. \qquad (5)$$

To facilitate comparison with our simulation results, we adopt the peak virial mass $M_{\text{peak}}$, defined as the largest virial mass a subhalo attains along its main branch in our zoom-in simulations (see Section 4.1), as the initial mass $M_i$. For simplicity and because $M_{\text{peak}}$ is achieved at fairly high redshift for typical subhalos in our simulations, we adopt the age of the universe $t_H$ as $\Delta t$. Neither assumption has an important effect on the qualitative results of our mass-loss calculation, which is shown in the left panel of Figure 1. The majority of observed MW satellites that enter our analysis inhabit subhalos in the peak virial mass range of $10^8 M_{\odot} \lesssim M_{\text{peak}} \lesssim 10^{10} M_{\odot}$, where $M_{\text{peak}}$ is the largest virial mass each subhalo ever achieves over its entire main-branch history (Nadler et al. 2020b). In this calculation and in our comparison to the MW satellite population, we additionally assume that halos that fall below a minimum mass threshold at $z=0$ cannot host observed MW satellites. In our cosmological DDM simulations, subhalos are no longer identified when their mass falls below a resolution limit of $\approx 2 \times 10^7 M_{\odot}$, corresponding to roughly 100 particles (see Section 4.1); we adopt this as our minimum subhalo mass threshold for MW satellites, $M_{\text{thresh}}$. This mass threshold is reasonable given the dynamical properties of observed MW satellites (e.g., Strigari et al. 2008).

Imposing our minimum $z=0$ mass threshold yields a condition for the maximum $M_{\text{peak}}$ that is significantly impacted by decays:

$$M_{\text{sub}}(t_H) \le M_{\text{thresh}} \Rightarrow \exp\left[-\frac{t}{\tau}\right] M_{\text{peak}} \le M_{\text{thresh}}. \qquad (6)$$

Figure 1 shows the relation between final halo mass (in units of $M_{\text{peak}}$) and decay lifetime, while indicating our mass threshold over the range of peak subhalo masses relevant for our analysis. Thus, for decay lifetimes $\tau \lesssim 10$ Gyr, we expect low-mass subhalos ($M_{\text{peak}} \sim 10^8 M_{\odot}$) to approach the threshold mass and (by construction) the resolution limit of our simulations due to mass loss from decays alone. This foreshadows the suppression of subhalo abundances at these mass scales derived from our simulations in Section 4.2. Again, we note that our mass-loss calculation applies equally to isolated halos, and thus predicts that the abundance of isolated halos with similar peak masses that remain above our threshold mass at $z=0$ is suppressed by decays. Appendix C examines this effect in our cosmological zoom-in simulations.

We reiterate that the results in Figure 1 represent order-of-magnitude estimates of DDM effects on low-mass subhalos, and that our constraints are instead based on detailed cosmological DDM simulations (see Section 4). Nonetheless, the toy model presented above demonstrates that the abundance and structure of subhalos expected to host MW satellite galaxies are sensitive to DM particle lifetimes of $\mathcal{O}(10 \text{ Gyr})$ and to the microphysics of these decays encapsulated by $V_{\text{kick}}$.

## 4. Decaying Dark Matter Subhalo Populations

The DDM models we consider differ from many popular alternatives to CDM in that they feature late-time suppression of small subhalos, rather than a suppression in the linear matter power spectrum, which eventually manifests as a reduction of subhalo abundances (i.e., smaller subhalos simply do not form in, for example, WDM, whereas in DDM they can form initially but are disrupted by late times). Importantly, this precludes a first-principles mapping between the impact of DDM and other well-studied CDM alternatives like WDM on small-scale structure formation. Instead, we study the effects of DM decays at late times and their impact on subhalo populations using cosmological zoom-in simulations of MW-mass systems to model the nonlinear impact of late-time DM decays. We use these simulations to fit the suppression of the subhalo mass function as a function of $\tau$ for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$ separately, and we use these fits in our likelihood framework to constrain $\tau$ as described in Section 5.2.<fnref n="55" />

Section 4.1 gives an overview of the DDM and CDM simulations used in this work, Section 4.2 describes the general characteristics of the DDM subhalo populations, and Section 4.3 studies the suppression of subhalo abundances relative to CDM and derives an analytic model for this effect to enable the statistical inference performed in Section 5.

### 4.1. Cosmological Zoom-in Simulations

To study the impact of DDM physics on low-mass subhalos and to derive predictions for the suppression of the subhalo mass function (SHMF), we use an expanded set of cosmological zoom-in simulations of MW-mass halos in DDM models based on those presented in Wang et al. (2014). In particular, we study six simulations with DDM parameters ($\tau/\text{Gyr}$, $V_{\text{kick}}/\text{kms}^{-1}) \in \{(10, 20), (20, 20), (20, 40), (40, 20), (40, 40), (80, 40)\}$ and a corresponding CDM simulation, all with

<footnote n="55">55 We focus on $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$ following Peter (2010) and Wang et al. (2014), who found that these kick velocities are sufficient to significantly impact subhalo abundances and internal densities.</footnote>
identical initial conditions. These simulations were performed using a version of the GADGET-2 (Springel 2005) and GADGET-3 N-body codes as modified by Peter et al. (2010b). The simulations are run in a box of length $50 h Mpc^{-1}$ per side and zoom in on a halo mass of $M \approx 10^{12} M_{\odot}$ (Wang et al. 2014). The highest-resolution region is simulated with a Plummer-equivalent force softening of 143 pc and a particle mass of $1.92 \times 10^5 M_{\odot}$.<fnref n="56" /> Cosmological parameters are set to $\Omega_M = 0.266$, $\Omega_{\Lambda} = 0.734$, $n_s = 0.963$, $h = 0.71$, and $\sigma_8 = 0.801$ (Jarosik et al. 2011), and we analyze these simulations using the ROCKSTAR (Behroozi et al. 2013a) halo finder and CONSISTENT-TREES (Behroozi et al. 2013b) merger tree code.

We emphasize that the Wang et al. (2014) simulations are not specifically selected to resemble the MW in properties beyond its host halo mass. To derive the DDM constraints in Section 5, we therefore analyze the two additional MW-like N-body zoom-in simulations originally presented in Mao et al. (2015) and studied in Nadler et al. (2020b, 2021c) to model the MW satellite galaxy population. The properties of these host halos are consistent with the mass and concentration of the MW halo and include both realistic Large Magellanic Cloud (LMC) analog systems on recent infall and Gaia-Sausage-Enceladus-like merger events at early times (Belokurov et al. 2018; Helmi et al. 2018). We describe these simulations in detail in Appendix A. We also perform several DDM resimulations of these systems, which we use in Appendix B to validate that the impact of DDM on the subhalo populations in our expanded suite of Wang et al. (2014) simulations is consistent with its impact on MW-like systems. Throughout this paper, the expanded suite of simulations from Wang et al. (2014) are represented with a yellow-green-blue-purple color scheme. Plots using the MW-like resimulations only appear in Appendix B & C and are represented with a red-purple color scheme.

## 4.2. Decaying Dark Matter Subhalo Populations

We begin by summarizing the main differences between the DDM and CDM subhalo populations in the Wang et al. (2014) simulations. Figure 2 shows the SHMFs in the Wang et al. (2014) simulations introduced in Section 4.1, defined as the cumulative abundance of subhalos as a function of peak subhalo virial mass. The SHMFs for DDM models with a range of lifetimes for each $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (40 kms$^{-1}$) are shown on the upper left (right) panels. Lowering the decay lifetime and increasing the kick velocity both systematically decrease the abundance of surviving subhalos at low masses. Note that the ratio of the DDM SHMF relative to CDM at peak subhalo masses above $\sim 10^9 M_{\odot}$ is consistent with unity within the statistical precision of our simulations; this behavior is conservatively reflected in our fitting function predictions (lower panels of Figure 2; see Section 4.3). We leave a detailed study of dark matter decays on the abundance and density profiles of more massive subhalos to future work. Importantly, our main results are driven by the abundance of subhalos with masses below this threshold (see Nadler et al. 2020b, 2021c).

Figure 3 shows the distribution of subhalo present-day maximum circular velocity, $V_{\text{max}}$, divided by its peak value, $V_{\text{peak}}$, again for DDM models with $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (40 kms$^{-1}$)

on the left (right). There is a pronounced decrement in this ratio for each of our DDM simulations relative to their CDM counterparts, implying that subhalos' central densities are reduced in DDM, consistent with the findings of Peter et al. (2010b) and Wang et al. (2014). This effect could be leveraged to improve DDM constraints by incorporating MW satellites' stellar velocity dispersion measurements in the inference (e.g., Drlica-Wagner et al. 2019; Kim & Peter 2021), rather than relying solely on their surface brightness function and radial distribution.

Figure 4 shows the $z=0$ radial distribution of surviving subhalos in the DDM and CDM simulations. Note that the suppression of subhalo abundances in these DDM simulations is roughly radially independent, which is consistent with the intuition developed in Section 3: because decays alone significantly deplete the mass contained in small subhalos, the suppression of subhalo abundance is mainly determined by $M_{\text{peak}}$, which is not strongly correlated with present-day radial distance (e.g., Springel et al. 2008), even in the presence of the LMC (Nadler et al. 2021a).<fnref n="57" /> On the other hand, note that suppression of subhalo abundances due to the Galactic disk is a strong function of radius and only weakly depends on subhalo mass (Garrison-Kimmel et al. 2017; Kelley et al. 2019). We will exploit the fact that the shape of the radial distribution of surviving subhalos is approximately unchanged in DDM relative to CDM when deriving our constraints in Section 5.

In Appendix C, we study the evolution of DDM subhalo abundances in our MW-like resimulations described in Appendix A. We find that the suppression of the DDM SHMF relative to CDM sets in at late times ($z \lesssim 2$), which is expected based on the long timescale (of the order of the Hubble time) of the decays in the models we consider. This is consistent with the intuition that small-scale structure is only suppressed in DDM at late times. This implies that probes of small-scale structure in the low-redshift universe, including MW satellite galaxies, stellar streams, and strong gravitational lenses, are uniquely suited to constrain DDM models. On the other hand, probes of the same comoving scales at earlier times, including the Ly$\alpha$ forest (Viel et al. 2013; Iršič et al. 2017; Palanque-Delabrouille et al. 2020) and the high-redshift galaxy luminosity function (Schultz et al. 2014; Menci et al. 2016; Corasaniti et al. 2017; Rudakovskiy et al. 2021), offer relatively less constraining power for these models.

## 4.3. Subhalo Mass Function Suppression

Having explored the main features of the subhalo populations in our DDM simulations, we now derive a fitting function for the suppression of the SHMF relative to CDM. In particular, we express the DDM SHMF as

$$ \left( \frac{dN_{\text{sub}}}{dM} \right)_{\text{DDM}} \equiv f_{\text{DDM}}(M, \tau, V_{\text{kick}}) \left( \frac{dN_{\text{sub}}}{dM} \right)_{\text{CDM}} , \quad (7) $$

where $f_{\text{DDM}}(M, \tau, V_{\text{kick}})$ is the suppression of the DDM SHMF relative to that in CDM as a function of subhalo peak virial mass, decay lifetime, and recoil kick velocity. Note that, due to the resolution limit of our simulations, $f_{\text{DDM}}$ describes the suppression of subhalo abundances above a present-day mass threshold of $\sim 2 \times 10^7 M_{\odot}$. To measure this quantity, we fit the

<footnote n="56">⁵⁶ We define virial quantities according to the Bryan & Norman (1998) virial definition, with overdensities $\Delta_{\text{vir}}$ set according to the cosmological parameters in our two simulation suites.</footnote>

<footnote n="57">⁵⁷ This is also supported by Figures 11–12, which show that the suppression of isolated halo and subhalo abundances in our DDM resimulations of MW-like systems is comparable.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top panels: subhalo mass functions at $z=0$ for the MW-mass host halo in the expanded suite of zoom-in simulations from Wang et al. (2014). SHMFs are shown at $z=0$ in CDM (dashed gray) and in DDM models with $\tau = 10, 20, 40$, and $80$ Gyr (from yellow to purple), with $V_{\text{kick}} = 20 \text{ km s}^{-1}$ (left panel) and $40 \text{ km s}^{-1}$ (right panel) as a function of peak subhalo virial mass $M_{\text{peak}}$. All SHMFs are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ km s}^{-1}$ and $V_{\text{max}} > 9 \text{ km s}^{-1}$. Shorter DDM decay lifetimes result in fewer surviving subhalos at low peak masses relative to CDM, and this effect is more pronounced for models with higher kick velocities. Bottom panels: suppression of the subhalo mass function relative to CDM at $z=0$ for the same DDM models shown in the top panels. Solid lines show the SHMF suppression measured directly from our expanded suite of zoom-in simulations based on Wang et al. (2014), and dashed lines show the best-fit function derived in Section 4.3 for each DDM model. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements, which are consistent with no suppression (i.e., $N_{\text{DDM}}/N_{\text{CDM}} = 1$ at peak masses $M_{\text{peak}} \gtrsim 10^9 M_{\odot}$). Thus, our fitting functions approach unity in the high-mass regime.)</figcaption></figure>

output of the DDM simulations from Wang et al. (2014) described above using a functional form of the SHMF similar to that proposed for warm dark matter (Benito et al. 2020; Lovell 2020):

$$f_{\text{DDM}}(M, \tau, V_{\text{kick}}) = \left(1 + \left(\frac{M_0}{M}\right)^{\beta}\right)^{-\gamma(\tau, V_{\text{kick}})}, \quad (8)$$

where

$$\gamma(\tau, V_{\text{kick}}) \equiv \gamma_0 \frac{V_{\text{kick}}}{\tau}. \qquad (9)$$

To determine the parameters $M_0$, $\beta$, and $\gamma_0$, we perform least-squares fitting over a grid of SHMF values as a function of $M_{\text{peak}}$ with varying $\tau$ and $V_{\text{kick}}$. We perform this fit over a mass range of $3.8 \times 10^7 - 9.4 \times 10^9 M_{\odot}$ and derive best-fit values of

$$\log(M_0/M_{\odot}) = 8.5, \qquad (10)$$

$$\beta = 0.6, \qquad (11)$$

$$\gamma_0 = 0.8, \qquad (12)$$

with a reduced $\chi^2$ goodness of fit of 0.3 (which reflects the relatively large Poisson uncertainties on the SHMFs).<fnref n="58" /> The fit of Equation (8) with these best-fit values is shown in the lower panels of Figure 2. Although this fit slightly overpredicts SHMF suppression at peak masses below $\sim 5 \times 10^8 M_{\odot}$, we estimate that this only influences our limits at the $\sim 5$ Gyr level because most of the constraining power results from satellites associated with higher-mass halos (Section 5.2).

It is important to note that, while Equation (8) empirically describes the suppression of the DDM SHMF relative to CDM, it is neither derived from first principles nor physically motivated. Moreover, it only accurately describes SHMF suppression for the range of $(\tau, V_{\text{kick}})$ sampled in our expanded suite of simulations based on Wang et al. (2014), which is

<footnote n="58">⁵⁸ The suppression of the DDM SHMF has a slightly weaker dependence on subhalo mass than in WDM, for which $\beta \approx 1.0$ and $\gamma \approx -0.99$ (Lovell et al. 2014).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Cumulative distribution function of subhalo $z = 0$ maximum circular velocity, $V_{max}$, relative to its peak value, $V_{peak}$, for subhalos of the MW-mass host halo in the expanded suite of zoom-in simulations from Wang et al. (2014). Distributions are shown for CDM (gray) and in DDM models with $\tau = 10, 20, 40, \text{and } 80$ Gyr (from yellow to purple), with $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (left column) and $40 \text{ kms}^{-1}$ (right column). All distributions are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ kms}^{-1}$ and $V_{\text{max}} > 9 \text{ kms}^{-1}$. Subhalos in our DDM simulations exhibit systematically lower values of $V_{\text{max}}/V_{\text{peak}}$ relative to CDM, and this effect is more pronounced for shorter decay lifetimes and higher kick velocities. This indicates that the central densities of DDM subhalos are reduced, consistent with previous studies (e.g., Peter et al. 2010b; Wang et al. 2014).</figcaption></figure>

bounded below (above) by $\tau/\text{Gyr}$, $V_{\text{kick}}/\text{kms}^{-1} = 10, 20 (80, 40)$, with $z=0$.

# 5. Constraints from Milky Way Satellite Galaxies

## 5.1. Forward Model and Fitting Procedure

To constrain the effects of DDM physics using observations of MW satellite galaxies, we incorporate the suppression of the SHMF derived from our DDM simulations into a forward-modeling framework that generates realizations of the satellite populations observed by DES and PS1. In particular, we apply the galaxy-halo connection model presented in Nadler et al. (2020b) to the subhalo populations in the two CDM MW-like simulations from the Mao et al. (2015) suite in order to predict the absolute magnitude, half-light radius, and Galactocentric distance distributions of satellites corresponding to subhalos in each simulation. We model the suppression of the SHMF in DDM by applying a fitting function to these CDM simulations because they have realistic LMC analogs and to enable a more direct comparison with Nadler et al. (2021c).

Our galaxy-halo connection model includes eight free parameters that control the abundance-matching model that relates satellite luminosity to subhalo peak maximum circular velocity, the size model that relates satellite half-light radius to subhalo size at accretion, the efficiency of subhalo disruption due to the Galactic disk, and the minimum peak halo mass and scatter of the galaxy occupation fraction. These parameters are defined in Appendix D and Table 2; we refer the reader to Nadler et al. (2020b) for a comprehensive description of the galaxy-halo connection model. Following Nadler et al. (2021c), we add one free parameter to this model, $\tau$, which controls the suppression of the DDM SHMF at a given $V_{\text{kick}}$ according to Equation (8), and we perform the analysis for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$ separately. The contribution of each satellite to the mock observed number count is then weighted according to the probability that its corresponding subhalo survives for a given set of DDM parameters.

As in Nadler et al. (2021c), this procedure assumes that the shape of the subhalo radial distribution is unchanged in DDM

relative to CDM, which is demonstrated in Figure 4 for our MW-mass simulations and Figure 10 for the MW-like simulations used in the inference. Furthermore, we do not modify the fiducial subhalo disruption probabilities predicted by the Nadler et al. (2018) algorithm, which was calibrated on CDM hydrodynamic simulations. This is a conservative assumption because DDM reduces the central densities of surviving subhalos (see Figure 3), making them more susceptible to tidal disruption; however, because these disruption probabilities are marginalized over in our fitting procedure, this assumption is not expected to significantly impact our constraints. Note that we do not account for adiabatic expansion of satellite sizes due to decays, which is also a conservative strategy because this effect would push some predicted satellites below the detectability threshold, thereby forcing even lower-mass subhalos to contribute and leading to more stringent DDM constraints.

For a given set of galaxy-halo connection and DDM model parameters, we perform mock observations of the DES and PS1 satellite populations using the observational selection functions presented in Drlica-Wagner et al. (2020) by self-consistently orienting the survey footprints relative to the LMC analogs in our MW-like simulations. Thus, our procedure explicitly incorporates inhomogeneities in the spatial distribution and detectability of MW satellites and yields realizations of the observed DES and PS1 satellite populations that are compared to the data from Drlica-Wagner et al. (2020) by assuming that satellite surface brightness is distributed according to a Poisson point process in each survey footprint.

Following Nadler et al. (2020b, 2021c), we use the sample of kinematically confirmed and probable dwarfs from Drlica-Wagner et al. (2020). We remind the reader that, unlike in the WDM-like analyses from Nadler et al. (2021c), we additionally assume that observed MW satellites occupy subhalos above a minimum present-day mass threshold of $\sim 2 \times 10^7 \text{M}_\odot$, which was imposed when deriving our DDM SHMF suppression predictions. Although subhalos stripped below this resolution threshold are technically included in our analysis through our orphan satellite model following Nadler et al. (2020b), the abundance of these systems is assumed to follow an extrapolation of the DDM SHMF
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Subhalo radial distributions at $z=0$ as a function of distance from the center of the host halo in the expanded suite of zoom-in simulations from Wang et al. (2014) for CDM (dashed gray) and for DDM models with $\tau = 10, 20, 40$, and 80 Gyr (from yellow to purple), with $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (left panel) and $40 \text{ kms}^{-1}$ (right panel). The bottom panels show the ratio of the radial distribution in each DDM model relative to CDM. All results are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ kms}^{-1}$ and $V_{\text{max}} > 9 \text{ kms}^{-1}$. Shorter DDM decay lifetimes result in fewer surviving subhalos at all radii relative to CDM, and this effect is more pronounced for models with higher kick velocities. The virial radius is taken to be 300 kpc. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements.</figcaption></figure>

suppression derived in Section 4.2, and their contribution to our results is therefore negligible.

We use the Markov Chain Monte Carlo code emcee (Foreman-Mackey et al. 2013) to simultaneously fit for eight parameters governing the galaxy–halo connection and the efficiency of subhalo disruption due to the Galactic disk, and one parameter governing the impact of DDM. In particular, we fit for the power-law slope of the satellite luminosity function ($\alpha$), the scatter in luminosity at fixed $V_{\text{peak}}$ ($\sigma_M$), the mass at which 50% of halos host galaxies ($\log(M_{50}/M_{\odot})$), the strength of subhalo disruption due to baryons ($B$), the scatter in the galaxy occupation fraction ($\sigma_{\text{gal}}$), the amplitude of the galaxy–halo size relation ($A$), the scatter in half-light radius at fixed halo size ($\sigma_{\log R}$), the power-law index of the galaxy–halo size relation ($n$), and the lifetime of the DM particle ($\log(\tau/\text{Gyr})$). Note that we fit for $\log \tau$ (where $\tau$ is measured in gigayears), rather than $\tau$ itself, using a uniform prior on the logarithmic quantity in the range $\log \tau \in [1.3, 1.9]$; the remaining prior distributions are identical to those adopted in Nadler et al. (2020b).<fnref n="59" /> We perform two separate fits at fixed $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$, each of which uses 36 walkers, discards a burn-in period of ~20 autocorrelation lengths, and retains ~10⁵ samples corresponding to ~100 autocorrelation lengths.

## 5.2. Results

The posterior distributions from our $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $V_{\text{kick}} = 40 \text{ kms}^{-1}$ fits are summarized in Table 1 and shown in Figure 5. The marginalized posterior distributions for the eight

galaxy–halo connection parameters are nearly identical in both cases, and are consistent with the CDM fit in Nadler et al. (2020b). For both values of $V_{\text{kick}}$, we obtain a lower limit on $\log \tau$, which is expected because the CDM model (i.e., the large $\tau$ limit) is consistent with the data. At 95% confidence, we obtain $\log(\tau/\text{Gyr}) > 1.46$ ($\tau > 29 \text{ Gyr}$) for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $\log(\tau/\text{Gyr}) > 1.63$ ($\tau > 43 \text{ Gyr}$) for $V_{\text{kick}} = 40 \text{ kms}^{-1}$ (Table 1).

Following Nadler et al. (2021c), we scale these constraints to conservatively account for uncertainty in the MW host halo mass, which is known to impact limits on non-CDM models derived from the MW satellite population (e.g., Newton et al. 2021). In particular, for each $V_{\text{kick}}$, we compute the value of $\tau$ that decreases $f_{\text{DDM}}$ by 27% relative to its value at the original $\tau$ constraint when evaluated at a peak subhalo mass of $3.2 \times 10^8 M_{\odot}$. This value corresponds to the minimum halo mass probed by the DES and PS1 data in CDM (Nadler et al. 2020b), and therefore represents an upper limit on the minimum halo mass in our DDM inference. The 27% uncertainty corresponds to the ratio of the maximum allowed MW halo mass from Callingham et al. (2019), appropriately converted to our virial mass definition, relative to the average host halo mass from our MW-like simulations used to perform the inference. We assume that this uncertainty affects the allowed amount of SHMF suppression linearly when deriving our conservative constraints on $\tau$. For lower values of $\tau$, the number of MW satellites with $L \lesssim 10^4 L_{\odot}$ predicted to be observed in the DES and PS1 footprints is significantly lower than the data, similar to the alternative DM models shown in Figure 1 of Nadler et al. (2021c).

This procedure yields our fiducial and conservative 95% confidence constraints of $\tau > 18 \text{ Gyr}$ for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $\tau > 29 \text{ Gyr}$ for $V_{\text{kick}} = 40 \text{ kms}^{-1}$. The SHMF suppression

<footnote n="59">⁵⁹ Although our DDM simulations sample down to $\tau = 10 \text{ Gyr}$, we use a slightly higher lower-limit of the $\log \tau$ prior, below which the marginalized posterior is flat and nearly zero; this is a conservative choice.</footnote>
<table><caption>**Table 1**
Prior Distributions and 95% Credible Intervals Derived from the Marginalized Posterior for Each Parameter from Our DDM Fits to the MW Satellite Population for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $V_{\text{kick}} = 40 \text{ kms}^{-1}$ (Figure 5)</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Prior Distribution</th><th colspan="2">95% Credible Interval</th></tr><tr><th><i>V</i><sub>kick</sub> = 20 kms<sup>-1</sup></th><th><i>V</i><sub>kick</sub> = 40 kms<sup>-1</sup></th></tr></thead><tbody><tr><td>α</td><td>arctan α ~ unif(-1.1, -0.9)</td><td>[-1.452, -1.378]</td><td>[-1.454, -1.381]</td></tr><tr><td>σ<sub>M</sub></td><td>σ<sub>M</sub> ~ unif(0, 2) dex</td><td>[0.00*, 0.36] dex</td><td>[0.00*, 0.37] dex</td></tr><tr><td>log(<i>M</i><sub>50</sub>/<i>M</i><sub>☉</sub>)</td><td>log(<i>M</i><sub>50</sub>/<i>M</i><sub>☉</sub>) ~ unif(7.5, 11.0)</td><td>[7.50*, 8.01]</td><td>[7.50*, 7.99]</td></tr><tr><td><i>B</i></td><td>ln <i>B</i> ~ <i>N</i>(μ = 1.0, σ = 0.5)</td><td>[0.33, 1.84]</td><td>[0.33, 1.75]</td></tr><tr><td>σ<sub>gal</sub></td><td>σ<sub>gal</sub> ~ unif(0, 1) dex</td><td>[0.03*, 0.79] dex</td><td>[0.03*, 0.76] dex</td></tr><tr><td><i>A</i></td><td><i>A</i> ~ unif(0.0, 0.5)kpc</td><td>[15, 93] pc</td><td>[16, 76] pc</td></tr><tr><td>σ<sub>log R</sub></td><td>σ<sub>log R</sub> ~ unif(0, 2) dex</td><td>[0.01*, 0.86] dex</td><td>[0.00*, 0.75] dex</td></tr><tr><td>n</td><td>n ~ <i>N</i>(μ = 1.0, σ = 0.5)</td><td>[0.50, 2.00<sup>*</sup>]</td><td>[0.56, 2.00<sup>*</sup>]</td></tr><tr><td>log(τ/ Gyr)</td><td>log(τ/ Gyr) ~ unif(1.3, 1.9)</td><td>[1.46, 1.90<sup>*</sup>]</td><td>[1.63, 1.90<sup>*</sup>]</td></tr></tbody></table>

**Note.** The first eight parameters describe the galaxy-halo connection model used to associate satellite galaxies with subhalos in our MW-like zoom-in simulations (Nadler et al. 2020b), and the final parameter corresponds to the DDM particle lifetime. Definitions for each parameter are given in Appendix D. Asterisks denote prior-constrained limits. See Table 2 of Nadler et al. (2020b) for motivations for the prior distributions of the eight parameters describing the galaxy-halo connection model; note that the prior on *n* is bounded between 0.0 and 2.0 for convergence.

relative to CDM and the predicted luminosity function of DES and PS1 satellites for each of these models are shown in the left and right panels of Figure 6, respectively. The right panel of Figure 6 also compares these predictions to the observed luminosity function, demonstrating that the DDM models that our analysis rules out yield significantly fewer ultra-faint satellites than observed by DES and PS1, after accounting for observational selection effects and conservatively marginalizing over modeling uncertainties.

Figure 7 shows these lower limits on log(τ/ Gyr) for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$ alongside the preferred region of DDM parameter space that potentially alleviates the $S_8$ tension (Abellán et al. 2022). As discussed in detail in Section 6.1, these constraints very conservatively exclude roughly half of the DDM parameter space favored to resolve $S_8$ tension (Abellán et al. 2022) and the $H_0$ tension (Vattis et al. 2019). Our results are conservative in this context because our fits are performed at extremely low $V_{\text{kick}}$ relative to the typical values used to alleviate these cosmological tensions. Moreover, our fit only directly incorporates the effects of DDM physics on subhalo and satellite abundances, and therefore does not leverage the reduced central densities of DDM subhalos, which may yield additional constraining power.

# 6. Discussion

In this section, we place our DDM constraints in the context of previous studies (Section 6.1) and compare our constraints to those derived for other dark matter models using a similar analysis of the MW satellite population (Section 6.2).

## 6.1. Comparison to Previous Studies

We now place our results in the context of previous DDM studies, including (but not limited to) studies of MW satellite galaxies. We reiterate that our analysis makes conservative assumptions regarding both the microphysics of DM decays and their impact on structure formation. In particular, our constraints only directly apply to two-body decays that yield a cold, stable DM daughter particle, and to DDM models with sufficiently long lifetimes such that small-scale structure is only significantly affected relative to CDM at late times. We therefore caution that it is not straightforward to compare our

results with limits from DDM models with different decay mechanisms or to limits on short-lifetime decays that affect the pre-recombination universe. As a result, the following discussion focuses on limits derived for the same family of DDM models considered in our analysis. Beyond these cases, mapping the suppression of the subhalo mass function for our fiducial two-body decay model to that in a single-body decay would allow for a more direct comparison with many large-scale structure analyses (e.g., Chen et al. 2021; Hubert et al. 2021); this is left to a future work.

### 6.1.1. Limits from Milky Way Satellites

Several authors have studied the impact of decays on MW satellite galaxies for the same class of DDM models we consider. In particular, Peter & Benson (2010) compared the population of classical and Sloan Digital Sky Survey-discovered MW satellite galaxies to DDM predictions. These authors found that $\tau \lesssim 30$ Gyr is ruled out for $20 \text{ kms}^{-1} \lesssim V_{\text{kick}} \lesssim 200 \text{ kms}^{-1}$, where specific results depend on assumptions regarding the star formation histories of satellites and the evolution of subhalos and satellites in their semi-analytic model in detail.<fnref n="60" /> These limits are consistent with and slightly weaker than our fiducial constraints.

The most stringent limits from Peter & Benson (2010) are driven by the inferred mass enclosed within 300 pc of each MW satellite rather than the overall abundance of these systems. Thus, these limits are mainly set by the mass loss and reduction in central densities of DDM subhalos, rather than their enhanced disruption relative to CDM. From an observational standpoint, the enclosed mass depends on the measured stellar velocity dispersion for each satellite and a mapping between this quantity and a halo mass proxy (e.g., $V_{\text{max}}$), both of which are accompanied by significant systematic uncertainties.

On the other hand, our DDM constraints are driven by the abundance of confirmed and candidate dwarf galaxies in the DES and PS1 data sets as a function of absolute magnitude, half-light radius, and heliocentric distance, which we leverage in a statistical forward model. Integrating comparisons of predicted and inferred central dynamical masses into our model

<footnote n="60">⁶⁰ Peter et al. (2010a) synthesized the results of Peter et al. (2010b) and Peter & Benson (2010), reporting that $\tau \lesssim 40$ Gyr is ruled out for $V_{\text{kick}} > 20 \text{ kms}^{-1}$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Posterior distributions from our DDM fits to the DES and PS1 satellite populations for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (green) and $V_{\text{kick}} = 40 \text{ kms}^{-1}$ (purple). The effects of dark matter decays are more pronounced due to the greater kick velocity in the $V_{\text{kick}} = 40 \text{ kms}^{-1}$ case, raising the lower bound on the decay lifetime relative to the $V_{\text{kick}} = 20 \text{ kms}^{-1}$ case. Constraints on the eight galaxy–halo connection parameters are consistent for both values of $V_{\text{kick}}$. These parameters govern the faint-end slope ($\alpha$) and scatter ($\sigma_M$) of the satellite abundance-matching relation, the peak subhalo mass at which 50% of subhalos host galaxies ($\log M_{50}$, in units of $M_\odot$), the efficiency of subhalo disruption due to the Galactic disk ($B$), the scatter in the galaxy occupation fraction ($\sigma_{\text{gal}}$), the amplitude ($A$), scatter ($\sigma_{\log R}$), and power-law index ($n$) of the satellite–subhalo size relation, and the dark matter particle lifetime ($\log \tau$, in units of gigayears). Note that $\sigma_M$, $\sigma_{\text{gal}}$, and $\sigma_{\log R}$ are reported in dex and $A$ is reported in parsecs. Definitions for each parameter are given in Appendix D.</figcaption></figure>

is an important area for future work that will likely yield even more stringent constraints on DDM.

Several other authors have studied whether the family of DDM models that we analyze can reconcile the apparent tension between the predicted and inferred density profiles of dwarf galaxy halos. In particular, Wang et al. (2014) found that lifetimes $\tau \sim \mathcal{O}(10 \text{ Gyr})$ and kick velocities $V_{\text{kick}} \sim 20 \text{ kms}^{-1}$ are consistent with classical MW satellite galaxies. Meanwhile,

Chen & Chu (2021) studied the density profile of isolated dwarf galaxies, finding that $\tau \lesssim 7.0 \text{ Gyr}$ is needed to explain these measurements for $V_{\text{kick}} = 20 \text{ kms}^{-1}$; these results are consistent with those in Sánchez-Salcedo (2003) and Abdelqa-der & Melia (2008).

Comparing our constraints to these previous results implies that DDM models that significantly alter dwarf galaxy central densities are not viable because they simultaneously reduce
<figure><img src="image_8.png" /></figure>

**Figure 6.** Predictions for the DDM models disfavored at the 95% confidence level; these models represent the parameters for DDM sufficient to predict significant observational differences relative to observations. Left panel: SHMF suppression relative to CDM for the disfavored DDM models. The vertical dashed line indicates the 95% confidence upper limit on the lowest-mass halo inferred to host MW satellite galaxies (Nadler et al. 2020b). Right panel: predicted MW satellite galaxy luminosity functions for the disfavored DDM models compared to DES and PS1 observations, evaluated at the best-fit MW satellite model parameters from Nadler et al. (2020b). The shaded band illustrates the uncertainty of the DDM prediction due to the stochasticity of our galaxy-halo connection model and the limited number of simulations used in our analysis; the size of the uncertainty is very similar to that in CDM. Note that this panel is a simple one-dimensional representation of our MW satellite and DM model fit to the luminosity, size, and spatial distribution of satellites in the DES and PS1 survey footprints. The comparison of our CDM model to data is described in Nadler et al. (2020b).

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

**Figure 7.** Lower limits from our fits for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ and $40 \text{ kms}^{-1}$ are shown as points with black arrows indicating allowed parameter space. The corresponding excluded region is identified with black crosshatching. The region of DDM parameter space that potentially alleviates the $S_8$ tension (Abellán et al. 2022) is shown by the lavender contour, where dark (light) shading shows 68% (95%) confidence intervals. Note that the posterior distribution from Abellán et al. (2022) is limited by their prior range and does not extend to values of $V_{\text{kick}}$ as low as those studied in this paper. The typical magnitudes of peak velocities for MW satellite halos and host halo are indicated by vertical dashed gray lines.

MW satellite galaxy abundances to unacceptable levels. Similar conclusions have been drawn for warm and fuzzy DM, which drastically suppress the abundance of ultra-faint dwarf galaxies before they yield an observable impact on the density profiles of brighter dwarfs (Macciò et al. 2012; Safarzadeh & Spergel 2020; Nadler et al. 2021c).

### 6.1.2. Limits from Galaxy Clusters

Peter et al. (2010b) used the halo mass-concentration relation and mass function derived from galaxy clusters to estimate limits on the DDM kick velocity and lifetime, finding

that galaxy cluster observations rule out decay times less than a few times the age of the universe for kick velocities greater than ~100 kms⁻¹. These results are qualitatively consistent with our constraints in Figure 7, in the sense that larger values of $V_{\text{kick}}$ result in more stringent constraints on $\tau$, and they disfavor regions of parameter space that lie along this degeneracy at larger $V_{\text{kick}}$ than sampled by our analysis. As noted by Peter et al. (2010b), these results are approximate and warrant a detailed statistical analysis based on cosmological DDM simulations.

### 6.1.3. Limits from the Ly$\alpha$ Forest

Aside from the impact of DDM microphysics on low-mass subhalos at late times, DM decays can leave a substantial imprint on small-scale structure throughout cosmic history. For example, Wang et al. (2013) used observations of the Ly$\alpha$ forest to exclude $\tau \lesssim 10$ Gyr and $V_{\text{kick}} \gtrsim 30–70 \text{ kms}^{-1}$ for the same class of DDM models we consider. These constraints are again consistent with, and weaker than, our fiducial results, which is reasonable given that recent MW satellite constraints on WDM are more stringent than WDM constraints derived using the Ly$\alpha$ forest data considered in Wang et al. (2013; e.g., see Nadler et al. 2021c). As discussed in Wang et al. (2013), more precise limits likely entail a joint inference of intergalactic medium (IGM) properties and small-scale clustering based on hydrodynamic DDM simulations because IGM properties are partially degenerate with DM properties (e.g., Garzilli et al. 2021).

We note that low-redshift tracers of small-scale structure, including observations of ultra-faint dwarf galaxies, stellar streams, and strong gravitational lenses, are well suited to test DDM physics because the impact of DM decays on small-scale structure becomes more severe at late times. Thus, a detailed study of the synergies between small-scale structure probes that sample distinct scales and redshifts, including ultra-faint dwarf galaxies, stellar streams, strong gravitational lensing, and the Ly$\alpha$ forest, is particularly interesting in the context of DDM and promises to yield precise joint measurements (see
Enzi et al. 2021; Nadler et al. 2021b for examples of joint WDM constraints).

### 6.1.4. Limits from Other Cosmological Probes

DDM has recently gained interest as a potential solution to the $H_0$ and $\sigma_8$ tensions because it can potentially reduce the late-time expansion rate and matter power spectrum without strongly affecting early universe observables including the CMB. In particular, Vattis et al. (2019) found that DDM can resolve the $H_0$ tension for $\tau=35$ Gyr and $\epsilon=0.16$ ($V_{\text{kick}} \approx 10^4$ kms$^{-1}$) by combining late-time measurements of the expansion rate including distance-ladder, baryonic acoustic oscillation (BAO), quasar, and Ly$\alpha$ auto- and cross-correlation data. However, subsequent analyses have shown that this model is excluded by additional data sets. Specifically, Haridasu & Viel (2020) used Type Ia supernovae, BAO, and time delay measurements of gravitationally lensed quasars with priors set by the CMB to derive a limit of $\tau > 9$ Gyr for two-body decays at the maximum allowed value of $\epsilon = 0.5$ ($V_{\text{kick}} \approx 10^5$ kms$^{-1}$) in their analysis. Similarly, Clark et al. (2021b) used Planck CMB power spectra measurements to derive $\tau > 1000$ Gyr for $\epsilon = 0.5$. Neither of these analyses found that viable DDM models (at such high mass splitting ratios) can significantly affect the present-day expansion rate.

Meanwhile, Abellán et al. (2022) found that DDM can resolve the $S_8$ tension for $\tau = 56$ Gyr and $\epsilon \approx 0.007$ ($V_{\text{kick}} \approx 10^3$ kms$^{-1}$) by combining Planck CMB lensing and power spectra, BAO data, and Type Ia supernovae with KIDS1000+BOSS+2dfLenS measurements of $S_8$. Abellán et al. (2021) found that this conclusion is robust to the inclusion of additional cosmological and experimental constraints, and that it can potentially explain the anomalous Xenon1T electron recoil excess (Aprile et al. 2020; Kannike et al. 2020).

Our analysis excludes $\tau < 18$ Gyr for $V_{\text{kick}} = 20$ kms$^{-1}$ ($\epsilon \approx 10^{-4}$); this value of $V_{\text{kick}}$ is significantly lower than those preferred by the $H_0$ and $S_8$ analyses described above. Because the impact of DDM on small-scale structure becomes more severe for larger values of $V_{\text{kick}}$, we can extremely conservatively interpret this fiducial constraint on $\tau$ as a limit on the DDM models considered by, e.g., Vattis et al. (2019) and Abellán et al. (2022). As shown by the crosshatched regions in Figure 7, this disfavors the preferred DDM parameter space reported by Abellán et al. (2022) for low values of $\tau$ and $V_{\text{kick}}$.<fnref n="61" /> In practice, we expect our constraints on $\tau$ to become much more stringent for these models. For example, because halos with characteristic virial velocities of $\mathcal{O}(V_{\text{kick}})$ are affected by DM decays, we expect that the $V_{\text{kick}} \approx 10^3$ kms$^{-1}$ model considered by Abellán et al. (2022) would significantly affect the structure and abundance of halos that host MW-mass and even larger galaxies. In turn, structure on dwarf galaxy scales is likely completely different than in CDM (and therefore incompatible with MW satellite data) for all values of $\tau$ preferred by these analyses, though dedicated future work is necessary to quantify this claim.

## 6.2. Comparison to Constraints on Other Dark Matter Models from Milky Way Satellites

We now place our DDM results in the context of constraints on other DM models derived from the MW satellite population

and discuss implications for small-scale structure observables. We limit our quantitative comparison to the WDM constraints from Nadler et al. (2021c) because this study used an identical modeling framework with the exception of the assumed SHMF suppression; however, we note that many other studies have used the MW satellite population to derive WDM constraints (e.g., Macciò & Fontanot 2010; Polisensky & Ricotti 2011; Anderhalden et al. 2013; Kennedy et al. 2014; Kim et al. 2018; Nadler et al. 2019a; Dekker et al. 2021; Newton et al. 2021).<fnref n="62" />

Nadler et al. (2021c) found that thermal relic WDM masses below 6.5 keV are excluded by the same MW satellite census that we consider at 95% confidence after marginalizing over an identical set of galaxy-halo connection and MW halo mass parameters. This WDM model suppresses subhalo abundances by ~25% relative to CDM at the minimum observed halo mass scale of $3.2 \times 10^8 M_{\odot}$; the SHMF declines rapidly at smaller masses, reaching ~75% suppression relative to CDM at its half-mode mass of $3.8 \times 10^7 M_{\odot}$. For comparison, the DDM models we rule out at 95% confidence suppress the subhalo mass function by ~50% at the minimum observed halo mass scale, and their SHMFs decline less rapidly than for the ruled-out WDM model at lower peak subhalo masses. We reiterate that our fit to the DDM SHMF suppression is only valid for subhalos above a present-day mass threshold of ~$2 \times 10^7 M_{\odot}$, and that robust estimates for subhalo abundances at lower masses require higher-resolution simulations.

Although regions of parameter space for both WDM-like and DDM models are excluded by the abundance of known MW satellites, these scenarios make distinct predictions for other dwarf galaxy and small-scale structure observables. In particular, decays can significantly deplete both low-mass isolated halos and subhalos of dark matter at late times, lowering predicted mass-to-light ratios for both satellite and field dwarf galaxies relative to CDM. Similarly, mass loss and momentum transfer due to decays reduce halos' central densities, which future observations of dynamical tracers in MW satellites will better inform (Simon et al. 2019). On the other hand, WDM halos with masses well above the half-mode scale do not significantly differ in present-day mass relative to CDM, though delayed formation lowers their concentration (Bose et al. 2016; Stücker et al. 2022). Combining our constraints with small-scale structure observables that are sensitive to the abundance and internal structure of halos and subhalos at late times, including strong gravitational lensing (e.g., Minor et al. 2017; Gilman et al. 2020a, 2020b; Hsueh et al. 2020), will therefore help differentiate these classes of models.

In addition to the WDM-like models discussed above, it is also interesting to contrast the effects of DDM with those of self-interacting DM (SIDM), which can also suppress the abundance of low-mass subhalos while altering their density profiles (Vogelsberger et al. 2012; Zavala et al. 2013; Tulin & Yu 2018; Robles et al. 2019; Nadler et al. 2020a, 2021a). Unlike DDM, which (to first order) equally depletes both isolated halos and subhalos of dark matter relative to CDM, subhalos' mass loss and disruption in SIDM are closely tied to their orbital histories (e.g., Dooley et al. 2016; Jiang et al. 2021).

<footnote n="61">⁶¹ Note that the posterior from Abellán et al. (2022) is prior-limited and not strongly constrained as a function of $\tau$.</footnote>

<footnote n="62">⁶² Nadler et al. (2021c) also derived constraints on DM-baryon interactions and fuzzy DM based on the suppression of the linear matter power spectrum and low-mass halo abundances in these models. We consider these “WDM-like” models for this discussion, although they may have distinct effects on the MW satellite population in detail.</footnote>
Thus, comparing the abundance and internal structure of isolated halos and subhalos—for example, through joint analyses of field and satellite dwarf galaxy populations or of line-of-sight and subhalo perturbations in strong lensing data—will help disentangle these forms of dynamical DM microphysics. These differences can also potentially be tested by comparing the abundance and properties of LMC-associated MW satellites with the remainder of the MW satellite population (Nadler et al. 2021a).

## 7. Conclusions

We have used a state-of-the-art census of the MW satellite galaxy population to set robust and stringent constraints on the DM particle lifetime that are among the most robust and stringent to date while making conservative assumptions about the decay mechanism (i.e., late-time decays that include a stable CDM-like daughter product). In particular, we combined cosmological zoom-in simulations of DDM with a forward model of the MW satellite population to jointly infer the connection between these galaxies and their DM subhalos and the potential impact of DM decays on the abundance of these systems.

For DM that undergoes late-time two-body decays to a massless dark radiation species and a cold, stable daughter DM particle, we find that:

1. For DM particle lifetimes of $O(t_H)$ and recoil kick velocities of $O(10 \text{ kms}^{-1})$, the smallest subhalos that host ultra-faint dwarf galaxies ($M_{\text{peak}} \sim 10^8 M_{\odot}$) can lose a significant fraction of their peak mass due to decays alone;

2. DM decays can suppress the abundance of surviving subhalos above a minimum mass threshold of $\approx 4 \times 10^7 M_{\odot}$ at the $\sim 50\%$ level relative to CDM. This suppression is approximately independent of Galactocentric radius, and Equation (8) provides a fitting function for this effect derived from cosmological zoom-in simulations;

3. The population of MW satellite galaxies observed by DES and PS1 excludes DDM models with decay lifetime $\tau < 18 \text{ Gyr}$ (29 Gyr) for $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (40 kms$^{-1}$) at 95% confidence;

4. These constraints can be conservatively extrapolated to higher $V_{\text{kick}}$ values to exclude approximately half of the DDM parameter space preferred to alleviate the $H_0$ and $S_8$ tensions. These constraints are expected to become more stringent for the $V_{\text{kick}}$ values considered in those analyses and with the inclusion of MW satellite stellar velocity dispersion measurements;

5. Combining our DDM constraints based on MW satellites with complementary small-scale structure probes at low and high redshifts—including field dwarf galaxy luminosity functions, strong gravitational lensing, and the Ly$\alpha$ forest—will help differentiate the effects of decays from other DM microphysics.

Our analysis only directly leverages the reduction of DDM subhalo abundances and its impact on MW satellite abundances, rather than the (potentially observable) effects on the internal dynamics of satellite galaxies. Nonetheless, our results are consistent with and more stringent than previous limits driven by MW satellite stellar velocity dispersion measurements. Future work that combines our approach with the inferred dynamical masses and density profiles of MW

satellites promises to further improve DDM constraints, as do future detections of dwarf galaxies within and beyond the MW. These observational advances will, respectively, be enabled by forthcoming spectroscopic facilities and giant segmented mirror telescopes (Simon et al. 2019) and observational facilities including the Vera C. Rubin Observatory (Drlica-Wagner et al. 2019; Mutlu-Pakdil et al. 2021) and the Nancy Grace Roman Space Telescope.

Although other cosmological observables disfavor DDM as a solution to the $H_0$ and $S_8$ tensions (e.g., Haridasu & Viel 2020; Clark et al. 2021b),<fnref n="63" /> we emphasize that our results inform DDM physics as a solution to these tensions in a way that is complementary to expansion history and large-scale structure probes. In particular, large-scale structure probes are primarily sensitive to the DM lifetime, regardless of the microphysical decay mechanism, while the abundance of low-mass DM subhalos traced by MW satellite galaxies is sensitive to both the lifetime of the DM particle and the decay mechanism as encapsulated by $V_{\text{kick}}$ in our model. Thus, to more fully inform DDM physics, it is crucial to combine observables that cover a wide range of cosmological epochs and scales, including probes of small-scale structure.

The authors thank Guillermo Abellán for sharing their results from Abellán et al. (2022) and Andrew Benson for comments on the manuscript.

This research received support from the National Science Foundation (NSF) under grant No. NSF DGE-1656518 through the NSF Graduate Research Fellowship received by S.M. and E.O.N. and from the U.S. Department of Energy under contract No. DE-AC02-76SF00515 to SLAC National Accelerator Laboratory. Y.-Y.M. is supported by NASA through the 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, Inc., under NASA contract NAS5-26555.

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 arXiv.org (https://arXiv.org) and NASA's Astrophysics Data System for bibliographic information. This research made use of adstex (https://github.com/yymao/adstex). This research made use of the cubehelix color scheme (Green 2011).

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

<footnote n="63">⁶³ We note that Clark et al. (2021b) approximated the perturbations of the massive daughter particles, which were found to be essential to account for the effects on the $S_8$ parameter by Abellán et al. (2022).</footnote>
à 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 Labora-
tory, the University of California at Santa Cruz, the University
of Cambridge, Centro de Investigaciones Energéticas, Med-
ioambientales y Tecnológicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Con-
sortium, the University of Edinburgh, the Eidgenössische
Technische Hochschule (ETH) Zürich, Fermi National Accel-
erator 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, NSF's NOIRLab, 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.

Based in part on observations at Cerro Tololo Inter-
American Observatory at NSF's NOIRLab (NOIRLab Prop.
ID 2012B-0001; PI: J. Frieman), which is managed 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 No. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MICINN under grants ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. 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) do e-Universo (CNPq grant 465376/2014-2).

This manuscript has been authored by 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 NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), seaborn (Waskom 2021), emcee (Foreman-Mackey et al. 2013), ChainConsumer,<fnref n="64" /> incredible,<fnref n="65" /> GADGET (Springel 2005), ROCKSTAR (Behroozi et al. 2013a), CONSISTENT-TREES (Behroozi et al. 2013b).

Appendix A
Milky Way–like Simulations

To perform the inference in Section 5, we use cosmological
zoom-in simulations of two MW-like halos originally presented

in Mao et al. (2015) and studied in Nadler et al. (2020b, 2021c) to analyze the MW satellite population. Furthermore, to validate our DDM subhalo population predictions derived from the expanded Wang et al. (2014) simulation suite, we perform DDM resimulations of these MW-like systems. These host halos are selected based on mass and concentration estimates for the MW and due to the presence of realistic LMC analogs and early Gaia-Sausage-Enceladus-like merger events with the properties described in Nadler et al. (2020b). We simulate both of these systems in DDM models with ($\tau$/ Gyr, $V_{\text{kick}}$/ kms$^{-1}$) $\in$ {(20, 20, (40, 30)}), roughly corresponding to the 95% confidence level constraints determined in Section 5.

To perform the DDM resimulations, we use the same modified version of the GADGET-2 and GADGET-3 N-body codes from Peter et al. (2010b). The original CDM MW-like simulations and the DDM resimulations are run with $\Omega_M = 0.286$, $\Omega_\Lambda = 0.714$, $n_s = 1$, $h = 0.7$, and $\sigma_8 = 0.82$ (Hinshaw et al. 2013).<sup>66</sup> The highest-resolution region is simulated with a Plummer-equivalent force softening of 170 pch<sup>-1</sup> and a particle mass of 3.0 × 10<sup>5</sup> M<sub>☉</sub>h<sup>-1</sup>. We analyze these simulations using the ROCKSTAR (Behroozi et al. 2013a) halo finder and CONSISTENT-TREES (Behroozi et al. 2013b) merger tree code.

We note that the differences in low-mass subhalo abundances introduced by changes to the numerical and cosmological parameters with respect to the Wang et al. (2014) simulations are minor compared to the differences between the DDM and CDM simulations within each suite (see, e.g., Dooley et al. 2014 for a study of the impact of cosmological parameters on subhalo statistics). The largest difference in cosmological parameters is the value of $n_s$; however, once the orbital phase of the LMC is fixed, the SHMF in the $n_s = 1$ simulations we use is enhanced by ~10% at all subhalo masses relevant for our study compared to simulations with $n_s = 0.96$. This allows for more severe SHMF suppression due to DDM when comparing to the data, meaning that the constraints we derive are conservative.

Appendix B
Subhalo Mass Function Suppression Validation

We use our MW-like resimulations to test that the impact of
DDM on low-mass subhalos derived from the Wang et al.
(2014) simulations—and particularly the suppression of the
DDM SHMF—is applicable to hosts of similar masses that
specifically resemble the MW system. To do so, we compare
predictions from (i) Equation (8) and (ii) an interpolation of the
SHMF suppression from the Wang et al. (2014) simulations to
the SHMF suppression determined directly from our MW-like
resimulations. For the second comparison, we use piecewise
linear interpolation to smoothly connect the SHMF suppres-
sion, N(>M<sub>peak</sub>)<sub>DDM</sub>/N(>M<sub>peak</sub>)<sub>CDM</sub>, between the points in
DDM parameter space simulated by Wang et al. (2014):
(τ/ Gyr, V<sub>kick</sub>/ kms<sup>-1</sup>) ∈ {(10, 20, (20, 20), (20, 40), (40, 20),
(40, 40), (80, 40)}}. This provides an alternative estimate of the
SHMF suppression in regions of DDM parameter space that
were not directly simulated. The interpolating function is
compared to the Wang et al. (2014) simulation results in
Figure 8, which demonstrates agreement at the ~1σ level.

As shown in Figure 9, the DDM SHMF provided by
the interpolation function also matches that derived from our

<footnote n="64">⁶⁴ https://github.com/Samreay/ChainConsumer</footnote>

<footnote n="65">⁶⁵ https://github.com/abmantz/incredible</footnote>

<footnote n="66">⁶⁶ Note that the $n_s$ value is in fact 1 instead of 0.96 as stated in previous studies (e.g., Mao et al. 2015).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** SHMF suppression at $z=0$ for the MW-mass host halo in our expanded suite of zoom-in simulations based on Wang et al. (2014) measured directly from the simulations (solid lines) and using the interpolating procedure described in Appendix B (dotted lines). Results are shown for DDM models with $\tau = 10, 20, \text{ and } 40$ Gyr (from yellow to purple), with $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (left panel) and $40 \text{ kms}^{-1}$ (right panel) as a function of peak subhalo virial mass $M_{\text{peak}}$. All SHMFs are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ kms}^{-1}$ and $V_{\text{max}} > 9 \text{ kms}^{-1}$. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements; the SHMF suppression predicted by our interpolating procedure is consistent with that measured directly from the simulations at the ~1$\sigma$ level.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Average SHMF suppression at $z=0$ for the MW-like host halos in our MW-like resimulations (described in Appendix A) measured directly from the resimulations (solid lines) and using the interpolating procedure described in Appendix B. Results are shown for the resimulated DDM models with $\tau = 20$ Gyr, $V_{\text{kick}} = 20 \text{ kms}^{-1}$ (left panel) and $\tau = 30$ Gyr, $V_{\text{kick}} = 40 \text{ kms}^{-1}$ (right panel) as a function of peak subhalo virial mass $M_{\text{peak}}$. All SHMFs are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ kms}^{-1}$ and $V_{\text{max}} > 9 \text{ kms}^{-1}$. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements. The SHMF suppression predicted by our interpolating procedure based on the Wang et al. (2014) simulations is consistent with that measured from our MW-like resimulations at the ~1$\sigma$ level.</figcaption></figure>

MW-like resimulations at the ~1$\sigma$ level. We note that, while one of these resimulations covers the same point in $\tau$, $V_{\text{kick}}$ space as one of the Wang et al. (2014) simulations, the second does not; thus, the interpolating function is used to compare with our MW-like resimulations in both cases.

Combining the results from Figure 8 and Figure 9 demonstrates that the DDM SHMF suppression derived and

adopted in our fiducial analysis based on the Wang et al. (2014) simulations is consistent with that from simulations of systems specifically chosen to resemble the MW, lending confidence to our constraints. Finally, we validate the behavior of the subhalo radial distribution of the Wang et al. (2014) simulations (Figure 4) by comparing to the distribution for our MW-like simulations in Figure 10.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Average subhalo radial distributions at $z=0$ as a function of distance from the center of the host halo for the CDM (dashed) and DDM (solid) resimulations of MW-like halos described in Appendix A. Results are shown for the resimulated DDM models with $\tau = 20$ Gyr, $V_{\text{kick}} = 20 \text{ km s}^{-1}$ (red) and $\tau = 30$ Gyr, $V_{\text{kick}} = 40 \text{ km s}^{-1}$ (purple). The bottom panels show the ratio of the radial distribution in each DDM model relative to CDM. All results are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ km s}^{-1}$ and $V_{\text{max}} > 9 \text{ km s}^{-1}$. Shorter DDM decay lifetimes result in fewer surviving subhalos relative to CDM. This effect is more pronounced for models with higher kick velocities and does not strongly depend on radius from the center of the host halo. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements.</figcaption></figure>

# Appendix C

## Evolution of the DDM Subhalo and Halo Mass Functions

We use the DDM resimulations of MW-like systems presented in Appendix A to study the evolution of the DDM subhalo and halo mass functions. In particular, Figure 11 shows the evolution of the average SHMF for subhalos of the two MW-like hosts at $z=0, 1$, and $2$. The DDM SHMFs are consistent with the corresponding CDM SHMFs for $z \gtrsim 2$, accounting for Poisson uncertainties, and only become significantly suppressed at later times. This reflects the combined effects of mass loss due to DM decays, which can push subhalos below the mass resolution limit of our simulations (Section 3), and the enhanced tidal disruption of these systems with reduced central densities relative to CDM.

This late-time suppression differentiates DDM from WDM-like models that suppress the linear matter power spectrum.

Finally, Figure 12 shows peak velocity functions at $z=0, 1$, and $2$ for *isolated* halos surrounding our resimulated MW-like systems, most of which lie within $a \sim 3$ Mpc radius from the center of the host halo that contains $\sim 90\%$ of the highest-resolution particles (Wang et al. 2021).<fnref n="67" /> Like the SHMF, the suppression of the isolated DDM halo mass function only sets in significantly at late times. The suppression of isolated halo abundances is only slightly less severe than that for subhalos, consistent with mass loss due to DM decays driving the disruption. This effect—i.e., severe mass loss for isolated halos —differentiates DDM from models like self-interacting DM in which late-time physics preferentially disrupts subhalos at late times (e.g., Tulin & Yu 2018; Nadler et al. 2020a).

<footnote n="67">⁶⁷ In particular, we analyze isolated halos by choosing systems with ROCK-STAR *upid* equal to -1.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the subhalo $M_{\text{peak}}$ functions for the CDM and DDM resimulations of MW-like halos described in Appendix A. Results are shown at $z=0$ (left panel), $z=1$ (middle panel), and $z=2$ (right panel), and the bottom panels show the corresponding SHMF suppression. Results are shown for the resimulated DDM models with $\tau = 20$ Gyr, $V_{\text{kick}} = 20 \text{ km s}^{-1}$ (red) and $\tau = 30$ Gyr, $V_{\text{kick}} = 40 \text{ km s}^{-1}$ (purple). All SHMFs are restricted to subhalos above a conservative resolution threshold of $V_{\text{peak}} > 10 \text{ km s}^{-1}$ and $V_{\text{max}} > 9 \text{ km s}^{-1}$; note that $V_{\text{peak}}$ is typically achieved at much earlier times than shown here ($z \sim 4$). The suppression of the DDM SHMF only sets in significantly at late times, consistent with the intuition developed in Section 3. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of the isolated halo $M_{\text{peak}}$ functions for the CDM and DDM resimulations of MW-like halos described in Appendix A. These measurements are presented analogously to Figure 11. The abundance of isolated halos is significantly suppressed at late times in DDM and is slightly less severe than the suppression of subhalo abundances, consistent with mass loss due to DM decays driving this effect. Shaded bands indicate 68% confidence interval Poisson uncertainties on the simulation measurements.</figcaption></figure>
<table><caption>**Table 2**
Galaxy–Halo Connection Model Parameters</caption><thead><tr><th>Parameter</th><th>Physical Interpretation</th><th>Symbol</th><th>Units</th></tr></thead><tbody><tr><td>Faint-end slope</td><td>Power-law slope of satellite luminosity function</td><td>α</td><td>none</td></tr><tr><td>Luminosity scatter</td><td>Scatter in luminosity at fixed V<sub>peak</sub></td><td>σ<sub>M</sub></td><td>dex</td></tr><tr><td>50% occupation mass</td><td>Mass at which 50% of halos host galaxies</td><td>M<sub>50</sub></td><td>M<sub>☉</sub></td></tr><tr><td>Baryonic effects</td><td>Strength of subhalo disruption due to baryons</td><td>B</td><td>none</td></tr><tr><td>Occupation scatter</td><td>Scatter in galaxy occupation fraction</td><td>σ<sub>gal</sub></td><td>dex</td></tr><tr><td>Size amplitude</td><td>Amplitude of galaxy–halo size relation</td><td>A</td><td>pc</td></tr><tr><td>Size scatter</td><td>Scatter in half-light radius at fixed halo size</td><td>σ<sub>log R</sub></td><td>dex</td></tr><tr><td>Size power-law index</td><td>Power-law index of galaxy–halo size relation</td><td>n</td><td>none</td></tr><tr><td>Decay lifetime</td><td>Inverse of the DM particle decay rate</td><td>τ</td><td>Gyr</td></tr></tbody></table>

## Appendix D
Galaxy–Halo Connection Model and Parameters

We follow Nadler et al. (2020b) in modeling the galaxy–halo connection used in our forward model (Section 5.1) with eight free parameters and introduce one additional parameter for DDM particle lifetime. These parameters are summarized in Table 2.

### D.1. Satellite Luminosities

We employ an abundance-matching procedure that relates the absolute V-band magnitude of satellites, $M_V$, to the peak circular velocity of subhalos, $V_{peak}$ (Nadler et al. 2019b). This relation is extended to dim satellites by allowing the faint-end slope of the satellite luminosity function, $\alpha$, and the lognormal scatter in luminosity at fixed $V_{peak}$, $\sigma_M$, to be free parameters. While this abundance-matching model does not capture the entire star formation histories of ultra-faint dwarf galaxies, it is consistent with current MW satellite data (Drlica-Wagner et al. 2020).

### D.2. Satellite Sizes

The mean predicted size of each satellite at accretion is set according to

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

where $A$ and $n$ are free parameters corresponding to the amplitude and power-law index of the galaxy–halo size relation, respectively, and $R_{vir}$ denotes the subhalo virial radius as measured at accretion; $R_0 = 10$ kpc is a normalization constant.

In our inference, satellite sizes are drawn from a lognormal distribution with mean given by Equation (D1) and standard deviation $\sigma_{\log R}$, a further free parameter. While post-infall effects—including adiabatic decays—can shrink or enlarge satellites, Nadler et al. (2020b) found that our results are not sensitive to these effects while using a model for satellite size evolution due to tidal stripping; thus, these effects are not modeled here.

### D.3. Subhalo Disruption due to Baryonic Effects

We incorporate the effects of baryonic physics—particularly the tidal influence on the Galactic disk—on our simulated subhalo populations following Garrison-Kimmel et al. (2017) and Nadler et al. (2018). The strength of the disruption is modeled using the

free parameter $B$ for which $B = 1$ corresponds to fiducial hydrodynamical predictions (Nadler et al. 2018) and larger (smaller) values of $B$ correspond to more (less) effective subhalo disruption. For each subhalo, we set

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

where $p_{\text{disrupt},0}$ is the fiducial disruption probability given by the random forest algorithm of Nadler et al. (2018).

### D.4. Galaxy Formation Efficiency

We parameterize the fraction of halos that host galaxies of any mass—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_{50}}{\sqrt{2} \sigma_{\text{gal}}} \right) \right], \qquad (\text{D3})$$

where $M_{peak}$ is the largest virial mass a subhalo ever attains, $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 inference, $M_{50}$ and $\sigma_{gal}$ are free parameters.

### D.5. Orphan Satellites

We account for orphan satellites—subhalos that have been artificially disrupted by approaching the resolution limit of our simulations—by following the prescription of Nadler et al. (2019b), which identifies disrupted subhalos in each simulation, interpolations 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. The effective abundance of orphans is parameterized by setting their disruption probabilities equal to

$$p_{\text{disrupt}} \equiv (1 - a_{\text{acc}})^{\mathcal{O}}, \qquad (\text{D4})$$

where $a_{acc}$ is the final scale factor at which each subhalo enters the virial radius of the MW analog in the simulations, and $\mathcal{O}$ captures deviations from disruption probabilities in hydrodynamic simulations. We follow Nadler et al. (2019b) by fixing $\mathcal{O} = 1$; note that Nadler et al. (2020b) found that the results of a CDM fit to the observed MW satellite population are insensitive to the value of $\mathcal{O}$.

## ORCID iDs

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

E. O. Nadler https://orcid.org/0000-0002-1182-3825
R. H. Wechsler https://orcid.org/0000-0003-2229-011X

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

K. Bechtol https://orcid.org/0000-0001-8156-0429

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

D. Huterer https://orcid.org/0000-0001-6558-0112

T. S. Li https://orcid.org/0000-0002-9110-6163

Y.-Y. Mao https://orcid.org/0000-0002-1200-0820

C. E. Martínez-Vázquez https://orcid.org/0000-0002-9144-7726

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

B. Mutlu-Pakdil https://orcid.org/0000-0001-9649-4815

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

A. Peter https://orcid.org/0000-0002-8040-6785

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

L. Strigari https://orcid.org/0000-0001-5672-6079

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

M. Aguena https://orcid.org/0000-0001-5679-6747

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

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

D. Bacon https://orcid.org/0000-0002-2562-8537

E. Bertin https://orcid.org/0000-0002-3602-3664

S. Bocquet https://orcid.org/0000-0002-4900-805X

D. Brooks https://orcid.org/0000-0002-8458-5047

D. L. Burke https://orcid.org/0000-0003-1866-1950

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

M. Costanzi https://orcid.org/0000-0001-8158-1449

M. Crocce https://orcid.org/0000-0002-9745-6228

T. M. Davis https://orcid.org/0000-0002-4213-8783

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

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

I. Ferrero https://orcid.org/0000-0002-1295-1132

B. Flaugher https://orcid.org/0000-0002-2367-5049

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

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

G. Giannini https://orcid.org/0000-0002-3730-1750

D. Gruen https://orcid.org/0000-0003-3270-7644

R. A. Gruendl https://orcid.org/0000-0002-4588-6517

J. Gschwend https://orcid.org/0000-0003-3023-8362

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

S. R. Hinton https://orcid.org/0000-0003-2071-9349

D. L. Hollowood https://orcid.org/0000-0002-9369-4157

K. Honscheid https://orcid.org/0000-0002-6550-2023

D. J. James https://orcid.org/0000-0001-5160-4486

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

O. Lahav https://orcid.org/0000-0002-1134-9035

M. A. G. Maia https://orcid.org/0000-0001-9856-9307

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

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

R. Morgan https://orcid.org/0000-0002-7016-5471

R. L. C. Ogando https://orcid.org/0000-0003-2120-1154

F. Paz-Chinchón https://orcid.org/0000-0003-1339-2683

A. Pieres https://orcid.org/0000-0001-9186-6042

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

S. Serrano https://orcid.org/0000-0002-0211-2861

I. Sevilla-Noarbe https://orcid.org/0000-0002-1831-1953

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

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

C. To https://orcid.org/0000-0001-7836-2261

D. L. Tucker https://orcid.org/0000-0001-7211-5729

J. Weller https://orcid.org/0000-0002-8282-2010

References

Abdelqader, M., & Melia, F. 2008, MNRAS, 388, 1869

Abellán, G. F., Murgia, R., & Poulin, V. 2021, PhRvD, 104, 123533

Abellán, G. F., Murgia, R., Poulin, V., & Lavalle, J. 2022, PhRvD, 105, 063525

Ackermann, M., Ajello, M., Albert, A., et al. 2015, PhRvD, 91, 122002

Anderhalden, D., Schneider, A., Macciò, A. V., Diemand, J., & Bertone, G. 2013, JCAP, 2013, 014

Ando, S., Barik, S. K., Feng, Z., et al. 2021, PhRvD, 104, 023022

Aprile, E., Aalbers, J., Agostini, F., et al. 2020, PhRvD, 102, 072004

ATLAS Collaboration 2021, JHEP, 2021, 226

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

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

Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

Benito, M., Criado, J. C., Hütsi, G., Raidal, M., & Veermäe, H. 2020, PhRvD, 101, 103023

Bhattacharya, A., & Slatyer, T. R. 2019, JCAP, 2019, 029

Blackadder, G., & Koushiappas, S. M. 2014, PhRvD, 90, 103527

Blackadder, G., & Koushiappas, S. M. 2016, PhRvD, 93, 023510

Blinov, N., Keith, C., & Hooper, D. 2020, JCAP, 2020, 005

Bose, S., Hellwing, W. A., Frenk, C. S., et al. 2016, MNRAS, 455, 318

Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40

Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2012, MNRAS, 422, 1203

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

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

Cen, R. 2001, ApJ, 546, L77

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Chen, A., Huterer, D., Lee, S., et al. 2021, PhRvD, 103, 123528

Chen, J., & Chu, M. C. 2021, ApJ, 906, 86

Clark, S. J., Vattis, K., Fan, J., & Koushiappas, S. M. 2021a, arXiv:2110.09562

Clark, S. J., Vattis, K., & Koushiappas, S. M. 2021b, PhRvD, 103, 043014

Corasaniti, P. S., Agarwal, S., Marsh, D. J. E., & Das, S. 2017, PhRvD, 95, 083512

Dekker, A., Ando, S., Correa, C. A., & Ng, K. C. Y. 2021, arXiv:2111.13137

DES Collaboration 2005, IJMPA, 20, 3121

DES Collaboration 2016, MNRAS, 460, 1270

DES Collaboration 2018, PhRvD, 98, 043526

Dooley, G. A., Griffen, B. F., Zukin, P., et al. 2014, ApJ, 786, 50

Dooley, G. A., Peter, A. H. G., Vogelsberger, M., Zavala, J., & Frebel, A. 2016, MNRAS, 461, 710

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

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

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2015, JCAP, 2015, 067

Enzi, W., Murgia, R., Newton, O., et al. 2021, MNRAS, 506, 5848

Feng, J. L., Rajaraman, A., & Takayama, F. 2003a, PhRvL, 91, 011302

Feng, J. L., Rajaraman, A., & Takayama, F. 2003b, PhRvD, 68, 011302

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

Garrison-Kimmel, S., Boylan-Kolchin, M., Bullock, J. S., & Kirby, E. N. 2014, MNRAS, 444, 222

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

Garzilli, A., Magalich, A., Ruchayskiy, O., & Boyarsky, A. 2021, MNRAS, 502, 2356

Gaskins, J. M. 2016, ConPh, 57, 496

Gilman, D., Birrer, S., Nierenberg, A., et al. 2020a, MNRAS, 491, 6077

Gilman, D., Du, X., Benson, A., et al. 2020b, MNRAS, 492, L12

González, M. C., & Toro, N. 2021, arXiv:2108.13422

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

Green, D. A. 2011, BASI, 39, 289

Haridasu, B. S., & Viel, M. 2020, MNRAS, 497, 1757

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hsueh, J. W., Enzi, W., Vegetti, S., et al. 2020, MNRAS, 492, 3047

Hubert, J., Schneider, A., Potter, D., Stadel, J., & Giri, S. K. 2021, JCAP, 2021, 040

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

Ichiki, K., Oguri, M., & Takahashi, K. 2004, PhRvL, 93, 071302

Iršič, V., Viel, M., Haehnelt, M. G., Bolton, J. S., & Becker, G. D. 2017,
PhRvL, 119, 031302

Jarosik, N., Bennett, C. L., Dunkley, J., et al. 2011, ApJS, 192, 14

Jiang, F., Kaplinghat, M., Lisanti, M., & Slone, O. 2021, arXiv:2108.03243

Jordan, J. R., Kahn, Y., Krnjaic, G., Moschella, M., & Spitz, J. 2018, PhRvD, 98, 075020

Kannike, K., Raidal, M., Veermäe, H., Strumia, A., & Teresi, D. 2020, PhRvD, 102, 095002

Kaplinghat, M. 2005, PhRvD, 72, 063510

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, S. Y., & Peter, A. H. G. 2021, arXiv:2106.09050

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

Lovell, M. R. 2020, ApJ, 897, 147

Lovell, M. R., Frenk, C. S., Eke, V. R., et al. 2014, MNRAS, 439, 300

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

Macciò, A. V., Paduroiu, S., Anderhalden, D., Schneider, A., & Moore, B. 2012, MNRAS, 424, 1105

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

McDonald, P., Seljak, U., Burles, S., et al. 2006, ApJS, 163, 80

Menci, N., Sanchez, N. G., Castellano, M., & Grazian, A. 2016, ApJ, 818, 90

Minor, Q. E., Kaplinghat, M., & Li, N. 2017, ApJ, 845, 118

Mutlu-Pakdil, B., Sand, D. J., Crnojević, D., et al. 2021, ApJ, 918, 88

Nadler, E. O., Banerjee, A., Adhikari, S., Mao, Y.-Y., & Wechsler, R. H. 2020a, ApJ, 896, 112

Nadler, E. O., Banerjee, A., Adhikari, S., Mao, Y.-Y., & Wechsler, R. H. 2021a, ApJL, 920, L11

Nadler, E. O., Birrer, S., Gilman, D., et al. 2021b, ApJ, 917, 7

Nadler, E. O., Drlica-Wagner, A., Bechtol, K., et al. 2021c, PhRvL, 126, 091101

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

Nadler, E. O., Wechsler, R. H., Bechtol, K., et al. 2020b, ApJ, 893, 48

Newton, O., Leo, M., Cautun, M., et al. 2021, JCAP, 2021, 062

Palanque-Delabrouille, N., Yèche, C., Schöneberg, N., et al. 2020, JCAP, 2020, 038

Pandey, K. L., Karwal, T., & Das, S. 2020, JCAP, 2020, 026

Papastergis, E., Giovanelli, R., Haynes, M. P., & Shankar, F. 2015, A&A, 574, A113

Peter, A. H. G. 2010, PhRvD, 81, 083511

Peter, A. H. G., & Benson, A. J. 2010, PhRvD, 82, 123521

Peter, A. H. G., Moody, C. E., Benson, A. J., & Kamionkowski, M. 2010a, arXiv:1011.4970

Peter, A. H. G., Moody, C. E., & Kamionkowski, M. 2010b, PhRvD, 81, 103501

Planck Collaboration 2020, A&A, 641, A6

Polisensky, E., & Ricotti, M. 2011, PhRvD, 83, 043506

Poulin, V., Serpico, P. D., & Lesgourgues, J. 2016, JCAP, 2016, 036

Robles, V. H., Kelley, T., Bullock, J. S., & Kaplinghat, M. 2019, MNRAS, 490, 2117

Rudakovskiyi, A., Mesinger, A., Savchenko, D., & Gillet, N. 2021, MNRAS, 507, 3046

Safarzadeh, M., & Spergel, D. N. 2020, ApJ, 893, 21

Salucci, P. 2019, A&ARv, 27, 2

Sánchez-Salcedo, F. J. 2003, ApJL, 591, L107

Schultz, C., Óñorbe, J., Abazajian, K. N., & Bullock, J. S. 2014, MNRAS, 442, 1597

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

Simon, J. D. 2019, ARA&A, 57, 375

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685

Strigari, L. E., Bullock, J. S., Kaplinghat, M., et al. 2008, Natur, 454, 1096

Stücker, J., Angulo, R. E., Hahn, O., & White, S. D. M. 2022, MNRAS, 509, 1703

Tulin, S., & Yu, H.-B. 2018, PhR, 730, 1

Vattis, K., Koushiappas, S. M., & Loeb, A. 2019, PhRvD, 99, 121302

Verde, L., Treu, T., & Riess, A. G. 2019, NatAs, 3, 891

Viel, M., Becker, G. D., Bolton, J. S., & Haehnelt, M. G. 2013, PhRvD, 88, 043502

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vogelsberger, M., Zavala, J., Cyr-Racine, F.-Y., et al. 2016, MNRAS, 460, 1399

Vogelsberger, M., Zavala, J., & Loeb, A. 2012, MNRAS, 423, 3740

Wang, M.-Y., Croft, R. A. C., Peter, A. H. G., Zentner, A. R., & Purcell, C. W. 2013, PhRvD, 88, 123515

Wang, M.-Y., Peter, A. H. G., Strigari, L. E., et al. 2014, MNRAS, 445, 614

Wang, M.-Y., & Zentner, A. R. 2012, PhRvD, 85, 043514

Wang, Y., Nadler, E. O., Mao, Y.-Y., et al. 2021, ApJ, 915, 116

Waskom, M. L. 2021, JOSS, 6, 3021

Zavala, J., Vogelsberger, M., & Walker, M. G. 2013, MNRAS, 431, L20