$\gamma$-ray and $\nu$ searches for dark matter subhalos in the Milky Way with a baryonic potential

TL;DR

This study investigates how the baryonic potential in the Milky Way affects dark matter subhalo distributions and their gamma-ray and neutrino signals, revealing significant impacts on detection prospects and subhalo characteristics.

Contribution

It introduces updated subhalo spatial distributions considering baryonic effects into the CLUMPY code, providing new insights into subhalo signals and properties.

Findings

Brightest subhalo flux decreases by a factor of 2 to 7 for annihilating DM.

No impact on decay signals from subhalos.

Mass distribution of brightest subhalo suggests different natures depending on tidal effects.

Abstract

The distribution of dark matter (DM) subhalos in our Galaxy remains disputed, leading to varying $\gamma$-ray and $\nu$ flux predictions from their annihilation or decay. In this work, we study how, in the inner Galaxy, subhalo tidal disruption from the Galactic baryonic potential impacts these signals. Based on state-of-the art modeling of this effect from numerical simulations and semi-analytical results, updated subhalo spatial distributions are derived and included in the CLUMPY code. The latter is used to produce a thousand realizations of the $\gamma$-ray and $\nu$ sky. Compared to predictions based on DM only, we conclude a decrease of the flux of the brightest subhalo by a factor 2 to 7 for annihilating DM and no impact on decaying DM: the discovery prospects or limits subhalos can set on DM candidates are affected by the same factor. This study also provides probability density functions for the distance, mass, and angular distribution of the brightest subhalo, among which the mass may hint at its nature: it is most likely a dwarf spheroidal galaxy in the case of strong tidal effects from the baryonic potential, whereas it is lighter and possibly a dark halo for DM only or less pronounced tidal effects.

Figures

Tables

Table 1. Table 1: Subhalo parameters for the models investigated in this study, with model #1 based on results of DM-only numerical simulations, while models #2 to #4 are different implementations of DM subhalos post-processed in the Milky Way halo and baryonic disk potential. For models #2, #3, and #4, we also show the number of surviving subhalos with tidal masses between 10 8 ​ M ⊙ superscript 10 8 subscript M direct-product 10^{8}\,\rm M_{\odot} and 10 10 ​ M ⊙ superscript 10 10 subscript M direct-product 10^{10}\,\rm M_{\odot} . See Sect.  3 for details and parameters common to all subhalo configurations.

	Model #1	Model #2	Model #3	Model #4
$\frac{d𝒫}{dV}$	Aquarius Springel et al. (2008)	Phat-ELVIS Kelley et al. (2018)	SL17 Stref and Lavalle (2017) with ${ϵ}_t={10}^{-2}$​	​ SL 17 Stref and Lavalle (2017) with ${ϵ}_t=1$​
	Einasto	​​​Sigmoid-Einasto Eq. (5)​​​	$\propto {\rho }_{\mathrm{sm}}$	$\propto {\rho }_{\mathrm{sm}}$
	${\alpha }_E=0.68$	${\alpha }_E=0.68$	NFW ⋆	NFW ⋆
	$r_{-2}=199$ kpc	$r_{-2}=128$ kpc	$r_{-2}=r_s=19.6$ kpc ⋆	$r_{-2}=r_s=19.6$ kpc ⋆
	-	$r_0=29.2$ kpc	-	-
	-	$r_c=4.24$ kpc	-	-
$N_{\text{calib}}$	300	-	276 ⋆	276 ⋆
$N_{\mathrm{surviving}}$	-	90	$114\pm 11$	$112\pm 10$
$c\left(m\right)$	​Moliné et al. (2017)​	Moliné et al. (2017)	​Sánchez-Conde & Prada Sánchez-Conde and Prada (2014)​	​Sánchez-Conde & Prada Sánchez-Conde and Prada (2014)​

Table 2. Table 2: Properties of the brightest subhalos in the random realizations displayed in Fig.  4 .

	Position in	Distance	Mass ${𝒎}^{\star }$	${𝑱}_{0.5{}^{\circ }}$	${𝑱}_{\mathrm{𝐭𝐨𝐭}}$	${𝑫}_{0.5{}^{\circ }}$	${𝑫}_{\mathrm{𝐭𝐨𝐭}}$
	Map $(𝒍,𝒃)$	${𝒍}^{\star }$ $\left[\mathrm{𝐤𝐩𝐜}\right]$	$\left[{𝐌}_{\odot }\right]$	$\left[{\mathrm{𝐆𝐞𝐕}}^{\mathrm{𝟐}}{\mathrm{𝐜𝐦}}^{-\mathrm{𝟓}}\right]$		$\left[\mathrm{𝐆𝐞𝐕}{\mathrm{𝐜𝐦}}^{-\mathrm{𝟐}}\right]$
DM only	$(-54{}_{}{}^{\circ },-13{}_{}{}^{\circ })$	$16.3$	$4.0\times {10}^9$	$1.8\times {10}^{20}$	$5.5\times {10}^{20}$	$1.1\times {10}^{19}$	$1.8\times {10}^{21}$
Phat-ELVIS	$(-141{}_{}{}^{\circ },+31{}_{}{}^{\circ })$	$42.2$	$2.0\times {10}^9$	$1.4\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$	$2.3\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$	$3.6\times {10}^{18}$	$1.3\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟐𝟎}}$
	$(+37{}_{}{}^{\circ },+8{}_{}{}^{\circ })$	$69.5$	$3.8\times {10}^9$	$1.4\times {10}^{19}$	$2.0\times {10}^{19}$	$4.1\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟖}}$	$9.2\times {10}^{19}$
SL17, ${ϵ}_t={10}^{-2}$	$(-19{}_{}{}^{\circ },-73{}_{}{}^{\circ })$	$7.8$	$3.7\times {10}^6$	$3.3\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$	$6.1\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$	$2.3\times {10}^{18}$	$7.1\times {10}^{18}$
	$(-73{}_{}{}^{\circ },+20{}_{}{}^{\circ })$	$54.8$	$2.3\times {10}^9$	$1.0\times {10}^{19}$	$2.6\times {10}^{19}$	$4.1\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟖}}$	$9.0\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$
SL17, ${ϵ}_t=1$	$(+71{}_{}{}^{\circ },-25{}_{}{}^{\circ })$	$30.1$	$1.3\times {10}^8$	$9.0\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟖}}$	$1.5\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$	$2.3\times {10}^{18}$	$1.6\times {10}^{19}$
	$(+93{}_{}{}^{\circ },-50{}_{}{}^{\circ })$	$68.5$	$3.4\times {10}^9$	$5.2\times {10}^{18}$	$1.4\times {10}^{19}$	$3.4\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟖}}$	$8.4\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟗}}$

Equations

