Bounds on Velocity-dependent Dark Matter–Proton Scattering from Milky Way Satellite Abundance

After decades of versatile experimental searches, observations of the universe remain the sole source of evidence for dark matter (DM). Identifying its nature amounts to understanding a major constituent of matter and thus inspires investigations across different fields of physical science. As the laboratory bounds on the most popular theoretical candidate models grow stronger, cosmological and astrophysical observations have emerged as an alternative and complementary probe of DM microphysics (for reviews, see Buckley & Peter 2018; Drlica-Wagner et al. 2019; Gluscevic et al. 2019; Grin et al. 2019).

The standard model of cosmology assumes cold dark matter (CDM) whose nongravitational interactions are observationally insignificant. Nearly all deviations from pure CDM considered in the current literature affect the way matter is distributed in the universe, including warm DM (WDM; Schneider 2016; Abazajian 2017; Adhikari et al. 2017), fuzzy DM (FDM; Hu et al. 2000; Hui et al. 2017), self-interacting DM (SIDM; Tulin & Yu 2018), DM interacting with dark radiation (Cyr-Racine et al. 2016), and with Standard Model particles (IDM; Dvorkin et al. 2014; Boddy & Gluscevic 2018; Escudero et al. 2018; Gluscevic & Boddy 2018; Nadler et al. 2019a, 2020a). One of the original incentives to consider beyond-CDM models was the perceived "missing satellites problem"—the apparent mismatch between the observed satellite galaxies orbiting the Milky Way and their predicted population from CDM simulations of structure formation (Klypin et al. 1999; Moore et al. 1999). Different properties of DM were invoked to account for the apparent mismatch, including appreciable free streaming (in WDM models), a macroscopic de Broglie wavelength (in FDM models), and particle interactions (in SIDM and IDM models). Virtually all of them suppress the formation of low-mass DM halos and reduce the abundance of galaxies that would inhabit them. However, a more recent census of faint galaxies in our galactic neighborhood, combined with advanced modeling of the galaxy–halo connection, has shown consistency between the CDM predictions and the observed satellite abundance, down to a halo mass of ∼10⁸ M_⊙ (Jethwa et al. 2018; Kim et al. 2018; Newton et al. 2018; Nadler et al. 2019b, 2020b). With this new development, 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; Nadler et al. 2020a, 2020b; Drlica-Wagner et al. 2020) 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. 2020a). 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 DM–baryon 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., Emken & Kouvaris 2018; Agnese et al. 2019; The XENON collaboration et al. 2020). 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² = 0.1153, baryon density Ω_b h² = 0.02223, radiation density Ω_rad ≈ 9.23 × 10⁻⁵, the Hubble constant h = 0.6932, optical depth to reionization τ_reio = 0.081, the amplitude of the scalar perturbations A_s = 2.464 × 10⁻⁹, and the scalar spectral index n_s = 0.9608; we set c = k_B = 1. ⁵

We consider elastic scattering between DM particles and protons that predominantly takes place in the early universe. ⁶ We consider any scattering process with a momentum-transfer cross section of the form σ_MT = σ₀ vⁿ , where v is the relative particle velocity and σ₀ is a free parameter of the model; we focus on power-law index values n ∈ {0, 2, 4, 6} and consider a range of DM particle masses ${m}_{\chi }\in \left[15\,\mathrm{keV},100\,\mathrm{GeV}\right]$. ⁷ 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 low-energy effective field theory of DM scattering with nucleons, broadly considered in DM searches (Fitzpatrick & Zurek 2010; Anand et al. 2014). 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; Boddy & Gluscevic 2018; Boddy et al. 2018; Gluscevic & Boddy 2018; Slatyer & Wu 2018; Xu 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}_{\chi }^{{\prime} }$ are proportional to σ₀, 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⁴), the associated R_χ evolves monotonically with redshift z, as the universe cools (Figure 1). For models with n ≥ 0, DM decouples from protons well before cosmic recombination and deep into this regime (Dvorkin et al. 2014; Boddy & Gluscevic 2018). In such cases, the interactions affect structure today primarily by means of suppressing the linear matter power spectrum P(k) at small scales (large wavenumbers k), early on in cosmic history; the square of the transfer function T²(k) ≡ P(k)/P_CDM(k) (the ratio between the IDM and the CDM power spectrum) features a cutoff, shown in Figure 2 in colored lines. Models with n < 0, on the other hand, feature scattering at late times, after structure formation commences. Their effects are more challenging to compute and we leave their consideration for a future study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For n ≥ 0, DM scattering affects physical scales that enter the particle horizon prior to DM–baryon decoupling. For the n = 0 case of a velocity-independent interaction, the resulting cutoff in T²(k) is the main signature of IDM physics. However, for the velocity-dependent interactions with n > 0, there are also prominent "dark acoustic oscillations" (DAO; see Cyr-Racine et al. 2016) that appear at scales below the cutoff, due to the tight coupling between the photon–baryon fluid and the DM fluid at early times (Figure 2). In both cases, the IDM-induced suppression of small-scale density perturbations ultimately leads to a decrement in the abundance of low-mass halos that host dwarf galaxies, as compared to the CDM cosmology.

