Remnants of Galactic subhalos and their impact on indirect dark matter searches

TL;DR

This paper reevaluates the abundance of dark matter subhalos in the Milky Way, accounting for simulation artifacts, and finds that the expected gamma-ray and cosmic-ray signals are roughly doubled, impacting indirect dark matter detection.

Contribution

It introduces a semi-analytical method to accurately estimate subhalo abundance, correcting for over-efficient tidal disruption seen in simulations.

Findings

Boost factor for gamma rays increased by ~2

Boost factor for cosmic-ray antiprotons increased by ~2

Revised subhalo abundance impacts dark matter search strategies

Abstract

Dark matter subhalos, predicted in large numbers in the cold dark matter scenario, should have an impact on particle dark matter searches. Recent results show that tidal disruption of these objects in computer simulations is over-efficient due to numerical artifacts and resolution effects. Accounting for these results, we re-estimate the subhalo abundance in the Milky Way using semi-analytical techniques. In particular, we show that the boost factor for gamma rays and cosmic-ray antiprotons is increased by roughly a factor of two.

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}m\mathrm{d}c}(\vec{r},m,c)\propto \frac{\mathrm{d}{\mathcal{P}}_v}{\mathrm{d}V}(\vec{r})\times \frac{\mathrm{d}{\mathcal{P}}_m}{\mathrm{d}m}(m)\times \frac{\mathrm{d}{\mathcal{P}}_c}{\mathrm{d}c}(c,m),$$

$${\rho }_{\mathrm{sm}}(\vec{r})={\rho }_{\mathrm{DM}}(\vec{r})-⟨{\rho }_{\mathrm{cl}}⟩(\vec{r}),$$

$$\frac{\mathrm{d}{\mathcal{P}}_m}{\mathrm{d}m}(m)\propto m^{-{\alpha }_{\mathrm{m}}}\Theta (m-m_{\min })\Theta (m_{\max }-m),$$

$$m_{\mathrm{t}}(\vec{r},m,c)=4\pi {\int }_0^{r_{\mathrm{t}}(\vec{r},m,c)}\mathrm{d}x{x}^2{\rho }_{\mathrm{sub}}(x)\le m,$$

$$\frac{r_{\mathrm{t}}(\vec{r},m,c)}{r_{\mathrm{s}}(m,c)}<{ϵ}_{\mathrm{t}}\iff \mathrm{subhalo}\mathrm{is}\mathrm{disrupted}$$

$$\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}m\mathrm{d}c}(\vec{r},m,c)=\frac{N_{\mathrm{tot}}}{K_{\mathrm{tot}}}\frac{\mathrm{d}{\mathcal{P}}_v}{\mathrm{d}V}(\vec{r})\times \frac{\mathrm{d}{\mathcal{P}}_m}{\mathrm{d}m}(m)\times \frac{\mathrm{d}{\mathcal{P}}_c}{\mathrm{d}c}(c,m)\Theta \left(\frac{r_{\mathrm{t}}(\vec{r},m,c)}{r_{\mathrm{s}}(m,c)}-{ϵ}_{\mathrm{t}}\right),$$

$$K_{\mathrm{tot}}=\int \mathrm{d}V\frac{\mathrm{d}{\mathcal{P}}_v}{\mathrm{d}V}(\vec{r})\int \mathrm{d}m\frac{\mathrm{d}{\mathcal{P}}_m}{\mathrm{d}m}(m)\int \mathrm{d}c\frac{\mathrm{d}{\mathcal{P}}_c}{\mathrm{d}c}(c,m)\Theta \left(\frac{r_{\mathrm{t}}(\vec{r},m,c)}{r_{\mathrm{s}}(m,c)}-{ϵ}_{\mathrm{t}}\right).$$

$$f_{\mathrm{sub}}(m_1,m_2)=\frac{1}{M_{\mathrm{DM}}}\int \mathrm{d}V{\int }_{m_1}^{m_2}\mathrm{d}m\int \mathrm{d}c\times m\times {\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}m\mathrm{d}c}|}_{\mathrm{DMO},{ϵ}_{\mathrm{t}}=1}.$$

$$⟨{\rho }_{\mathrm{cl}}⟩(\vec{r})={\int }_{m_{\min }}^{m_{\max }}\mathrm{d}m{\int }_1^{\infty }\mathrm{d}c\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}m\mathrm{d}c}{m}_{\mathrm{t}}(\vec{r},m,c),$$

$$\frac{N_{\mathrm{DMO},{ϵ}_{\mathrm{t}}=1}}{K_{\mathrm{DMO},{ϵ}_{\mathrm{t}}=1}}=\frac{N_{\mathrm{tot}}}{K_{\mathrm{tot}}}.$$

$$\frac{\mathrm{d}N}{\mathrm{d}V}(\vec{r})$$

$$=$$

$${\mathcal{L}}_0(\vec{r})={\rho }_{\mathrm{DM}}^2(\vec{r}).$$

$$L_{\mathrm{t}}(\vec{r},m,c)={\int }_{V_{\mathrm{t}}}{\mathrm{d}}^3\vec{r}{\rho }_{\mathrm{sub}}^2(\vec{r}),$$

$${\mathcal{L}}_{\mathrm{cl}}(\vec{r})$$