$$\frac{\mathrm{d}\Phi_{\gamma,\nu}}{\mathrm{d}E}\biggl{(}E,\psi,\theta,\Delta\Omega\biggl{)}^{\rm Annihil.}_{\rm Decay}=\underbrace{\frac{1}{4\pi}\,\sum_{f}\frac{dN^{f}_{\gamma,\nu}}{dE}\,B_{f}\times\left.\begin{cases}\displaystyle\frac{\langle\sigma v\rangle}{m_{\rm DM}^{2}\,\delta}\\ \vspace{-0.4cm}\\ \displaystyle\,\,\frac{1}{\tau_{\rm DM}\,m_{\rm DM}}\end{cases}\right\}}_{\text{Particle physics term: $\displaystyle\frac{\mathrm{d}\Phi^{PP}_{\gamma,\nu}}{\mathrm{d}E}(E)$}}\;\times\;\underbrace{\int_{0}^{\Delta\Omega}\!\!\int_{\rm{l.o.s.}}\!\!\mathrm{d}l\,\mathrm{d}\Omega\times\left.\begin{cases}\displaystyle\,\rho^{2}(\vec{r})\\ \vspace{0.2cm}\\ \displaystyle\,\rho(\vec{r})\end{cases}\right\}}_{\text{Astrophysics term: $J$- or $D$-factor}}\,,$$

$$\frac{\mathrm{d}{J}_{\mathrm{tot}}}{\mathrm{d}\Omega }={\int }_{l.o.s.}{\left({\rho }_{\mathrm{sm}}+\sum _i^{N_{\mathrm{subs}}}{\rho }_i\right)}^2\mathrm{d}l=\underset{\mathrm{d}{J}_{\mathrm{sm}}/\mathrm{d}\Omega }{\underbrace{\int \int {\rho }_{\mathrm{sm}}^2\mathrm{d}l}}+\underset{\mathrm{d}{J}_{\mathrm{subs}}/\mathrm{d}\Omega }{\underbrace{\int \int {\left(\sum _i^{N_{\mathrm{subs}}}{\rho }_i\right)}^2\mathrm{d}l}}+\underset{\mathrm{d}{J}_{\mathrm{cross}-\mathrm{prod}}/\mathrm{d}\Omega }{\underbrace{2\int \int {\rho }_{\mathrm{sm}}\sum _i^{N_{\mathrm{subs}}}{\rho }_i\mathrm{d}l}}.$$

$$\frac{\mathrm{d}{J}_{\mathrm{tot}}}{\mathrm{d}\Omega }=\frac{\mathrm{d}{J}_{\mathrm{sm}}}{\mathrm{d}\Omega }+⟨\frac{\mathrm{d}{J}_{\mathrm{subs}}}{\mathrm{d}\Omega }⟩+⟨\frac{\mathrm{d}{J}_{\mathrm{cross}-\mathrm{prod}}}{\mathrm{d}\Omega }⟩+\sum _{k=1}^{N_{\mathrm{subs}}^{\mathrm{drawn}}}\frac{\mathrm{d}{J}_{\mathrm{drawn}}^k}{\mathrm{d}\Omega }.$$

$$\frac{{\mathrm{d}}^3\mathcal{P}}{\mathrm{d}V\mathrm{d}m\mathrm{d}c}=\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}V}(\vec{r})\times \frac{\mathrm{d}\mathcal{P}}{\mathrm{d}m}(m)\times \frac{\mathrm{d}\mathcal{P}}{\mathrm{d}c}(c,m)\text{.}$$

$$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}V}(r)=\frac{A}{1+e^{-(r-r_0)/r_c}}\times \exp \left\{-\frac{2}{{\alpha }_E}\left[{\left(\frac{r}{r_{-2}}\right)}^{\alpha }-1\right]\right\},$$

$$⟨{\rho }_{\mathrm{tot}}(r)⟩={\rho }_{\mathrm{s}}{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{-\gamma }{\left[1+{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{\alpha }\right]}^{(\gamma -\beta )/\alpha }.$$

$$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}V}(r,\mathrm{initial})\propto ⟨{\rho }_{\mathrm{tot}}(r)⟩.$$

$$\frac{r_{\mathrm{t}}}{r_{\mathrm{s}}}<{ϵ}_{\mathrm{t}}\iff \mathrm{subhalo}\mathrm{is}\mathrm{disrupted}.$$

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\history \Title$\gamma$-ray and $\nu$ Searches for Dark-Matter Subhalos in the Milky Way with a Baryonic Potential

 \AuthorMoritz Hütten *1,∗*\orcidA, Martin Stref *2,∗*\orcidE, Céline Combet *3*\orcidC, Julien Lavalle *2*\orcidD and David Maurin *3*\orcidB

\AuthorNamesMoritz Hütten , Martin Stref, Céline Combet, Julien Lavalle, David Maurin

\corresCorrespondence: mhuetten@mpp.mpg.de (M.H.); martin.stref@umontpellier.fr (M.S.)

1 Introduction

In this contribution to The Role of Halo Substructure in Gamma-Ray Dark-Matter Searches, we revisit a previous study on the detectability of galactic dark clumps in $\gamma$-rays Hütten et al. (2016). The latter relied on the best knowledge we had a few years ago of the properties of dark-matter (DM) clumps in the Milky Way. These properties (e.g., the mass and spatial distributions of galactic subhalos) were inferred from numerical simulations with a typical mass resolution of a few $10^{5}$ $\rm M_{\odot}$, and extrapolated down ten orders of magnitude to the model-dependent minimal masses of subhalos Profumo et al. (2006); Bringmann (2009). Functional parametrizations were incorporated in the CLUMPY code Charbonnier et al. (2012); Bonnivard et al. (2016); Hütten et al. (2019) to generate $\gamma$-ray skymaps, accounting for the whole population of subhalos. For each combination of subhalo properties we explored, hundreds of skymap realizations were drawn, allowing us to calculate the average properties of the brightest clump. In the context of the future CTA $\gamma$-ray observatory Acharya et al. (2013) and its foreseen extragalactic survey, we concluded that the limits on DM set from this brightest clump should be “competitive and complementary to those based on long observation of a single bright dwarf spheroidal galaxy”.

In the recent years, numerical simulations D’Onghia et al. (2010); Zhu et al. (2016); Kelley et al. (2018) and semi-analytical studies Green and Goodwin (2007); Goerdt et al. (2007); Berezinsky et al. (2008, 2014); Stref and Lavalle (2017) have investigated the impact of the baryonic components of disk galaxies on their subhalo population by tidal stripping and disruption. These works reached the generic conclusion of a strong depletion of subhalos in the disk regions (i.e., also the Solar neighborhood), though with different quantitative estimates. Such a difference is expected due to the diversity of assumptions and methods used by different groups. In any case, this immediately questions the conclusion reached in Hütten et al. (2016), where the brightest subhalo was found, on average, at $\sim$10 kpc from the Galactic center and at similar distance from Earth. Experimental limits on DM from galactic subhalos, derived from Fermi-LAT Berlin and Hooper (2014); Bertoni et al. (2015); Schoonenberg et al. (2016); Mirabal et al. (2016); Bertoni et al. (2016); Wang et al. (2016); Liang et al. (2017); Calore et al. (2017); Hooper and Witte (2017) or expected from prolonged operation Charles et al. (2016), from HAWC Harding and Dingus (2015); Abeysekara et al. (2018), or future instruments such as GAMMA 400 Egorov et al. (2018) and CTA Brun et al. (2011); Hütten et al. (2016), should also be impacted by this result.