To forward-model a population of galaxies in a beyond-CDM cosmology and confront it with observations in principle requires a suite of fully consistent cosmological simulations, including beyond-CDM physics. Such simulations are computationally expensive and were only performed for certain sets of beyond-CDM scenarios that feature a suppression of P(k), most notably WDM, SIDM, and ETHOS models (Schneider et al. 2012; Angulo et al. 2013; Lovell et al. 2014; Cyr-Racine et al. 2016; Bose et al. 2017; Murgia et al. 2017). In Nadler et al. (2019a), we utilized the fact that the velocity-independent DM–proton scattering (with n = 0) suppresses P(k) in a way that resembles WDM, and used the mapping between the parameters of the two models to derive bounds on this specific IDM case. However, the mapping between WDM and IDM breaks down for velocity-dependent (n > 0) scattering we focus on here because of the large DAO features. In fact, the IDM transfer function does not straightforwardly map onto any other previously explored scenario with a power cutoff. For this reason, we develop and apply a new method that relates the bounds on the matter transfer function, inferred from DES and PS1 measurements, to the most robust features of T²(k) within IDM with velocity-dependent scattering.

We use the recent measurements of the Milky Way satellite population by DES and PS1 (Drlica-Wagner et al. 2020), and their inferred bounds on the matter transfer function T(k). The bounds are based on a probabilistic inference that combines (i) models for satellite detectability in the relevant survey footprints (Drlica-Wagner et al. 2020), (ii) high-resolution DM-only simulations (Mao et al. 2015) chosen to match the observed characteristics of the Milky Way system, and (iii) an empirical model of the galaxy–halo connection (Nadler et al. 2019b, 2020b). By performing mock observations of the satellite populations and statistically comparing them to the luminosity, size, and radial distributions of the observed satellites, the model—including suppression of the subhalo mass function—is fit to the data using a Markov Chain Monte Carlo approach (Nadler et al. 2020b, 2020a).

The results are cast in terms of the bounds on models that suppress T²(k), notably a lower limit on the thermal-relic WDM mass of 6.5 keV at 95% confidence (Nadler et al. 2020a). The corresponding WDM T²(k) is shown as the solid black line in Figure 2. The same results were also cast as an upper limit on the minimum mass of halos that host observed satellite galaxies, ${M}_{\min }\lt 3.2\times {10}^{8}$ M_⊙. The corresponding comoving wavenumber k_crit = 33.2 h Mpc⁻¹, given by

$$\begin{eqnarray}&&{M}_{\min }=\displaystyle \frac{4\pi }{3}{{\rm{\Omega }}}_{\mathrm{dm}}\bar{\rho }{\left(\displaystyle \frac{\pi }{{k}_{\mathrm{crit}}}\right)}^{3},\end{eqnarray} \tag{ 1 }$$

represents the largest k effectively probed by these data, shown as the vertical dotted line in Figure 2. $\bar{\rho }=4.75\times {10}^{-6}$ is the mean density of the universe today and λ_crit = 2π/k_crit is the comoving size of perturbations giving rise to halos of average mass M_min.

The WDM transfer function corresponding to the WDM mass limit, together with the minimum-halo-mass limit, in Figure 2, delineates the allowed region for T²(k): a viable beyond-CDM model must not suppress T²(k) more than the WDM model in this figure, unless the suppression occurs beyond k_crit, where the data has no constraining power. The implications to IDM are as follows:

• 1.  
  Since T²(k) at the thermal-relic-WDM mass limit delineates the maximum suppression tolerated by current data, IDM models that produce a more suppressed T²(k) are inconsistent with the data.
• 2.  
  M_min is the minimum mass of halos whose abundance is demonstrably consistent with CDM; IDM models that exclusively alter the abundance of lower-mass halos are consistent with the data.

We incorporate these points within a numerical approach and an analytic estimate, described below, to translate the M_min and WDM mass bounds described in the previous section into a bound on IDM. Before proceeding, we highlight an important distinction between them: the numerical approach yields an upper bound on σ₀, while the analytic estimate does not provide upper bounds by itself; rather, it roughly estimates a potential maximal improvement over the numerical limit, if full forward modeling of the satellite population is applied to the same data.