$${\mathcal{L}}_{\mathrm{sm}}(\vec{r})={\rho }_{\mathrm{sm}}^2(\vec{r}),$$

$${\mathcal{L}}_{\mathrm{cs}}(\vec{r})$$

$$=$$

$$\mathcal{L}(\vec{r})={\mathcal{L}}_{\mathrm{cl}}(\vec{r})+{\mathcal{L}}_{\mathrm{sm}}(\vec{r})+{\mathcal{L}}_{\mathrm{cs}}(\vec{r}).$$

$$1+\mathcal{B}(\vec{r})=\frac{\mathcal{L}(\vec{r})}{{\mathcal{L}}_0(\vec{r})}.$$

$$\frac{\mathrm{d}{\Phi }_{\gamma }}{\mathrm{d}E\mathrm{d}\Omega }=\frac{1}{4\pi }\frac{⟨\sigma v⟩}{2m_{\chi }^2}\frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}E}\int \mathrm{d}s{\rho }^2,$$

$$J(\psi )=\int \mathrm{d}s{\rho }^2(s,\psi ),$$

$$1+{\mathcal{B}}_{\gamma }(\psi )=\frac{J(\psi )}{J_0(\psi )},$$

$$1+{\mathcal{B}}_{\overline{p}}(T)=\frac{{\Phi }_{\overline{p}}(T)}{{\Phi }_{\overline{p},0}(T)},$$

$$-K\Delta \Psi +\vec{\nabla }.({\vec{V}}_{\mathrm{c}}\Psi )+{\partial }_E\left[b\Psi -K_{EE}{\partial }_E\Psi \right]+2h\delta (z){\Gamma }_{\mathrm{ann}}\Psi =Q_{\mathrm{DM}},$$

$$Q_{\mathrm{DM}}(E,\vec{r})=\frac{⟨\sigma v⟩}{2}\frac{\mathrm{d}{N}_{\overline{p}}}{\mathrm{d}E}{\left(\frac{\rho (\vec{r})}{m_{\chi }}\right)}^2,$$

$$\frac{\mathrm{d}{\Phi }_{\overline{p}}}{\mathrm{d}T\mathrm{d}\Omega }=\frac{v_{\overline{p}}}{4\pi }\frac{⟨\sigma v⟩}{2m_{\chi }^2}\int \mathrm{d}{E}_{\mathrm{s}}\int {\mathrm{d}}^3{\vec{r}}_{\mathrm{s}}G(E\leftarrow E_{\mathrm{s}};{\vec{r}}_{\odot }\leftarrow {\vec{r}}_{\mathrm{s}})\frac{\mathrm{d}{N}_{\overline{p}}}{\mathrm{d}E}(E_{\mathrm{s}}){\rho }^2({\vec{r}}_{\mathrm{s}}).$$

$$G(E\leftarrow E_{\mathrm{s}};{\vec{r}}_{\odot }\leftarrow {\vec{r}}_{\mathrm{s}})=\delta (E-E_{\mathrm{s}})\times \overline{G}({\vec{r}}_{\odot }\leftarrow {\vec{r}}_{\mathrm{s}}),$$

$$1+{\mathcal{B}}_{\overline{p}}(T)=\frac{\int {\mathrm{d}}^3{\vec{r}}_{\mathrm{s}}\overline{G}({\vec{r}}_{\odot }\leftarrow {\vec{r}}_{\mathrm{s}})\mathcal{L}({\vec{r}}_{\mathrm{s}})}{\int {\mathrm{d}}^3{\vec{r}}_{\mathrm{s}}\overline{G}({\vec{r}}_{\odot }\leftarrow {\vec{r}}_{\mathrm{s}}){\mathcal{L}}_0({\vec{r}}_{\mathrm{s}})}.$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date \continuouspagesyes \TitleRemnants of Galactic subhalos and their impact on indirect dark matter searches \AuthorMartin Stref *1,†**\orcidA, Thomas Lacroix 2,‡\orcidB and Julien Lavalle 1,§\orcidC \AuthorNamesMartin Stref, Thomas Lacroix and Julien Lavalle

\firstnotemartin.stref@umontpellier.fr \secondnotethomas.lacroix@uam.es \thirdnotelavalle@in2p3.fr

1 Introduction

There is overwhelming evidence that most of the matter in the Universe is non-baryonic Aghanim et al. (2018). An exciting possibility to account for these puzzling observations is that the Universe is filled with exotic particles which interacts only very weakly with ordinary matter Bertone et al. (2005); Feng (2010). One of the most elegant and popular particle dark matter (DM) candidate is the WIMP (Weakly Interacting Massive Particle). These hypothetical particles are being looked for in particle colliders Fairbairn et al. (2007); Goodman et al. (2010); Kahlhoefer (2017), in direct detection experiments Cerdeno and Green (2010); Baudis (2012); Freese et al. (2013) and in cosmic radiation Lavalle and Salati (2012); Gaskins (2016), so far without success. Although one of the motivations for WIMPs is related to the fact that they emerge naturally in particle theories addressing the hierarchy problem Jungman et al. (1996); Servant and Tait (2003); Agashe and Servant (2004), WIMPs are also attractive stemming from their very simple thermal production mechanism in the early universe. Moreover, a large fraction of the available parameter space is still unconstrained and currently actively explored Leane et al. (2018). For completeness, it is worth recalling that many alternatives to WIMPs exist which we will not discuss here, like axions Marsh (2016), sterile neutrinos Boyarsky et al. (2009), primordial black holes Carr et al. (2016) and extended dark sectors Alexander et al. (2016). The cosmological paradigm best supported by current probes is that DM is cold, i.e. collisionless and non-relativistic. This implies a structuring of matter on scales smaller than typical galaxies, with a model-dependent cutoff Blumenthal et al. (1984); Davis et al. (1985). Interestingly, subgalactic scales are those where there could be departures from the predictions of the cold DM paradigm, because of some observational issues Bullock and Boylan-Kolchin (2017). This might sign new specific properties of the dark matter (e.g. Tulin and Yu (2018)), or it could be due to baryonic effects (e.g. Pontzen and Governato (2012)). This motivates a detailed inspection of the impact of DM properties on the smallest scales, irrespective of the underlying scenario.

