Remnants of Galactic Subhalos and Their Impact on Indirect Dark-Matter Searches

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


Introduction
There is overwhelming evidence that most of the matter in the Universe is non-baryonic [1]. An exciting possibility to account for these puzzling observations is that the Universe is filled with a new type of particles which interacts only very weakly with ordinary matter [2,3]. These hypothetical dark matter (DM) particles are being looked for in particle colliders [4][5][6], in direct detection experiments [7][8][9] and in cosmic radiation [10,11], so far without success. The cosmological paradigm is that DM is cold, meaning collisionless and non-relativistic. 1 This implies the structuring of DM on scales smaller than typical galaxies [13,14] which translates into a large population of subhalos within galactic halos [15][16][17]. Modeling these subhalos is crucial if one is to make accurate predictions for direct and indirect DM searches. This is a difficult task as numerical simulations are far from resolving the smallest structures predicted by the cold DM paradigm. To incorporate the smallest structures, one can extrapolate the results of simulations over orders of magnitude in scales, see e.g. [18], but this represents a leap of faith. On the other hand, one can employ semi-analytical models, see e.g. [19][20][21][22]. The difficulty with the latter is to account for the tidal effects experienced by subhalos within the host galaxy. These models can be calibrated on cosmological simulations, which are supposed to consistently describe the tidal stripping of subhalos in their host halo. However, it was recently pointed out by van den Bosch and collaborators [23,24] that simulations are plagued with numerical artifacts that lead to a significant overestimate of the tidal stripping efficiency, and therefore to an underestimate of the actual subhalo population even within the numerical resolution limit. An alternative and complementary way to study the tidal stripping of subhalos is to rely on analytical or semi-analytical methods, which are based on first principles and allow to deal with subhalo mass scales 1 Although the compatibility of the cold DM paradigm with observations on sub-galactic scales is still under debate [12]. arXiv:1905.02008v1 [astro-ph.CO] 6 May 2019 down to the free-streaming scale. Here, we review the semi-analytical model developed by Stref and Lavalle [25] (SL17 hereafter), which incorporates a realistic and kinematically constrained Milky Way mass model (including baryons) and predicts the Galactic subhalo abundance. This model accounts for different sources of tidal effects, and can easily accommodate to different prescriptions for the tidal disruption efficiency.
This paper is structured as follows. In Sec. 2, we briefly review the SL17 model and discuss the resilience of subhalos to tidal effects in light of recent analyses of simulation results [23,24]. In Sec. 3, we compute the DM mass density within subhalos as well as the number density of these objects in the Milky Way. Finally, in Sec. 4, we look at the impact of our results on indirect searches for annihilating DM focusing on gamma rays and cosmic-ray antiprotons.

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

