measurements of the Milky Way satellite population can be reinterpreted to place stringent bounds on the microphysics of DM.

In this study, we focus on a scenario in which DM elastically scatters with normal matter (baryons), altering matter perturbations in the early Universe and consequently reducing the present-day population of small galaxies. We further rely on the concordance between the CDM predictions and measurements of the Milky Way satellite abundance from the Dark Energy Survey (DES) and Pan-STARRS1 (PS1; Drlica-Wagner et al. 2020;Nadler et al. 2020a,b) to place the most stringent astrophysical bounds on a variety of velocity-dependent interactions between DM particles and protons.

In a previous pilot study (Nadler et al. 2019a), we developed a method to constrain velocity-independent scattering only, which we later applied to DES data (Nadler et al. 2020b). Here, we generalize our analysis to include a whole class of velocity-dependent interaction models; this generalization has required a new approach to quantifying the impact of DM interactions on the satellite population. This method relies on predicting the most robust features of the matter transfer function in IDM cosmology and relating those features to the present-day abundance of low-mass halos. It does not necessitate precise modeling of the intricacies of galaxy formation within IDM. As such, it is only suited for placement of conservative upper bounds on the interaction cross section. Even so, the upper bounds we obtain are 3-5 orders of magnitude more stringent than the previous observational limits. They are also the first near-field cosmological limit on velocity-dependent DMbaryon interactions.

We address the same low-energy physics-and the same DM parameter space-as direct detection experiments. However, as with most other observational approaches, it is particularly well-suited for probing relatively large interaction cross sections and sub-GeV particle masses, outside the target sensitivity of most nuclear-recoil-based underground experiments (e. g. Agnese et al. 2019;The XENON collaboration et al. 2020;Emken & Kouvaris 2018). It is thus directly complementary to laboratory searches for DM interactions with Standard Model particles and substantially reduces the allowed parameter space for IDM models.

The paper is organized as follows. In Section 2, we review the theoretical models of DM-proton scattering and their effects on observations. In Section 3, we review the observational constraints on the Milky Way satellite galaxy population. In Section 4, we describe our approach to inferring upper limits on the interaction cross section from the measured abundance of the Milky Way satellites. In Section 5, we present our results. We discuss our findings and conclude in Section 6. Throughout, we adopt the following cosmological parameters: DM density Ω dm h 2 = 0.1153, baryon density Ω b h 2 = 0.02223, radiation density Ω rad ≈ 9.23 × 10 -5 , the Hubble constant h = 0.6932, optical depth to reionization τ reio = 0.081, the amplitude of the scalar perturbations A s = 2.464× 10 -9 , and the scalar spectral index n s = 0.9608; we set c = k B = 1.foot_0 2. THEORY We consider elastic scattering between DM particles and protons that predominantly takes place in the early Universe. 2 We consider any scattering process with a momentum-transfer cross section of the form σ MT = σ 0 v n , where v is the relative particle velocity and σ 0 is a free parameter of the model; we focus on powerlaw index values n ∈ {0, 2, 4, 6} and consider a range of DM particle masses m χ ∈ [15 keV, 100 GeV]. 3 We choose this empirical parameterization and values of n because they are representative of wide variety of relativistic DM models which can be described by a lowenergy effective field theory of DM scattering with nucleons, broadly considered in DM searches (Anand et al. 2014;Fitzpatrick & Zurek 2010). The models are represented here by an appropriate choice of n (Boddy & Gluscevic 2018). For example, n = 0 represents a cross section with no velocity dependence and corresponds to a spin-independent or spin-dependent contact interaction, well-studied in context of direct detection; n = 2 arises at leading order from DM with an electric dipole moment, induced by a heavy mediator that kinetically mixes with the photon (Fitzpatrick et al. 2013). Aside from its connection to particle theory, the power-law parameterization is sufficient to fully capture the effects of scattering on structure formation and thermal history of the Universe, and is thus adopted as a standard approach in observational searches for DM interactions (e. g. Dvorkin et al. 2014;Gluscevic & Boddy 2018;Boddy & Gluscevic 2018;Xu et al. 2018;Slatyer & Wu 2018;Boddy et al. 2018).

In an IDM cosmology, DM-baryon scattering leads to heat and momentum transfer between the cosmological fluids, smoothing out small-scale density perturbations through collisional damping. The momentum-transfer rate R χ and the heat-transfer rate R χ are proportional to σ 0 , and their redshift evolution is largely dictated by the evolution of the relative particle velocities (Dvorkin et al. 2014;Gluscevic & Boddy 2018). Since particle velocities are primarily sourced by thermal motions in the early Universe (z 10 4 ), the associated R χ evolves of Figure 2), that the data is likely even more constraining. However, non-linear evolution of structure at such small scales makes it difficult to improve the limit on the basis of linear-theory considerations, without running IDM cosmological simulations. Nonetheless, our numerical limit presents a tremendous improvement over the other observational bounds on IDM, as quantified in Section 5.

As noted before, the analytic-limit prescription of Nadler et al. (2019a) does not strictly apply in generic IDM models. We consider it here only as a rough estimate for how much our numerical limits could be improved in principle. We start by noting that IDM scattering affects matter perturbations until DM and baryons decouple, at z dec . Following Nadler et al. (2019a), we find z dec by setting

where (Boddy et al. 2018) 3) into the Eq. ( 2), we get z dec (σ 0 |m χ , n).