The small-scale structuring of DM, as treated for instance in the WIMP scenario, translates into a large population of subhalos within galactic halos Moore et al. (1999); Diemand et al. (2008); Springel et al. (2008). Modeling these subhalos is crucial if one is to make accurate predictions for direct and indirect DM searches. This is a difficult task as numerical simulations are far from resolving the smallest structures predicted by the cold DM paradigm. To incorporate the smallest structures, one can extrapolate the results of simulations over orders of magnitude in scales, see e.g. Springel et al. (2008), but this represents a leap of faith. On the other hand, one can employ semi-analytical models, see e.g. Berezinsky et al. (2003); Penarrubia and Benson (2005); Lavalle et al. (2008); Bartels and Ando (2015). The difficulty with the latter is to account for the tidal effects experienced by subhalos within the host galaxy. These models can be calibrated on cosmological simulations, which are supposed to consistently describe the tidal stripping of subhalos in their host halo. However, it was recently pointed out by van den Bosch and collaborators van den Bosch et al. (2018); van den Bosch and Ogiya (2018) that simulations are plagued with numerical artifacts that lead to a significant overestimate of the tidal stripping efficiency, and therefore to an underestimate of the actual subhalo population even within the numerical resolution limit. An alternative and complementary way to study the tidal stripping of subhalos is to rely on analytical or semi-analytical methods, which are based on first principles and allow to deal with subhalo mass scales down to the free-streaming scale. Here, we review the semi-analytical model developed by Stref and Lavalle (2017) (SL17 hereafter), which incorporates a realistic and kinematically constrained Milky Way mass model (including baryons) and predicts the Galactic subhalo abundance.111We refer to this model as semi-analytical because it involves integrals that must be computed numerically. The model does not rely on numerical simulations except at the level of a calibration described in Sec. 2.3 This model accounts for different sources of tidal effects, and can easily accommodate to different prescriptions for the tidal disruption efficiency.

This paper is structured as follows. In Sec. 2, we briefly review the SL17 model and discuss the resilience of subhalos to tidal effects in light of recent analyses of simulation results van den Bosch et al. (2018); van den Bosch and Ogiya (2018). In Sec. 3, we compute the DM mass density within subhalos as well as the number density of these objects in the Milky Way. Finally, in Sec. 4, we look at the impact of our results on indirect searches for annihilating DM focusing on gamma rays and cosmic-ray antiprotons.

2 A semi-analytical model of Galactic subhalos

In this section, we review the SL17 Galactic subhalo population model and discuss the tidal effects experienced by subhalos. We then propose a way of incorporating the recent results of van den Bosch and collaborators in the model in a consistent calibration procedure.

2.1 Review of the Stref & Lavalle model

SL17 is a semi-analytical model of Galactic subhalos which is built upon dynamical constraints and cosmological considerations. The main input of the model is the initial subhalo phase-space density

[TABLE]

where phase space refers to the position-mass-concentration space. The functions $\mathrm{d}\mathcal{P}_{v}/\mathrm{d}V$, $\mathrm{d}\mathcal{P}_{m}/\mathrm{d}m$ and $\mathrm{d}\mathcal{P}_{c}/\mathrm{d}c$ are the spatial, mass and concentration distributions, respectively. It is assumed that, should subhalos behave as hard spheres (as is the case for single DM "particles" in a cosmological simulation), they would be spatially distributed as $\mathrm{d}\mathcal{P}_{v}/\mathrm{d}V\propto\rho_{\rm DM}$ where $\rho_{\rm DM}$ is the total DM density profile of the Galaxy. This sets our initial conditions before tidal disruption. The smooth DM mass density is computed through

[TABLE]

where $\left<\rho_{\rm cl}\right>$ is the average DM mass density inside clumps (this quantity is explicitly computed in Sec. 3). In the following, we use the Galactic mass models constrained by McMillan (2017) on pre-Gaia data for the DM and baryonic mass distributions. In this framework, Eq. (2) ensures the compatibility of our subhalo model with the constrained DM profile $\rho_{\rm DM}$. Although this work is devoted to the study of Milky Way subhalos, SL17 can in principle be used to study the substructure population in any virialized DM system. One only needs a mass model for the system in question and a proper calibration of the subhalo mass fraction through the procedure outlined in Sec. 2.3.

