γ-ray and ν Searches for Dark-Matter Subhalos in the Milky Way with a Baryonic Potential

The distribution of dark-matter (DM) subhalos in our galaxy remains disputed, leading to varying γ-ray and ν 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 γ-ray and ν sky. Compared to predictions based on DM only, we conclude a decrease of the flux of the brightest subhalo by a factor of 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.


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 γ-rays [1]. 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 M , and extrapolated down ten orders of magnitude to the model-dependent minimal masses of subhalos [2,3]. Functional parametrizations were incorporated in the CLUMPY code [4][5][6] to generate γ-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 γ-ray observatory [7] 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 [8][9][10] and semi-analytical studies [11][12][13][14][15] 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 [1], where the brightest subhalo was found, on average, at ∼ 10 kpc from the Galactic center and at similar distance from Earth. Experimental limits on DM from Galactic subhalos, derived from Fermi-LAT [16][17][18][19][20][21][22][23][24] or expected from prolonged operation [25], from HAWC [26], or future instruments like GAMMA 400 [27] and CTA [1,28], should also be impacted by this result.
The paper is organized as follows: Sect. 2 presents the overall methodology and recalls how all relevant subhalos are efficiently accounted for with CLUMPY. Sect. 3 lists the subhalo critical parameters, highlighting the very different spatial distributions considered in this analysis. Sect. 4 presents an 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 (via the so-called J-factors), or decaying DM (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 Sect. 5.

Important quantities and methodology
2.1. γ-ray and ν fluxes from dark matter The γ-ray or ν flux from annihilating/decaying DM particles, at energy E, along the line-of-sight (l.o.s.) in the direction (ψ, θ), and integrated over the solid angle ∆Ω = 2π (1 − cos α int ), is given by where m χ is the mass of the DM particle, σv is the velocity-averaged annihilation cross-section 1 (resp. τ DM is the DM particle lifetime), dN f γ /dE and B f correspond to the spectrum and branching ratio of annihilation or decay channel f , l is the distance from the observer, and ρ( r) is the Galactic DM distribution. The latter can be cast formally as the sum of a smooth distribution ρ sm of unclustered DM particles, and a collection of i = 1 . . . N subs subhalos, each with a density ρ i ( r). The astrophysical term for annihilating DM 2 then reads, [4,5,29] dJ tot dΩ = l.o.s.
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. [30]) or extracted from numerical simulations, as discussed below.