4.1. Numerical Limits

4.2. Analytic Estimates

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

$$\begin{eqnarray}&&{\left.{aH}={R}_{\chi }\right|}_{{z}_{\mathrm{dec}}},\end{eqnarray} \tag{ 2 }$$

where (Boddy et al. 2018)

$$\begin{eqnarray}&&{R}_{\chi }=\displaystyle \frac{{N}_{n}a{\rho }_{b}{Y}_{p}{\sigma }_{0}}{{m}_{\chi }+{m}_{p}}{\left(\displaystyle \frac{{T}_{b}}{{m}_{p}}+\displaystyle \frac{{T}_{\chi }}{{m}_{\chi }}\right)}^{\tfrac{1+n}{2}},\end{eqnarray} \tag{ 3 }$$

${N}_{n}={2}^{5+n/2}{\rm{\Gamma }}(3+n/2)/(3\sqrt{\pi })$, a is the scale factor, Y_p = 0.75 is the proton mass fraction, ρ_b is the baryon energy density, and m_p is the proton mass. During radiation domination, the temperature of baryons evolves as T_b = T₀(1 + z), where T₀ = 2.73 K. The temperature of the DM fluid T_χ is strongly coupled to T_b until thermal decoupling at z_th, and afterwards evolves adiabatically, T_χ = T₀(1 + z)²/(1 + z_th). Thermal decoupling occurs when the heat-transfer rate, ${R}_{\chi }^{{\prime} }\equiv {R}_{\chi }{m}_{\chi }/({m}_{\chi }+{m}_{p})$, matches the Hubble rate, ${aH}={R}_{\chi }^{{\prime} }{| }_{{z}_{\mathrm{th}}}$. Substituting Equation (3) into the Equation (2), we get z_dec(σ₀∣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

$$\begin{eqnarray}&&{k}_{\mathrm{crit}}\equiv 2{aH}{| }_{{z}_{\mathrm{dec}}}\approx 2{H}_{0}{z}_{\mathrm{dec}}\sqrt{{{\rm{\Omega }}}_{\mathrm{rad}}}.\end{eqnarray} \tag{ 4 }$$

Substituting z_dec(σ₀∣m_χ , n) in Equation (4), we obtain k_crit(σ₀∣m_χ , n). Finally, we use Equation (1), to relate σ₀ 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²(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\times {10}^{8}$ M_⊙ roughly corresponds to z_dec ≈ 4 × 10⁷ (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 σ₀. In other words, the values of σ₀ that satisfy ${M}_{\mathrm{crit}}({\sigma }_{0}| {m}_{\chi },n)\lt {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 σ₀ is set using the analytic estimate. The corresponding T²(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 σ₀ 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 σ₀. 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 Lyα 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. (2019a, 2020a).

Zoom In Zoom Out Reset image size

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 σ₀ 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 σ₀ 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.

We use the latest measurements of the Milky Way satellite population from DES and Pan-STARRS1 to infer the most stringent astrophysical bound to date on velocity-dependent interactions between dark matter particles and protons. We generalize methods we previously developed for velocity-independent scattering and apply them to any velocity-dependent interaction that dominates over Hubble expansion in the early universe. We do not assume any specific high-energy behavior of dark matter, and thus probe the parameter space for sub-GeV dark matter in a generic way, complementary to laboratory experiments, providing an important guide for identifying viable candidate models. Our results exclude interaction cross sections that can be reached with the next-generation cosmic microwave background experiments (Li et al. 2018). The methods we developed here are applicable to future data from facilities such as the Vera C. Rubin Observatory (LSST Science Collaboration et al. 2009) and from line-intensity tomography surveys (Muñoz et al. 2020).

We identify a few promising directions for follow-up studies. First, we note that our numerical results use linear perturbation theory and rely only on the most robust features of the matter transfer function that arise as a consequence of dark matter scattering. A simulation-based approach to consistently and fully forward-model satellite populations within an interacting dark matter cosmology can improve upon our results. Furthermore, the behavior of the dark acoustic oscillations we observe in Figure 2 may be possible to model semianalytically to understand their effects on dark matter substructure 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 CDM 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 (Bœhm 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.

K.M. acknowledges the support through the Provost's Research Fellowship for undergraduate students at USC. V.G. is supported by the National Science Foundation under grant No. PHY-2013951. E.O.N. is supported by the National Science Foundation (NSF) under grant No. NSF DGE-1656518 through the NSF Graduate Research Fellowship. The authors thank Yacine Ali-Haïmoud and Keith Bechtol for useful comments.