The mass $m$ and concentration $c$ refer to the cosmological mass $m_{200}$ and concentration $c_{200}$ (defined with respect to the critical density) where we have dropped the 200 index for convenience. The subhalo mass function measured in simulations is consistent with a power-law Diemand et al. (2008); Springel et al. (2008)

[TABLE]

where $\Theta$ is the Heaviside step function, and the power-law index is $\alpha_{\rm m}=1.9$ or $\alpha_{\rm m}=2$. These values of $\alpha_{\rm m}$ encompass the Press and Schechter (1974) mass function and the Sheth and Tormen (1999) mass function, as illustrated in Fig. 1. These functions can be computed directly from the matter power spectrum in the framework of the excursion set theory Bond et al. (1991), for the spherical collapse (Press-Schechter) and the ellipsoidal collapse (Sheth-Tormen). Thus the two power-law indices we consider bracket the theoretical uncertainties on the small-scale mass function. If the DM is made of weakly interacting massive particles (WIMPs), the mass cutoff $m_{\rm min}$ can be related to the kinetic decoupling of the DM particle and is found to lie between $10^{-4}\,\rm M_{\odot}$ and $10^{-10}\,\rm M_{\odot}$ Hofmann et al. (2001); Berezinsky et al. (2003); Loeb and Zaldarriaga (2005); Green et al. (2005); Bertschinger (2006); Bringmann and Hofmann (2007); Bringmann (2009). The maximal mass $m_{\rm max}$ is set to $0.01\times M_{\rm DM}$ where $M_{\rm DM}$ is the total DM mass in the Milky Way. The concentration distribution $\mathrm{d}\mathcal{P}_{c}/\mathrm{d}c$ is generically found to exhibit a log-normal distribution for field halos Bullock et al. (2001); Macciò et al. (2008), which will define our initial concentration distribution (before tidal stripping). We adopt the peak value and variance fit in Sánchez-Conde and Prada (2014), which was shown to provide a good description of cosmological simulations run independently by several groups. Subhalos are assumed to have a Navarro-Frenk-White (NFW) profile Navarro et al. (1997) with parameters set by $m$ and $c$ (the impact of choosing an Einasto profile Einasto (1965) instead of NFW was investigated in Stref and Lavalle (2017)).

The subhalo population is strongly affected by tidal interactions with the potential of the host galaxy (Binney and Tremaine, 1987). This is accounted for in the model through the calculation of a tidal radius $r_{\rm t}$ for each subhalo. The tidal radius should be interpreted as the physical extension of a subhalo, which is in general smaller than the extension it would have on a flat background. The physical mass of a subhalo is then

[TABLE]

where $\rho_{\rm sub}$ is the subhalo inner mass density profile. In our modeling, tidal stripping only removes the outer layers of subhalos while leaving the inner parts unchanged. In reality, DM should rearrange itself into a new equilibrium state. The central density however, should be left essentially unchanged Drakos et al. (2017). Since we are interested in indirect searches for annihilating DM, and the central density gives the dominant contribution to the annihilation rate, our modeling should lead to a reasonable approximation. Two important contributions are accounted for in SL17: the effect of the smooth Galactic potential (including both DM and baryons), and the gravitational shocking induced by the baryonic disk Ostriker et al. (1972); Gnedin and Ostriker (1999). The latter effect turns out to be very efficient at stripping subhalos in the inner 20 kpc of the Galaxy, a result also found in numerical studies D’Onghia et al. (2010); Errani et al. (2017); Kelley et al. (2018). The strength of SL17 over simulations is that it accounts for the constrained potential of the MW, with a detailed description of baryons.

2.2 Subhalo disruption?

Whether a subhalo can be completely disrupted by tidal effects is an open question. A number of numerical studies have found that a subhalo is completely disrupted when the total energy gained through tidal-stripping or disk-shocking effects is comparable to the binding energy (Hayashi et al., 2003; D’Onghia et al., 2010). On the other hand, some studies (Goerdt et al., 2007; Berezinsky et al., 2008; Peñarrubia et al., 2010) have found that cuspy subhalos almost always survive mass loss, leaving a small bound remnant behind even after gaining an energy far greater than their binding energy. These contradictory results may have been reconciled in a recent series of papers by van den Bosch and collaborators (van den Bosch, 2017; van den Bosch et al., 2018; van den Bosch and Ogiya, 2018). In these studies, it is shown that subhalo disruption in N-body simulations can actually be entirely explained by numerical artifacts. In particular, disruption is shown to be highly sensitive to the value of the force-softening length. If this length is taken sufficiently small, the authors show that subhalos survive tidal mass loss in the form of a small bound remnant. We aim at quantifying the impact of these results on the whole subhalo population. Tidal disruption is modeled in a very simple way in SL17: given a subhalo with scale radius $r_{\rm s}$ and tidal radius $r_{\rm t}$, we assume

[TABLE]

In Eq. (5), $\epsilon_{\rm t}$ is a dimensionless free parameter assumed universal i.e. independent of the subhalo’s mass, concentration or position. In SL17, the value of the disruption parameter was set to $\epsilon_{\rm t}=1$ in agreement with numerical results, see e.g. Hayashi et al. (2003). The results of van den Bosch and collaborators point toward a much lower value for $\epsilon_{\rm t}$. In this work, we consider two extreme values: $\epsilon_{\rm t}=1$ and $\epsilon_{\rm t}=0.01$. The latter means a subhalo is disrupted when it has lost around 99.99% of its mass. In the following, we refer to these two configurations as "fragile subhalos" ($\epsilon_{\rm t}=1$) and "resilient subhalos" ($\epsilon_{\rm t}=0.01$). The final subhalo phase-space density can now be written