Review of the Stref & Lavalle model
SL17 is a semi-analytical model of Galactic subhalos which is built upon dynamical constraints and cosmological considerations. The main input of the model is the initial subhalo phase-space density where phase space refers to the position-mass-concentration space. The functions dP v /dV, dP m /dm and dP c /dc are the spatial, mass and concentration distributions, respectively. It is assumed that, should subhalos behave as hard spheres (as is the case for single DM "particles" in a cosmological simulation), they would be spatially distributed as dP v /dV ∝ ρ DM where ρ DM is the total DM density profile of the Galaxy. This sets our initial conditions before tidal disruption. The smooth DM mass density is computed through where ρ cl is the average DM mass density inside clumps (this quantity is explicitly computed in Sec. 3).
In the following, we use the Galactic mass models constrained by McMillan [26] on pre-Gaia data for the DM and baryonic mass distributions. In this framework, Eq. (2) ensures the compatibility of our subhalo model with the constrained DM profile ρ DM . The mass m and concentration c refer to the cosmological mass m 200 and concentration c 200 (defined with respect to the critical density) where we have dropped the 200 index for convenience. The subhalo mass function measured in simulations is consistent with a power-law [16,17] where Θ is the Heaviside step function, and the power-law index is α m = 1.9 or α m = 2. These values of α m encompass the Press and Schechter [27] mass function and the Sheth and Tormen [28] mass function, as illustrated in Fig. 1. These functions can be computed directly from the matter power spectrum in the framework of the excursion set theory [29], for the spherical collapse (Press-Schechter) and the ellipsoidal collapse (Sheth-Tormen). Thus the two power-law indices we consider bracket the theoretical uncertainties on the small-scale mass function. If the DM is made of weakly interacting massive particles (WIMPs), the mass cutoff m min can be related to the kinetic decoupling of the DM particle and is found to lie between 10 −4 M and 10 −10 M [19, [30][31][32][33][34][35]. The maximal mass m max is set to 0.01 × M DM where M DM is the total DM mass in the Milky Way. The concentration distribution dP c /dc is generically found to exhibit a log-normal distribution for field halos [36,37], which will define our initial concentration distribution (before tidal stripping). We adopt the peak value and variance fit in Sánchez-Conde and Prada [38], which was shown to provide a good description of cosmological simulations run independently by several groups. Subhalos are assumed to have a Navarro-Frenk-White (NFW) profile [39] with parameters set by m and c (the impact of choosing an Einasto profile [40] instead of NFW was investigated in [25]).  The subhalo population is strongly affected by tidal interactions with the potential of the host galaxy [42]. This is accounted for in the model through the calculation of a tidal radius r t for each subhalo. The tidal radius should be interpreted as the physical extension of a subhalo, which is in general smaller than the extension it would have on a flat background. The physical mass of a subhalo is then where ρ sub is the subhalo inner mass density profile. In our modeling, tidal stripping only removes the outer layers of subhalos while leaving the inner parts unchanged. In reality, DM should rearrange itself into a new equilibrium state. The central density however, should be left essentially unchanged [43]. Since we are interested in indirect searches for annihilating DM, and the central density gives the dominant contribution to the annihilation rate, our modeling should lead to a reasonable approximation. Two important contributions are accounted for in SL17: the effect of the smooth Galactic potential (including both DM and baryons), and the gravitational shocking induced by the baryonic disk [44,45]. The latter effect turns out to be very efficient at stripping subhalos in the inner 20 kpc of the Galaxy, a result also found in numerical studies [46][47][48]. The strength of SL17 over simulations is that it accounts for the constrained potential of the MW, with a detailed description of baryons.

Subhalo disruption?
Whether a subhalo can be completely disrupted by tidal effects is an open question. A number of numerical studies have found that a subhalo is completely disrupted when the total energy gained through tidal-stripping or disk-shocking effects is comparable to the binding energy [46,49]. On the other hand, some studies [50][51][52] have found that cuspy subhalos almost always survive mass loss, leaving a small bound remnant behind even after gaining an energy far greater than their binding energy. These contradictory results may have been reconciled in a recent series of papers by van den Bosch and collaborators [23,24,53]. In these studies, it is shown that subhalo disruption in N-body simulations can actually be entirely explained by numerical artifacts. In particular, disruption is shown to be highly sensitive to the value of the force-softening length. If this length is taken sufficiently small, the authors show that subhalos survive tidal mass loss in the form of a small bound remnant. We aim at quantifying the impact of these results on the whole subhalo population. Tidal disruption is modeled in a very simple way in SL17: given a subhalo with scale radius r s and tidal radius r t , we assume In Eq. (5), t is a dimensionless free parameter assumed universal i.e. independent of the subhalo's mass, concentration or position. In SL17, the value of the disruption parameter was set to t = 1 in agreement with numerical results, see e.g. [49]. The results of van den Bosch and collaborators point toward a much lower value for t . In this work, we consider two extreme values: t = 1 and t = 0.01. The latter means a subhalo is disrupted when it has lost around 99.99% of its mass. In the following, we refer to these two configurations as "fragile subhalos" ( t = 1) and "resilient subhalos" ( t = 0.01). The final subhalo phase-space density can now be written where N tot is the total number of substructures within the virial radius of the Milky Way and K tot is a normalization factor:

Calibration procedure
In its current version, the SL17 model requires a calibration of the subhalo abundance in a given mass range (this will change in future versions). To be consistent with results from the highest-resolution simulations available, the calibration is done by demanding that the subhalo mass fraction is similar to what is found in the dark matter-only Via Lactea II simulation [16]. This amounts to 11% of the total dark halo mass in the form of subhalos in the virial mass range [m 1 , where M DM is the total DM mass of the Galaxy. We stress that these numbers are expressed in terms of virial masses, not tidal masses (see [25] for further details). To reproduce the (likely overestimated) tidal disruption efficiency in simulations, we have to set t = 1 at the calibration stage, and we also neglect the impact of baryons. The disruption efficiency parameter t can safely be changed after the calibration has been completed. It is much safer to perform this calibration on dark matter-only simulations because the tidal stripping induced by baryons strongly depends on the details of the stellar distribution, which is acutely constrained in the Milky Way. More formally, the normalization procedure reads: Fixing f sub (m 1 , m 2 ) = 0.11 leads to the total number of clumps N DMO, t =1 in the simulation-like configuration. Note that this value assumes that m is really m 200 in the equation above, not the tidal mass. Now that the model is properly calibrated, we incorporate all the effects that are not included in the calibration i.e. the tidal effects due to the baryons and possibly t < 1. This is done by assuming that subhalos in the outskirts of the Galaxy are not affected by baryonic tides or the value of t . This is motivated by the observation that tidal effects are inefficient far from the center of the Galaxy, and subhalos almost behave like isolated halos. The DM mass within clumps per unit of volume can be expressed as where m t is the tidal subhalo mass introduced in Eq. (4). Equating ρ cl (r 200 ), where r 200 is the virial radius of the Galaxy, in the DM-only+ t = 1 configuration, to the same quantity in the realistic configuration (including baryons and t ≤ 1) leads to the simple relation The two normalization factors K can be computed using Eq. (7) and we get the value of N tot . The number of subhalos within the Solar radius r = 8.21 kpc is shown in Tab. 1. This number is highly sensitive to the parameters of the mass function α m and m min as already shown in [25]. Furthermore, it is quite sensitive to t : going from t = 1 to t = 0.01, the number of subhalos increases by at least an order of magnitude. The impact of t on the subhalo mass and number density is investigated in the next section.

Mass and number densities of subhalos
In this section, we compute the mass density within subhalos as well as the subhalo number density. The subhalo mass density is defined in Eq. (9). Once the subhalo density is known, it is used to determine the amount of DM smoothly distributed across the Galaxy through Eq. (2). The DM mass inside subhalos is compared to the total DM density ρ DM in Fig. 2. The mass density in the form of subhalos is predicted to be much higher, by orders of magnitude, for resilient subhalos than for fragile subhalos. The former case is also more justified theoretically, although the latter one allows us to compare with very conservative assumptions. At the position of the Solar System, the impact is around one order of magnitude. Although these are large differences, we note the subhalo mass density is still far below the total DM density. This means that most of the DM mass within the orbit of the Sun is smoothly distributed rather than clumpy, irrespective of the efficiency of the tidal disruption set by t . The subhalo number density in SL17 can be formally written dN dV ( r) = dm dc dN dV dm dc (11) The results obtained are shown in Fig. 3. The number density, just like the mass density, is highly sensitive to α m and m min , as well as the disruption parameter. Interestingly, the values we get in the Solar neighborhood are comparable to the local number density of stars n * ∼ 1 pc −3 . For a low value of the minimal mass m min , the subhalo number density can even be much higher, possibly going as high as 10 5 pc −3 . This could have a number of interesting implications for the interactions between subhalos and stars. The tidal heating of subhalos by stars has been investigated in a number of studies [19,50,[54][55][56][57][58], with different conclusions. In the next section, we look at the impact of our results on indirect DM searches.

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

Annihilation profiles and local boost factors
The number of self-annihilation of DM particles at position r is proportional to ρ 2 ( r) where ρ is the DM mass density. If subhalos are discarded, the Galactic annihilation profile is Let us now consistently include the contribution of subhalos. The luminosity of a single clump is where V t ( r, m, c) is the volume of the clump within its tidal radius. The annihilation of the full subhalo population is simply obtained by integrating the luminosity of a single object over the subhalo phase-space number density The full annihilation profile must also incorporate the annihilation in the smooth halo (different from L 0 which is the density assuming all the DM is smoothly distributed), as well as the annihilation of subhalo particles onto smooth halo particles. The first contribution can be written where ρ sm ( r) = ρ DM ( r) − ρ cl ( r) is the smooth DM density. The clump-smooth contribution is The total annihilation profile is simply the sum of all contributions This should be compared to L 0 ( r) to evaluate the impact of clustering on the annihilation rate. This is usually done in terms of a boost factor which we define as This is not quite the boost factor used in indirect searches, which is defined through a ratio of fluxes (see Eq. (23) and Eq. 29)). The boost in Eq. (20) is rather the local increase in the annihilation rate due to clustering. According to this definition, the boost is zero if L = L 0 i.e. substructures are not included. 2 The annihilation profiles are shown on the top panels in Fig. 4, and the associated boost factors are shown on the bottom panels. As already shown a long time ago, see e.g. [19,60], the boost is an increasing function of the galactocentric radius r = | r|. This is due to the morphology of the annihilation profiles which is modified by the inclusion of clumps: we have L cl ∝ ρ DM while L 0 ∝ ρ 2 DM . Also noticeable is the high sensitivity of the annihilation profile to the mass function index α m which is by far the largest source of uncertainty on the clump contribution. The value of the disruption parameter t has almost no impact on the boost above 20 kpc due to the ineffectiveness of tidal effects far from the center of the Galaxy. Below 20 kpc, the boost is strongly sensitive to the disruption parameter. In the inner few kiloparsecs, fixing t = 0.01 leads to a boost orders of magnitude larger than in the t = 1 configuration. We have however B( r) 1 below 3 kpc regardless of the value of t , meaning L L 0 and subhalos do not have any impact on the annihilation rate in that region. The region where the impact of t on the annihilation profile is the most important is located between 3 kpc and 10 kpc. This region coincidentally includes the Solar system, located at r 8 kpc. This motivates a more detailed investigation of two standard annihilation channels: gamma rays and cosmic-ray antiprotons, which are sensitive to different annihilation regions.

Application to gamma rays
The energy-differential flux of gamma rays originating from DM self-annihilation is, on a given line of sight where σv is the thermally-averaged annihilation cross-section, m χ is the DM mass, dN γ /dE is the gamma-ray spectrum at annihilation and s the distance coordinate along the line of sight. 3 If the 2 This differs by one unit from the definition used in [25]. 3 The expression of dΦ γ /dE should be multiplied by 1/2 if the DM particle is not its own antiparticle.  annihilation cross-section is velocity-independent, astrophysical ingredients enter only through the J-factor defined as where ψ is the angle between the direction of the Galactic center and the line of sight (spherical symmetry of the dark halo is assumed). The impact of small-scale clustering on this J-factor has been considered in a number of studies [22,38,[60][61][62][63][64][65][66][67][68]. We define the gamma-ray boost factor as where J 0 = ds L 0 is the J-factor without subhalos. Unlike the local boost B in Eq. (20), the gamma-ray boost B γ depends on the line of sight rather than the position in the Galaxy [60]. The boost is shown as a function of ψ in Fig. 5. The growth of the local boost B( r) as a function of r translates into a growth of B γ (ψ) as a function of ψ i.e. the maximal gamma-ray boost is reached at the anti-center. This maximal boost ranges from 0.5 to 9, depending on the values of α m and m min . The survival of clumps noticeably increases the boost at all latitudes. The gain is greater at small latitudes where substructures are more impacted by tidal effects. Below ψ 40 deg, the boost is increased by a factor of at least two in all configurations. This should have important consequences for indirect searches using gamma rays, especially at high latitudes. Interestingly, high latitudes have been shown to be a very sensitive probe of DM annihilation even without the inclusion of clumps [69].

Application to cosmic-ray antiprotons
Charged cosmic rays constitute an indirect detection channel complementary to gamma rays [10,70]. Since their original proposal as a probe of DM annihilation [71], cosmic-ray antiprotons have been shown to be especially sensitive, see e.g. [72][73][74][75][76][77][78]. Antiprotons have been the subject of much scrutiny since the latest measurement of the antiproton flux performed by the AMS-2 collaboration [79]. A number of studies [80][81][82][83][84][85] have found a discrepancy between the measured flux and a purely secondary origin of antiprotons. This discrepancy could be interpreted as evidence for annihilating DM, although the significance of the excess is debated as it depends on the propagation model used and the modeling of systematic uncertainties. In this context, it is worth evaluating systematic uncertainties coming from the small-scale structuring.
Since antiprotons have a random motion due to their diffusion on the inhomogeneities of the magnetic halo, their detection gives little information on their source. This implies that the antiproton boost factor, and the boost for charged cosmic rays in general, is not direction-dependent unlike for gamma rays. Instead, this boost is energy-dependent [86] and has been shown to be mild, at most a factor of two [21,64]. Although smaller than the gamma-ray boost, this can still be larger than the systematic uncertainties on cosmic-ray propagation. This motivates a new computation of the boost, which we perform here. The antiproton boost factor is defined as where T is the antiproton kinetic energy, Φ p is the DM-induced antiproton flux including the subhalo contribution and Φ p,0 the same flux assuming all the DM is smoothly distributed. To get the flux, one must solve the cosmic-ray steady-state propagation equation which accounts for spatial diffusion, convection, energy losses, diffusive reacceleration, and spallation processes in the disk (taken as infinitely thin). In Eq. (25), Ψ is the antiproton number density per unit energy which is related to the flux through Φ p = v p /(4π) × Ψ where v p is the antiproton speed. Antiprotons are sourced by DM annihilation where dN p /dE is the antiproton spectrum at annihilation. 4 Several unknown propagation parameters enter Eq. (25). These can be constrained using the measured boron to carbon ratio (B/C) [87]. We use the best-fit model derived by Reinert and Winkler [82], which includes an energy-break in the diffusion coefficient. The B/C ratio can only constrain K 0 /L where K 0 is the normalization of the diffusion coefficient and L is the half-height of the magnetic halo. As shown in [73], the DM-induced antiproton flux depends crucially on L hence we consider two extremal values in this work. A lower bound on L can be obtained from low-energy positron data [88] and the authors of [82] find L = 4.1 kpc. For the largest value, we take L = 15 kpc. According to the analysis of [82], the B/C data are consistent with negligible reacceleration. Furthermore we neglect energy losses which are unimportant for high-energy antiprotons. The resulting transport equation can be solved semi-analytically using the Green's function formalism (see [21] for the solution) and the differential flux can be written Since all the energy-dependent terms have been neglected in Eq. (25), the energy part of the Green's function is trivial and the boost factor can be simply written The boost factor is shown as a function of the antiproton kinetic energy in Fig. 6. We first note that the boost is roughly energy-independent. This is because antiprotons probe the entire volume of the magnetic halo during their lifetime, independently of their energy. This is not true at low energies, below a few GeVs, where energy losses become relevant. The half-height of the magnetic halo has a small impact on the boost, with L = 15 kpc leading to a slightly larger B p than L = 4.1 kpc. The main source of uncertainties are coming from the subhalo parameters α m , m min and t . The survival parameter t has a significant impact on the result, with a small value t = 0.01 leading to a boost roughly twice as large as in the t = 1 case. As for gamma rays, the most critical parameter is α m . For a low value α m = 1.9, the boost never exceeds 10% while it is always higher for α m = 2. Overall, playing with the propagation and subhalo parameters, we find that the antiproton boost can conservatively range from 2% to 140%. These values are in agreement with earlier results. Although it is conservative to ignore the small-scale clustering when deriving limits on the annihilation cross-section using data, the boost should be included when interpreting an excess as a signature of DM annihilation. Indeed, a factor of two in the DM contribution would change the inferred mass and cross-section of the hypothetical DM particle.

Conclusion
Subhalos suffer mass loss due to their interaction with the tidal field of the Galaxy, which makes their modeling very challenging. Consequently, most subhalo models rely at least partly on numerical simulations to calibrate their predictions. However, it was recently shown that numerical simulations might not properly account for the tidal disruption of subhalos, as artificial effects lead to a serious overestimation of the efficiency of these processes [23,24]. We note that the resistance of subhalos to tidal stripping is further supported by theoretical arguments, like adiabatic invariance that should prevail in their inner parts [45], as already emphasized in [25]. We have derived some of the consequences of these results using the semi-analytical Galactic subhalo population model of Stref and Lavalle [25], assuming a tidal disruption efficiency t = 0.01. We predict the spatial dependence of the subhalo properties due to tides induced both by the global gravitational potential and baryonic disk shocking. We remind the reader that this model is built from constrained mass models for the Milky Way and is therefore consistent with current kinematic constraints, which is usually not the case in extrapolations from "Milky-Way-like" simulations. We find that the local mass density is still dominated by the smooth component of the dark halo. The local number density of subhalos is increased by roughly one order of magnitude with respect to estimates based on a tidal disruption efficiency similar to that inferred from simulations ( t = 1). This makes the subhalo number density comparable, for a broad range of minimal subhalo mass m min , to the local star number density. The resilience of subhalos increases the local DM annihilation rate, which in turn affects the predictions for indirect searches. For gamma rays, we find that the boost factor is increased by at least a factor of two for ψ < 40 deg, and slightly less for higher values of ψ. The boost factor for antiprotons is also increased by a rough factor of two if subhalos are resilient to tidal disruption. For a complemetary study comparing the SL17 model to simulation results regarding indirect searches in gamma rays and neutrinos, we refer the reader to [89].
For future work, we plan on including a more detailed mass function directly deriving from the primordial power spectrum, as well as the tidal heating of subhalos due to individual stars and study the consequences of having a large population of small objects in the Solar neighborhood.