The paper is organized as follows: Section 2 presents the overall methodology and recalls how all relevant subhalos are efficiently accounted for with CLUMPY. Section 3 lists the subhalo critical parameters, highlighting the very different spatial distributions considered in this analysis. Section 4 presents updated statistics of the subhalo population and provides probability density functions (PDFs) of the brightest subhalo’s properties (distance to the observer, mass, brightness, etc.). The analysis is performed for both annihilating DM Gunn et al. (1978); Lake (1990); Silk and Stebbins (1993) via the so-called $J$-factors, or decaying DM Berezinsky et al. (1991); Bertone et al. (2007) ($D$-factors). We also show one realization of a subhalo skymap for all configurations considered. We conclude and briefly comment on the consequences for DM indirect detection limits in Section 5.

2 Important Quantities and Methodology

2.1 $\gamma$-ray and $\nu$ Fluxes from Dark Matter

The $\gamma$-ray or $\nu$ flux from annihilating/decaying DM particles, at energy $E$, along the line-of-sight (l.o.s.) in the direction $(\psi,\theta)$, and integrated over the solid angle $\Delta\Omega=2\pi\,(1-\cos\,\alpha_{\rm int})$, is given by

[TABLE]

where we made the distinction between the cases of DM self-annihilation (top) and decay (bottom). In both cases, $m_{\chi}$ is the mass of the DM particle, $\mathrm{d}N_{\gamma}^{f}/\mathrm{d}E$ and $B_{f}$ correspond to the spectrum per interaction and branching ratio of annihilation or decay channel $f$, and $l$ is the distance from the observer. In case of annihilation, $\langle\sigma v\rangle$ is the velocity-averaged cross-section111We assume here the DM particle to be a Majorana particle, so that $\delta=2$ (for a Dirac, $\delta=4$)., while $\tau_{\rm DM}$ is the DM particle lifetime in the decay scenario. Finally, $\rho(\vec{r})$ is the overall Galactic DM density distribution. The latter can be cast formally as the sum of a smooth distribution $\rho_{\rm sm}$ of unclustered DM particles, and a collection of $i=1\dots N_{\rm subs}$ subhalos, each with a density $\rho_{i}(\vec{r})$. The astrophysical term for annihilating DM222For conciseness, we present in Equations (2) and (3) formulae for annihilating DM only. Analogous formulae without a cross-product term and linear in the DM density can also be written for the $D$-factor of decaying DM. then reads Bergström et al. (1998); Charbonnier et al. (2012); Bonnivard et al. (2016),

[TABLE]

The above formula corresponds to a single realization of the underlying distribution of subhalos in the galaxy. The statistical properties of this distribution can be partly obtained from the formalism of hierarchical structure formation (e.g., Press and Schechter (1974)) or extracted from numerical simulations, as discussed below.

2.2 Generating Skymaps with CLUMPY v3.0

For all our calculations, we rely on the public CLUMPY code described in Charbonnier et al. (2012); Bonnivard et al. (2016); Hütten et al. (2019). It is a flexible code that efficiently emulates the end-product of numerical simulations in terms of $\gamma$-ray and neutrino signals for DM annihilation or decay. It allows easy exploration of how results are affected when changing the properties of the DM halos. CLUMPY v2.0 Bonnivard et al. (2016) was used for this purpose, to estimate the sensitivity of the CTA Hütten et al. (2016) and HAWC Abeysekara et al. (2018) $\gamma$-ray telescopes to galactic DM subhalos. Aside from galactic subhalo studies, CLUMPY has also been used by several teams to model DM annihilation or decay in $\gamma$-rays or $\nu$ in many targets: dwarf spheroidal galaxies Walker et al. (2011); Charbonnier et al. (2011); Bonnivard et al. (2015a, b, c, 2016); Genina and Fairbairn (2016); Walker et al. (2016); Campos et al. (2017); Albert et al. (2018), the galactic halo The ANTARES Collaboration (2015); Albert et al. (2017); Balázs et al. (2017), the Smith HI cloud Nichols et al. (2014), nearby galaxies Albert et al. (2018), galaxy clusters Combet et al. (2012); Nezri et al. (2012); Maurin et al. (2012), and also for the extragalactic diffuse emission Hütten et al. (2018).

The present analysis is performed with CLUMPY v3.0 Hütten et al. (2019).333When releasing CLUMPY v3.0, we corrected a misprint that was present in v2.0, related to our implementation of the virial overdensity from Bryan and Norman (1998). This issue was responsible for obtaining in Hütten et al. (2016) about a factor 3 more subhalos than expected per flux decade (see full details in the CLUMPY documentation). We recall that in Hütten et al. (2016), we found that galactic variance is responsible for a factor $\lesssim 10$ uncertainty on the value of the brightest subhalo, and that other subhalo properties were responsible for another factor $\sim 6$. Given these very large uncertainties, the conclusions on DM limits set from dark clumps with CLUMPY v2.0 are not qualitatively changed, but we urge users to rely on CLUMPY v3.0 for future studies. For completeness, we recap below the main steps of the CLUMPY calculation used for this work:

• •