[TABLE]

where $N_{\rm tot}$ is the total number of substructures within the virial radius of the Milky Way and $K_{\rm tot}$ is a normalization factor:

[TABLE]

2.3 Calibration procedure

In its current version, the SL17 model requires a calibration of the subhalo abundance in a given mass range (this will change in future versions). To be consistent with results from the highest-resolution simulations available, the calibration is done by demanding that the subhalo mass fraction is similar to what is found in the dark matter-only Via Lactea II simulation Diemand et al. (2008). This amounts to 11% of the total dark halo mass in the form of subhalos in the virial mass range $[m_{1},m_{2}]=[2.2\times 10^{-6}{M_{\rm DM}},8.8\times 10^{-4}{M_{\rm DM}}]$ where $M_{\rm DM}$ is the total DM mass of the Galaxy. We stress that these numbers are expressed in terms of virial masses, not tidal masses (see Stref and Lavalle (2017) for further details). To reproduce the (likely overestimated) tidal disruption efficiency in simulations, we have to set $\epsilon_{\rm t}=1$ at the calibration stage, and we also neglect the impact of baryons. The disruption efficiency parameter $\epsilon_{\rm t}$ can safely be changed after the calibration has been completed. It is much safer to perform this calibration on dark matter-only simulations because the tidal stripping induced by baryons strongly depends on the details of the stellar distribution, which is acutely constrained in the Milky Way. More formally, the normalization procedure reads:

[TABLE]

Fixing $f_{\rm sub}(m_{1},m_{2})=0.11$ leads to the total number of clumps $N_{{\rm DMO},\,\epsilon_{\rm t}=1}$ in the simulation-like configuration. Note that this value assumes that $m$ is really $m_{200}$ in the equation above, not the tidal mass.

Now that the model is properly calibrated, we incorporate all the effects that are not included in the calibration i.e. the tidal effects due to the baryons and possibly $\epsilon_{\rm t}<1$. This is done by assuming that subhalos in the outskirts of the Galaxy are not affected by baryonic tides or the value of $\epsilon_{\rm t}$. This is motivated by the observation that tidal effects are inefficient far from the center of the Galaxy, and subhalos almost behave like isolated halos. The DM mass within clumps per unit of volume can be expressed as

[TABLE]

where $m_{\rm t}$ is the tidal subhalo mass introduced in Eq. (4). Equating $\left<\rho_{\rm cl}\right>(r_{200})$, where $r_{200}$ is the virial radius of the Galaxy, in the DM-only+$\epsilon_{\rm t}=1$ configuration, to the same quantity in the realistic configuration (including baryons and $\epsilon_{\rm t}\leq 1$) leads to the simple relation

[TABLE]

The two normalization factors $K$ can be computed using Eq. (7) and we get the value of $N_{\rm tot}$. The number of subhalos within the Solar radius $r_{\odot}=8.21$ kpc is shown in Tab. 1. This number is highly sensitive to the parameters of the mass function $\alpha_{\rm m}$ and $m_{\rm min}$ as already shown in Stref and Lavalle (2017). Furthermore, it is quite sensitive to $\epsilon_{\rm t}$: going from $\epsilon_{\rm t}=1$ to $\epsilon_{\rm t}=0.01$, the number of subhalos increases by at least an order of magnitude. The impact of $\epsilon_{\rm t}$ on the subhalo mass and number density is investigated in the next section.

3 Mass and number densities of subhalos

In this section, we compute the mass density within subhalos as well as the subhalo number density. The subhalo mass density is defined in Eq. (9). Once the subhalo density is known, it is used to determine the amount of DM smoothly distributed across the Galaxy through Eq. (2). The DM mass inside subhalos is compared to the total DM density $\rho_{\rm DM}$ in Fig. 2. The mass density in the form of subhalos is predicted to be much higher, by orders of magnitude, for resilient subhalos than for fragile subhalos. The former case is also more justified theoretically, although the latter one allows us to compare with very conservative assumptions. At the position of the Solar System, the impact is around one order of magnitude. Although these are large differences, we note the subhalo mass density is still far below the total DM density. This means that most of the DM mass within the orbit of the Sun is smoothly distributed rather than clumpy, irrespective of the efficiency of the tidal disruption set by $\epsilon_{\rm t}$.

The subhalo number density in SL17 can be formally written

[TABLE]

The results obtained are shown in Fig. 3. The number density, just like the mass density, is highly sensitive to $\alpha_{\rm m}$ and $m_{\rm min}$, as well as the disruption parameter. Interestingly, the values we get in the Solar neighborhood are comparable to the local number density of stars $n_{*}\sim 1\,\rm pc^{-3}$. For a low value of the minimal mass $m_{\rm min}$, the subhalo number density can even be much higher, possibly going as high as $10^{5}\,\rm pc^{-3}$. This could have a number of interesting implications for the interactions between subhalos and stars. The tidal heating of subhalos by stars has been investigated in a number of studies Berezinsky et al. (2003, 2006); Zhao et al. (2007); Goerdt et al. (2007); Angus and Zhao (2007); Green and Goodwin (2007); Schneider et al. (2010), with different conclusions. In the next section, we look at the impact of our results on indirect DM searches.

