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

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

   The lifetime of the decay;

2. 2.

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

3. 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, b; 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}(\rm 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 mater (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).

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^{\prime}$ and a massless dark radiation species $\gamma^{\prime}$,

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

These decays impart a recoil velocity

assuming $V_{\rm 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\mathcal{O}(\mathrm{\,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^{\prime}$, or $\gamma^{\prime}$) 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.

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/\mathrm{\,Gyr},V_{\rm kick}/\mathrm{\,km}\mathrm{\,s}^{-1})\in\{(10,20),\allowbreak(20,20),\allowbreak(20,40),\allowbreak(40,20),\allowbreak(40,40),\allowbreak(80,40)\}$ and a corresponding CDM simulation, all with 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\mathrm{\,Mpc}^{-1}$ per side and zoom in on a halo mass of $M\approx 10^{12}\mathrm{\,M_{\odot}}$ (Wang et al., 2014). The highest-resolution region is simulated with a Plummer-equivalent force softening of $143\mathrm{\,pc}$ and a particle mass of $1.92\times 10^{5}\mathrm{\,M_{\odot}}$.³³3We define virial quantities according to the Bryan & Norman (1998) virial definition, with overdensities $\Delta_{\rm vir}$ set according to the cosmological parameters in our two simulation suites. 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 the appendix and are represented with a red–purple color scheme.

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_{\rm kick}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ ($40\mathrm{\,km}\mathrm{\,s}^{-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}\mathrm{\,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_{\rm max}$, divided by its peak value, $V_{\rm peak}$, again for DDM models with $V_{\rm kick}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ ($40\mathrm{\,km}\mathrm{\,s}^{-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_{\mathrm{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).⁴⁴4This 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. 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 Lyman-$\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; Rudakovskyi et al., 2021), offer relatively less constraining power for these models.

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

where $f_{\rm DDM}(M,\tau,V_{\rm 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_{\mathrm{DDM}}$ describes the suppression of subhalo abundances above a present-day mass threshold of $\sim 2\times 10^{7}\mathrm{\,M_{\odot}}$. To measure this quantity, we fit the 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):

where

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_{\rm peak}$ with varying $\tau$ and $V_{\rm kick}$. We perform this fit over a mass range of $3.8\times 10^{7}$–$9.4\times 10^{9}\mathrm{\,M_{\odot}}$ and derive best-fit values of

with a reduced $\chi^{2}$ goodness of fit of 0.3 (which reflects the relatively large Poisson uncertainties on the SHMFs).⁵⁵5The 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). 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}\mathrm{\,M_{\odot}}$, we estimate that this only influences our limits at the $\sim 5\mathrm{\,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_{\rm kick})$ sampled in our expanded suite of simulations based on Wang et al. (2014), which is bounded below (above) by $\tau/\mathrm{\,Gyr},V_{\rm kick}/\mathrm{\,km}\mathrm{\,s}^{-1}=10,20$ ($80,40$), with $z=0$.

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_{\mathrm{kick}}$ according to Equation (8), and we perform the analysis for $V_{\mathrm{kick}}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ and $40\mathrm{\,km}\mathrm{\,s}^{-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}\mathrm{\,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 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_{\rm peak}$ ($\sigma_{M}$), the mass at which 50% of halos host galaxies ($\log(\mathcal{M}_{50}/\mathrm{\,M_{\odot}})$), the strength of subhalo disruption due to baryons ($\mathcal{B}$), the scatter in the galaxy occupation fraction ($\sigma_{\mathrm{gal}}$), the amplitude of the galaxy–halo size relation ($\mathcal{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/\mathrm{\,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).⁶⁶6Although our DDM simulations sample down to $\tau=10\mathrm{\,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. We perform two separate fits at fixed $V_{\mathrm{kick}}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ and $40\mathrm{\,km}\mathrm{\,s}^{-1}$, each of which uses $36$ walkers, discards a burn-in period of $\sim 20$ autocorrelation lengths, and retains $\sim 10^{5}$ samples corresponding to $\sim 100$ autocorrelation lengths.

The posterior distributions from our $V_{\rm kick}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ and $V_{\rm kick}=40\mathrm{\,km}\mathrm{\,s}^{-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_{\mathrm{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/\mathrm{\,Gyr})>1.46$ ($\tau>29\mathrm{\,Gyr}$) for $V_{\mathrm{kick}}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ and $\log(\tau/\mathrm{\,Gyr})>1.63$ ($\tau>43\mathrm{\,Gyr}$) for $V_{\mathrm{kick}}=40\mathrm{\,km}\mathrm{\,s}^{-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_{\mathrm{kick}}$, we compute the value of $\tau$ that decreases $f_{\mathrm{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}\mathrm{\,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_{\mathrm{\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\mathrm{\,Gyr}$ for $V_{\mathrm{kick}}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ and $\tau>29\mathrm{\,Gyr}$ for $V_{\mathrm{kick}}=40\mathrm{\,km}\mathrm{\,s}^{-1}$. The SHMF suppression 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(\tau/\mathrm{\,Gyr})$ for $V_{\rm kick}=20\mathrm{\,km}\mathrm{\,s}^{-1}$ and $40\mathrm{\,km}\mathrm{\,s}^{-1}$ alongside the preferred region of DDM parameter space that potentially alleviates the $S_{8}$ tension (Abellán et al., 2020). 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., 2020) 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_{\mathrm{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.

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.

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; Newton et al. 2021; Dekker et al. 2021).⁹⁹9Nadler 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.

Nadler et al. (2021c) found that thermal relic WDM masses below $6.5\ \mathrm{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 $\sim 25\%$ relative to CDM at the minimum observed halo mass scale of $3.2\times 10^{8}\mathrm{\,M_{\odot}}$; the SHMF declines rapidly at smaller masses, reaching $\sim 75\%$ suppression relative to CDM at its half-mode mass of $3.8\times 10^{7}\mathrm{\,M_{\odot}}$. For comparison, the DDM models we rule out at $95\%$ confidence suppress the subhalo mass function by $\sim 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 $\sim 2\times 10^{7}\mathrm{\,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; Hsueh et al. 2020; Gilman et al. 2020a, b), 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). 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).