CLUMPY* and the particle physics term:* Equation (1) shows that the particle physics term and the astrophysical terms are decoupled.444Strictly speaking, this factorization holds true only for DM candidates for which $\langle\sigma v\rangle$ is independent of the velocity and consideration of small redshift cells, $\Delta z/z\ll 1$. As the flux depends on the specific DM candidate chosen, we provide results in terms of $J$- and $D$-factors only; CLUMPY can easily be used to transform those into $\gamma$-ray or $\nu$ fluxes for any user-defined DM candidate (see CLUMPY’s online documentation555https://lpsc.in2p3.fr/clumpy).

• •

CLUMPY* and the astrophysics term:* to calculate skymaps of $\mathrm{d}J/\mathrm{d}\Omega$, one should rely in principle on Equation (2). However, this is impractical in terms of computing time, as \sim$$10^{14} subhalos are expected in a Milky Way-sized DM halo. This problem can be overcome by formally separating Equation (2) in an average and “resolved” component,

[TABLE]

With this ansatz, only a limited number $N_{\rm subs}^{\rm{drawn}}$ of subhalos need to have their $J$-factor profiles calculated individually, while an average description is sought for the remaining “unresolved” DM. The criterion to discriminate between resolved and unresolved components often relies on a simple subhalo mass threshold, e.g., as done in works directly relying on numerical simulations Springel et al. (2008) or their subhalo catalogs Lange and Chu (2015). CLUMPY has been developed to treat this problem in a more efficient way, acknowledging the fact that rather light, but close-by subhalos may show $J$-factors comparable to heavier, more distant objects. The CLUMPY approach relies on the notion that the overall DM signal fluctuates around an average description, $\langle J_{\rm tot}\rangle\pm\sigma_{J_{\rm subs}{}}$, and we refer to Charbonnier et al. (2012) for a detailed description of our criterion to accordingly discriminate between unresolved and resolved halos. For the purpose of this work (and also the previous Hütten et al. (2016)), this approach allows us to preselect halos likely to shine bright at Earth and to consider all decades down to the smallest subhalo masses in the calculation.

3 Modeling the Galactic Subhalo Distribution

In this study, we focus on the impact of tidal disruption of subhalos in interaction with the baryonic components of the Milky Way, and compare four parametric models of the resulting galactic subhalo abundance and their $J$/$D$-factors. We consider three quantities to be most sensitive to the $J$/$D$-factor distribution: (i) the spatial PDF of subhalos in the Milky Way, $\mathrm{d}{\cal P}/\mathrm{d}V$, (ii) the mass-concentration666The concentration is defined to be $c=r_{\Delta}/r_{-2}$, with $r_{\Delta}$, taken to be the subhalo boundary, is the radius at which the mean subhalo density is $\Delta$ times the critical density (see, e.g., Hütten et al. (2019)), and $r_{-2}$ is the position in the subhalo for which the slope of the density is $-2$. We use $\Delta=200$ in this work. relationship, $c(m,r)$, and (iii) the calibration of the total number of subhalos in the Milky Way, $N_{\rm calib}$. The latter number is determined from numerical simulations, in a range where subhalos are resolved. Here, $N_{\rm calib}$ is defined for the mass range $10^{8}$–$10^{10}~{}\rm M_{\odot}$, and it typically falls in the range 100–300.777This range was recently shown to be in agreement with the observed number of dwarf spheroidal galaxies SDSS corrected by the detection efficiency Kim et al. (2018), alleviating the tension caused by the so-called missing satellite problem in CDM scenarios Klypin et al. (1999); Strigari et al. (2007). Given the minimal mass of subhalos, $m_{\rm min}$, $N_{\rm calib}$ can be used to calculate the mass fraction, $f_{\rm DM}$, of DM in subhalos.

For modeling the subhalo distribution with these parameters, we start from an “unevolved” distribution, where we assume the position and mass of a subhalo to be uncorrelated,

[TABLE]

Here, ${\mathrm{d}{\cal P}}/{\mathrm{d}m}$ and ${\mathrm{d}{\cal P}}/{\mathrm{d}c}$ describe the PDFs for a subhalo to have a given mass and a given concentration $c$. In reality, the factorization in Equation (4) may break down when subhalos gravitationally interact with the DM and baryonic potentials of their host halo Han et al. (2016); Stref and Lavalle (2017), entangling their mass and positional distributions. We will consider this effect by “evolving” the distribution of Equation (4) in the model presented in Section 3.2.3.

3.1 Fixed Subhalo-Related Quantities

This work focuses on the impact of a baryonic disk potential on the subhalo population, mainly through the spatial PDF $\mathrm{d}{\cal P}/\mathrm{d}V$ (see Section 3.2). We keep several other subhalo-related quantities fixed to isolate this effect. For details on these parameters and how they affect the subhalo emission, we refer to our earlier work in Hütten et al. (2016) and only provide a brief summary below:

• •

Index $\alpha_{m}$ of the power-law subhalo mass PDF $\mathrm{d}{\cal P}/\mathrm{d}m\propto m^{-\alpha_{m}{}}$ and subhalo mass range: We choose $\alpha_{m}{}=1.9$, $m_{\rm min}=10^{-6}\,\rm M_{\odot}$, and $m_{\rm max}=1.3\times 10^{10}\,\rm M_{\odot}$. The maximum clump mass for all models is set to $10^{-2}\times M_{200}$ of the NFW halo from Section 3.2.3. This is motivated by the fact that we do not consider the possibility of any subhalos heavier than the Magellanic clouds, the heaviest satellites of our galaxy. The minimal clump mass and $\alpha_{m}$ mostly affect the diffuse emission boost from unresolved halos. For a fixed normalization $N_{\rm calib}$, a steeper mass function ($\alpha_{m}=2.0$) decreases the number of bright halos ($J\gtrsim 10^{20}~{}{\rm GeV}^{2}~{}{\rm cm}^{-3}$) by not more than $\sim$30%.

• •

Width of $\mathrm{d}{\cal P}/\mathrm{d}c$: We set $\sigma_{c}=0.14$ Sánchez-Conde and Prada (2014). Using a larger scatter $\sigma_{c}=0.24$ Bullock et al. (2001) increases only by a few percent the number of subhalos per flux decade. Reducing the scatter or no scatter has the opposite effect.

• •

Subhalo density profile: We model all subhalos with a spherically symmetric NFW profile Navarro et al. (1996). Using an Einasto profile Navarro et al. (1996); Einasto and Haud (1989) instead amounts to a global increase $\sim$2 of the number of subhalos per flux decade within the considered integration regions $\Delta\Omega$. Please note that micro-halos with $m\ll\rm M_{\odot}$ may show steeper inner slopes (Ishiyama et al., 2010; Ishiyama, 2014; Angulo et al., 2017); however, we have found that these micro-halos do not provide new bright, resolved subhalo candidates Hütten et al. (2016).

• •

Level of sub-substructures: We do not consider an emission boost from substructure within subhalos. Such a boost from additional levels of substructure888As shown in Bonnivard et al. (2016) (see their Figure 1), only the first level of substructure significantly boosts the halo luminosity, the next levels bringing a few extra percent at most. increases the number of subhalos per flux decade, with the largest increase of almost a factor 2 for the largest luminosities. Sub-substructures actually increase the signal in the outskirts of halos (see Figure 4 of Hütten et al. (2016)), the impact of which depends on the instrument angular resolution or containment angle used in the analysis. For instance, in Hütten et al. (2016), no impact was found for dark clumps within the angular resolution of CTA.

Please note that our choice for these constant parameters will lead to a rather conservative number of detectable subhalos and $J$-factor of the brightest ‘resolved’ subhalo—hence conservative limits for DM indirect detection—compared to other choices. Pushing all the parameters to get the most optimistic case would lead to a factor $\sim$2–3 increase of the $J$-factor of the brightest subhalo Hütten et al. (2016).

3.2 The Spatial Distribution ${\mathrm{d}\cal~{}P}/{\mathrm{d}V}$ of Subhalos

We consider four configurations of ${\mathrm{d}\cal P}/{\mathrm{d}V}$ in this work, which are described below and summarized in Tab. 1. The first configuration (model #1) is close to one of our 2016 study Hütten et al. (2016), i.e., based on results from DM-only simulations. It is used as a reference to which the other configurations describing interaction with the baryonic disk are compared to. We consider a spherically symmetric Galactic DM halo and correspondingly, also $\mathrm{d}{\cal P}/\mathrm{d}V$ distribution. The maximum distance of any subhalo from the Galactic center is set to $R_{200}=231.7\,\rm{kpc}$ for all configurations, inspired by the NFW halo from Sect. 3.2.3. We show later that the brightest subhalo is found only with negligible probability at larger Galactocentric radii for all models. Please note that despite the common value for $R_{200}$, the total Galactic DM profile is different from one configuration to another. We however focus here on the clumpy part of the halo only and do not consider the smoothly distributed DM.999We still provide later the total DM profile for the semi-analytical configurations (model #3 and model #4) because it is one of the building blocks of the model. While this article is dedicated to the derivation and comparison of statistical properties of the subhalo population brightness for these different models, we still emphasize that when going to firm predictions and limits, the overall DM profile should matter. Indeed, the Milky Way is a strongly constrained system Piffl et al. (2014); McMillan (2017), which must be taken into account when extrapolating simulation results. This aspect goes beyond the scope of this work though, but is worth mentioning as the Gaia mission is currently boosting our handle on Milky Way dynamics Binney (2017); Gaia Collaboration et al. (2018).

3.2.1 Model #1: DM only (as Implemented in Hütten et al., 2016)

This first configuration uses the position-dependent concentration parametrization of Moliné et al. (2017). The latter is based on the analysis of the DM-only simulations VL-II Diemand et al. (2008) and ELVIS Garrison-Kimmel et al. (2014), and it predicts that subhalos of a given mass are more concentrated close to the galactic center than in the outer parts of their host halos. This effect is related to tidal disruption in the DM potential of the host halo. In the outer parts, when tidal disruption in the DM potential becomes negligible, the concentration is very close to the concentration found for field halos Sánchez-Conde and Prada (2014).101010The difference between using the space-dependent or field halo concentration was found to be a factor $\sim 2$ larger on the brightest subhalo in Hütten et al. (2016). DM-only tidal effects also affect the spatial PDF of subhalos Han et al. (2016), and all recent DM-only simulations found it to be flatter than the DM distribution in the smooth halo. We use an Einasto profile for $\mathrm{d}{\cal P}/\mathrm{d}V$ with $\alpha_{E}=0.68$ and $r_{-2}=199$ kpc, following the results of the Aquarius A-1 halo Springel et al. (2008), and we fix the number of subhalos above $10^{8}~{}\rm M_{\odot}$ to $N_{\rm calib}=300$, as an upper bound to what was found in the Aquarius simulations Springel et al. (2008). Subhalo outer radii and corresponding masses in this model are kept to be cosmological radii and masses, and subhalos are truncated where their mean density reaches $200$ times the critical density of the Universe (see, e.g., Hütten et al. (2019) for details). This model is contained in the parameter space already explored in our previous study Hütten et al. (2016).

3.2.2 Model #2: DM + Galactic Disk Potential (Numerical, Phat-ELVIS)

In the recent Phat-ELVIS simulations, Kelley et al. (2018) have accounted for the effect of the Milky Way baryon potential (including stellar and gas disk, and bulge). They found that subhalos were strongly destroyed in the inner part of the halo, leaving basically no subhalo with mass $m\lesssim 5\times 10^{6}\;{\rm M}_{\odot}$ within the inner $\sim$30 kpc of the host DM halo. This is at odds with the predicted $\mathrm{d}{\cal P}/\mathrm{d}V$ of previous simulation sets (e.g., Aquarius Springel et al. (2008) that was previously used as reference in Hütten et al. (2016)).

Using the subhalo catalog from the simulations provided by the authors, we compute the normalized subhalo PDF per unit volume, $\mathrm{d}{\cal P}/\mathrm{d}V$, at $z=0$, averaged over the 12 Milky-Way-like halos available. The dispersion in each bin provides the error bar. The data are fitted using the following parametrization:

[TABLE]

and results are shown in the left panel of Fig. 1. The first term in Eq. (5) is a Sigmoid function, centered on $r_{0}$ and increasing from 0 to $A$ given $r_{c}$, that allows us to capture the sharp decrease of the number of subhalos in the inner regions of the parent halo. The Sigmoid then transitions to an Einasto profile with characteristic scale $r_{-2}$ and index $\alpha_{E}$. In order to have an Einasto profile close to the DM-only case in the outer parts (see previous paragraph), we fix $\alpha_{E}=0.68$, and the best-fit parameters, obtained on all subhalos in the catalog are $r_{0}=29.2$ kpc, $r_{c}=4.24$, and $r_{-2}=128$ kpc.111111The Milky-Way-like halos in the Phat-ELVIS simulations have masses ranging from $7\times 10^{11}\,\rm M_{\odot}$ to $1.9\times 10^{12}\,\rm M_{\odot}$, and virial radii from 235 kpc to 329 kpc respectively. The fit above was performed on the radial range common to all host halos, namely from 0 to 235 kpc. The value of the parameters remain compatible at the one-sigma level when increasing the radial range to 329 kpc, or when cutting on masses above $5\times 10^{6}~{}\rm M_{\odot}$. The constant $A$ is set to ensure $\mathrm{d}{\cal P}/\mathrm{d}V$ is a PDF.

Finally, among the 12 Milky-Way-like host halos of the Phat-ELVIS simulations, we select halos with masses similar to the NFW halo introduced in the following Sect. 3.2.3, in the range of $1.1\times 10^{12}-1.4\times 10^{12}\;\rm M_{\odot}$; averaging the number of subhalos between $10^{8}$ and $10^{10}\;\rm M_{\odot}$ that have survived interaction with the baryonic potential, we fix in CLUMPY the normalization of the number of subhalos to $N_{\rm calib}=N_{\rm surviving}=90$. In the same way as in the DM-only model # 1, we define and calculate subhalo masses based on their cosmological radii.

3.2.3 Models #3 and #4: DM + Disk Potential (Semi-Analytical, SL17

Complementary to numerical approaches, semi-analytical models have also been considered by several authors, e.g., Green and Goodwin (2007); Goerdt et al. (2007); Berezinsky et al. (2008, 2014); Stref and Lavalle (2017). In this work, we use the study by Stref and Lavalle (2017) (SL17 hereafter) to capture the effects of tidal stripping from the smooth galactic potential of DM and baryons and disk shocking by the baryonic disk. The model in SL17 is built on the dynamically constrained mass models from McMillan (2017), where the latter used kinematic data (including maser observations, the Solar velocity, terminal velocity curves, the vertical force and the mass within large radii) to determine the Milky Way’s DM distribution following the parametrization of Zhao (1996)

[TABLE]

In this work, we only consider the results based on the NFW parametrization ($\alpha,\beta,\gamma=1,3,1$), for which the best-fit parameters are $r_{\rm s}=19.6$ kpc, $\rho_{\rm s}=8.517\times 10^{6}~{}\rm M_{\odot}{}~{}{\rm kpc}^{-3}$, resulting in $R_{200}=231.7$ kpc and $M_{200}=1.31\times 10^{12}\;\rm M_{\odot}$. We checked that our conclusions are left unchanged if considering a cored profile ($\alpha,\beta,\gamma=1,3,0$) instead. In SL17, the initial population of subhalos traces the above Galactic DM halo mass PDF and there is initially no correlation between a subhalo’s mass and its position,

[TABLE]

The addition of tidal interactions modifies this picture because tidal effects select subhalos based on their mass and concentration to produce an ’evolved’ subhalo population. This evolution manifests itself in two aspects: (i) subhalos with a given mass at a given position with too small a concentration are disrupted, while (ii) subhalos that survive are stripped from a large fraction of their mass. Mass stripping is encoded into the subhalo tidal radius $r_{\rm t}$, which is computed as described in (Stref and Lavalle, 2017). Please note that tidal effects in SL17 are computed assuming circular orbits for the clumps. Cosmological simulations show that subhalos are efficiently disrupted by tidal interactions. For instance, Hayashi et al. (2003) find that a subhalo with tidal radius $r_{\rm t}$ and scale radius $r_{\rm s}$ is disrupted if $r_{\rm t}\lesssim 0.77\,r_{\rm s}$. It has however been recently pointed out by van den Bosch et al. van den Bosch et al. (2018); van den Bosch and Ogiya (2018) that disruption within simulations might be largely explained due to a lack of numerical resolution, Poisson noise, or runaway instabilities. According to these authors, subhalos are far more resilient to tidal disruption than numerical simulations tend to show, implying that a subhalo could survive even if $r_{\rm t}\ll r_{\rm s}$ (this result is expected from theoretical grounds Weinberg (1994); Gnedin et al. (1999) and in agreement with earlier findings by Peñarrubia et al. (2010)). The SL17 model includes a free parameter $\epsilon_{\rm t}$ that allows us to simply investigate both possibilities. The $\epsilon_{\rm t}$ parameter is defined such that

[TABLE]

This disruption criterion can also be expressed in terms of the concentration: the subhalo is disrupted if $c<c_{\rm min}$ where $c_{\rm min}$ is referred to as the minimal concentration (which depends a priori on the subhalo’s position and mass). A “cosmological-simulation-like” configuration, where subhalos are efficiently disrupted, corresponds to $\epsilon_{\rm t}\sim 1$. Conversely, a model of very resilient subhalos is obtained by setting $\epsilon_{\rm t}\ll 1$. In the following, we will consider two extreme values : $\epsilon_{\rm t}=0.01$ (model #3) and $\epsilon_{\rm t}=1$ (model #4). The former value implies a subhalo is disrupted when it has lost around 99.99% of its cosmological mass. Please note that the value of $r_{\rm t}$ does depend on the choice of $\epsilon_{\mathrm{t}}$.

The behavior of the SL17 model is illustrated in Fig. 1 (right panel) where we show the spatial distributions of surviving subhalos for different mass decades. We see that the distribution of lighter objects extends to lower radii than the distribution of heavier ones. This is because, as already pointed out in Stref and Lavalle (2017), smaller objects are more concentrated on average and therefore more resilient to tidal disruption. Interestingly, if simulations overestimate tidal disruption as pointed out in van den Bosch et al. van den Bosch et al. (2018); van den Bosch and Ogiya (2018), i.e., $\epsilon_{\rm t}=0.01$ with our parametrization, SL17 predicts a large population of very light subhalos in the innermost regions of the Milky Way.

The number of subhalos in SL17 is calibrated onto the Via Lactea II cosmological simulations Diemand et al. (2008), by the mass fraction in resolved subhalos in these simulations. Since Via Lactea II is a DM-only simulation, the calibration is performed without the baryonic tides and setting $\epsilon_{\rm t}=1$ for consistency. This leads to $N_{\rm calib}$ = 276 for initial cosmological subhalo masses, $m_{200}$, between $10^{8}\,\rm M_{\odot}$ and $10^{10}\,\rm M_{\odot}$ as quoted in Tab. 1. Adding a baryonic potential and considering tidally stripped masses leads to $\sim$110 surviving subhalos in the same mass range for both choices of $\epsilon_{\mathrm{t}}$ (see also Tab. 1).

3.2.4 ${\mathrm{d}\cal P}/{\mathrm{d}V}$ Model Comparison

The four subhalo PDFs considered in this study are compared in Fig. 2. The configurations based on the Phat-ELVIS and the Aquarius simulations are shown in red and blue, respectively, while the SL17 models are shown in yellow ($\epsilon_{\rm t}=0.01$) and green ($\epsilon_{\rm t}=1$). In order to make a meaningful comparison with the simulation results, we only show the SL17 prediction for large subhalo tidal masses ($m>10^{6}\,\rm M_{\odot}$), comparable to the mass of the smallest objects identified in Phat-ELVIS.

The effect of a galactic disk as implemented in the Phat-ELVIS simulation is to disrupt most subhalos in the inner 30 kpc of the galactic halo, as opposed to a DM-only simulation such as Aquarius where subhalos can survive down to much lower radii. The subhalo distribution predicted by SL17, which accounts for the effect of the disk, resembles the one of Phat-ELVIS in that massive subhalos are disrupted at the center. We note that $\mathrm{d}\mathcal{P}/\mathrm{d}V$ peaks at a lower radius in SL17 compared to Phat-ELVIS. Setting $\epsilon_{\rm t}=0.01$ pushes the peak radius to an even lower value with respect to the $\epsilon_{\rm t}=1$ case, as expected since subhalos are then more resilient to disruption.

A more detailed understanding of the difference between these models is beyond the scope of this paper, since the SL17 models are semi-analytically constructed from the kinematically constrained mass model of the Milky Way from McMillan (2017) taking into account the mass dependence of the subhalo spatial distribution, while the DM-only and Phat-ELVIS configurations are based on approximate fits to Milky-Way-like halos from numerical simulations. We still note a similar trend between the SL17 ($\epsilon_{\mathrm{t}}=1$) and Phat-ELVIS configurations, which may indicate that the semi-analytical method associated with the former (pending some simplifying assumptions) somewhat capture the main features of the latter (pending numerical effects likely to dominate in the central regions).

4 Results

Using CLUMPY, we generate fullsky subhalo populations and corresponding $J$- and $D$-factors according to the models from Tab. 1. For all configurations, the distance between the Sun and the Galactic center is set to $R_{\odot}=8.21$ kpc McMillan (2017). We consider two estimations of the $J$-/$D$-factors, integrating either over the full angular extent of a subhalo or up to a radius of $\alpha_{\rm int}=0.5^{\circ}$. Averages and PDFs are then obtained from a statistical sample of 1000 realizations of the subhalo population for each model.

4.1 Subhalo Source Count Distributions

Figure 3 presents the source count distribution for all models and $J$-factors $J>2\times 10^{16}\,\rm{GeV^{2}\,cm^{-5}}$ ($D>2\times 10^{15}\,\rm{GeV\,cm^{-2}}$),121212We checked for convergence in order to obtain a complete ensemble of objects above these thresholds. similar to Figure 3 of our earlier work Hütten et al. (2016). Solid lines show $J$-/$D$-factors within $\alpha_{\rm int}=0.5^{\circ}$ ($\Delta\Omega=2.5\times 10^{-4}\,\rm{sr}$), dashed lines the full signal. For two configurations, we also show the variance bands of the distributions. For the first time, we also show in this work the $D$-factor distribution for decaying DM. Comparing the solid and dashed curves, it can be seen that in the case of annihilation, most of the fainter halos’ emission is contained within the innermost $0.5^{\circ}$ except for the $\sim$100 brightest halos (left panel). This is different in the case of DM decay, for which the emission profile is much more extended (right panel). Comparing the models of this work, we reach the following conclusions for annihilating DM (left panel):

• •

Model #2 (Phat-ELVIS, red lines) predicts about a factor 5 less halos per flux decade than the Aquarius-like DM-only reference model #1. The average brightest halo (within $0.5^{\circ}$) is about a factor 4 fainter than expected for the DM-only case. This drastic decrease of bright objects is both attributed to the fact that the Phat-ELVIS simulations Kelley et al. (2018) find (i) overall less subhalos in Milky Way-like galactic halos ($N_{\rm calib}=90$ vs. $N_{\rm calib}=300$) and (ii), no subhalos are found close to Earth in the innermost 30 kpc of the galactic halo.

• •

Model #3 (SL17, yellow lines) predicts almost the same abundance of bright halos as the DM-only model #1. Model #3 starts from an initial subhalo distribution biased towards the Galactic center following the overall DM distribution, and accounts for tidal subhalo disruption and stripping afterwards according to the semi-analytical model of Stref and Lavalle (2017). With $\epsilon_{\rm t}=0.01$ few subhalos are affected. In turn, the DM-only model #1 already includes a subhalo distribution anti-biased towards the Galactic center in an evolved galactic halo according to the Aquarius simulations (although the considered fitting to the Aquarius simulations Springel et al. (2008) does not account for a mass dependence of the halo depletion).

• •

Model #4 (SL17, green lines) applies a much stronger condition on tidal stripping and total depletion than the model #3 configuration within the semi-analytical approach of Stref and Lavalle (2017). Illustratively, we calculate a total of $1.41\times 10^{6}$ initial subhalos (in the full range between $m_{\rm min}=10^{-6}\,\rm M_{\odot}$ and $m_{\rm max}=1.3\times 10^{10}\,\rm M_{\odot}$) for the subhalo models #3 and #4, out of which 20,000 are completely disrupted for the model #3 ($\epsilon_{\rm t}=10^{-2}$). In contrast, 530,000 halos are disrupted for $\epsilon_{\rm t}=1$ in the model #4.131313Please note that most drawn subhalos are disrupted at masses below a tidal mass of $10^{4}\,\rm M_{\odot}$, the scale above which subhalos are shown in the later Fig. 4. In result, a factor 2 less halos are present above the lower end of the displayed brightness distributions, the ratio increasing for the brightest decades. Recall that surviving halos are truncated at the same tidal radius in models #3 and #4.

For unstable, decaying DM, the flux is proportional to the distance-scaled DM column density (Fig. 3, right panel). Contrary to the case of annihilation, we find:

• •

Changing the signal integration region $\Delta\Omega$ drastically impacts the collected signal, as the emission shows a much broader profile than for annihilation. This loss is most drastic for the brightest halos.

• •

For an integration angle of $\alpha_{\rm int}=0.5^{\circ}$, all models are in remarkable agreement at the brightest end. For fainter flux decades and considering the signal over the full halos extent, models differ by a factor $\sim 5$ (however, with a rather large spread in the $D$-factor PDF of the brightest halo in the individual models, see the later Fig. 5). This suggests that predictions for the largest subhalo flux from decaying DM should be rather model-independent.

For illustrative purposes, Fig. 4 displays subhalo skymaps of a random realization of each of the four models under scrutiny. For each model and to ease comparison, the same DM subhalo sky is used in case of annihilation ($J$-factors, left) or decay ($D$-factors, right). The sky realization varies of course from one model to the other. We do not include the average and smooth DM distributions here. The maps include all subhalos with masses above $10^{4}\,\rm M_{\odot}$ (cosmological masses in the models #1 & #2, tidal masses in the models #3 & #4). For example, in the DM-only case, this corresponds to 1,214,313 halos included in the map, and for a HEALPix resolution of ${N_{\rm side}}=1024$, CLUMPY requires \sim$$30\,\rm CPUh for its computation in case of annihilation (\sim$$20\,\rm CPUh in case of decay). Please note that we did not select the shown random sky realizations to reflect some particular average or extreme case. In Appendix A, we list some properties of the brightest objects in these maps, which can be compared to the average properties derived in the remainder of this section.

4.2 Statistical Properties of the Brightest Halo

Finally, we focus on the statistics linked to the halo with the largest $J$ or $D$-factor (its properties are marked with a “$\star$” symbol in the following). For cold DM particles structuring on small scales, fully dark subhalos represent interesting targets for indirect searches. Conversely, not detecting the brightest of them can be used to set limits. We remark that focusing on the brightest halo alone for setting limits is a simplistic assumption in some circumstances. For example, there are numerous yet unidentified sources detected by the Fermi-LAT which could include a (i.e., the brightest) DM subhalo Berlin and Hooper (2014); Bertoni et al. (2015); Schoonenberg et al. (2016); Mirabal et al. (2016); Bertoni et al. (2016); Wang et al. (2016); Ackermann et al. (2018); The Fermi-LAT collaboration (2019), so constraints set should be correspondingly weaker in this scenario.

The bottom row of Fig. 5 shows the PDFs of $J^{\star}$ and $D^{\star}$, distilled from 1000 realizations of a DM subhalo sky for each model.141414Probability distributions were derived using a kernel density estimation (KDE) with an adaptive Gaussian kernel according to Scott (2015) (except $\mathrm{d}\mathcal{P}/\mathrm{d}(\cos\theta^{\star})$, which was obtained from a histogram). To handle the boundary conditions of $\mathrm{d}\mathcal{P}/\mathrm{d}m(m>m_{\rm{max}}=1.3\times 10^{10}\,\rm M_{\odot})=0$ and $\mathrm{d}\mathcal{P}/\mathrm{d}V(R>R_{200}=231.7\,\rm kpc)=0$, we use the PyQt-Fit KDE implementation by P. Barbier de Reuille, https://pyqt-fit.readthedocs.io (not anymore maintained as of submission of the manuscript), which accounts for a renormalization algorithm at the boundary. Please note that for a precise power-law source count distribution, the  PDF of $J^{\star}$/$D^{\star}$ follows a Fréchet distribution, see App. B of Hütten et al. (2016). From this follows that the brightest expected signal from subhalos has in median a $J$-factor of $\widetilde{J}_{0.5^{\circ}}^{\star}=8.8^{+11}_{-4.0}\times 10^{19}\,\rm{GeV^{2}\,cm^{-5}}$ ($\widetilde{J}_{\rm{full}}^{\star}=1.6^{+2.9}_{-0.9}\times 10^{20}\,\rm{GeV^{2}\,cm^{-5}}$) for the optimistic DM-only model. The signal is expected a factor 7 lower, at $\widetilde{J}_{0.5^{\circ}}^{\star}=1.3^{+0.8}_{-0.4}\times 10^{19}\,\rm{GeV^{2}\,cm^{-5}}$ ($\widetilde{J}_{\rm{full}}^{\star}=3.2^{+2.5}_{-1.3}\times 10^{19}\,\rm{GeV^{2}\,cm^{-5}}$; factor 5 lower) in the case of a largely depleted inner galactic halo (model #4). For decaying DM, all models produce remarkably similar fluxes within $\alpha_{\rm int}=0.5^{\circ}$ (lower left panel of Fig. 5) of $D_{0.5^{\circ}}^{\star}\sim 4^{+3}_{-1}\times 10^{18}\,\rm{GeV\,cm^{-2}}$. Over the full extent of the halo, however, the width of the $D$-factor distributions is much larger, and $D$-factors between $D_{\rm{full}}^{\star}\gtrsim 10^{20}\,\rm{GeV\,cm^{-2}}$ up to $D_{\rm{full}}^{\star}\lesssim 10^{22}\,\rm{GeV\,cm^{-2}}$ are obtained.

The PDFs of the brightest subhalo’s properties may also be derived. The top row of Fig. 5 (left) shows that for annihilating DM, the brightest halo is found at a distance of \sim$$10\,\rm{kpc} from Earth for the models #1 (see also Hütten et al. (2016)) and #3, and also on similar Galactocentric radii (left panel of second row). For the models #2 and #4, which reflect strong tidal disruption of halos in the inner galactic region, bright halos are found at about $\sim$30–40$\,\rm{kpc}$ distance from the Galactic center, and equivalent distances from Earth. For decaying DM, models #2 and #4 give similar predictions, while models #1 and #3 tend to predict the brightest halo at larger Galactocentric and observer distances (right panels of the top rows).

More importantly, the third row of Fig. 5 sheds light on the question whether the brightest DM halo is likely to be a dark halo or a satellite galaxy. For annihilating DM, models #2 and #4 predict subhalos with masses $m^{\star}\gtrsim 10^{8}\,\rm M_{\odot}$ to shine brightest in $\gamma$-rays or $\nu$; these objects are most probably associated with a (dwarf) satellite hosted in their center Strigari (2018). For the DM-only case #1 this is not anymore so obvious, as discussed in Hütten et al. (2016), as lighter objects become probable candidates to provide the highest fluxes. Finally, model #3 predicts very light and close-by halos to shine the brightest: If this model reflects the true nature of DM in our galaxy, the highest $\gamma$-rays or $\nu$ signal from Galactic DM substructure will likely arise from a dark spot in the sky. For decaying DM, the brightest Galactic DM subhalo is likely a satellite galaxy with mass $m^{\star}\gtrsim 10^{8}\,\rm M_{\odot}$.

For blind searches of this brightest halo, it is finally useful to check whether there is a preferred direction in the sky to search for the brightest halo. As all our configurations are symmetrical around the direction towards the Galactic center, we present in the fourth row of Fig. 5 the probability per unit area, $\mathrm{d}\mathcal{P}/\mathrm{d}\cos(\theta^{\star})$, to find the brightest halo at the angular distance $\theta^{\star}$ from the Galactic center (GC). As found in Hütten et al. (2016), in the DM-only case #1, the brightest halo is more probably found in a direction close to the GC. In contrast, model #3 suggests to preferentially search in a ring-like region at \sim$$90^{\circ} distance from the GC; while models #2 and #4 predict the highest probability to find the brightest subhalo towards the galactic anticenter. While these trends also apply for searches for DM decay in subhalos (right panel), they are much more pronounced for signals from annihilating DM.

5 Conclusions

In this work, we have studied the impact of tidal disruption of subhalos by Milky Way’s baryonic potential on the properties of the $\gamma$-ray and $\nu$ signals from subhalos. This effect is mostly encoded in the spatial distribution of subhalos, and four models have been considered. The first model serves as reference and is based on ‘DM-only’ simulations, not including tidal disruption in the baryonic potential. Similar models based on this assumption have been applied by many authors in the past, e.g., Brun et al. (2011); Berlin and Hooper (2014); Lange and Chu (2015); Schoonenberg et al. (2016); Hütten et al. (2016). A second model was obtained by implementing the recent results from the Phat-ELVIS numerical simulation (resolution-limited to a few $10^{6}\,\rm M_{\odot}$) where the inner 30 kpc of the Milky Way are depleted of subhalos. The last two models relied on semi-analytical calculations applicable to the whole subhalo mass range: these calculations find a Galactocentric radius $\lesssim$10 kpc, slightly dependent on the subhalo mass, below which subhalos are stripped of their outer parts or even disrupted, depending on a disruption criterion $\epsilon_{t}$, taken to be to 1 or $10^{-2}$ in this study.

These models lead to significantly different brightness populations of DM annihilation and decay signals. To quantify the difference, we have simulated 1,000 realizations of the subhalo population for each model. Focusing on the brightest subhalo, whose properties can be used to study DM detectability Hütten et al. (2016), we find (for an integration angle of $0.5^{\circ}$) a factor 2 to 7 less signal compared to the case of subhalos in the DM-only configuration, but no significant difference for the decay signal. Our large statistical sample also allowed us to reconstruct the PDF of several properties of this brightest subhalo (mass, distance to the observer, angular distribution in the sky). In particular, the mass information indicates that in models without or little disruption in the disk potential, the brightest subhalo can be close by, and with a mass range below that of known dwarf spheroidal galaxies, i.e., a dark clump. On the other hand, in models with strong tidal disruption, the brightest subhalo is farther away, and its preferred mass is shifted to values similar to those of satellite galaxies, i.e., it could be a known dSph. The latter situation would worsen the prospects of blind galactic dark clump searches with background-dominated instruments that were discussed in Hütten et al. (2016).

In any case, our results highlight the importance of better characterizing the spatial PDF of the subhalo population, in particular by constraining further the level of tidal disruption. Although both numerical and semi-analytical approaches show the same trend in reducing the number and brightness of subhalos, there remain serious quantitative differences. We recall that it is still debated whether or not tidal disruption could be significantly amplified by numerical artifacts in simulations van den Bosch et al. (2018); van den Bosch and Ogiya (2018). On the other hand, present-day semi-analytical methods currently rely on simplifying assumptions, some of which should be relaxed, e.g., include a more realistic distribution of orbits—see complementary studies in Hiroshima et al. (2018); Ando et al. (2019). A further question is related to the possible redistribution of DM inside stripped subhalos, see, e.g., Goerdt et al. (2007); Peñarrubia et al. (2008); Drakos et al. (2017). Until a clearer picture surfaces, all possibilities from weak to strong tidal effects must be equally considered for indirect DM searches. Finally, we stress again that any complete DM halo model comprising a subhalo population should be checked against kinematic constraints, which are increasingly stringent for the Milky Way in the context of the Gaia results Binney (2017); Gaia Collaboration et al. (2018). This is particularly important for tidal disruption as it strongly depends on the detailed description of the baryonic components of the galaxy. Predictions or limits based on galactic halo models which do not account for these constraints should be taken with caution.

All our calculations were performed with the public code CLUMPY, and our results illustrate how this code can quickly be used to incorporate and exploit any progress made by numerical simulations and/or semi-analytical calculations. All computations and drawing of random realizations of the discussed models can be repeated at one’s own account. Also, the subhalo skymaps and catalogs shown here as illustration for the various models are available upon request.

\authorcontributions

Conceptualization, all authors; Formal analysis, M.H., M.S., C.C.; Software, all authors; Data curation, M.H.; Writing—original draft preparation, D.M. and C.C.; Writing—review & editing, all authors.

\funding

This work was supported by the “Investissements d’avenir, Labex ENIGMASS” and by the Max Planck society (MPG). We also acknowledge support from the ANR project GaDaMa (ANR-18-CE31-0006), the OCEVU Labex (ANR-11-LABX-0060), the CNRS IN2P3-Theory/INSU-PNHE-PNCG project “Galactic Dark Matter”, and  European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreements N∘ 690575 and N∘ 674896.

Acknowledgements.

We thank T. Kelley and J. S. Bullock for providing us, prior to the public release, with the Phat-ELVIS subhalo catalog. We also thank M. Doro and M. A. Sánchez-Conde for inviting us on this special issue, and for organizing the Halo Substructure and Dark-Matter Searches workshop held in Madrid in 2018, where discussions prompted this work. Calculations were performed at the Max Planck Computing & Data Facility at Forschungszentrum Garching. We finally thank the anonymous referees for their reviews and the valuable comments which helped to improve the quality of the manuscript. \conflictsofinterestThe authors declare no conflict of interest. The founding sponsors 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. \appendixtitlesyes \appendixsectionsmultiple

Appendix A Properties of the Brightest Subhalos in the Example Maps

In Tab. 2, we list the properties of the brightest subhalos in the random realizations displayed in Fig. 4. Please note that different objects may provide the largest flux for annihilation or decay and for different integration regions $\Delta\Omega$, and the largest values are marked in boldface (except for the DM-only halo, where the very same object is the brightest in all considered scenarios). Further properties of the objects (structural parameters, brightness of lower ranked subhalos, etc.) can be retrieved from the full subhalo catalogs which can be provided to the reader upon request.

\reftitle

References \externalbibliographyyes