4 Impact on indirect searches for annihilating dark matter

In this section, we quantify the impact of Galactic clumps on indirect searches for self-annihilating DM. Inhomogeneities are known to enhance the DM annihilation rate in Galactic halos (Silk and Stebbins, 1993). We compute the local DM self-annihilation rate and evaluate the enhancement due to the survival of clump remnants, referred to as the boost factor. Two complementary channels are then investigated: gamma rays and antiproton cosmic rays.

4.1 Annihilation profiles and local boost factors

The number of self-annihilation of DM particles at position $\vec{r}$ is proportional to $\rho^{2}(\vec{r})$ where $\rho$ is the DM mass density. If subhalos are discarded, the Galactic annihilation profile is

[TABLE]

Let us now consistently include the contribution of subhalos. The luminosity of a single clump is

[TABLE]

where $V_{\rm t}(\vec{r},m,c)$ is the volume of the clump within its tidal radius. The annihilation of the full subhalo population is simply obtained by integrating the luminosity of a single object over the subhalo phase-space number density

[TABLE]

The full annihilation profile must also incorporate the annihilation in the smooth halo (different from $\mathcal{L}_{0}$ which is the density assuming all the DM is smoothly distributed), as well as the annihilation of subhalo particles onto smooth halo particles. The first contribution can be written

[TABLE]

where $\rho_{\rm sm}(\vec{r})=\rho_{\rm DM}(\vec{r})-\left<\rho_{\rm cl}\right>(\vec{r})$ is the smooth DM density. The clump-smooth contribution is

[TABLE]

The total annihilation profile is simply the sum of all contributions

[TABLE]

This should be compared to $\mathcal{L}_{0}(\vec{r})$ to evaluate the impact of clustering on the annihilation rate. This is usually done in terms of a boost factor which we define as

[TABLE]

This is not quite the boost factor used in indirect searches, which is defined through a ratio of fluxes (see Eq. (23) and Eq. 29)). The boost in Eq. (20) is rather the local increase in the annihilation rate due to clustering. According to this definition, the boost is zero if $\mathcal{L}=\mathcal{L}_{0}$ i.e. substructures are not included.222This differs by one unit from the definition used in Stref and Lavalle (2017).

The annihilation profiles are shown on the top panels in Fig. 4, and the associated boost factors are shown on the bottom panels. As already shown a long time ago, see e.g. Bergstrom et al. (1999); Berezinsky et al. (2003), the boost is an increasing function of the galactocentric radius $r=|\vec{r}|$. This is due to the morphology of the annihilation profiles which is modified by the inclusion of clumps: we have $\mathcal{L}_{\rm cl}\propto\rho_{\rm DM}$ while $\mathcal{L}_{0}\propto\rho_{\rm DM}^{2}$. Also noticeable is the high sensitivity of the annihilation profile to the mass function index $\alpha_{\rm m}$ which is by far the largest source of uncertainty on the clump contribution. The value of the disruption parameter $\epsilon_{\rm t}$ has almost no impact on the boost above 20 kpc due to the ineffectiveness of tidal effects far from the center of the Galaxy. Below 20 kpc, the boost is strongly sensitive to the disruption parameter. In the inner few kiloparsecs, fixing $\epsilon_{\rm t}=0.01$ leads to a boost orders of magnitude larger than in the $\epsilon_{\rm t}=1$ configuration. We have however $\mathcal{B}(\vec{r})\ll 1$ below 3 kpc regardless of the value of $\epsilon_{\rm t}$, meaning $\mathcal{L}\simeq\mathcal{L}_{0}$ and subhalos do not have any impact on the annihilation rate in that region. The region where the impact of $\epsilon_{\rm t}$ on the annihilation profile is the most important is located between 3 kpc and 10 kpc. This region coincidentally includes the Solar system, located at $r\simeq 8$ kpc. This motivates a more detailed investigation of two standard annihilation channels: gamma rays and cosmic-ray antiprotons, which are sensitive to different annihilation regions.

4.2 Application to gamma rays

The energy-differential flux of gamma rays originating from DM self-annihilation is, on a given line of sight

[TABLE]

where $\left<\sigma v\right>$ is the thermally-averaged annihilation cross-section, $m_{\chi}$ is the DM mass, $\mathrm{d}N_{\gamma}/\mathrm{d}E$ is the gamma-ray spectrum at annihilation and $s$ the distance coordinate along the line of sight.333The expression of $\mathrm{d}\Phi_{\gamma}/\mathrm{d}E$ should be multiplied by $1/2$ if the DM particle is not its own antiparticle. If the annihilation cross-section is velocity-independent, astrophysical ingredients enter only through the J-factor defined as

[TABLE]

where $\psi$ is the angle between the direction of the Galactic center and the line of sight (spherical symmetry of the dark halo is assumed). The impact of small-scale clustering on this J-factor has been considered in a number of studies (Bergstrom et al., 1999; Diemand et al., 2007; Strigari et al., 2007; Pieri et al., 2008, 2011; Blanchet and Lavalle, 2012; Sánchez-Conde and Prada, 2014; Ishiyama, 2014; Bartels and Ando, 2015; Moliné et al., 2017; Hiroshima et al., 2018). We define the gamma-ray boost factor as