Generating skymaps with CLUMPY v3.0
For all our calculations, we rely on the public CLUMPY code described in [4][5][6]. It is a flexible code that efficiently emulates the end-product of numerical simulations in terms of γ-ray and neutrino signals for DM annihilation or decay. It allows to easily explore how results are affected when changing the properties of the DM halos. CLUMPY v2.0 [5] was used to this purpose, to estimate the sensitivity of the CTA [1] and HAWC [26] γ-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 γ-rays or ν in many targets: dwarf spheroidal galaxies [31][32][33][34][35][36][37][38][39][40], the Galactic halo [41][42][43], the Smith HI cloud [44], nearby galaxies [45], galaxy clusters [46][47][48], and also for the extragalactic diffuse emission [49]. 1 We assume here the DM particle to be a Majorana particle, so that δ = 2 (for a Dirac, δ = 4). 2 For conciseness, we present in Eqs. (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. The present analysis is performed with CLUMPY v3.0 [6]. 3 For completeness, we recap below the main steps of the CLUMPY calculation used for this work: • CLUMPY and the particle physics term: Eq. (1) shows that the particle physics term and the astrophysical terms are decoupled. 4 As the flux depends on the specific DM candidate chosen, we provide results in terms of Jand D-factors only; CLUMPY can easily be used to transform those into γ-ray or ν fluxes for any user-defined DM candidate (see CLUMPY's online documentation 5 ). • CLUMPY and the astrophysics term: to calculate skymaps of dJ/dΩ, one should rely in principle on Eq.
(2). However, this is impractical in terms of computing time, as ∼ 10 14 subhalos are expected in a Milky Way-sized DM halo. This problem can be overcome by formally separating Eq. (2) in an average and "resolved" component, With this ansatz, only a limited number N drawn subs 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 [51] or their subhalo catalogs [52]. 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, J tot ± σ J subs , and we refer to [4] 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 [1]), this approach allows us to pre-select halos likely to shine bright at Earth and to consider all decades down to the smallest subhalo masses in the calculation.

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, dP /dV, (ii) the mass-concentration 6 relationship, c(m, r), and (iii) the calibration of the total number of subhalos in the Milky Way, N calib . The latter number is determined from numerical simulations, in a range where subhalos are resolved. Here, N calib is defined for the mass range 10 8 − 10 10 M , and it typically falls in the range 100 − 300. 7 3 When releasing CLUMPY v3.0, we corrected a misprint that was present in v2.0, related to our implementation of the virial overdensity from [50]. This issue was responsible for obtaining in [1] about a factor 3 more subhalos than expected per flux decade (see full details in the CLUMPY documentation). We recall that in [1], we found that galactic variance is responsible for a factor 10 uncertainty on the value of the brightest subhalo, and that other subhalo properties were responsible for another factor ∼ 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. 4 Strictly speaking, this factorization holds true only for DM candidates for which σv is independent of the velocity and consideration of small redshift cells, ∆z/z 1. 5 https://lpsc.in2p3.fr/clumpy 6 The concentration is defined to be c = r ∆ /r −2 , with r ∆ , taken to be the subhalo boundary, is the radius at which the mean subhalo density is ∆ times the critical density (see, e.g., [6]), and r −2 is the position in the subhalo for which the slope of the density is -2. We use ∆ = 200 in this work. 7 This range was recently shown to be in agreement with the observed number of dwarf spheroidal galaxies SDSS corrected by the detection efficiency [53], alleviating the tension caused by the so-called missing satellite problem in CDM scenarios [54,55]. Given the minimal mass of subhalos, m min , N calib can be used to calculate the mass fraction, f 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, Here, dP /dm and dP /dc describe the PDFs for a subhalo to have a given mass and a given concentration c. In reality, the factorization in Eq. (4) may break down when subhalos gravitationally interact with the DM and baryonic potentials of their host halo [15,56], entangling their mass and positional distributions. We will consider this effect by "evolving" the distribution of Eq. (4) in the model presented in Sect. 3.2.3.

Fixed subhalo-related quantities
This work focuses on the impact of a baryonic disk potential on the subhalo population, mainly through the spatial PDF dP /dV (see Sect. 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 [1] and only provide a brief summary below: • Index α m of the power-law subhalo mass PDF dP /dm ∝ m −α m and subhalo mass range: We choose α m = 1.9, m min = 10 −6 M , and m max = 1.3 × 10 10 M . The maximum clump mass for all models is set to 10 −2 × M 200 of the NFW halo from Sect. 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 α m mostly affect the diffuse emission boost from unresolved halos. For a fixed normalization N calib , a steeper mass function (α m = 2.0) decreases the number of bright halos (J 10 20 GeV 2 cm −3 ) by not more than ∼ 30%.
• Width of dP /dc: We set σ c = 0.14 [57]. Using a larger scatter σ c = 0.24 [58] 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 [59].
Using an Einasto profile [59,60] instead amounts to a global increase ∼ 2 of the number of subhalos per flux decade within the considered integration regions ∆Ω. Note that micro-halos with m M may show steeper inner slopes [61][62][63], however, we have found that these micro-halos do not provide new bright, resolved subhalo candidates [1].
• Level of sub-substructures: We do not consider an emission boost from substructure within subhalos. Such a boost from additional levels of substructure 8 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 Fig. 4 of [1]), the impact of which depends on the instrument angular resolution or containment angle used in the analysis. For instance, in [1], no impact was found for dark clumps within the angular resolution of CTA.
Note that our choice for the 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 ∼ 2 − 3 increase of the J-factor of the brightest subhalo [1]. 8 As shown in [5] (see their Fig. 1), only the first level of substructure significantly boosts the halo luminosity, the next levels bringing a few extra percent at most.  Properties of the initial subhalos from which the surviving ones are obtained after interaction with the baryonic potential.

The spatial distribution dP /dV of subhalos
We consider four configurations of dP /dV 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 [1], 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 dP /dV distribution. The maximum distance of any subhalo from the Galactic center is set to R 200 = 231.7 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. 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. 9 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 dark matter profile should matter. Indeed, the Milky Way is a strongly constrained system [64,65], which has to 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 [66,67] This first configuration uses the position-dependent concentration parametrization of Moliné et al. [69]. The latter is based on the analysis of the DM-only simulations VL-II [70] and ELVIS [71], 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 [57]. 10 DM-only tidal effects also affect the spatial PDF of subhalos [56], 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 dP /dV with α E = 0.68 and r −2 = 199 kpc, following the results of the Aquarius A-1 halo [68], and we fix the number of subhalos above 10 8 M to N calib = 300, as an upper bound to what was found in the Aquarius simulations [68]. 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., [6] for details).
3.2.2. Model #2: DM + Galactic disk potential (numerical, Phat-ELVIS [10]) In the recent Phat-ELVIS simulations, Kelley et al. [10] 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 5 × 10 6 M within the inner ∼ 30 kpc of the host DM halo. This is at odds with the predicted dP /dV of previous simulation sets (e.g., Aquarius [68] that was previously used as reference in [1]).
Using the subhalo catalog from the simulations provided by the authors, we compute the normalized subhalo PDF per unit volume, dP /dV, at z = 0, averaged over the twelve Milky-Way-like halos available. The dispersion in each bin provides the error bar. The data are fitted using the following parametrization: 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 α E . In order to have an Einasto profile close to the DM-only case in the outer parts (see previous paragraph), we fix α 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. 11 The constant A is set to ensure dP /dV is a PDF. Finally, among the twelve Milky-Way-like host halos of the Phat-ELVIS simulations, we select halos with masses similar to the SL17 NFW halo (see below), in the 1.1 × 10 12 − 1.4 × 10 12 M range; averaging the number of subhalos between 10 8 and 10 10 M that have survived interaction with the baryonic potential, we fix in CLUMPY the normalization of the number of subhalos to N calib = N 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 [15]) Complementary to numerical approaches, semi-analytical models have also been considered by several authors, e.g., [11][12][13][14][15]. In this work, we use the study by Stref and Lavalle [15] (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 [65], 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 [72]  In this work, we only consider the results based on the NFW parametrization (α, β, γ = 1, 3, 1), for which the best-fit parameters are r s = 19.6 kpc, ρ s = 8.517 × 10 6 M kpc −3 , resulting in R 200 = 231.7 kpc and M 200 = 1.31 × 10 12 M . We checked that our conclusions are left unchanged if considering a cored profile (α, β, γ = 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, dP dV (r, initial) ∝ ρ tot (r) .
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 having 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 t , which is computed as described in [15]. 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. [73] find that a subhalo with tidal radius r t and scale radius r s is disrupted if r t 0.77 r s . It has however been recently pointed out by van den Bosch et al. [74,75] 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 t r s (this result is expected from theoretical grounds [76,77] and in agreement with earlier findings by Peñarrubia et al. [78]). The SL17 model includes a free parameter t that allows us to simply investigate both possibilities. The t parameter is defined such that r t r s < t ⇔ subhalo is disrupted. The four spatial PDFs of subhalos considered in this work: SL17 models for subhalos with m > 10 6 M (yellow and green), PDF based on the Phat-ELVIS simulation (red) [10], and on the Aquarius simulation (blue) [68]. To highlight the behavior at low radii where tidal effects are the most relevant, the curves are shifted to match the value of dP /dV(231.7 kpc) in the Phat-ELVIS configuration.
This disruption criterion can also be expressed in terms of the concentration: the subhalo is disrupted if c < c min where c 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 t ∼ 1. Conversely, a model of very resilient subhalos is obtained by setting t 1. In the following, we will consider two extreme values : t = 0.01 (model #3) and t = 1 (model #4). The former value implies a subhalo is disrupted when it has lost around 99.99% of its cosmological mass. Note that the value of r t does depend on the choice of 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 [15], 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. [74,75], i.e. 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 [70], 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 t = 1 for consistency. This leads to N calib =276 for initial cosmological subhalo masses, m 200 , between 10 8 M and 10 10 M as quoted in Tab. 1. Adding a baryonic potential and considering tidally stripped masses leads to ∼ 110 surviving subhalos in the same mass range for both choices of t (see also Tab. 1).

dP /dV 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 ( t = 0.01) and green ( 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 M ), 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 like 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 dP /dV peaks at a lower radius in SL17 compared to Phat-ELVIS. Setting t = 0.01 pushes the peak radius to an even lower value with respect to the 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 [65] 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 ( 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).

Results
Using CLUMPY, we generate fullsky subhalo populations and corresponding Jand 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 = 8.21 kpc [65]. 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 α int = 0.5 • . Averages and PDFs are then obtained from a statistical sample of 1,000 realizations of the subhalo population for each model. • 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 • ) 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 [10] find (i) overall less subhalos in Milky Way-like galactic halos (N calib = 90 vs. N 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 [15]. With 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 [68] does not account for a mass-dependence of the halo depletion).

Subhalo source count distributions
• 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 [15]. Illustratively, we calculate a total of 1.41 × 10 6 initial subhalos (in the full range between m min = 10 −6 M and m max = 1.3 × 10 10 M ) for the subhalo models #3 and #4, out of which 20,000 are completely disrupted for the model #3 ( t = 10 −2 ). In contrast, 530,000 halos are disrupted for t = 1 in the model #4. 13 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 ∆Ω 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 α int = 0.5 • , 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 ∼ 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.  drawn halos); model #2 emulating the Phat-ELVIS simulations [10] (364,064 drawn halos); and the semi-analytical models #3 (subhalos more resilient against tidal disruption, 549,572 surviving halos) and #4 (less subhalos surviving tidal destruction, 546,096 surviving halos) according to SL17 [15]. The displayed maps (50MB in file size) include all subhalos above a mass of 10 4 M (stripped mass for models #3 & #4), and can be provided in the fits format upon request.

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 " " symbol in the following). For cold dark matter 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 [16][17][18][19][20][21]79,80], so constraints set should be correspondingly weaker in this scenario.
The bottom row of Fig. 5 shows the PDFs of J and D , distilled from 1,000 realizations of a DM subhalo sky for each model. 14 From this follows that the brightest expected signal from subhalos has in median a J-factor of J 0.5 • = 8.8 +11 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 ∼ 10 kpc from Earth for the models #1 (see also [1]) 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 ∼ 30 − 40 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 10 8 M to shine brightest in γ-rays or ν; these objects are most probably associated with a (dwarf) satellite hosted in their center [82]. For the DM-only case #1 this is not anymore so obvious, as discussed in [1], 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 γ-rays or ν 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 10 8 M . 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, dP /d cos(θ ), to find the brightest halo at the angular distance θ from the Galactic center (GC). As found in [1], 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 ∼ 90 • 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.

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 γ-ray and ν 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. [1,16,18,28,52]. A second model was obtained by implementing the recent results from the Phat-ELVIS numerical simulation (resolution-limited to a few 10 6 M ) 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 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 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 [1], we find (for an integration angle of 0.5 • ) 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 [1].
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 [74,75]. 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 [83,84]. A further question is related to the possible redistribution of DM inside stripped subhalos, see, e.g., [12,85,86]. Until a clearer picture surfaces, all possibilities from weak to strong tidal effects have to 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 more and more stringent for the Milky Way in the context of the Gaia results [66,67]. 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. The subhalo skymaps shown as illustration for the various models are available upon request.