We can further estimate a critical comoving scale below which collisional damping suppresses the matter transfer function; this scale corresponds to the size of the particle horizon at z dec , given by

Substituting z dec (σ 0 |m χ , n) in Eq. ( 4), we obtain k crit (σ 0 |m χ , n). Finally, we use Eq. ( 1), to relate σ 0 to the mean mass of the smallest halos affected by IDM physics, for a given n and m χ . We note that a particular value of k crit translates to a different amount of suppression in T 2 (k), depending on the interaction model; however, this analytic prescription does not predict the amount of suppression. We also note that M min = 3.2 × 10 8 M roughly corresponds to z dec ≈ 4 × 10 7 (the intersection point in Figure 1). The benefit of the analytic calculation for n > 0 is that it provides a rough estimate of the largest mass at which halo abundances are affected by IDM, for a given σ 0 . In other words, the values of σ 0 that satisfy M crit (σ 0 |m χ , n) < M min only affect halos of masses below the current detection threshold. As such, they are largely allowed by the current data. For illustration, the right panel of Figure 2 shows transfer functions for all our IDM models, where σ 0 is set using the analytic estimate. The corresponding T 2 (k) curves present outer envelopes of the "disallowed" (shaded) region for most values of k. We thus expect that the analytically-estimated bounds roughly capture maximal improvement that can be obtained with detailed forward modeling of the same data; however, this is a rough estimate that only holds true for some DM masses, as we show in the following.

Our numerical bounds on σ 0 as a function of m χ are presented in Figure 3 and Table 1 for n ∈ {0, 2, 4, 6}. In the same Figure, we present the results of our analytic estimates, cast as an equivalent limit on σ 0 . We also show the previous limits from Planck measurements of the CMB temperature and polarization anisotropy (Boddy & Gluscevic 2018), the limits from FIRAS spectral-distortion bounds (Ali-Haïmoud et al. 2015), and the limits from Lyman-α forest analysis (Xu et al. 2018). Our numerical limits are orders of magnitude more constraining than those in previous studies and currently present the most stringent astrophysical bounds on these IDM models. Comparing to the Planck limits, we report an improvement of approximately 3 and 5 orders of magnitude for n = 2 and n = 4, respectively. For n = 0, our findings are consistent with Nadler et al. (2019aNadler et al. ( , 2020b)).

Our numerical limits are the most conservative upper bounds on the momentum-transfer cross section from linear perturbation theory, in the sense that larger values of σ 0 lead to dramatic decrements in power on scales that are measured to be consistent with CDM. The analytic estimates, on the other hand, roughly identify values of σ 0 below which current data has a limited constraining power. However, the analytic prescription is a poor predictor of the bound at low DM masses for n > 0 models, as the the analytically-estimated cross sections fall into the excluded regions of the parameter space.

We note that the mass dependence of the numerical bound shown in Figure 3 differs from the dependence of the analytic estimate. While the analytic estimate directly inherits its mass dependence from the DM-baryon momentum transfer rate R χ , the numerical bound is additionally modulated by the requirement that the DAO features fall strictly beneath the WDM transfer function. The size of DAO features as a function of m χ is not straightforwardly modeled, but it affects the mass dependence of the numerical limit. 1 GeV 5.5 × 10 -16 1.9 × 10 -18

10 GeV 1.4 × 10 -14 3.9 × 10 -17 100 GeV 1.4 × 10 -13 4.4 × 10 -16

15 keV 1.9 × 10 -29 6.6 × 10 -28 100 keV 7.0 × 10 -24 2.1 × 10 -25 6 10 MeV 1.0 × 10 -12 3.0 × 10 -22

1 GeV 6.2 × 10 -10 7.8 × 10 -14 10 GeV 7.8 × 10 -9 2.5 × 10 -12 100 GeV 1.0 × 10 -7 3.2 × 10 -11

Table 1. Bounds on the normalization σ0 of the momentumtransfer cross section, σMT = σ0v n , obtained via the analytic and numerical approaches, for a set of DM masses mχ and power-law dependencies on particle velocity v, with an index n. Table entries correspond to the limits shown in Figure 3.

the behavior of the dark acoustic oscillations we observe in Figure 2 may be possible to model semi-analytically to understand their effects on dark matter substruc-ture in galaxies like the Milky Way. Indeed, such approaches will be necessary to move beyond limits and toward a discovery of new dark matter physics, should inconsistencies with the cold dark matter paradigm arise in future measurements. Simulations that include dark matter-baryon scattering could also uncover other potentially observable signatures of the interactions, such as impacts on halo density profiles. A combined analysis of all available observational probes is perhaps the most robust way to search for new physics of dark matter with upcoming surveys. Finally, we note that the validity of our results does not explicitly require thermal production of dark matter. However, we do assume that dark matter follows a Maxwell-Boltzmann distribution, achieved by the sufficiently strong coupling with baryons, through any one of the interactions we considered. Deviations from this assumption may occur because thermal decoupling takes place before the interactions themselves decouple (Ali-Haïmoud 2019). Furthermore, thermally-produced dark matter can alter primordial element abundances (Boehm et al. 2013;Nollett & Steigman 2015;Krnjaic & McDermott 2020) and there are corresponding limits on thermal dark matter candidate mass, complementary to our results. However, these limits rely on details of a specific model of dark matter, such as its spin statistics and the high-energy behavior of its interactions; in contrast, we need not make any assumptions about these details. We thus leave detailed comparisons of these bounds and considerations related to the velocity distributions for the future.