[TABLE]

where $J_{0}=\int\mathrm{d}s\,\mathcal{L}_{0}$ is the J-factor without subhalos. Unlike the local boost $\mathcal{B}$ in Eq. (20), the gamma-ray boost $\mathcal{B}_{\gamma}$ depends on the line of sight rather than the position in the Galaxy Bergstrom et al. (1999). The boost is shown as a function of $\psi$ in Fig. 5. The growth of the local boost $\mathcal{B}(\vec{r})$ as a function of $r$ translates into a growth of $\mathcal{B}_{\gamma}(\psi)$ as a function of $\psi$ i.e. the maximal gamma-ray boost is reached at the anti-center. This maximal boost ranges from 0.5 to 9, depending on the values of $\alpha_{\rm m}$ and $m_{\rm min}$. The survival of clumps noticeably increases the boost at all latitudes. The gain is greater at small latitudes where substructures are more impacted by tidal effects. Below $\psi\simeq 40$ deg, the boost is increased by a factor of at least two in all configurations. This should have important consequences for indirect searches using gamma rays, especially at high latitudes. Interestingly, high latitudes have been shown to be a very sensitive probe of DM annihilation even without the inclusion of clumps Chang et al. (2018).

4.3 Application to cosmic-ray antiprotons

Charged cosmic rays constitute an indirect detection channel complementary to gamma rays (Maurin et al., 2002; Lavalle and Salati, 2012). Since their original proposal as a probe of DM annihilation (Silk and Srednicki, 1984), cosmic-ray antiprotons have been shown to be especially sensitive, see e.g. Bergström et al. (1999); Donato et al. (2004); Barrau et al. (2005); Bringmann and Salati (2007); Cirelli et al. (2009); Boudaud et al. (2015); Giesen et al. (2015). Antiprotons have been the subject of much scrutiny since the latest measurement of the antiproton flux performed by the AMS-2 collaboration (Aguilar et al., 2016). A number of studies (Cuoco et al., 2017; Cui et al., 2017; Reinert and Winkler, 2018; Cui et al., 2018; Cuoco et al., 2019; Cholis et al., 2019) have found a discrepancy between the measured flux and a purely secondary origin of antiprotons. This discrepancy could be interpreted as evidence for annihilating DM, although the significance of the excess is debated as it depends on the propagation model used and the modeling of systematic uncertainties. In this context, it is worth evaluating systematic uncertainties coming from the small-scale structuring.

Since antiprotons have a random motion due to their diffusion on the inhomogeneities of the magnetic halo, their detection gives little information on their source. This implies that the antiproton boost factor, and the boost for charged cosmic rays in general, is not direction-dependent unlike for gamma rays. Instead, this boost is energy-dependent (Lavalle et al., 2007) and has been shown to be mild, at most a factor of two (Lavalle et al., 2008; Pieri et al., 2011). Although smaller than the gamma-ray boost, this can still be larger than the systematic uncertainties on cosmic-ray propagation. This motivates a new computation of the boost, which we perform here. The antiproton boost factor is defined as

[TABLE]

where $T$ is the antiproton kinetic energy, $\Phi_{\overline{p}}$ is the DM-induced antiproton flux including the subhalo contribution and $\Phi_{\overline{p},0}$ the same flux assuming all the DM is smoothly distributed. To get the flux, one must solve the cosmic-ray steady-state propagation equation

[TABLE]

which accounts for spatial diffusion, convection, energy losses, diffusive reacceleration, and spallation processes in the disk (taken as infinitely thin). In Eq. (25), $\Psi$ is the antiproton number density per unit energy which is related to the flux through $\Phi_{\overline{p}}=v_{\overline{p}}/(4\pi)\times\Psi$ where $v_{\overline{p}}$ is the antiproton speed. Antiprotons are sourced by DM annihilation

[TABLE]

where $\mathrm{d}N_{\overline{p}}/\mathrm{d}E$ is the antiproton spectrum at annihilation.444If DM is not its own antiparticle, $Q$ should be divided by 2. Several unknown propagation parameters enter Eq. (25). These can be constrained using the measured boron to carbon ratio (B/C) (Strong and Moskalenko, 1998). We use the best-fit model derived by Reinert and Winkler (2018), which includes an energy-break in the diffusion coefficient. The B/C ratio can only constrain $K_{0}/L$ where $K_{0}$ is the normalization of the diffusion coefficient and $L$ is the half-height of the magnetic halo. As shown in Donato et al. (2004), the DM-induced antiproton flux depends crucially on $L$ hence we consider two extremal values in this work. A lower bound on $L$ can be obtained from low-energy positron data (Lavalle et al., 2014) and the authors of Reinert and Winkler (2018) find $L=4.1\,\rm kpc$. For the largest value, we take $L=15\,\rm kpc$. According to the analysis of Reinert and Winkler (2018), the B/C data are consistent with negligible reacceleration. Furthermore we neglect energy losses which are unimportant for high-energy antiprotons. The resulting transport equation can be solved semi-analytically using the Green’s function formalism (see (Lavalle et al., 2008) for the solution) and the differential flux can be written

[TABLE]

Since all the energy-dependent terms have been neglected in Eq. (25), the energy part of the Green’s function is trivial

[TABLE]

and the boost factor can be simply written

[TABLE]

The boost factor is shown as a function of the antiproton kinetic energy in Fig. 6. We first note that the boost is roughly energy-independent. This is because antiprotons probe the entire volume of the magnetic halo during their lifetime, independently of their energy. This is not true at low energies, below a few GeVs, where energy losses become relevant. The half-height of the magnetic halo has a small impact on the boost, with $L=15\,\rm kpc$ leading to a slightly larger $\mathcal{B}_{\overline{p}}$ than $L=4.1\,\rm kpc$. The main source of uncertainties are coming from the subhalo parameters $\alpha_{\rm m}$, $m_{\rm min}$ and $\epsilon_{\rm t}$. The survival parameter $\epsilon_{\rm t}$ has a significant impact on the result, with a small value $\epsilon_{\rm t}=0.01$ leading to a boost roughly twice as large as in the $\epsilon_{\rm t}=1$ case. As for gamma rays, the most critical parameter is $\alpha_{\rm m}$. For a low value $\alpha_{\rm m}=1.9$, the boost never exceeds 10% while it is always higher for $\alpha_{\rm m}=2$. Overall, playing with the propagation and subhalo parameters, we find that the antiproton boost can conservatively range from $2\%$ to $140\%$. These values are in agreement with earlier results. Although it is conservative to ignore the small-scale clustering when deriving limits on the annihilation cross-section using data, the boost should be included when interpreting an excess as a signature of DM annihilation. Indeed, a factor of two in the DM contribution would change the inferred mass and cross-section of the hypothetical DM particle.

5 Conclusion

Subhalos suffer mass loss due to their interaction with the tidal field of the Galaxy, which makes their modeling very challenging. Consequently, most subhalo models rely at least partly on numerical simulations to calibrate their predictions. However, it was recently shown that numerical simulations might not properly account for the tidal disruption of subhalos, as artificial effects lead to a serious overestimation of the efficiency of these processes van den Bosch et al. (2018); van den Bosch and Ogiya (2018). We note that the resistance of subhalos to tidal stripping is further supported by theoretical arguments, like adiabatic invariance that should prevail in their inner parts Gnedin and Ostriker (1999), as already emphasized in Stref and Lavalle (2017). We have derived some of the consequences of these results using the semi-analytical Galactic subhalo population model of Stref and Lavalle (2017), assuming a tidal disruption efficiency $\epsilon_{\rm t}=0.01$. We predict the spatial dependence of the subhalo properties due to tides induced both by the global gravitational potential and baryonic disk shocking. We remind the reader that this model is built from constrained mass models for the Milky Way and is therefore consistent with current kinematic constraints, which is usually not the case in extrapolations from "Milky-Way-like" simulations. We find that the local mass density is still dominated by the smooth component of the dark halo. The local number density of subhalos is increased by roughly one order of magnitude with respect to estimates based on a tidal disruption efficiency similar to that inferred from simulations ($\epsilon_{\rm t}=1$). This makes the subhalo number density comparable, for a broad range of minimal subhalo mass $m_{\rm min}$, to the local star number density. Since our description of the subhalo population relies solely on gravitational principles, it should be valid for a wide range of cold DM candidate such as WIMPs, axions or primordial black holes. One only needs to modify the mass function, in particular $m_{\rm min}$, to explore alternatives to the WIMP scenario. The resilience of subhalos increases the local WIMP annihilation rate, which in turn affects the predictions for indirect searches. For gamma rays, we find that the boost factor is increased by at least a factor of two for $\psi<40$ deg, and slightly less for higher values of $\psi$. The boost factor for antiprotons is also increased by a rough factor of two if subhalos are resilient to tidal disruption. For a complemetary study comparing the SL17 model to simulation results regarding indirect searches in gamma rays and neutrinos, we refer the reader to Hütten et al. (2019).

For future work, we plan on including a more detailed mass function directly deriving from the primordial power spectrum, as well as the tidal heating of subhalos due to individual stars and study the consequences of having a large population of small objects in the Solar neighborhood.

\authorcontributions

conceptualization, all authors; software, M.S; investigation, M.S.; writing–original draft preparation, M.S.; writing–review and editing, T.L. and J.L.; supervision, J.L.;

Acknowledgements.

We wish to thank M. A. Sánchez-Conde and M. Doro for inviting us to contribute to this topical review. \fundingJ. L. and M. S. are partly supported by the Agence Nationale pour la Recherche (ANR) Project No. ANR-18-CE31-0006, the Origines, Constituants, et EVolution de l’Univers (OCEVU) Labex (No. ANR-11-LABX-0060), the CNRS IN2P3-Theory/INSU-PNHE-PNCG project “Galactic Dark Matter,” and the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreements No. 690575 and No. 674896, in addition to recurrent funding by Centre National de la Recherche Scientifique (CNRS) and the University of Montpellier. T. L. is supported by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 713366. The work of TL was also supported by the Spanish Agencia Estatal de Investigación through the grants PGC2018-095161-B-I00, IFT Centro de Excelencia Severo Ochoa SEV-2016-0597, and Red Consolider MultiDark FPA2017-90566-REDC. \conflictsofinterestThe authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. \externalbibliographyyes