# Classification of gamma-ray targets for velocity-dependent and subhalo-boosted dark-matter annihilation

Thomas Lacroix, Gaétan Facchinetti, Judit Pérez-Romero, Martin Stref, Julien Lavalle, David Maurin, Miguel A. Sánchez-Conde

► To cite this version:

Thomas Lacroix, Gaétan Facchinetti, Judit Pérez-Romero, Martin Stref, Julien Lavalle, et al.. Classification of gamma-ray targets for velocity-dependent and subhalo-boosted dark-matter annihilation. Journal of Cosmology and Astroparticle Physics, 2022, 10, pp.021. 10.1088/1475-7516/2022/10/021 hal-03625365

HAL Id: hal-03625365

https://hal.science/hal-03625365v1

Submitted on 20 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Classification of gamma-ray targets for
velocity-dependent and
subhalo-boosted dark-matter
annihilation

Thomas Lacroix,ª Gaétan Facchinetti,ᵇ Judit Pérez-Romero,ª
Martin Stref,ᶜ Julien Lavalle,ᵈ David Maurin,ᵉ and Miguel
A. Sánchez-Condeª

<sup>a</sup>Instituto de Física Teórica (IFT) UAM-CSIC, Departamento de Física Teórica (DFT), Universidad Autónoma de Madrid, Calle Nicolás Cabrera, 13-15, 28049 Madrid, Spain

<sup>b</sup>Service de Physique Théorique, CP225, Université Libre de Bruxelles, Boulevard du Triomphe, B-1050 Brussels, Belgium

<sup>c</sup>Laboratoire d'Annecy-le-Vieux de Physique Théorique (LAPTh), Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, F-74000 Annecy, France

<sup>d</sup>Laboratoire Univers et Particules de Montpellier (LUPM), Université de Montpellier & CNRS, Place Eugène Bataillon, 34095 Montpellier Cedex 05, France

<sup>e</sup>Laboratoire de Physique Subatomique & Cosmologie (LPSC), Université Grenoble Alpes, CNRS/IN2P3, 53 avenue des Martyrs, 38026 Grenoble, France
E-mail: thomas.lacroix@uam.es, gaetan.facchinetti@ulb.be, judit.perez@uam.es, martin.stref@univ-smb.fr, lavalle@in2p3.fr, david.maurin@lpsc.in2p3.fr, miguel.sanchezconde@uam.es

**Abstract.** Gamma-ray observations have long been used to constrain the properties of dark matter (DM), with a strong focus on weakly interacting massive particles annihilating through velocity-independent processes. However, in the absence of clear-cut observational evidence for the simplest candidates, the interest of the community in more complex DM scenarios involving a velocity-dependent cross-section has been growing steadily over the past few years. We present the first systematic study of velocity-dependent DM annihilation (in particular $p$-wave annihilation and Sommerfeld enhancement) in a variety of astrophysical objects, not only including the well-studied Milky Way dwarf satellite galaxies, but nearby dwarf irregular galaxies and local galaxy clusters as well. Particular attention is given to the interplay between velocity dependence and DM halo substructure. Uncertainties related to halo mass, phase-space and substructure modelling are also discussed in this velocity-dependent context. We show that, for $s$-wave annihilation, extremely large subhalo boost factors are to be expected, up to $10^{11}$ in clusters and up to $10^6 - 10^7$ in dwarf galaxies where subhalos are usually assumed not to play an important role. Boost factors for $p$-wave annihilation are smaller but can still reach $10^3$ in clusters. The angular extension of the DM signal is also significantly impacted, with e.g. the cluster typical emission radius increasing by a factor of order 10 in the $s$-wave case. We also compute the signal contrast of the objects in our sample with respect to annihilation happening in the Milky Way halo. Overall, we find that the hierarchy between the brightest considered targets depends on the specific details of the assumed particle-physics model.

**Keywords:** Dark matter searches, Gamma rays, Galactic dynamics
# Contents

1 Introduction 1

2 Velocity-dependent annihilation: theoretical ingredients 3

2.1 Dark matter annihilation and self-interaction: Sommerfeld effect 3

2.2 Generalized J-factors and phase-space modelling 6

2.3 Host halo and subhalos: generalized boost factor 7

3 Selected targets and mass modelling 9

3.1 Dwarf spheroidal galaxies 9

3.2 Dwarf irregular galaxies 11

3.3 Galaxy clusters 12

4 Generalised J-factors for host halos without substructures 15

4.1 General features for s-wave and p-wave annihilations 15

4.2 Uncertainties from the mass modelling 18

4.3 Impact of uncertainties on the ranking of targets 19

5 Generalised J-factors with substructure boost 20

5.1 Impact of subhalos: boost factors for generalised J-factors 20

5.2 Ranking of target classes 22

5.3 Comparison to previous works and prospects for γ-ray DM searches 24

6 Summary and conclusions 31

A Uncertainties related to phase-space modelling 33

B The SL17 subhalo population model 34

C Numerical calculation of the subhalo boost 36

D Dependence of $J_S^{\text{target}}/J_S^{\text{MW}}$ on $\theta_{\text{int}}$ and $\epsilon_{\phi}$ 37

## 1 Introduction

With the absence of robust non-gravitational evidence for dark matter (DM), astrophysical observations remain a prime avenue to find DM; the latter can be made of exotic particles, macroscopic compact objects, or be a result of an incomplete understanding of gravity. The crucial role of astrophysics has already proven essential in the quest for particle self-annihilating DM candidates, thanks to multi-wavelength and multi-messenger constraints on the properties of the underlying DM particle candidates [1–4]. In particular, γ-ray searches, traditionally used to look for signatures of the annihilation of weakly interacting massive particles (WIMPs), are a very powerful tool to explore broader classes of particle DM candidates.

In the past decades, significant efforts have been devoted to the “vanilla” WIMP DM scenario, in which the DM particles annihilate through an s-wave process with a thermal cross
section of ~3×10⁻²⁶ cm³ s⁻<fnref n="1" />. The absence of a firm detection in the 10-100 GeV range by the *Fermi* Large Area Telescope (*Fermi* LAT) on board the NASA Fermi satellite [5, 6] has led to a broadening of search strategies in recent years. First of all, the thermal *s*-wave scenario can only be robustly excluded in a small range of masses [7], while the (multi-)TeV range is still essentially unconstrained despite the close scrutiny of Cherenkov γ-ray telescopes from the ground [8].<fnref n="1" /> Moreover, the bulk of the effort in γ-ray DM searches has been focused on canonical velocity-independent (*s*-wave) DM annihilations, while velocity dependence of the annihilation cross section can appear in many theoretical contexts [11]. This case leads to a much richer phenomenology, upon which we build in this work.

Present-day WIMP annihilations may be naturally suppressed if the *s*-wave partial wave contribution is negligible and annihilations result from a *p*-wave process instead, i.e., $\sigma v_{rel} \propto v_{rel}^2$. This arises, for instance, if the interaction between WIMPs and Standard Model particles is mediated by a scalar. Due to the small DM-induced γ-ray fluxes expected for *p*-wave annihilations, constraints focusing on dwarf spheroidal satellite galaxies (dSphs) are fairly weak and do not reach the level of the thermal *p*-wave annihilation cross section [12–14]. Moreover, the standard paradigm has moved from scenarios involving one DM candidate to more complex dark sectors in which DM annihilation can be mediated, for instance, by new light states [15, 16]. In these scenarios, the exchange of light mediators induces a long-range interaction between the DM particles. This modifies the short-range annihilation cross section, enhancing it for an attractive interaction [17–22]. This effect is referred to as the Sommerfeld enhancement, and leads to a specific velocity dependence of the annihilation cross section; the latter boils down, in some regions of parameter space, to inverse powers of the relative velocity (i.e., enhanced annihilation at low velocities).

Another crucial aspect of the calculation is the DM distribution in subhalos, which has long been recognized to play a very important part in predictions of annihilation signatures [23–29]. The presence of subhalos enhances γ-ray fluxes in the outskirts of their host halos, boosting the total signal. However, this boost is strongly dependent on the host halo mass and the structural properties of its subhalos; see for instance [30–35] for boost calculated for velocity-independent annihilation fluxes. The calculation of the boost is further complicated for velocity-dependent annihilations, as discussed in [15, 36–39]. In this work, we go beyond previous results and provide a systematic study of the substructure boosts resulting from the interplay of both host halo mass and subhalo structural properties.

Overall, many studies have focused on refining the modelling of velocity-independent annihilations for γ-ray DM searches, considering and ranking potential astrophysical targets inside the Galaxy (Galactic centre, subhalos, and dSphs) and outside (nearby galaxies and galaxy clusters, diffuse extragalactic emission). These studies have helped to refine pointing and analysis strategies for ground-based Cherenkov telescopes and spaceborne γ-ray instruments. In contrast, velocity-dependent annihilations [15, 36, 40–42] have been less systematically considered. Recent studies have mostly focused on dSphs [12–14, 41, 43–47], with fewer studies on other targets (Galactic centre [48, 49], Milky-Way [50, 51], subhalos [52, 53], galaxy clusters [54], and diffuse extragalactic emission [55, 56]).<fnref n="2" /> For this reason, it is timely to perform a more systematic study and comparison for a wider variety of target objects and particle physics models, in particular for *p*-wave annihilation.

<footnote n="1">¹The upcoming Cherenkov Telescope Array (CTA), that will lead the field of γ-ray DM searches in the near future, is expected to provide competitive DM constraints in the TeV regime soon [9, 10].</footnote>

<footnote n="2">²Complementary to γ-ray searches, velocity-dependent DM annihilation has also been studied in a cosmological context, see e.g. [57–59].</footnote>
This paper is dedicated to detailed predictions of the astrophysical factors that determine expected DM-induced γ-ray fluxes, as well as a thorough discussion of the associated theoretical uncertainties; connecting these predictions to actual γ-ray data analysis and discussing the full implications in terms of DM models is left for a follow-up study. For our purpose, we consider several targets, avoiding those with too diffuse or too extended expected signals. These targets are taken among three families, namely dSphs, dwarf irregular galaxies (dIrrs) – these ones included for the first time in a velocity-dependent study – and galaxy clusters. Our goal is to address how the intra- and inter-family ranking is impacted by the particle-physics model considered. In the process, we make several improvements with respect to previous calculations. The main novelty is that, for all targets, we consider the subhalo boost obtained from a self-consistent semi-analytic model, i.e., reconstructing the velocity dependence at all scales. This is done via the reconstruction of the phase-space distribution function from the mass modelling. We study in particular the complex interplay between velocity dependence and substructure boost, and the implications in terms of ranking of the various targets. A more in-depth and analytical study of the boost factor is presented in a companion paper [60].

The paper is organized as follows. In Sec. 2, we introduce the various ingredients needed to perform the calculation of the generalised astrophysical factor of the flux, the so-called *J*-factor (velocity-weighted DM squared density integrated over the phase-space distribution): (i) we recall the regimes where the Sommerfeld enhancement can occur, in particular we focus on the *s*- and *p*-wave cases; (ii) we describe how we invert the mass model to obtain the phase-space DM distribution entering the *J*-factor calculations; (iii) we discuss the calculation of the boost factors in the context of generalised *J*-factors. In Sec. 3, we detail the mass modelling of specific targets, selected among the sample of known dSph galaxies, dIrr galaxies, and galaxy clusters. In Sec. 4, we present our *J*-factor results for the host halos alone (i.e., no substructures), highlighting the two main uncertainties at this stage of the calculation, namely the DM density profile and phase-space modelling. Then, in Sec. 5, we detail how the targets of interest are boosted by the presence of subhalos and how this may affect the intensity and the spatial morphology of the DM signal. We also provide a ranking of these targets (in terms of their generalised *J*-factor) and highlight the potential of galaxy clusters for *p*-wave annihilation. Finally, we conclude and discuss the next steps of our work in Sec. 6. To ease the reading, we postpone to the appendix the discussion of uncertainties related to phase-space modelling (App. A), the details of the subhalo model (App. B), the expressions used to perform the numerical computation of the subhalo boost factor (App. C), and a discussion of the signal contrast between our targets and the Milky Way (MW) DM annihilation foreground (App. D).

## 2 Velocity-dependent annihilation: theoretical ingredients

### 2.1 Dark matter annihilation and self-interaction: Sommerfeld effect

In this section, we provide a brief review of the impact of DM self-interaction on the physics relevant to γ-ray searches. We consider the phenomenological scenario in which DM particles self-interact through the exchange of a light mediator. In the absence of such interactions, the annihilation cross section can be computed perturbatively from the physics of the short-range annihilation processes. However, a light mediator leads to a long-range interaction which can distort the wave function of the corresponding two-body system in a non-perturbative way,
leading to Sommerfeld enhancement<fnref n="3" /> of the annihilation cross section in the non-relativistic regime [17–19]. Since the effect appears in the non-relativistic limit, the enhancement can be computed by solving the Schrödinger equation for the scattering of two DM particles. The radial part of the wave function $R_ℓ(r) = χ_ℓ(r)/r$, for the partial wave with angular momentum $ℓ$, solves

$$ \left(-\frac{\hbar^2}{m_\chi} \partial_r^2 + \frac{\hbar^2 \ell(\ell+1)}{m_\chi r^2} + V(r)\right) \chi_\ell(r) = \frac{(\hbar k)^2}{2m_\chi} \chi_\ell(r), \qquad (2.1) $$

where $m_χ$ is the DM mass, $V$ is the interaction potential and $k = m_χv/ħ$ is the wave vector at infinity of one of the incoming DM particles in the centre-of-mass frame.

Equation (2.1) is solved with the boundary conditions that the interaction only leads to outgoing spherical waves at infinity, and with $R_ℓ(r) \propto r^ℓ$ as $r \to 0$. The Sommerfeld enhancement factor reads [61–63]

$$ S_{\ell} = \left| \frac{(2\ell + 1)!! \chi_{\ell}^{(\ell+1)}(0)}{(\ell + 1)! k^{\ell+1}} \right|^2 . \qquad (2.2) $$

The factor $S_ℓ$ multiplies the corresponding term in the partial wave expansion of the annihilation cross-section

$$ \sigma v_{\text{rel}} = S_s \sigma_0 c + S_p \sigma_1 c \left(\frac{v_{\text{rel}}}{c}\right)^2 + \mathcal{O}\left(\left(\frac{v_{\text{rel}}}{c}\right)^4\right), \qquad (2.3) $$

where $S_s$ is the s-wave factor ($ℓ = 0$), $S_p$ is the p-wave factor ($ℓ = 1$), etc. To make an explicit computation, we must fix the interaction potential. The relevant one for an attractive interaction through a massive mediator is the Yukawa potential

$$ V_Y(r) = -\alpha_D \frac{e^{-m_\phi r}}{r}, \qquad (2.4) $$

with $m_φ$ the mediator mass and $α_D$ the dark fine-structure constant. Unfortunately, solutions of the Schrödinger equation for the Yukawa potential are only known numerically. Luckily, analytical expressions are known for the closely related Hulthén potential

$$ V_H(r) = -\alpha_D \frac{m_* e^{-m_* r}}{1 - e^{-m_* r}}. \qquad (2.5) $$

Fixing $m_* = (\pi^2/6)m_φ$ makes the Hulthén analytical solution close to the Yukawa numerical solution [61], therefore we use the former in the following. From the Schrödinger equation, one sees that the enhancement factor is only a function of two dimensionless parameters

$$ \epsilon_v \equiv \frac{v}{\alpha_D c} \quad \text{and} \quad \epsilon_\phi \equiv \frac{m_\phi}{\alpha_D m_\chi}. \qquad (2.6) $$

In this paper, we focus on s-wave and p-wave annihilation processes, corresponding to $ℓ = 0$ and $ℓ = 1$, respectively. This already covers a broad variety of underlying particle-physics models. For an s-wave annihilation process, the enhancement factor can be written

<footnote n="3">³We restrict ourselves to symmetric DM with attractive interactions, for which the Sommerfeld factor is effectively an enhancement factor.</footnote>
as [63]

$$ S_s(v) \approx \begin{cases} \frac{\sinh\left(\frac{2\pi\epsilon_v}{\epsilon_\phi^*}\right)}{\epsilon_v \cosh\left(\frac{2\pi\epsilon_v}{\epsilon_\phi^*}\right) - \cos\left(2\pi\sqrt{\frac{1}{\epsilon_\phi^*} - \frac{\epsilon_v^2}{\epsilon_\phi^{*2}}}\right)} & \text{if } \epsilon_v \le \sqrt{\epsilon_\phi^*} \\ \frac{\pi/\epsilon_v}{1 - e^{-\pi/\epsilon_v}} & \text{otherwise,} \end{cases} \quad (2.7) $$

where $\epsilon_\phi^* \equiv (\pi^2/6)\epsilon_\phi$. It should be noted that the expression of the Sommerfeld factor in the first line of Eq. (2.7) — which is the standard result for the Hulthén potential in the literature — is only valid when both $\epsilon_\phi, \epsilon_v \ll 1$, i.e., in the regime of large enhancements, as discussed in Ref. [63]. When $\epsilon_v \gg \sqrt{\epsilon_\phi}$, we recover the standard solution for the Coulomb potential corresponding to a massless mediator. For p-wave annihilation, the enhancement factor is

$$ S_p(v) = \frac{\left(1 - \epsilon_{\phi}^{*2}\right) + 4 \epsilon_{v}^{2}}{\epsilon_{\phi}^{*2} + 4 \epsilon_{v}^{2}} \times S_s(v). \quad (2.8) $$

Different regimes arise according to the values of $\epsilon_v$ and $\epsilon_\phi$, which encode the dependence of the Sommerfeld effect on the relative velocity of the DM particles and the masses of the DM candidate and light mediator:

* At large velocities, $\epsilon_v \gg 1$, there is no enhancement: $S_s \approx 1$ and $S_p \approx 1$;

* In the intermediate regime, for which $\epsilon_\phi \ll \epsilon_v \ll 1$, we have $S_s \approx \pi/\epsilon_v \propto 1/v$ and $S_p \approx \pi/(4\epsilon_v^3) \propto 1/v^3$. This contains the regime in which the interaction potential tends to a Coulomb potential (for $\epsilon_v \gg \sqrt{\epsilon_\phi}$) but spans a broader range of values of $\epsilon_v$;

* The regime of small velocities, i.e., $\epsilon_v \ll \epsilon_\phi \ll 1$, corresponds to the saturation regime of the Sommerfeld effect which is almost independent of the velocity of the DM particles, except at a series of resonances — namely $\epsilon_\phi = 6/(\pi^2 n^2)$ with n is an integer — for which $S_s \approx 1/(n^2\epsilon_v^2) \propto 1/v^2$ and $S_p \approx (n^2 - 1)^2/(n^2\epsilon_v^2) \propto 1/v^2$;

* Finally when the mediator is heavy, i.e., $\epsilon_\phi \gg 1$, there is again no enhancement, so $S_s \approx 1$ and $S_p \approx 1$.

The analytic solution typically reproduces the numerical result within 10%, except close to resonances where larger differences arise since the analytic resonances for the Hulthén potential are slightly offset from the ones obtained for the Yukawa potential, as discussed for instance in Refs. [43, 64]. However, for the purpose of this work the features of the solution, especially the resonances, are sufficiently well accounted for by the analytic solution. For numerical calculations, we consider a benchmark value of $\alpha_D = 10^{-2}$, but generalised J-factors can be easily rescaled, and we provide the scalings whenever relevant.

It should be noted that long-range interaction can lead to the formation of unstable bound states which modify the annihilation cross-section and are not taken into account in the equations above. In the resonant regime, we regularise the resonances by performing the replacement $v \to v + \alpha_D^4$ which accounts for the finite lifetime of these intermediate states [19, 64, 65]. Bound states can also form in the Coulomb regime where they might significantly change the overall cross-section, however the velocity dependence is left unchanged [66, 67] hence we choose to ignore this effect.
## 2.2 Generalized J-factors and phase-space modelling

The DM-induced γ-ray flux integrated over a sky region of solid angle ΔΩ reads<fnref n="4" />

$$ \frac{d\Phi_{\gamma}}{dE_{\gamma}} = \frac{1}{4\pi} \frac{(\sigma v_{\text{rel}})_{0}}{\eta m_{\chi}^{2}} \frac{dN}{dE_{\gamma}} J_{\text{S}}(\Delta\Omega), \qquad (2.9) $$

where $dN/dE_γ$ is the γ-ray spectrum per annihilation, $η = 2$ for self-conjugate DM ($η = 4$ for non-self-conjugate DM), and the astrophysical factor $J_S$ encodes the information on the DM spatial and velocity distribution.

**Generalised J-factors.** We introduce the following notation

$$ \bar{S} = \begin{cases} S_s \left( \frac{v_{\text{rel}}}{2} \right) & (\text{for s-wave annihilation}) \\ \left( \frac{v_{\text{rel}}}{c} \right)^2 S_p \left( \frac{v_{\text{rel}}}{2} \right) & (\text{for p-wave annihilation}) \end{cases} \qquad (2.10) $$

to treat the s- and p-wave annihilations on equal footing. This allows to write

$$ J_S(\Delta\Omega) = \int_{\Delta\Omega} d\Omega \int_{\text{l.o.s.}} ds \int d^3\vec{v}_1 \int d^3\vec{v}_2 f(r(s, \Omega), \vec{v}_1) f(r(s, \Omega), \vec{v}_2) \bar{S}\left(\frac{v_{\text{rel}}}{2}\right), \quad (2.11) $$

where s the line-of-sight (l.o.s.) coordinate, Ω the solid angle, $\vec{v}_{rel} = \vec{v}_2 - \vec{v}_1$ is the relative velocity with $v_{rel} = |\vec{v}_{rel}|$ and $f(r, \vec{v})$ is the phase-space distribution function (PSDF) of the DM (assuming spherical symmetry), normalised to the total mass of the gravitational system of interest, such that the DM density $ρ_χ$ at galactocentric radius r is

$$ \rho_{\chi}(r) = \int f(r, \vec{v}) d^3\vec{v}. \qquad (2.12) $$

Equation (2.11) is referred to as the generalised J-factor. As the name indicates, it is a generalisation of the standard J-factor relevant for s-wave and p-wave annihilation without Sommerfeld enhancement (recovered for $S_s = S_p = 1$). For the sake of clarity, we will denote $J_{S,s}$ the generalised J-factor associated to the s-wave and $J_{S,p}$ the one associated to the p-wave.

Assuming spherical symmetry of the DM halo, the integral over the solid angle becomes an integral over the angular distance $ψ$ from the centre of the object, with $dΩ = 2π sinψ dψ$ and $r(s, Ω) = \sqrt{s^2 + D^2 - 2sD cosψ}$, where D is the distance from the observer to the centre of the object. In the following, we perform the integral over an angular size $θ_{int}$ that depends on the target and can also depend on the γ-ray detection technique. For instance, most dSphs are observed as point-like by Fermi-LAT, whereas galaxy clusters are extended targets, so we take $θ_{int} = 0.5°$ for dSphs and $θ_{int} ≈ R_{200}/D$ for clusters, with $R_{200}$ the virial radius.

In practice, Eq. (2.11) can be rewritten in terms of a J-factor for an effective squared density profile as

$$ J_S(\theta_{\text{int}}) = 2\pi \int_0^{\theta_{\text{int}}} d\psi \sin\psi \int ds \left\langle \bar{S}\left(\frac{v_{\text{rel}}}{2}\right) \right\rangle(r(s, \psi)) \rho_{\chi}^2(r(s, \psi)), \qquad (2.13) $$

<footnote n="4">⁴To ease the comparison with the majority of previous works in the literature, we do not include in the definition of the J-factor the $1/(4π)$ pre-factor. The latter appears in the derivation of an intensity from a volume emissivity, which is accounted for in Eq. (2.9) here. As a result, the J-factors given in this work are expressed in GeV$^2$ cm$^{-5}$ sr.</footnote>
where the average of an observable $O(v_{rel})$ that depends on the relative velocity is given by

$$ \langle O(v_{rel}) \rangle(r) = \frac{1}{\rho_{\chi}^2(r)} \int d^3\vec{v}_{rel} O(v_{rel}) \int d^3\vec{v}_c f(r, \vec{v}_1) f(r, \vec{v}_2), \quad (2.14) $$

with $\vec{v}_c = (\vec{v}_1 + \vec{v}_2)/2$ the centre-of-mass velocity. It should be noted that the profile is truncated at the tidal radius for dSphs, or (conventionally) at the virial radius for the other objects, with no contribution to the line-of-sight integral outside that radius.

**Phase-space modelling.** The main results of this work are based on the Eddington formalism, which provides the full PSDF $f(\vec{r}, \vec{v})$ of a given component of a system in dynamical equilibrium associated with a given density-potential pair. More specifically, under the assumptions of maximal symmetries, namely spherical symmetry of the system, and an isotropic velocity tensor, the PSDF $f(r, \vec{v})$ can be written as a function of the relative energy $\mathcal{E} = \Psi(r) - v^2/2$ only, $f(r, \vec{v}) \equiv f(\mathcal{E})$, where $\Psi$ is the total gravitational potential of the system. In that case, Eq. (2.12) can be uniquely inverted, leading to the well-known Eddington formula [68, 69]:

$$ f(\mathcal{E}) = \frac{1}{\sqrt{8\pi^2}} \left[ \frac{1}{\sqrt{\mathcal{E}}} \left( \frac{d\rho_{\chi}}{d\Psi} \right)_{\Psi=0} + \int_0^{\mathcal{E}} \frac{d^2\rho_{\chi}}{d\Psi^2} \frac{d\Psi}{\sqrt{\mathcal{E} - \Psi}} \right], \quad (2.15) $$

where the first term between brackets is related to the radial boundary of the system. We disregard this term in the following, since we model the system as infinite as far as the PSDF is concerned when computing $J$-factors, considering that the latter are not sensitive to the very outer parts of the system.

In App. A we provide a brief overview of the prediction methods for the PSDF $f(r, \vec{v})$ of DM particles from first principles that we use in this study. In particular, the Eddington formalism can be extended to anisotropic PSDFs under some specific assumptions. We use this extended formalism to quantify the uncertainty on generalised $J$-factors from the modelling of the PSDF itself, which we discuss in App. A. In the main text, we restrict the presentation to the Eddington method which provides a very good approximation. Technical details regarding the semi-analytic derivation of averages over the relative velocity distribution, Eq. (2.14), for various assumptions on the anisotropy of the velocity distribution, can be found in Ref. [70].

## 2.3 Host halo and subhalos: generalized boost factor

DM subhalos, which are characteristic of any self-annihilating CDM particle scenario, are expected to boost the gamma-ray signals that would be predicted assuming DM is smoothly distributed in target halos [23, 24]. We account for this boost factor by means of the analytical subhalo population model developed in Ref. [35] (SL17 henceforth) — see also Refs. [71–74]. Considering subhalos is particularly important when annihilation is velocity-dependent because the internal velocity dispersion in these objects is much smaller than that of the host halo. Indeed, from dimensional arguments and assuming virial equilibrium, the velocity dispersion should scale as $m^{1/3}$ where $m$ is the virial halo mass. This is quite relevant for Sommerfeld-enhanced processes, which depend on inverse powers of the velocity, leading to a potentially strong enhancement when the subhalo mass range extends down to very small masses. Interestingly, heavy WIMPs beyond ~10 TeV are naturally subject to Sommerfeld effects [18, 21], and also lead to subhalo virial masses as small as ~$10^{-12} M_{\odot}$ [75].
**Modelling the host halo and subhalos.** The SL17 model assumes that the host halo can be described by a spherically symmetric and smoothed total DM density profile, $\rho_{\text{host}}(r)$, comprising a genuine smooth component $\rho_{\text{sm}}(r)$ and a subhalo component made of individual objects (DM inhomogeneities), but globally described by an average density profile $\rho_{\text{sub}}(r)$. These components are simply bound to obey the relation

$$ \rho_{\text{host}}(r) = \rho_{\text{sm}}(r) + \rho_{\text{sub}}(r) \ge 0, \qquad (2.16) $$

where $r$ is the radial distance to the host's centre. The density profile $\rho_{\text{host}}$ can in principle be constrained both theoretically and observationally, in particular its shape (cuspy or cored, external tail, etc.), and both its global and internal properties (virial and/or tidal mass, concentration, etc.). For all of the host halos studied in this paper, we do use observationally constrained density profiles for $\rho_{\text{host}}$ (see next section), and the smooth halo component $\rho_{\text{sm}}$ is obtained by subtracting $\rho_{\text{sub}}$ from $\rho_{\text{host}}$, see Eq. (2.16). The average density profile of subhalos, $\rho_{\text{sub}}$, is calculated from the SL17 model. It can be expressed in terms of a continuous number density of subhalos, $n_{\text{sub}}$, depending on the virial$^5$ mass $m = M_{200}$ and concentration $c = c_{200}$ of the subhalos, and their radial position $r$ in the host halo:

$$ \begin{aligned} \rho_{\text{sub}}(r) &= \int dm \int dc \, m_t(m, c, r) \left\{ \frac{d^5 N_{\text{sub}}(m, c, r)}{d^3 \vec{r} \, dm \, dc} \right\} \equiv N_{\text{tot}} \frac{d^5 \mathcal{P}_{\text{sub}}(m, c, r)}{d^3 \vec{r} \, dm \, dc} \\ &= n_{\text{sub}}(r) \times \langle m_t(m, c) \rangle_{(m,c)}(r), \end{aligned} \qquad (2.17) $$

where $m_t \neq m$, the physical tidal (not virial) mass of subhalos, critically depends on the position $r$. Tidal stripping effects are sourced by all gravitational components of the host, leading in particular to a calculated total number $N_{\text{tot}}$ of surviving subhalos. Moreover, although the spatial dependence of the mass and concentration PDFs were initially set homogeneous (from cosmological considerations), the tidal effects make the calculated probabilistic parameter phase space $d\mathcal{P}_{\text{sub}}(m, c, r)$ fully intricate and non-separable; the latter is normalised to unity over the whole halo phase-space volume (position, mass, and concentration parameters). A short presentation of the SL17 model with more technical details is given in App. B.

**Generalized J-factor for subhalos.** From Eq. (2.17), we can write

$$ \begin{aligned} \underline{\rho}_{\text{sub}}^2(r) &= N_{\text{tot}} \rho_{\oplus}^2 \int dm \int dc \, \xi_t(m, c, r) \frac{d^5 \mathcal{P}_{\text{sub}}(m, c, r)}{d^3 \vec{r} \, dm \, dc} \\ &= n_{\text{sub}}(r) \rho_{\oplus}^2 \langle \xi_t(m, c) \rangle_{(m,c)}(r), \end{aligned} \qquad (2.18) $$

where we defined the subhalo effective tidal annihilation volume$^6$

$$ \xi_t(m, c, r) \equiv \int_{x \le x_t(m, c, r)} d^3 \vec{x} \left\{ \frac{\rho(x, m, c)}{\rho_{\oplus}} \right\}^2, \qquad (2.19) $$

with $\rho(x, m, c)$ the inner subhalo profile and $x_t$ its tidal extension. This tidal annihilation volume can be generalized to the velocity-dependent Sommerfeld enhancement case by writing

$$ \xi_{S,t}(m, c, r) \equiv \int_{x \le x_t(m, c, r)} d^3 \vec{x} \left\{ \frac{\rho(x, m, c)}{\rho_{\oplus}} \right\}^2 \langle \bar{S} \rangle(x), \qquad (2.20) $$

<footnote>$^5$We use the conventional definitions, where the index “200” indicates that quantities are defined with respect to some virial radius $R_{200}$ (or $x_{200}$ for subhalos) over which the average density of a halo is 200 times the critical density at redshift zero. The actual extension of a halo is not necessarily its virial radius; for subhalos, the physical extension is taken to be the tidal radius, $r_t \le R_{200}$.</footnote>

<footnote>$^6$This is the volume a subhalo would have in order to sustain its own annihilation rate if it had an arbitrary constant DM density of $\rho_{\oplus}$ (similar to an intrinsic annihilation luminosity except for physical dimensions).</footnote>
where $\langle \rangle$ denotes the velocity average over the 2-particle phase-space volume introduced in Eq. (2.14). We stress that this average is taken over the subhalo PSDF, which depends on c and m and is very different from the host PSDF. Eventually, the generalized J-factor associated with the total subhalo contribution reads

$$J_{\text{S,sub}}(\theta_{\text{int}}) = 2\pi \int_{0}^{\theta_{\text{int}}} d\psi \sin\psi \int ds \rho_{\text{S,sub}}^{2}(r(s, \psi)), \qquad (2.21)$$

with

$$\underline{\rho}_{\text{S,sub}}^{2}(r) = n_{\text{sub}}(r) \rho_{\oplus}^{2} \langle \xi_{\text{S,t}}(m, c, r) \rangle_{(m,c)}(r). \qquad (2.22)$$

**Total J-factor and boost.** The total generalized J-factor is obtained summing up all contributions (host and subhalos):

$$J_{\text{S,tot}}(\theta_{\text{int}}) = 2\pi \int_{0}^{\theta_{\text{int}}} d\psi \sin\psi \int ds \rho_{\text{S,tot}}^{2}(r(s, \psi)), \qquad (2.23)$$

where

$$\underline{\rho}_{\text{S,tot}}^{2} = \underline{\rho}_{\text{S,sub}}^{2} + \underline{\rho}_{\text{S,sm}}^{2} + 2 \underline{\rho}_{\text{S,sm}} \underline{\rho}_{\text{sub}} \neq \underline{\rho}_{\text{S,host}}^{2}. \qquad (2.24)$$

All these terms include a Sommerfeld-enhancement correction (subscript S). The last term before the inequality is the cross-product between the smooth DM component and subhalos, for which the relevant velocity field is that of the host halo; this term can actually safely be neglected [35]. For further technical details on how $J_{\text{S,tot}}$ is computed in practice, we refer the reader to App. C.

Finally, we can formally define the generalized subhalo boost factor as

$$\mathcal{B}_{\text{S}} \equiv \frac{J_{\text{S,tot}}}{J_{\text{S,host}}} \approx 1 + \frac{J_{\text{S,sub}}}{J_{\text{S,host}}}, \qquad (2.25)$$

where $J_{\text{S,host}}$ is evaluated from the (squared) smoothed host profile given in Eq. (2.16), also corrected for the Sommerfeld enhancement as above. The approximation on the right-hand-side is valid only when the host is distant enough so that most of the annihilation rate is contained within the angular resolution of the telescope; in that case, $J_{\text{S,sm}} \simeq J_{\text{S,host}}$. Note that a detailed analytical study of the subhalo boost factor in the context of the Sommerfeld enhancement has been carried out in a companion paper [60].

# 3 Selected targets and mass modelling

We list in this section the astrophysical targets considered for our work (dwarf spheroidal galaxies, dwarf irregular galaxies, and galaxy clusters). For each of these object classes, we motivate our specific selection and discuss the DM density profile used for our analyses.

## 3.1 Dwarf spheroidal galaxies

Owing to their close distance (tens of kpc), potentially high DM densities, and negligible astrophysical background, MW dSph satellites are among the most promising targets for indirect DM detection [76, 77]. In the absence of a clear signal, the best current limits on WIMP DM candidates in γ-rays were obtained from the combined analysis of Fermi-LAT
data on many dSphs [78–80]. Their DM content is inferred from the velocity dispersion of their stellar population (obtained from spectroscopic measurements), via moments of the Jeans equations [69, 81]. DSphs are typically separated in two categories: ‘classical’ and ‘ultra-faint’. The former are brighter, with hundreds to thousands of stars measured, while the latter are fainter, with only tens of known member stars. It ensues that the DM content of ultra-faint dSphs suffers larger uncertainties than that of the classical ones, which translates into less robust constraints on the DM particle properties in the former case. However, more and more ultra-faint dSphs are discovered thanks to optical surveys [82–91], or will be discovered in the next decade [92], and those potentially located just tens of kpc away from us could shine even brighter than the classical dSphs in terms of γ-rays from their annihilating DM halos.

In the last decade, many studies have refined and improved the calculation of J-factors, in order to rank the best targets [14, 93–102]. Although these studies are overall in broad agreement, the assumptions made on the underlying ingredients (light and DM profiles, anisotropy distribution, triaxiality), methodology (e.g., using higher moments of the Jeans equation), and statistical analysis framework and priors used (data-driven approach, DM-simulation or mock-data based priors, etc.) can lead to sizeable differences in the expected DM signal and also on the J-factor uncertainties of some dSphs (factor of a few). Discussing the relative merits of each approach to single out the best one goes beyond the scope of this paper, and is in any case a very difficult task: all studies consider slightly different but mostly relevant methodologies (with different limitations) for the reconstruction of DM density profiles. With the improvement on stellar structural parameters [103, 104] and new spectroscopic data [105, 106], predictions for the DM halo will hopefully become less uncertain, in particular for ultra-faint dSphs (see, e.g., [104] for a recent review).

For definiteness, we pick here two classical dSphs (Draco and Sculptor, respectively in the Northern and Southern sky) and one prototypical ultra-faint (Reticulum II), which were found to be among the best-ranked targets for DM annihilation in [97, 107]. In the latter studies, the DM profile parameters were reconstructed from a Markov Chain Monte Carlo (MCMC) engine coupled to a Jeans analysis with the CLUMPY code<fnref n="7" /> [33, 108]; see [109] for more details on the methodology. We use these chains to calculate the median profile, that we adopt as a reference for our analyses here<fnref n="8" />. We gather in Table 1 the position and DM profile parameters for the three selected dSphs, modelled following an Einasto profile:

$$ \rho_{\text{Ein}}(r) = \rho_{-2} \exp \left\{ -\frac{2}{\alpha} \left[ \left( \frac{r}{r_{-2}} \right)^{\alpha} - 1 \right] \right\}. \qquad (3.1) $$

The parameters $r_{-2}$, $\rho_{-2}$, and $\alpha$ are the radius for which the slope is -2, the DM density at this radius, and the slope of the Einasto profile, respectively.

The generalised J-factors for velocity-dependent cross-sections rely on the DM phase space distribution (see Sec. 2.2). To properly and fully propagate the DM profile uncertainties to the generalised J-factor, we start from 1000 profile parameter samples taken from the analysis of [97] ([107] for Reticulum II), apply the Eddington calculation to obtain the

<footnote n="7">⁷https://clumpy.gitlab.io/CLUMPY/</footnote>

<footnote n="8">⁸We do not use the best-fit profile parameters, because the scarcity of data in ultra-faint dSphs make them display an unphysical behaviour (e.g., a very flat and extended profile). Using ‘effective’ structural parameters matching the median profile cures this issue. The generalised J-factors calculated from these effective parameters are also found to be very close to the median generalised J-factors calculated over the MCMC values.</footnote>
<table><caption>**Table 1.** Relevant parameters of the three MW dSphs selected for this analysis. The first columns report the dSph position (Galactic longitude and latitude) and distance from the observer: for consistency with the analysis of Ref. [97], that we follow here, we take the distance from Ref. [110], although more recent estimates can slightly differ [111]. The last three columns list the Einasto profile parameters (normalisation $\rho_{-2}$, scale radius $r_{-2}$, and slope $\alpha$) corresponding to the median profile calculated over our MCMC sample.</caption><thead><tr><th rowspan="2">dSph</th><th>(l, b)</th><th>D</th><th rowspan="2">ρ<sub>-2</sub><br>[10<sup>7</sup>M<sub>☉</sub> kpc<sup>-3</sup>]</th><th rowspan="2">r<sub>-2</sub><br>[kpc]</th><th rowspan="2">α</th></tr><tr><th>[deg]</th><th>[kpc]</th></tr></thead><tbody><tr><td>Reticulum II (Ret2)</td><td>(266.3, -49.7)</td><td>30</td><td>2.53</td><td>0.92</td><td>0.46</td></tr><tr><td>Sculptor (scl)</td><td>(287.5, -83.2)</td><td>79</td><td>2.87</td><td>0.50</td><td>0.31</td></tr><tr><td>Draco (dra)</td><td>(86.4, +34.7)</td><td>82</td><td>1.06</td><td>2.09</td><td>0.46</td></tr></tbody></table>

phase-space associated to each profile (see also App. A), and then calculate the associated generalised *J*-factor. From the distribution of the 1000 calculated *J*-factors, we can calculate any quantile to derive the mean of the distribution and its uncertainties.

## 3.2 Dwarf irregular galaxies

Dwarf irregular (dIrrs) galaxies have recently entered in the list of prime targets for indirect γ-ray DM searches. Indeed, the existence of these isolated galaxies within the Local Group, at $O(1$ Mpc) distances, makes them interesting targets given both their proximity and typical masses $M_{200} \approx 10^7 - 10^{10}$M$_{☉}$. DIrrs are rotationally-supported objects, allowing to reconstruct the underlying DM density profiles from their measured rotation curves (RCs). Such RC studies show that dIrrs are DM-dominated objects at all radii [112–114]. Unlike dSphs, dIrrs are star-forming galaxies, yet the γ-ray emission associated to astrophysical processes has been estimated to be negligible compared to that expected from DM annihilation [115, 116]. One more reason that makes dIrrs promising targets for DM searches is the fact that, given their typical host halo masses, the so-called subhalo boost is expected to be significant in their case, reaching values up to ~5 [117], depending on the definition of the boost factor. This is in contrast to the case of dSphs, which are not only less massive than dIrrs but also tidally stripped objects, thus with expected subhalo boosts of the order of only a few percent [34]. Despite the above considerations, dIrrs have not been used for γ-ray DM searches up to just recently [116, 117].

With the current available observational data, the study of dIrrs RCs is not conclusive and, indeed, there is still a debate in the literature about the precise inner shape of the DM density profile in these objects. Fits to the RCs favor core-like profile [118], yet this conclusion is in contrast with that expected from N-body cosmological simulations, that point to a universal cuspy profile like NFW [119, 120] or Einasto [121]. Multiple studies have investigated the source of this apparent disagreement — not unique to this type of objects — between data and ΛCDM expectations, providing different solutions mainly based on the impact of baryonic feedback on the DM distribution and its ability to shallow the initial cusps in the innermost regions of the DM density profiles, especially at some particular mass scales [122–126]. As this issue is far from being solved, the authors in [117] adopted an agnostic path and decided to perform a DM modelling for dIrrs using the two different types
<table><thead><tr><th rowspan="2">dIrr</th><th rowspan="2">(l, b)<br>[deg]</th><th rowspan="2">D<br>[kpc]</th><th rowspan="2">M<sub>200</sub><br>[10<sup>10</sup> M<sub>☉</sub>]</th><th rowspan="2">Profile</th><th colspan="2">ρ</th><th colspan="2">R<sub>200</sub></th></tr><tr><th>[10<sup>7</sup> M<sub>☉</sub> kpc<sup>-3</sup>]</th><th>[kpc]</th><th>[kpc]</th></tr></thead><tbody><tr><td rowspan="2">NGC6822</td><td rowspan="2">(25.34, -18.40)</td><td rowspan="2">480</td><td rowspan="2">3.16</td><td>Burkert*</td><td>3.16</td><td>3.3</td><td>62.9</td></tr><tr><td>NFW</td><td>0.79</td><td>5.9</td><td>62.6</td></tr><tr><td rowspan="2">IC10</td><td rowspan="2">(118.96, -3.33)</td><td rowspan="2">790</td><td rowspan="2">3.98</td><td>Burkert*</td><td>15.85</td><td>2.0</td><td>71.3</td></tr><tr><td>NFW</td><td>0.63</td><td>6.8</td><td>70.3</td></tr><tr><td rowspan="2">WLM</td><td rowspan="2">(75.87, -73.86)</td><td rowspan="2">970</td><td rowspan="2">0.40</td><td>Burkert*</td><td>6.31</td><td>1.3</td><td>33.3</td></tr><tr><td>NFW</td><td>1.00</td><td>2.8</td><td>33.6</td></tr></tbody></table>

**Table 2.** Parameters of the three dIrrs in our sample. The first column reports the dIrrs position (galactic longitude and latitude) and distance from the observer (see [117] and references therein). We then list the best-fit profile parameters for Burkert and NFW (mass, normalisation and scale radius), as well as the virial radius assuming an overdensity of 200 times the critical density of the Universe. For columns 6 and 7, $\rho$ and $r$ stand for $\rho_c$ and $r_c$ in case of the Burkert profile, and $\rho_0$ and $r_s$ for NFW. Profiles marked with * are used as reference for all calculations and figures unless indicated otherwise.

of profiles, i.e. (i) a Burkert, core-like profile [127]

$$ \rho_{\text{Bur}}(r) = \frac{\rho_c r_c^3}{(r + r_c)(r^2 + r_c^2)}, \qquad (3.2) $$

where $r_c$ and $\rho_c$ are, respectively, a core radius and DM density, and (ii) an NFW cusp-like profile

$$ \rho_{\text{NFW}}(r) = \frac{\rho_0}{\left(\frac{r}{r_s}\right)\left(1 + \frac{r}{r_s}\right)^2}, \qquad (3.3) $$

where $r_s$ and $\rho_0$ are, respectively, a scale radius and a characteristic DM density.

In this work, for each of these profiles, we simply use the best-fit parameters obtained in [117], where authors analyze the RCs of 7 dIrrs and obtain a prediction of the J-factors for the two different models of the DM density profile under consideration here. According to the observed RCs, NGC6822, IC10 and WLM are the ones with more available data, thus in these cases the fits are more robust and stable than for the rest of objects in their sample. The mentioned three objects also yield the highest J-factor values independently of the selected DM profile or substructure boost values. Taking these findings in [117] into account, we thus decided to include NGC6822, IC10 and WLM in our sample, whose parameters are gathered in Table 2. As for the dSphs, the calculation of the generalised J-factors relies on the DM phase space distribution described in Sec. 2.2 (based on the inversion of the DM profile). However, at variance with the dSphs, the uncertainties for dIrrs are estimated from the comparison of the results obtained from the Burkert and NFW profiles. The modelling of subhalos in the context of Sommerfeld enhancement relies on the formalism described in Sec. 2.3.

## 3.3 Galaxy clusters

Galaxy clusters are the largest gravitationally-bound objects in the Universe. Their masses are between $M_{200} \approx 10^{14} - 10^{15} M_{\odot}$ and up to 80% of this mass is expected to be DM [128]. The rest is baryonic matter, in the form of galaxies, hot gas and dust in the intra-cluster medium (ICM). Even though clusters are supposedly stable and virialized objects at
present, the presence of hot gas, galaxies, and even Active Galactic Nuclei (AGNs), produces turbulence phenomena and complex baryonic feedback reactions in the ICM (where also significant high magnetic fields are involved). All these astrophysical processes end up acting as acceleration mechanisms, leading to the presence of cosmic rays (CRs), that have been confirmed through the observation of diffuse synchrotron emission produced by the leptonic CRs at different wavelengths [129]. Galaxy clusters have avoided detection in γ-rays so far [130, 131]<fnref n="9" />, but this high-energy emission is indeed expected from hadronic CRs [134, 136, 137].

Despite their expected CR-induced γ-ray emission, galaxy clusters are still considered excellent targets for γ-ray DM searches in the WIMP scenario (from DM annihilation or decay). The DM science case of galaxy clusters soon resulted in studies aimed at determining which galaxy clusters meet the most appropriate conditions to be searched in γ-rays [31, 132, 138–144] and at disentangling both the CR- and DM-induced γ-ray emissions from each other [138, 145]. First, there exists a significant number of local galaxy clusters ($z < 0.1$) for which substantial DM-induced fluxes are expected. Second, DM searches should focus on those with the lowest expected CR backgrounds [138]. In [31, 146], the authors studied the annihilation flux of the most promising galaxy clusters, once DM halo substructures — particularly relevant for clusters — were taken into account. It was found that the brightest galaxy clusters can yield total annihilation fluxes as large as some of the dSphs. Furthermore, for clusters, the annihilation flux profiles become comparatively more spatially extended, as most subhalos are located in the outer halo regions. Overall, this subhalo boost to the annihilation signal is expected to play a key role for clusters as compared to other targets, such as dSphs and dIrrs, for which the boost is negligible or much smaller, respectively [34, 95, 117, 146]. We note, however, that the inclusion of halo substructure, in the case of expanding the annihilation cross-section to p-waves and in the framework of Sommerfeld enhancement, becomes more complex and requires a specific approach that is addressed in Sec. 2.3.

**Halo mass modelling.** For this work, we follow Ref. [146] as a starting point to build our sample of most promising galaxy clusters for DM searches. Their sample was constituted by Virgo, Coma, Fornax, Ophiuchus and Perseus.<fnref n="10" /> Yet, some of these clusters present major observation inconveniences. While Virgo exhibits the highest *J*-factor, it is currently going through a major merger event with the neighbouring M49 galaxy cluster [143]. Also, its proximity to Earth results in an angular extension of several degrees. The observation of an object of this size is extremely challenging given the field of view of existing IACTs. On the other hand, the galactic diffuse emission should be ideally avoided as to simplify any potential DM analysis. This can be easily addressed by removing from our sample those objects located close to the Galactic plane and centre, where this emission is most extreme. This requirement leaves out Ophiuchus, less than 10 degrees far from the Galactic centre. Thus, in the following we will obtain predictions for Coma, Fornax and Perseus and will remove both Virgo and Ophiuchus from our list of clusters. Note that this number of targets is also similar to the numbers in our sample of dSphs and dIrrs. We build the DM density profile of galaxy clusters starting from their measured mass. For nearby galaxy clusters as the ones in our sample, $M_{200}$ can be obtained from X-ray observations of the surface brightness

<footnote n="9">⁹There is a growing evidence, though, for a potential detection in the vicinity of the Coma cluster [132–135].</footnote>

<footnote n="10">¹⁰A comprehensive and systematic ranking of galaxy clusters in terms of their expected annihilation signals can be found in [31], where other targets were also found at the level of those selected for this study.</footnote>
profiles. Indeed, these observations have been used to create catalogues containing the most relevant cluster parameters [147–150]. In our work, we adopt the mass estimates in [151] for Coma and Fornax, while for Perseus we use data from [147] (rescaled to our cosmology). First, we assume the NFW DM density profile given in Eq. (3.3). Assuming a spherical collapse model with an overdensity $\Delta = 200$ times the critical density of the Universe, we can obtain the corresponding virial radius $R_{200}$ than contains the mass $M_{200}$. Now, in order to obtain the two NFW profile parameters we need to assume a concentration-mass ($c-M$) relation. We adopt the parametrization proposed in [32] for main halos. From the value of the concentration and the already obtained $R_{200}$, we can then compute the NFW scale radius $r_s$ as well as the scale density $\rho_0$.

**Mass modelling uncertainties.** The main uncertainties in our DM modelling come from (i) the estimate of the mass as derived from X-rays data and (ii) the intrinsic scatter of the concentration-mass relation. Indeed, it is well known that different observational methods can yield different mass estimates for galaxy clusters. Deviation of mass estimates from surface brightness X-ray measurements with respect to the masses obtained by other observation methods is typically referred to as the hydrostatic bias. Yet, at present there is an on-going debate in the community about how to precisely quantify and treat its value [152]. A complementary approach is to compare X-ray masses, usually labelled as $M_{hydro}$, with the masses provided by other methods, whenever available. For example, the authors in [151] concluded that cluster masses in their catalogue showed a good agreement with the ones obtained from velocity dispersion measurements [153], while this was not the case for the objects in their sample for which SZ measurements [154] were also available. More precisely, for clusters with masses $M_{200} < 5 \times 10^{14} h^{-1} M_{\odot}$, they narrowed down the discrepancy to $M_{hydro}/M_{SZ} = 0.86 \pm 0.01$, and for clusters with larger masses to $M_{hydro}/M_{SZ} = 1.46 \pm 0.08$. From these results, we can conclude that the X-ray mass can be underestimated by $\sim 20\%$ in the case of less massive clusters, while $M_{200}$ can be overestimated by $\sim 50\%$ for the most massive ones.

Following these results, in this work we adopt two mass estimates for each galaxy cluster, that will translate into a bracketing of the $J$-factor uncertainties (due to the cluster mass uncertainty). Our default model is built starting from $M_{hydro}$ and, in addition, we assign each cluster a second mass depending on the above bias. For Fornax, a light cluster, we use $1.2 \times M_{hydro}$ as a second, upper bound mass estimate. In contrast, we adopt $0.5 \times M_{hydro}$ as a lower bound for Perseus and Coma, both massive clusters according to the classification scheme in [151]. As for the uncertainty associated to the scatter of concentrations values for a given mass, we adopt a value of 0.14 dex as suggested by the authors of Ref. [32]. In order to keep a limited number of models, we take advantage of the fact that the $J$-factor $\propto (M_{200}^2 c_{200}^3)/D^2$ to further increase the previous uncertainties by considering extreme values of the concentration scatter. To do so, we consider, for both the upper and lower mass bounds previously derived, $M_{200}^{(\text{min})}$ and $M_{200}^{(\text{max})}$, the concentrations $c_{200}(M_{200}^{(\text{min})}) \times 10^{-\sigma_c}$ and $c_{200}(M_{200}^{(\text{max})}) \times 10^{+\sigma_c}$. The obtained DM density profile parameters for our sample of galaxy clusters are given in Table 3.

**Impact of baryons.** As introduced before in this section, most of galaxy clusters' mass is in the form of non-visible DM, and the rest is accounted for baryonic matter. This baryonic content is mostly encoded in the form of super-heated ionized plasma, the so-called ICM, that accounts for $\sim 15\%$ of the cluster mass, while the remaining $\sim 5\%$ is in the form of galaxies. Because of this, the effect of these baryonic components on the DM modelling of the galaxy
<table><caption>**Table 3.** DM density profile parameters for the three galaxy clusters in our sample. For each target, we show two mass models in order to bracket uncertainties in the corresponding *J*-factors. See text for details. Profiles marked with * are used as reference for all calculations and figures unless indicated otherwise.</caption><thead><tr><th rowspan="2">Cluster</th><th>(l, b)</th><th>D</th><th rowspan="2">Mass estimate</th><th>M<sub>200</sub></th><th>R<sub>200</sub></th><th>ρ<sub>0</sub></th><th>r<sub>s</sub></th></tr><tr><th>[deg]</th><th>[Mpc]</th><th>[10<sup>14</sup> M<sub>☉</sub>]</th><th>[10<sup>2</sup> kpc]</th><th>[10<sup>6</sup> M<sub>☉</sub> kpc<sup>-3</sup>]</th><th>[10<sup>2</sup> kpc]</th></tr></thead><tbody><tr><td rowspan="2">Coma</td><td rowspan="2">(58.09, 87.96)</td><td rowspan="2">102.18</td><td>Hydrostatic</td><td>13.16</td><td>23.19</td><td>2.29</td><td>3.38</td></tr><tr><td>Lower*</td><td>8.77</td><td>20.26</td><td>5.37</td><td>5.58</td></tr><tr><td rowspan="2">Fornax</td><td rowspan="2">(236.72, -53.64)</td><td rowspan="2">20.35</td><td>Hydrostatic*</td><td>0.51</td><td>7.83</td><td>7.42</td><td>1.86</td></tr><tr><td>Upper</td><td>0.61</td><td>8.32</td><td>3.20</td><td>1.05</td></tr><tr><td rowspan="2">Perseus</td><td rowspan="2">(150.57, -13.26)</td><td rowspan="2">80.69</td><td>Hydrostatic</td><td>7.71</td><td>19.41</td><td>2.35</td><td>2.80</td></tr><tr><td>Lower*</td><td>5.14</td><td>16.96</td><td>5.57</td><td>4.59</td></tr></tbody></table>

cluster's main halo can be neglected (as done above), as their contribution to the total mass of the system is even smaller than the size of the uncertainty in the mass estimates themselves. However, the inclusion of the baryonic content in the mean gravitational potential may play a relevant role in the modelling of substructures and the computation of the boost factor since it directly impacts the tidal field experienced by these objects. Indeed, given the typical mass range of the substructures, this second-order effect could lead to different distributions and properties of the subhalo population, meaning that, ideally, we would need to obtain a density model for the baryonic matter.

We thus wanted to quantify this effect for our work, neglecting in a first approximation the galaxies and focusing on the ICM alone, e.g., [155]. Starting from standard X-ray gas density profiles, we built baryon density profiles that included not only electrons, but also protons and Helium following the methodology in [156]. The cluster X-ray parameters were taken from [155]. We found that including baryons in the modelling of the mean gravitational impacts the final boost factors at the level of one percent at most. Thus, in the following, we implicitly neglect the baryonic content in clusters and only show results related to their DM content.

# 4 Generalised *J*-factors for host halos without substructures

In this section, we describe, for the host DM halos of our selected targets, the salient features of $J_S$, the generalised *J*-factor, as a function of $\epsilon_\phi$ (Sec. 4.1). We then show how systematic errors on parameters of the smooth DM profiles translate into systematic uncertainties on $J_S$ (Sec. 4.2). Next, we discuss the ranking of our targets in the various regimes of the Sommerfeld enhancement (Sec. 4.3). We stress that this section only deals with $J_S$ from the smooth DM distribution in our targets; the full calculation of the generalised *J*-factors including the contribution of DM substructures is postponed to the next section (Sec. 5).

We emphasise that all our results, based on full numerical calculations, have been cross-validated (for all regimes) thanks to the analytical calculations presented in the companion paper [60]. This gives us a strong confidence in these results and the conclusions we draw.

## 4.1 General features for *s*-wave and *p*-wave annihilations

The generalised *J*-factors for host halos (i.e., smooth DM distribution) are shown in Fig. 1 for dSphs (top), dIrrs (middle) and galaxy clusters (bottom), for *s*-wave annihilations (left panels) and *p*-wave annihilations (right panels). The behavior of $J_S$ as a function of $\epsilon_\phi$
<table><caption>**Table 4.** Summary table of the characteristic values $\epsilon_{\phi}^*$ of the $\epsilon_{\phi}$ parameter, corresponding to the transition between the Coulomb and saturation regimes of Sommerfeld enhancement, associated with the typical velocity in the object of interest, for $\alpha_D = 10^{-2}$.</caption><thead><tr><th>Target class</th><th>Target</th><th>ε<sup>*</sup><sub>φ</sub></th></tr></thead><tbody><tr><td rowspan="3">dSphs</td><td>Draco</td><td>1.4 × 10<sup>-2</sup></td></tr><tr><td>Sculptor</td><td>5.7 × 10<sup>-3</sup></td></tr><tr><td>Reticulum II</td><td>9.5 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">dIrrs</td><td>IC10</td><td>2.2 × 10<sup>-2</sup></td></tr><tr><td>NGC6822</td><td>1.6 × 10<sup>-2</sup></td></tr><tr><td>WLM</td><td>9 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">Clusters</td><td>Fornax</td><td>1.8 × 10<sup>-1</sup></td></tr><tr><td>Coma</td><td>5.2 × 10<sup>-1</sup></td></tr><tr><td>Perseus</td><td>4.4 × 10<sup>-1</sup></td></tr></tbody></table>

results from the convolution of the Sommerfeld enhancement factor, $S(v)$, with the velocity distribution in each target. The results are directly associated with the various regimes of $S$ discussed in Sec. 2.1.

The main scale of the problem is the characteristic value $\epsilon_{\phi}^*$ at which the Sommerfeld enhancement saturates for a given object. In practice, this transition between the Coulomb and (resonant) saturation regimes can be well reproduced by

$$ \epsilon_{\phi}^* \sim \frac{\bar{v}}{\alpha_D c}, \qquad (4.1) $$

where $\bar{v}$ is the characteristic velocity of the object, for which a good order-of-magnitude estimate<fnref n="11" /> is given by the circular velocity at the scale radius of the DM profile,

$$ \bar{v} \sim \sqrt{\frac{G_N m(r_s)}{r_s}}. \qquad (4.2) $$

The corresponding values of $\epsilon_{\phi}^*$ are given in the last column of Table 4; these values are also relevant for the boost from DM substructure (see App. C). Depending on the ordering of $\epsilon_{\phi}$ and $\epsilon_{\phi}^*$, three different regimes for $J_S(\epsilon_{\phi})$ can be identified.

**s-wave annihilation (left panels of Fig. 1).**

* for $\epsilon_{\phi} \lesssim \epsilon_{\phi}^* \ll 1$, the Sommerfeld enhancement is in the Coulomb regime, i.e., $\bar{S} \propto 1/\epsilon_v$: as a result, $J_S$ does not depend on $\epsilon_{\phi}$—but is roughly proportional to $\bar{v}^{-1}$— and displays a plateau below $\epsilon_{\phi}^*$ (left-hand side of the curves);

* for $\epsilon_{\phi} \gg 1$, there is no enhancement, and $J_S$ boils down to the standard $J$-factor (right-hand side of the curves);

* for $\epsilon_{\phi}^* \lesssim \epsilon_{\phi} \lesssim 1$, this is the resonant (saturation) regime, where the behaviour depends whether $\epsilon_{\phi}$ falls at, or between, resonances (between the two plateaus in the curves):

<footnote n="11">¹¹Roughly speaking, the characteristic velocity is of order $\sqrt{4\pi G_N \rho_0 r_s^2}$, where $\rho_0$ and $r_s$ refer generically to the characteristic density and scale radius of the DM profile considered for each class of object (be it NFW, Einasto or Burkert).</footnote>
<figure><img src="image_2.png" /></figure>

**Figure 1.** Generalised J-factors $J_S$ for the host halos of the selected targets as a function of $\epsilon_\phi$, for s-wave (left panels) and p-wave (right panels) annihilations; all calculations rely on the Eddington (isotropic) PSDF. We recall that $\theta_{\text{int}} = 0.5^\circ$ for dSphs and dIrrs, and $\theta_{\text{int}} = R_{200}/D$ for clusters; **Top panels:** Draco (red), Sculptor (green) and Reticulum II (violet) dSphs, where solid colored lines and shaded/hatched bands represent the median and 68% confidence regions computed from 1000 samples of the DM profile parameters (see Sec. 3.1). **Middle panels:** IC10 (gray), NGC6822 (purple), and WLM (golden) dIrrs, for Burkert (solid lines) and NFW (dashed) DM density profiles (see Sec. 3.2). **Bottom panels:** Coma (cyan), Fornax (black) and Perseus (brown) galaxy clusters, for ‘hydro’ (solid lines) and ‘upper’ (dashed) DM density profiles (see Sec. 3.3).

at resonance, $J_S$ is roughly proportional to $\bar{v}^{-2}$ and follows a $1/n^2$ power law (where $n$ is the integer defining each resonance, see Sec. 2.1). Between resonances, $J_S \propto 1/\epsilon_\phi$
and does not depend on the velocity.

p-wave annihilation (right panels in Fig. 1).

* for $\epsilon_{\phi} \lesssim \epsilon_{\phi}^{*}(\ll 1)$, $\bar{S} \propto 1/\epsilon_{v}$ similar to the s-wave case and is also independent of $\epsilon_{\phi}$, so $J_{S,p}$ also features a plateau in this regime. The p-wave plateau is lower by a factor $\alpha_{D}^{2}$ compared to the s-wave case (left-hand side of the curves);

* for $\epsilon_{\phi} \gg 1$, as for the s-wave, there is no enhancement and $J_{S}$ boils down to the standard p-wave $J$ (right-hand side of the curves);

* for $\epsilon_{\phi}^{*} \lesssim \epsilon_{\phi} \lesssim 1$, this is also the resonant (saturation) regime: at resonances, $\bar{S}$ is independent of $\epsilon_{v}$ and $J_{S} \propto 1/\epsilon_{\phi}^{3}$; however, between resonances, the Sommerfeld factor is independent of $\epsilon_{v}$ so $J_{S}$ is shaped by the p-wave velocity-dependence $\langle v_{\text{rel}}^{2} \rangle \sim \bar{v}^{2}$. Because of this dependence, p-wave annihilation resonances are more clear-cut for objects with a low characteristic velocity, like dSphs (top right panel), compared to dIrrs and galaxy clusters (middle and bottom right panels).

4.2 Uncertainties from the mass modelling

Systematic errors on the parameters describing the smooth DM density profiles (see Sec. 3) translate into systematic uncertainties on the generalised *J*-factors, which we briefly discuss quantitatively in the following.

DSphs (top panels of Fig. 1). As described in Sec. 3.1, for dSphs we compute the uncertainty on the generalised *J*-factors from the posterior distribution on $J_{S}$ obtained from the kinematic analysis of [97]. For s-wave annihilation and in the absence of any Sommerfeld enhancement, the uncertainty on the reconstruction of the DM profile parameters leads to factors of a few for classical dSphs (Draco and Sculptor), and a factor $\sim 20$ for the Reticulum II ultra-faint dSph (top-left panel of Fig. 1). This is the same in the saturation regime, off resonance, where the Sommerfeld factor does not depend on the velocity. Yet, in the Coulomb regime, and at resonance peaks in the saturation regime, the additional velocity dependence goes in the opposite direction with respect to $\rho^{2}$. For instance, for a given value of the scale radius $r_{-2}$, a larger value of $\rho_{-2}$ gives larger $\rho^{2}(r)$ but at the same time a larger typical velocity which enters the Sommerfeld factor through $\bar{v}^{-1}$ or $\bar{v}^{-2}$, leading to a reduction of the generalised *J*-factors. This leads to 68% uncertainty bands that are typically smaller in the Coulomb regime — e.g., less than an order of magnitude for Reticulum II — than in the no-Sommerfeld case.

For p-wave annihilation (top-right panel of Fig. 1), the uncertainty band spans about two orders of magnitude for Draco in the no-Sommerfeld regime, but ‘only’ one for Sculptor. For Reticulum II, the uncertainty reaches almost four orders of magnitude, owing to the loose kinematic constraints that affect both $\rho^{2}$ and $\langle v_{\text{rel}}^{2} \rangle$. In the Coulomb regime, the uncertainty on $\rho^{2}$ is again balanced by the $\bar{v}^{-1}$ dependence, leading to small 68% bands for all dSphs in our sample.

DIrrs (middle panels of Fig. 1). For these objects, we bracket the systematic error on the generalised *J*-factors by considering the Burkert and NFW mass models obtained from fitting rotation curve data (see Sec. 3.2). This error is encoded in the ratio $R_{J_{S}}(\epsilon_{\phi}) = J_{S}^{(\text{Burkert})}/J_{S}^{(\text{NFW})}$, which differs for the s-wave and p-wave cases. First, for NGC6822 and WLM, the two DM profiles are almost degenerate for the rotation curve fits, so that all
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of generalised *J*-factors for the smooth (host) halo of selected targets as a function of $\epsilon_\phi$ for *s*-wave (left panel) and *p*-wave (right panel) annihilations. We reproduce here some of the curves shown in Fig. 1, highlighting as shaded/hatched areas the uncertainties estimated from mass modelling uncertainties. The selected objects are an ultra-faint dSph (Reticulum II), a classical dSph (Draco), a dIrr (IC10), and a galaxy cluster (Fornax). We recall that $\theta_{\text{int}} = 0.5^\circ$ for dSphs and dIrrs, and $\theta_{\text{int}} = R_{200}/D$ for clusters.</figcaption></figure>

associated ratios are close to one (for all $\epsilon_\phi$ values). However, for IC10, the best fit using the NFW profile differs more appreciably from the Burkert one: (i) in the Coulomb regime (small $\epsilon_\phi$ values), $R_{JS} \sim 3$ for both the *s*- and *p*-wave cases; (ii) in the saturation regime off-resonance (intermediate $\epsilon_\phi$), and in the standard regime with no enhancement (large $\epsilon_\phi$), $R_{JS} \sim 5$ for *s*-wave while $R_{JS} \sim 10$ for *p*-wave; (iii) on resonance, $R_{JS} \sim 1$ for *s*-wave while $R_{JS} \sim 5$ for *p*-wave. It should be noted that for *p*-wave on resonances, $J_S$ does not depend on the velocity, thus the difference between Burkert and NFW is the same as for *s*-wave without Sommerfeld enhancement.

**Galaxy clusters (bottom panels of Fig. 1).** In this case, the uncertainties on the generalised *J*-factors are related to uncertainties on the derived *X*-ray masses and the scatter on the mass-concentration relation. This allows to define a lower and upper bound on the modelling of the DM density profiles (see Sec. 3.3). The ratio $R_{JS}(\epsilon_\phi)$ of these two bounds is $J_S^{\text{(hydro)}}/J_S^{\text{(lower)}}$ for Coma and Perseus, and $J_S^{\text{(upper)}}/J_S^{\text{(hydro)}}$ for Fornax. We find, in the *s*-wave case, that $R_{JS} \sim 3-4$ in all regimes for all the clusters in our sample. We have otherwise in the *p*-wave case $R_{JS} \sim 3$ in the Coulomb regime (small $\epsilon_\phi$ values), and 5-6 in both the inter-resonance saturation regime (intermediate $\epsilon_\phi$ values) and the pure *p*-wave case with no enhancement (large $\epsilon_\phi$ values).

## 4.3 Impact of uncertainties on the ranking of targets

From the above discussion, we conclude that systematic errors — that stem from the data-driven modelling of the smooth DM density profile — have a strong impact on the generalised *J*-factors. As such, they can affect the hierarchy of targets according to their potential for γ-ray DM searches. In contrast, as discussed in App. A, the $\mathcal{O}(10\%)$ uncertainty on the PSDF itself — in particular how the anisotropy of the velocity distribution is accounted for — has little impact on the generalised *J*-factors, and does not affect the ranking of targets.

In Fig. 2, we compare the *J*-factors (as a function of $\epsilon_\phi$) for a selection of representative objects among the target classes considered in this work, namely Reticulum II (ultra-faint
dSph), Draco (classical dSph), IC10 (dIrr), and Fornax (galaxy cluster). To ease the comparison, we highlight our estimated uncertainties as shaded/hatched areas, although we remind that the different bands neither have the same origin nor the same statistical meaning (see discussion in Sec. 4.2). For s-wave annihilation (left panel), the uncertainty bands overlap in the regime with no Sommerfeld enhancement (large $\epsilon_\phi$ values) for Fornax and IC10, while this overlap disappears and is replaced by a gap for decreasing $\epsilon_\phi$; there is for instance almost a factor 10 difference between the Fornax lower edge (black-hatched band) and IC10 upper edge (gray-hatched band) in the Coulomb regime (small $\epsilon_\phi$ values). The situation is qualitatively similar for Draco and Reticulum II. For p-wave annihilation, the uncertainty bands overlap for the four representative targets in the regime of no enhancement (large $\epsilon_\phi$ values), whereas a clear hierarchy also appears when Sommerfeld enhancement becomes important, especially in the Coulomb regime (small $\epsilon_\phi$ values) due to the $1/v$ dependence.

As a conclusion from this section, which focused only on the signal from the host DM halo, we see that for s-wave annihilation, dSphs represent (in all regimes) the most promising targets for $\gamma$-ray searches in terms of generalised J-factors, even accounting for modelling uncertainties. Yet, the situation is less clear-cut for p-wave annihilation in the regime of no enhancement (large $\epsilon_\phi$ values). Indeed, in the latter case, within the uncertainties, some galaxy clusters can become the best targets. Nevertheless, as already highlighted in the literature (mostly for the standard J-factor calculations), accounting for DM substructures in all these different targets may change these conclusions. We discuss and detail in the next section how DM substructures are expected to boost the annihilation signal and impact the computation of the generalised J-factors.

# 5 Generalised J-factors with substructure boost

A fraction of the DM in halos is in the form of subhalos, which can boost the annihilation signal (compared to the case in which all the DM mass is smoothly distributed within the main halo). While the impact of these substructure boosts has been discussed extensively in the literature for the ‘classical’ J-factors (see, e.g. [30–35]), they have been discussed with lesser details in the context of generalised J-factors [15, 36–39].

The results derived in this section rely on the general formalism and methodology presented in Sec. 2.3, and our calculations are based on up-to-date models for both the properties of the subhalo population and the velocity distribution in each subhalo (determined by a phase space); the numerical calculations in this section (for subhalos) have also been cross-checked and validated with analytical approximations (see the companion paper, [60]). First, we discuss the boost factors obtained for our representative targets and highlight the differences observed between the s- and p-wave cases (Sec. 5.1). We then show the full calculation of the generalised J-factors for all our targets, and rank them according to their expected signals, also depending on the regime considered for s- or p-wave annihilations (Sec. 5.2). We finally show how these boosted signals compare to the ‘foreground’ DM annihilation signal coming from the smooth DM distribution in the MW, and briefly discuss the prospects for $\gamma$-ray searches (Sec. 5.3).

## 5.1 Impact of subhalos: boost factors for generalised J-factors

The generalised boost factor $B_S$, calculated for the generalised J-factor $J_S$, is given by

$$B_S = \frac{J_{S,tot}}{J_{S,host}}. \qquad (5.1)$$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Generalised boost factors $\mathcal{B}_S$ (see Eq. 5.1) as a function of $\epsilon_\phi$ for s-wave (left panel) and p-wave (right panel) annihilation. We illustrate the boosts obtained for each family of targets considered, i.e., galaxy clusters (Fornax, black lines), dIrrs (IC10, gray lines), and dSphs (Draco, red lines). The more massive the object, the more boosted the signal (see text for discussion). We recall that $\theta_{\text{int}} = 0.5^\circ$ for dSphs and dIrrs, and $\theta_{\text{int}} = R_{200}/D$ for clusters.</figcaption></figure>

In this definition, the denominator $J_{S,\text{host}}$ is the generalised J-factor for the host halo without substructures (i.e., assuming all the DM to be smoothly distributed), already calculated and presented in Sec. 4.1. The numerator $J_{S,\text{tot}}$ is the sum of the signals from the smooth halo of the host (which is now all the DM not in substructures), the population of subhalos and the cross-annihilation between host and subhalos (the latter is in general negligible compared to the sum of the other two contributions). In this definition, when the fraction of DM into substructures goes to zero, the boost goes to one, i.e., the overall signal is not boosted.

We show in Fig. 3 the generalised boost factors $\mathcal{B}_S$ as a function of $\epsilon_\phi$, for three representative objects from our three families of targets. Several striking features are observed. First, contrarily to the smooth host halo case where resonances disappear below $\epsilon_\phi \lesssim 10^{-2}$ (see, e.g., Fig. 2), resonances are present down to much smaller $\epsilon_\phi$ values here. This is because subhalos are less massive with smaller velocity dispersion, hence a smaller $\epsilon_\phi^*$ (see Eq. 4.1) below which the Coulomb regime is reached (compared to their host halo counterpart). Second, we see that a larger host is more boosted. This is a well-known feature of s-wave annihilation without Sommerfeld enhancement due to each decade in subhalo mass contributing to the annihilation at roughly the same level. When the Sommerfeld effect is included, this dependence on the host halo mass is preserved although the scaling is slightly modified. Moreover, it now extends to the p-wave case as well, and the scaling with the host mass is identical for both s-wave and p-wave. Third, a different scaling with $\epsilon_\phi$ is observed in the s-wave and p-wave case. This difference can be explained by considering which subhalos contribute most to the annihilation. We find that the s-wave signal is dominated by subhalos near the free-streaming cutoff, while the p-wave signal is dominated by subhalos near the mass scale which sets the transition between the Coulomb and saturation/resonant regime, which depends on $\epsilon_\phi$. Details and scaling relations are provided in the companion paper [60]. We stress that the resonances do not appear at arbitrarily low $\epsilon_\phi$, because subhalos cannot form below the free-streaming scale. We fixed this scale to $m_{\text{min}} = 10^{-6} \text{ M}_\odot$, which translates into a value $\epsilon_\phi^* \sim 10^{-7}$ below which all subhalos, and therefore all the DM in the object, are in the Coulomb regime. We chose to limit the x-axis to $\epsilon_\phi = 10^{-4}$ however, because lower
values have little motivations from the model-building point of view. The sensitivity of the s-wave annihilation to this minimum low-mass scale also explains why the s-wave boost is generically much larger than the p-wave boost. Indeed, for s-wave processes, the baseline and the resonant peaks have their amplitudes fixed by $\epsilon_{\phi}^{-1/1} m_{\text{min}}^{-\tilde{\alpha}_s}$, respectively,<fnref n="12" /> where $\tilde{\alpha}_s$ is some effective mass index that depends on the subhalo mass index, found positive here ($\sim 0.1$ on the baseline and $\sim 0.8$ on peaks—see [60] for details). This explains why the relative peak amplitude decreases as $\epsilon_{\phi}^2$ as $\epsilon_{\phi}$ decreases. In contrast, for p-wave processes, only the peaks have their amplitudes that scale like $\epsilon_{\phi}^{-1} m_{\text{min}}^{-\tilde{\alpha}_p}$, with $\tilde{\alpha}_p \sim 0.1$, while the baseline scales only $\propto \epsilon_{\phi}^{-1.3}$. In both the s- and p-wave cases, the overall amplitude of the boost, once the host smooth halo lies in the Coulomb regime, is further modulated by the host halo mass to some positive power, which explains the hierarchy between the different curves. All this allows to understand how changing $m_{\text{min}}$ may affect the final results.

To be more quantitative, the boost factors in the different regimes of the Sommerfeld enhancement are as follows: for large $\epsilon_{\phi}$, we have $\mathcal{B}_{S,s} \sim 1-20$ going from dSphs to galaxy clusters,<fnref n="13" /> while $\mathcal{B}_{S,p} = 1$ (no boost) for all targets. Moving down towards the saturation regime, and for the dSphs, dIrrs and galaxy clusters, respectively, we have for the first resonance $\mathcal{B}_{S,s} \sim 10^6, 10^8, 10^{10}$ and $\mathcal{B}_{S,p} \sim 1, 2, 20$.

## 5.2 Ranking of target classes

Now that we have detailed the behaviour of the generalised boost factors, we can go back to the generalised J-factors.

We show in Fig. 4 the $J_S$ values as a function of $\epsilon_{\phi}$ for the full calculation including the boost from substructures (solid lines). To our knowledge, these $J_S$ are the most complete and up-to-date estimates for such a variety of targets. For comparison purpose, we also reproduce some of the values shown in Fig. 2 for the case with no substructures (dashed lines). We note that the ratios between the pairs of solid and dashed lines in each panel are directly the boosts discussed in the previous section; we refer the reader to the details therein rather than repeating the discussion here. For brevity, it is enough to summarise the most salient features of the full calculation (solid lines): (i) $J_S$ values in the no-enhancement regime (large $\epsilon_{\phi}$) reach a plateau for both s- and p-waves, and these plateaus actually correspond to the minimum value of $J_S$ over $\epsilon_{\phi}$; (ii) the saturation regime at resonances gives the most favourable (and tremendous) signal for s-wave annihilation, but this is comparatively only mildly significant for p-wave annihilation; (iii) off-resonance and moving down towards the Coulomb regime (small $\epsilon_{\phi}$ values), the $J_S$ factor is increasing for both the s- and p-waves, but it increases less and converges faster towards a plateau in the former case.

With these results, we can now revisit our discussion on the ranking of the best targets (to either detect them or to set stringent constraints on DM particle candidates). We show in Fig. 5 a comparison between targets picked among each of the families considered in this study, namely dSphs (Draco and Reticulum II, with violet and red lines, respectively, and as

<footnote n="12">¹²Formally, as explained in the companion paper [60], the amplitudes of resonant peaks saturate at vanishingly small DM velocity, which translates into a universal unitarity cutoff mass $m_{\text{unit}} \neq m_{\text{min}}$, extremely sensitive to the DM fine structure constant, as it scales like $\alpha_D^{12}$. For $\alpha_D \sim 0.01$, we have $m_{\text{unit}} \sim 10^{-3} \text{M}_\odot > m_{\text{min}}$, which means that it is actually $m_{\text{unit}}$, still a universal parameter related to particle physics, that sets the peaks amplitudes in our calculations. Had we taken $\alpha_D \sim 10^{-3}$ instead, then $m_{\text{min}}$ would have been the peaks maker.</footnote>

<footnote n="13">¹³For s-wave, this is the regime where standard J-factor calculations are recovered, and the boost values obtained are in line with standard boost factors found in the literature, e.g. [33, 34].</footnote>
<figure><img src="image_5.png" /></figure>

**Figure 4.** Same as Fig. 1, i.e., generalised *J*-factors vs $\epsilon_\phi$ for our sample of dSphs (top panels), dIrrs (middle panels), and galaxy clusters (bottom panels), for both *s*-wave (left panels) and *p*-wave (right panels) annihilation. Two different calculations are shown: solid lines correspond to the full calculation accounting for the contribution of subhalos (see Sec. 5.1), whereas dashed lines (for comparison purpose) correspond to the ‘no-subhalos’ case already shown in Fig. 2 (we only show our ‘best’ mass modelling here).

representative examples of both classical and ultra-faint dSphs), dIrrs (IC10; gray lines), and galaxy clusters (Fornax; black lines). As before, solid lines correspond to our final results (with substructure boost), and dashed lines show their ‘no-substructures’ counterparts just for comparison purpose (taken from Fig. 2). From this Fig. 5, we notice that we can have
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as in Fig. 2, i.e., comparison of *J*-factors for *s*-wave (left panel) and *p*-wave (right panel) annihilation, for one object in each considered target class but dSphs, for which we include both a “classical” and a “ultra-faint” dwarf. Two different calculations are shown: solid lines correspond to the full calculation accounting for the contribution of subhalos (see Sec. 5.1), whereas dashed lines (for comparison purpose) correspond to the ‘no-subhalos’ case already shown in Fig. 2 (without uncertainties here).</figcaption></figure>

a complete inversion of the standard ranking for both the *s*- and *p*-wave cases. Indeed, instead of dSphs being the best targets (for standard *J*-factor calculation), galaxy clusters can now outrank dIrrs, which themselves outrank dSphs. This inversion can happen, e.g., at Sommerfeld resonances in the *s*-wave, and in both the Coulomb ($\epsilon_{\phi} \ll 1$) and no-Sommerfeld enhancement ($\epsilon_{\phi} \gg 1$) regimes for *p*-wave.

We recall that this inversion in the ranking of targets arises because of the role of substructures, that boost differently the different target classes. While the exact value of these boosts may vary by a factor of a few (due to uncertainties in the subhalo distribution, abundance and structural parameters), the trend of these boosts is not expected to change significantly. We also conclude that, for the generalised *J*-factors, the mass modelling uncertainties of the host halos play a subdominant part in almost all regimes: such uncertainties only impact the ranking in the no-enhancement regime ($\epsilon_{\phi} \gtrsim 1$) in *p*-wave annihilation (see discussion in Sec. 4.2 and also Fig. 2).

## 5.3 Comparison to previous works and prospects for γ-ray DM searches

With the results in previous sections, we can now draw some conclusions regarding the selection of the best targets for γ-ray DM searches, depending on the Sommerfeld regime considered (in *s*- or *p*-wave). Alternatively, given some γ-ray observations, we can also highlight the regimes where DM candidates are expected to be constrained the most.

**Comparison with previous studies.** It is interesting to compare our findings to what was previously obtained in other works. The largest body of results in the literature is for the standard *J*-factors — corresponding to $\epsilon_{\phi} \gg 1$ in the *s*-wave case. In this regime, which boils down to the calculation of the boost factors, our results agree with previous determinations; this is not a surprise since we recall that we rely, for the most part, on input ingredients taken from some of our previous works (e.g., [97] for dSphs, [117] for dIrrs, and [146] for galaxy clusters). As for the calculation of generalised *J*-factors, there is no study to compare to for dIrrs. For galaxy clusters, to our knowledge, the only previous study is that of [54], where the
authors do not directly calculate the generalised J-factors, but show limits on DM candidates from the observation of Fornax; hence it is difficult to make comparisons. There are several generalised J-factor calculations in the literature for dSphs, mostly ignoring substructures. Comparing the results obtained in the ‘no-substructure’ case (Fig. 2, top panels), we find our results to be comparable to those of [14, 43–46] for s-wave in all regimes.<fnref n="14" />

The interplay between subhalos and velocity-dependence has been investigated in several studies [36, 37, 53, 157–162]. Comparison with our results is difficult in most cases as alternatives targets and different regions of the parameter space are considered. The subhalo models are also quite different from the one we have used. In [37] the author performed a calculation of the subhalo boost factor in the presence of Sommerfeld enhancement, focusing on the s-wave case and dSph-sized hosts. The author found results that are qualitatively similar to ours: the resonant regime extends to very low values of $\epsilon_\phi$ when subhalos are considered, and the boost factor can reach extremely high values. Quantitatively, the boost factors in [37] seem to be higher than ours by one or two orders of magnitude. A possible reason for this discrepancy is the subhalo mass-concentration relation used in [37] which leads to subhalos that are much denser, and therefore over-annihilate, compared to what has been found in more recent numerical simulations (see [32] for the mass-concentration we have used instead).

**Angular extension of the signal.** The morphology or, for our targets, the radial dependence of the γ-ray signal is directly linked to the underlying emission processes and source spatial distribution (here, annihilations in the smooth halo and substructures). As advocated in past studies on standard J factors, the angular extension of the signal could be used to identify decaying from annihilating DM in dSphs [163] (objects in which boost factors are mostly irrelevant), or to disentangle CR-induced from DM-induced γ-ray signals in galaxy clusters [138, 145] (objects in which substructures both boost and enlarge the size of the object on the sky for annihilating DM; see also [146]). This reasoning has been further developed in the context of velocity-dependent annihilations in [161, 164], where the authors discuss how the radial dependence of the signal could help identifying the underlying particle physics model. However, these studies mostly focus on a single halo (although [161] briefly comments on the consequences for a distribution of subhalos), whereas we have shown that substructures may be important even for dSphs in some specific regimes. Another difference is that we perform a full numerical calculation while [161, 164] rely on analytical approximations.

Figure 6 shows the differential J-factor $\tilde{J}_S(\psi)$ defined as

$$J_S = \int_0^{\theta_{\text{int}}} d(\cos \psi) \tilde{J}_S(\psi) \qquad (5.2)$$

without (dashed lines) or accounting for (solid lines) substructures, as a function of the angle $\psi$ from the target centre. We recall that the relevant angular size of the γ-ray signal will be a combination of the physical size of the target, the distance to the observer, the steepness of the DM distribution, and, also for generalised J-factors, the velocity distribution profile. The latter is illustrated in Fig. 6 with the calculation of the signal without subhalos (dashed lines). Indeed, for standard J-factors, i.e., $\epsilon_\phi \gg 1$ in the s-wave case (top left panel), the typical radial extension of a given target boils down to $\psi_s = \arctan(r_s/D)$ (most of the emission is within $r_s$), i.e., 1.4° for Draco, 0.15° for IC10, and 0.5° for Fornax (using numbers taken from

<footnote n="14">¹⁴Note that, in some cases, a rescaling of $\alpha_D$ is needed to perform these comparisons.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Differential $J_S$-factor as a function of the scaled angle $\psi/\psi_s$ from the centre of the object, for three targets belonging to our three classes of interest. To ease the comparison, all profiles have been normalised to unity at their centre. The top row displays the results without Sommerfeld enhancement ($\epsilon_\phi \gg 1$), while the bottom row shows the results for $\epsilon_\phi \simeq 0.011$ (off resonance). The results for s-wave are shown in the left panels and those for p-wave are in the right panels. The reference angular size is $\psi_s = 1.4^\circ/0.15^\circ/0.5^\circ$ for Draco/IC10/Fornax.</figcaption></figure>

Tables 1, 2, and 3 respectively). These numbers compare well to the radial extensions seen in the top left panel of Figure 6, but we already see the impact of the velocity distribution function comparing the top left and bottom left dashed curves (no subhalos, no Sommerfeld enhancement) for s- and p-wave respectively. These radial extensions are also slightly different in the Sommerfeld regime (bottom panels).

The relatively small differences observed between different particle physics models for a smooth halo (without subhalos) are qualitatively similar to those highlighted in [164]. Yet, the presence of subhalos (solid lines) on $\tilde{J}_S$ is strongly model dependent. Subhalos have almost no impact for $\epsilon_\phi \gg 1$ in the p-wave case (right panels in Fig. 6) and a maximal impact in the s-wave resonant regime (bottom left panel), where a factor of ten increase of the typical angular size is observed for Fornax (compare the dashed and solid black lines). All in all, these results reinforce the case for the use of the angular dependence (or radial extension) of the $\gamma$-ray signal as a tool to discriminate between different particle models (if a DM signal is seen). Further study is necessary to decide/optimise which combinations of different targets are best to discriminate among particle physics models. Also, it is not clear whether degeneracies (and uncertainties) between both the DM and particle physics modellings would
prevent the use of this strategy in some cases (not to mention the observational challenge to characterise the spatial morphology of a putative γ-ray signal: current γ-ray telescopes possess typical angular resolutions of about one to few tenths of degree).

**Galactic DM foreground.** So far, we have considered γ-ray signals from isolated targets. However, the diffuse γ-ray emission originated from DM annihilations happening in the MW halo itself can be a sizable foreground when searching for γ-rays from various targets (the DM extragalactic diffuse signal also adds up to this foreground but will not be considered here; see, e.g., [165] for a comparison between both DM diffuse components in terms of their intensity flux). To optimise the detectability of such targets in γ-rays, their DM-induced signal should ideally lie above this Galactic DM foreground (typically, for a field of view corresponding to the angular resolution of the detector).

For completeness, we compute the astrophysical factor of the MW halo for a solid angle defined by $\Delta\Omega = 2\pi(1 - \cos\theta_{\text{int}})$ around the line of sight (l.o.s.) at an angle $\psi_{\text{GC}}$ from the Galactic centre, $J_S(\psi_{\text{GC}})$. We use Eq. (2.13) for the calculation of the annihilation. For the MW mass model, we use an NFW profile with $\rho_0 = 8.5 \times 10^6 \, \text{M}_\odot/\text{kpc}^3$ and $r_s = 19.6 \, \text{kpc}$ [166]. For the subhalos, we use again the SL17 model. An important difference when studying subhalos in the MW compared to the other targets is the role of baryons in shaping their distribution. In SL17, both the smooth tidal stripping induced by the baryonic potential and the gravitational shocking induced by the stellar disk are taken into account.<fnref n="15" /> The calculation for the contributions of the MW smooth halo and subhalos (see Sec. 2.3) is otherwise similar to that of the other targets, the main differences being that (i) the observer is now sitting inside the host halo, and (ii) the baryonic potential can no longer be neglected (see App. B on how it is accounted for in the modelling).

Figure 7 shows the resulting MW $J_S$ factors for s- and p-wave (left and right panels, respectively), without (dashed lines) or with (solid lines) subhalos. From top to bottom, going from a l.o.s. slightly offset from the Galactic centre ($\psi_{\text{GC}} = 0.1^\circ$ and $3^\circ$, top and middle panels) and moving towards the anticentre ($\psi_{\text{GC}} = 180^\circ$, bottom panel), we observe in all regimes, as expected, that the $J_S$ factors decrease. Also, as expected in the velocity-independent case ($\epsilon_\phi \gg 1$ in s-wave), subhalos only boost the signal towards the anticentre (the signal is dominated by the smooth halo towards the halo centre). Note that there is no boost from subhalos in the no-Sommerfeld enhancement regime ($\epsilon_\phi \gg 1$) for p-wave (right panels). This larger impact of subhalos, away from the Galactic centre, is also recovered in the s- (left panels) and p- (right panels) wave at intermediate and small $\epsilon_\phi$. Actually, the pattern is very similar to the one seen in Fig. 5, where the presence of subhalos, down to the cutoff mass, leads to larger peaks at resonances and growing boost with decreasing $\epsilon_\phi$ (see discussion in Sec. 5.1 for more details).

Based on numerical simulations, [51] recently showed that the predicted $J_S$ factor from the smooth halo of the MW is very sensitive to the DM velocity distribution function, comparing predictions from DM-only or hydrodynamical simulations. These authors also concluded that the impact of subhalos is subdominant in their work, given that their simulation only resolves the largest subhalos. Our results show that the distribution of subhalos down to the smallest masses is actually critical to correctly predict the MW signal on resonances and in the regime $\epsilon_\phi \ll 1$.

<footnote n="15">¹⁵We have checked that baryonic tidal effects are completely negligible in the other targets, even in galaxy clusters which have a sizeable baryonic content.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Generalised MW *J*-factors (integrated over $\theta_{\text{int}} = 0.5^{\circ}$) vs. $\epsilon_{\phi}$ for s-wave (left panels) and p-wave (right panels) annihilation, in various directions $\psi_{\text{GC}}$ from the Galactic centre: from top to bottom, $\psi_{\text{GC}} = 0.1^{\circ}, 3^{\circ}$, and $180^{\circ}$ (anti-centre). Solid lines correspond to the full calculation accounting for the contribution of subhalos (see Sec. 2.3), whereas dashed lines correspond to the ‘no-subhalos’ case.</figcaption></figure>

**Summary view of all targets against the Galactic DM foreground.** Thanks to the above calculation, we can now assess the contrast between the DM signal from all our targets and that from the MW DM foreground. This contrast is shown as a function of $\psi_{\text{GC}}$ (angle between the l.o.s. and the Galactic centre; the Galactic halo is spherically-symmetric) in Fig. 8. The panels in this figure are for two regimes of the s-wave (left panels) and p-wave
(right panels) annihilations, namely the no-Sommerfeld enhancement regime at $\epsilon_{\phi} \gg 1$ (top panels), the saturation regime off-resonance at $\epsilon_{\phi} \simeq 10^{-2}$ (middle panels) and the saturation regime on-resonance also at $\epsilon_{\phi} \simeq 10^{-2}$ (bottom panels). We recall that $\epsilon_{\phi} \gg 1$ corresponds to a mediator mass comparable or larger than the DM mass, while $\epsilon_{\phi} = 10^{-2}$ corresponds to a light mediator (e.g., $m_{\phi} = 1$ GeV for $m_{\chi} = 10$ TeV and our choice of $\alpha_{D} = 10^{-2}$). On each plot, the solid line corresponds to the MW $J_{S}$ values (see previous paragraph), and the symbols represent our various target values,$^{16}$ $J_{S_{\psi_{GC}}}$. Both the MW and the target values are calculated for an integration region $\theta_{\text{int}} = 0.5^{\circ}$.${}^{17}$

The various panels of Fig. 8 illustrate that: (i) in most regimes of s-wave annihilation (left panels), dSphs remain the most promising targets, whereas dIrr and galaxy cluster signals, $\sim 10$ times lower, remain interesting and complementary targets; (ii) in some regimes of p-wave annihilation (e.g., top right panel), an inversion of the ranking is observed and galaxy clusters become the most promising targets, while in some other cases (bottom right panel), dIrrs can shine as bright as some dSphs. In addition, it is interesting to note that for this value of $\theta_{\text{int}} = 0.5^{\circ}$, most targets in most regimes outshine the MW DM foreground. Yet, the signal contrast between the former and the latter is significantly smaller for the p-wave case: in the no-Sommerfeld enhancement regime (top right panel), the MW is even above all target signals but for the two brightest galaxy clusters. The non-trivial dependence of this contrast with the integration region ($\theta_{\text{int}}$) is further discussed in App. D. Overall, this variety of scenarios illustrate that devising the optimal signal region to search for DM signals (or to set limits on DM candidates) is not a trivial task and needs to be studied in detail in a case-by-case basis, as it is sensitive not only to both the individual target and target class but also to the specific particle physics model considered. Such a search strategy should also need to account for the Galactic DM foreground; for completeness, it should also account for the extragalactic DM background (see, e.g., [165]).

**Relevance to particle physics.** The results discussed here are relevant for particle-physics models involving mediators much lighter than the DM particle. Such a hierarchy appears in minimal setups such as Minimal Dark Matter [21, 167, 168] but also in models tailored to explain astrophysical observation such as the intense gamma-ray emission observed in the Galactic centre by Fermi-LAT [169] or the PAMELA/AMS-02 positron flux [170–172]. It is in general an expected feature of most models involving multi-TeV WIMPs [18, 19], for example in minimal models built from electroweak n-uplets [173]. In this case, DM annihilation can proceed through s-wave processes. There are also other models based on dark-sector extensions of the standard model (with a dark sector possibly secluded from the visible one), in which DM particles can in principle have masses down to tens of GeV [174, 175]. For example, for fermionic DM endowed with dark scalar (self-)interactions, p-wave annihilation is a natural outcome, providing a phenomenological setup that can therefore be tested through indirect searches from the targets we have investigated here. Our results, however, cannot always be *directly* applied to any specific particle-physics model, as one would first need to map that model onto our simplified parameter space (which only contains three parameters: $\alpha_{D}$, $m_{\chi}$ and $m_{\phi}$). Our results still provide decent estimates in particular when there is only one

<footnote>$^{16}$Each target position $\psi_{GC}$ is computed from its Galactic longitude $l$ and latitude $b$ as $\psi_{0,\text{target}} = \arccos(\cos l \cos b)$.</footnote>

<footnote>$^{17}$This value is motivated by the typical angular extensions of our targets. We note though that this integration angle may not enclose their total DM signals in some cases, especially for galaxy clusters. Yet, most of the signal will be still originated from this inner $0.5^{\circ}$. Also, from the data analysis point of view, more extended objects are more difficult to deal with, thus our compromise in terms of the chosen $\theta_{\text{int}}$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of generalised *J*-factors from all our targets (symbols) against the MW DM foreground (blue and red lines, depending on the integration angle) as a function of the angle $\Phi_0$ (away from the Galactic centre) for *s*-wave (left panels) and *p*-wave (right panels); the top panels are for $\epsilon_\phi \gg 1$, i.e., no Sommerfeld enhancement, the middle panels are for $\epsilon_\phi \simeq 10^{-2}$ in the saturation regime off-resonance and the bottom panels are on-resonance also at $\epsilon_\phi \simeq 10^{-2}$. $J_S$ values for all dSphs and dIrrs are calculated for an integration angle $\theta_{\text{int}} = 0.5^\circ$ while clusters are integrated up to the virial radius ($\theta_{\text{int}} = 2.2^\circ/1.1^\circ/1.2^\circ$ for Fornax/Coma/Perseus). See text for details.</figcaption></figure>

light mediator at play for the self-interactions. In turn, they only provide order-of-magnitude estimates in the presence of more complex dark sectors, involving for instance several mediators, provided one can identify a very few dominant interaction(s). Irrespective of any
mapping to specific models, we still expect our main result (i.e., very large subhalo boost factors and a model-dependent target hierarchy) to hold for any model featuring a similar velocity dependence of the annihilation cross-section.

# 6 Summary and conclusions

In this work, we have performed a comprehensive study of various classes of astrophysical targets (almost always discussed separately in the literature) for velocity-dependent WIMP annihilations. In particular, we computed their astrophysical *J*-factors entering in the computation of the annihilation flux and ranked them under different velocity-dependent scenarios. The main novelties of our generalised *J*-factor calculations are the following: (i) in addition to the ‘standard’ *s*-wave annihilation case, we also considered Sommerfeld-enhanced *p*-wave annihilation; (ii) for all classes of considered targets — dSphs, dIrrs, galaxy clusters, and even the MW —, we self-consistently derived the phase-space distribution function (PSDF) from host halo DM density profiles, accounting for annihilation boost factors from DM substructures (whose PSDFs were also self-consistently derived). It must be noted that very few studies carried out the calculation of the boost for dSphs in this velocity-dependent context, almost none did for galaxy clusters, and this is the first work where this was discussed for dIrrs.

The two most important and probably surprising results, obtained from a case study focusing on a few selected targets among the many available dSphs, dIrrs, and galaxy clusters, are the following. First, substructure boost factors can reach several orders of magnitude on-resonance for *s*-wave annihilation, and also in the Coulomb regime for both *s*- and *p*-wave annihilation; these large boost values are even present for dSphs. All these results are supported and cross-checked by analytical formulae derived (in various regime of the Sommerfeld enhancement) in a companion paper [60]. Second, the standard hierarchy of the most promising classes of targets for indirect DM searches (where, typically, dSphs rank first) can be drastically modified in the presence of both velocity-dependent annihilation and substructure boost. The most striking case is for *s*-wave on resonances and for *p*-wave in the no-Sommerfeld enhancement regime, where galaxy clusters can outshine all other classes of targets. This is a robust result: only in the no-Sommerfeld-enhancement regime we found that uncertainties in the modelling of the DM distribution may significantly impact the resulting ranking between dSphs, dIrrs, and galaxy clusters. We find that uncertainties related to velocity anisotropies are less important than those related to mass modelling. We stress that the modelling of tidal interactions experienced by subhalos, which is usually a source of large uncertainties, can be safely ignored here. This is because the boost factor is essentially set by subhalos in the outer regions of the targets which are not subject to strong tidal fields.

In this analysis, we have also inspected the spatial morphology of the velocity-dependent annihilation signal from various targets. We showed that subhalos could significantly en- large the angular size of the signal (up to 10 times for galaxy clusters in the Sommerfeld- enhancement regime). This provides both prospects to identify the underlying particle physics model (if DM-induced $\gamma$-rays are observed in several target classes), but it also brings complications in doing so owing to possible degeneracies between the particle physics parameter space and the many still uncertain DM distribution properties (subhalos, velocity). We have also investigated the signal contrast between the considered targets and the MW DM foreground that is present along the line of sight. The self-consistent calculation of the
MW signal showed that the presence of subhalos strongly boosts the signal in both s-wave
and p-wave annihilation (except in the no-Sommerfeld-enhancement regime): accounting for
these subhalos is critical on-resonance in all directions (Galactic centre or anticentre), and
is also critical especially towards the anticentre in the Coulomb regime. Though we found a
significant number of our targets to exhibit fluxes above the MW foreground, their specific
contrast possesses a non-trivial dependence on the signal integration angle and on the consid-
ered particle-physics model. Indeed, in some of the studied velocity-dependent regimes and
integration angles, some of our targets appear well below the level of the diffuse Galactic DM
signal, this way very likely complicating a potential γ-ray data analysis (typically focused
and optimized for point-like sources). On the other hand, a calculation of the extragalactic
diffuse DM component, similarly to what we did for the Galactic signal, may prove necessary
in some specific regimes for which targets outshine all the other ones and the MW foreground.
Overall, our results show that the analysis strategy (targets and signal regions) should prob-
ably be adapted specifically for each of these different particle-physics regimes, in order to
optimally search for (or to set constraints on) DM with γ-ray telescopes. These refinements
would complement the toolbox of existing strategies advocated and already explored in the
literature to best track down DM signals in γ rays (e.g., stacking of a large number of dSphs
or galaxy clusters; joint analyses between Fermi-LAT and ground-based instruments; opti-
mized search of extended sources; etc.). For velocity-dependent annihilation cross sections,
and given the hierarchy dependence of the ranking of targets on the particle physics model,
a combined analysis of DM signals from different classes of targets may actually be the most
optimal way to provide the most robust and consistent constraints on DM candidates.

This study is especially relevant for most models involving mediators much lighter than
the DM particle, from minimal extensions to the Standard Model to complex dark sectors.
This is especially important as indirect searches are moving toward multi-TeV DM masses
with the advent of CTA. However, rather than covering specific examples, we have studied a
simplified setup which allowed us to derive results that should qualitatively apply to a wide
range of models.

In a forthcoming effort, we plan to generalise our calculations to a larger list of astro-
physical targets and confront our predictions to existing γ-ray data, also focusing on devel-
oping and applying the best data analysis strategy. This should allow us to provide the most
robust, up-to-date, and competitive DM limits on generic velocity-dependent annihilation
models and associated DM particle candidates.

Acknowledgments

We acknowledge financial support by the CNRS-INSU programs PNHE and PNCG, the GaDaMa ANR project (ANR-18-CE31-0006), the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement N° 860881-HIDDen. MS acknowlegdes support from Université Savoie Mont Blanc and the PNHE through the AO INSU 2019, grant “DMSubG” (PI: F. Calore). GF acknowledges support of the ARC program of the Federation Wallonie-Bruxelles and of the Excellence of Science (EoS) project No. 30820817 - be.h “The H boson gateway to physics beyond the Standard Model”. JPR work is supported by grant SEV-2016-0597-17-2 funded by MCIN/AEI/10.13039/501100011033 and “ESF Investing in your future”. MASC was also supported by the *Atracción de Talento* contracts no. 2016-T1/TIC-1542 and 2020-5A/TIC-19725 granted by the Comunidad de Madrid in Spain. The work of JPR and MASC was
additionally supported by the grants PGC2018-095161-B-I00 and CEX2020-001007-S, both funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”.

## A Uncertainties related to phase-space modelling

We complement Sec. 4 with a discussion of the uncertainty on the generalised J-factors coming from the phase-space model itself. The phase space of the DM in our selected targets is modelled assuming equilibrium and spherical symmetry of the halo as well as isotropy of velocities. Under these assumptions, the PSDF of DM particles is given by the well-known Eddington formula [68, 69]:

$$f(\mathcal{E}) = \frac{1}{\sqrt{8\pi^2}} \left[ \frac{1}{\sqrt{\mathcal{E}}} \left( \frac{d\rho_{\chi}}{d\Psi} \right)_{\Psi=0} + \int_0^{\mathcal{E}} \frac{d^2\rho_{\chi}}{d\Psi^2} \frac{d\Psi}{\sqrt{\mathcal{E} - \Psi}} \right], \quad (A.1)$$

where $\Psi(r) = \Phi(R_{max}) - \Phi(r)$ is the (positive-defined) gravitational potential and $\mathcal{E} = \Psi - v^2/2$ the energy. The choice of the radial boundary $R_{max}$ can have a large impact on the phase space in the outer parts of the halo [70], however the annihilation signal is essentially set by the central region thus we can take $R_{max} \to \infty$.

The anisotropy of the DM velocity distribution is, for all intents and purposes, unconstrained in all the gravitational systems we consider. It is therefore important to go beyond the minimal assumption of isotropy and explore the associated uncertainty on predictions of DM-induced $\gamma$-ray fluxes. In general, the anisotropy of a given component in a gravitational system is quantified by the parameter [176]

$$\beta(r) = 1 - \frac{\sigma_{\theta}^{2} + \sigma_{\phi}^{2}}{2\sigma_{r}^{2}}, \qquad (A.2)$$

with $\sigma_r$, $\sigma_{\theta}$, and $\sigma_{\phi}$ the velocity dispersions in spherical coordinates. We consider the following ansatz for an anisotropic PSDF — this ansatz makes it possible to obtain a semi-analytic extension of the Eddington formula to an anisotropy profile defined by three parameters, namely the asymptotic values $\beta_0$ and $\beta_{\infty}$ at the centre and the outskirts of the galaxy, respectively, and a characteristic angular momentum $L_0$ that sets the transition radius between both regimes [177] —:

$$F(\mathcal{E}, L) = f_{\mathcal{E}}(\mathcal{E}) \left(1 + \frac{L^2}{2L_0^2}\right)^{-\beta_{\infty} + \beta_0} L^{-2\beta_0}. \qquad (A.3)$$

Constant anisotropy models, for which $\beta(r) = \beta_0 = \beta_{\infty}$ and

$$F(\mathcal{E}, L) \equiv f_{\beta_0}(\mathcal{E}, L) = f_{\mathcal{E}}(\mathcal{E})L^{-2\beta_0}, \qquad (A.4)$$

are a subset of the varying anisotropy models of Eq. (A.3), and have been extensively discussed in the literature (e.g. [69, 178–181]).

For all the targets in our sample, for s-wave annihilation the relative differences between anisotropic phase-space models with constant negative anisotropy or radius-dependent positive anisotropy and the isotropic (Eddington) result are shown in Table 5 for the Coulomb regime and the resonances. The value of $J_{S,s}$ for negative (positive) anisotropy is systematically larger (smaller) than the isotropic result by a few tens of %.
<table><caption>**Table 5.** Summary table of the relative difference between generalised *J*-factors obtained with anisotropic PSDFs and the benchmark Eddington (isotropic) case for each class of targets — both in the Coulomb regime and on top of resonances in the saturation regime. Sommerfeld enhancement factors used to derive the *J*-factors were computed for $\alpha_D = 10^{-2}$.</caption><thead><tr><th rowspan="2">Target class</th><th colspan="2">$\frac{J_{S,s}^{\beta=-0.5} - J_{S,s}^{\beta=0}}{J_{S,s}^{\beta=0}}$</th><th colspan="2">$\frac{J_{S,s}^{\beta(r)} - J_{S,s}^{\beta=0}}{J_{S,s}^{\beta=0}}$</th></tr><tr><th>Coulomb</th><th>Resonances</th><th>Coulomb</th><th>Resonances</th></tr></thead><tbody><tr><td>dSphs</td><td>+30%</td><td>+100%</td><td>-10%</td><td>-20%</td></tr><tr><td>dIrrs</td><td>+10%</td><td>+40%</td><td>-</td><td>-</td></tr><tr><td>Clusters</td><td>+20%</td><td>+50%</td><td>-10%</td><td>-30%</td></tr></tbody></table>

For *p*-wave annihilation, the dependence of $J_{S,p}$ on velocity is the same as for *s*-wave in the Coulomb regime, so the uncertainty from the unknown anisotropy is the same in both cases. However, in the saturation regime the dependence on velocity is different from the *s*-wave case. More specifically, there is no dependence on the velocity distribution at resonances, while in between resonances and in the high-$\epsilon_\phi$ regime in which there is no Sommerfeld enhancement — in both cases $J_{S,p} \propto \langle v_{rel}^2 \rangle$ — the relative difference between anisotropic and isotropic models is $\sim 10 - 20\%$.

Uncertainties from the phase-space model itself are therefore negligible with respect to uncertainties from mass modelling discussed in Sec. 4, and they do not affect in any sizable way either the values of the generalised *J*-factors or the ranking of targets. However, as discussed in Sec. 5, the DM substructure boost has a much more significant impact on generalised *J*-factors, and is the main source of theoretical uncertainties.

## B The SL17 subhalo population model

This section gives a short introduction on the SL17 subhalo population model [35], which was built to address DM searches in general (including the search for subhalos themselves), as flexible as possible to be applied to a diversity of DM candidates and to easily account for changes in relevant cosmological parameters—see also Refs. [71–74]. This model is analytical in its formulation, but semi-analytical in practice (numerical integrations, iterations, or interpolations are necessary). Assuming a global host halo profile $\rho_{host}$, the model predicts how DM distributes itself between a smooth component and a subhalo component, depending on the host halo profile and its baryonic component. To allow for fast semi-analytical calculations, the model is based on three main approximations: (i) spherically symmetric DM components, (ii) circular orbits for subhalos (i.e. positions are defined by radial distances to the host’s centre), (iii) subhalos are independent from each other. The model also includes gravitational tides as sourced by the different components of the host structure (including baryons), which induce a spatial dependence of the subhalo properties and makes the subhalo population model specific to all sorts of host halos (in particular to hosts constrained by observational data). The obtained subhalo distribution is not a fit extrapolated from cosmological simulations, and can consistently cover an arbitrary subhalo mass range. Still, the model allows to qualitatively recover and understand results of cosmological simulations, like the non-trivial spatial distribution of subhalos and the spatial evolution of their structural properties (e.g. “antibiasing” spatial distributions, spatial dependence of the concentration
and mass functions, etc. [30, 34, 182–185]). Other analytical or semi-analytical approaches to subhalo population models developed in a broader range of contexts can be found in e.g. Refs. [26, 186–189]. Monte-Carlo codes can also be used, e.g. Refs. [190, 191].

The SL17 model depicts the subhalo population of an arbitrary host halo through probability density functions (PDFs). The subhalo population is bound to be a component of the host halo, consistently with Eq. (2.16). Assuming that subhalos are all independent, the full population can be described with three parameters only: their virial (cosmological) mass $m = M_{200}$ and concentration $c = c_{200}$, and their distance to the centre of the host R. If subhalos were hard spheres, they would simply track the overall DM potential in the host, as “particles” do in N-body simulations. Thus, their population could be described by a separable parametric phase-space PDF, because the mass and concentration distributions would not depend on position. This is the starting point of the model, which assumes that the host halo builds up from the aggregation of hard spheres endowed with position, mass and concentration PDFs. The global initial PDF can be written as

$$ \frac{d^5 \mathcal{P}_{\text{sub}}^{\text{ini}}(m, c, R)}{d^3 \vec{R} \, dm \, dc} = \frac{d^3 \mathcal{P}_V(R)}{d^3 \vec{R}} \times \frac{d \mathcal{P}_m(m)}{dm} \times \frac{d \mathcal{P}_c(c, m)}{dc}, \quad (\text{B.1}) $$

i.e. a product of separable PDFs. The number of subhalos before tidal stripping is turned on can be predicted from first principles in a given cosmological framework, as will be explained below. Then, when tidal effects are plugged in, part of the DM initially confined into subhalos is redistributed as a smooth component. This tides intimately depend on the detailed distributions of the various components of the host halo. This induces a mass loss for subhalos, whose efficiency is position-dependent. The spatial dependence of tidal losses translates into a spatial dependence in the subhalo mass function. When tidal effects become disruptive, this turns into a selection in concentration space (more concentrated objects are more resilient to tides), which also implies a spatial dependence in the concentration function. Tidal effects are generically more efficient in the central parts of the host, where the subhalo number density strongly flattens up to almost full depletion close to the centre, depending on disruption criteria. Eventually, the final global PDF of subhalos is fully intricate and not separable anymore, due to mass losses and disruption. It can formally be written as

$$ \frac{d^5 \mathcal{P}_{\text{sub}}(m, c, R)}{d^3 \vec{R} \, dm \, dc} = \frac{1}{K_{\text{tidal}}} \frac{d^3 \mathcal{P}_V(R)}{d^3 \vec{R}} \times \frac{d \mathcal{P}_m(m)}{dm} \times \frac{d \mathcal{P}_c(c, m)}{dc} \times \mathcal{T}(m, c, R), \quad (\text{B.2}) $$

where $K_{\text{tidal}} \le 1$ ensures the correct normalisation of the global PDF, and $\mathcal{T}(m, c, R)$ symbolically encodes tidal stripping and disruption. It typically assigns a tidal radius $r_t(m, c, R)$ to a subhalo of virial mass $m$, concentration $c$, and position $R$, given the properties of the host halo components; or it moves the subhalo DM to the smooth component if disrupted.

There is no well-defined way to decide whether a subhalo should be tidally disrupted. It could actually be that a tiny core survives for ever if dense enough and no central collisions with stars occur, simply due to adiabatic protection [192, 193]. A practical criterion can still be inspired from studies of cosmological simulations [194], where it was found that fixing a lower threshold to $x_t = r_t/r_s$, where $r_s$ is the scale radius of the structure, was a way to efficiently capturing tidal disruption. We define this threshold as $\epsilon_t$. Initially found around $\sim 1$ [194], it was realized more recently that numerical artifacts could strongly bias these early estimates [195, 196], and that one could expect values for $\epsilon_t$ as low as 0.01 or even less, which is consistent with the argument given just above. For definiteness, we use two types
of tidal disruption criteria:

$$ \epsilon_t = \begin{cases} 1 & (\text{fragile subhalos}) \\ 0.01 & (\text{resilient subhalos}). \end{cases} \qquad (B.3) $$

We can now give a few details about the PDFs introduced above (see Refs. [35, 74] for an exhaustive presentation). For the initial spatial PDF, we simply assume $dP_V/dV = \rho_{\text{host}}(R)/M_{\text{host}}$, where $M_{\text{host}}$ is the total DM mass of the host. Note that the final spatial distribution strongly departs from the initial one after tidal effects are activated. For the concentration PDF, we use a log-normal distribution with $\sigma_c^{\text{dec}} = 0.14$ [197–200] centred about the mass–relation concentration $c(m)$ given in Ref. [32].

In order to avoid a calibration of the mass function on simulations, as was initially done in the SL17 model to determine the total number of substructures, we have implemented a cosmological mass function from first principles instead. This allows us to potentially change the cosmological parameters or the primordial power spectrum of density fluctuations. The procedure follows previous studies [185–187, 191, 201, 202], which extracted the subhalo cosmological mass function from merger-tree algorithms [203, 204], built upon the excursion set theory of structure formation [205–207]. We reproduce the same procedure with the merger tree introduced in Ref. [204] and we recover that the mass function can be well fitted by [186, 187]

$$ \frac{dN(m, M_{\text{host}})}{dm} = \frac{1}{M_{\text{host}}} \left[ \sum_{i=1,2} \gamma_i \left( \frac{m}{M_{\text{host}}} \right)^{-\alpha_i} \right] \exp \left\{ -\beta \left( \frac{m}{M_{\text{host}}} \right)^{\zeta} \right\}. \qquad (B.4) $$

We find $\gamma_1 = 0.014$, $\gamma_2 = 0.41$, $\alpha_1 = 1.965$, $\alpha_2 = 0.57$, $\beta = 20$, $\zeta = 3.4$, roughly independent of the cosmology and of the host mass—still, we used the cosmological parameters from the latest Planck analysis [208]. Contrarily to previous works (only interested in large masses) in this fit we also constrain the low mass part of the spectrum [74, 209]. Note that the above mass function is close to a power law in mass $\propto m^{-\alpha}$ with a spectral index $\alpha \simeq \alpha_1$.

The total number of subhalos before tidal stripping can be determined by integrating the mass function in the specified subhalo mass range, [$m_{\text{min}}, m_{\text{max}}$]

$$ N_{\text{tot}}^{\text{ini}} = \int_{m_{\text{min}}}^{m_{\text{max}}} dm \frac{dN(m, M_{\text{host}})}{dm}, \qquad (B.5) $$

such that the total number of subhalos, after tidal disruption effects are plugged, is given by

$$ N_{\text{tot}} = K_{\text{tidal}} \times N_{\text{tot}}^{\text{ini}}, \qquad (B.6) $$

where $K_{\text{tidal}} \le 1$ is the normalisation constant introduced in Eq. (B.2). In practice, we use $m_{\text{max}} = 10^{-2} M_{\text{host}}$, and $m_{\text{min}} = 10^{-10} M_{\odot}$, unless specified otherwise.

Finally, note that the SL17 model can be used with any assumption for the inner subhalo profiles. In this paper, we use NFW profiles for subhalos, whose properties are completely specified by the virial mass and concentration.

## C Numerical calculation of the subhalo boost

We present some technical details related to the calculation of the quantity

$$ \underline{\rho}_{\text{S,tot}}^2 = \underline{\rho}_{\text{S,sub}}^2 + \rho_{\text{S,sm}}^2 + 2 \rho_{\text{S,sm}} \rho_{\text{sub}}. \qquad (C.1) $$
First note that the smooth contribution in the presence of Sommerfeld enhancement is

$$ \rho_{S,sm}(r) = \rho_{sm}(r) \langle \bar{S} \rangle (r) \qquad (C.2) $$

where $\langle \rangle$ is the average over the relative velocity distribution of the host halo as defined in Eq. (2.14). This, combined with the average subhalo density defined in Eq. (2.17), also enables the calculation of the cross term $2 \rho_{S,sm} \rho_{sub}$.

To compute the contribution of subhalos, we need to evaluate

$$ \xi_{S,t}(m, c, r) = \int_{x \le x_t(r, m, c)} d^3 \vec{x} \left( \frac{\rho(x)}{\rho_{\oplus}} \right)^2 \langle \bar{S} \left( \frac{v_{\text{rel}}}{2} \right) \rangle (x) \qquad (C.3) $$

where the average is now taken over the velocity distribution of the subhalo with parameters $(m, c, R)$. We then need to perform the average over $m$ and $c$ to get $\rho_{S,sub}^2$ as shown in Eq. (2.18). This last step turns out to be very computationally expensive, because an integral over $\vec{v}_{rel}$ has to be performed for each subhalo mass $m$ and concentration $c$. To speed up the calculation, we instead rely on the following approximations:

$$ \langle \bar{S} \left( \frac{v_{\text{rel}}}{2} \right) \rangle (x) \approx \bar{S} \left( \frac{\langle v_{\text{rel}}^{-2} \rangle^{-1/2} (x)}{2} \right) \quad \text{for s-wave,} \qquad (C.4) $$

$$ \langle \bar{S} \left( \frac{v_{\text{rel}}}{2} \right) \rangle (x) \approx \bar{S} \left( \frac{\langle v_{\text{rel}}^2 \rangle^{1/2} (x)}{2} \right) \quad \text{for p-wave.} $$

These approximated expressions are much faster to compute, because the velocity moments have a simple scaling with the subhalo structural parameters $\langle v_{\text{rel}}^p \rangle \propto (\rho_0 r_s^2)^{p/2}$. We compute the PSDF $F_{\text{rel}}^{\text{sub}}$ and the velocity moments using the Eddington inversion method outlined in App. A. We have checked that our approximations introduce an error of at most 30% for some specific values of the subhalo mass and the $\epsilon_\phi$ parameter, and that the accuracy is better than 10% in most of the parameter space.

## D Dependence of $J_S^{\text{target}}/J_S^{\text{MW}}$ on $\theta_{\text{int}}$ and $\epsilon_\phi$

In this paper, we chose not to assess the impact of the MW signal on the detectability of the targets considered. Nevertheless, we briefly illustrate in this Appendix what would be the optimal integration angle $\theta_{\text{int}}$ so as to have $J_S^{\text{target}}/J_S^{\text{MW}} > 1$, that is a favorable contrast between the target and the MW signals. We also illustrate whether this optimal angle depends or not on the Sommerfeld parameter $\epsilon_\phi$.

We start in Fig. 9 with the contrast (between the target and the MW) as a function of the integration angle $\theta_{\text{int}}$. All curves for all configurations show the same behavior, i.e., a decreasing contrast with a growing $\theta_{\text{int}}$. The most favorable contrast is observed for small integration angles, because the diffuse MW DM foreground $J_S^{\text{MW}} \propto \theta_{\text{int}}^2$, while most of the signal remains located in the central regions of the target so that $J_S^{\text{target}}$ is independent of $\theta_{\text{int}}$ (if $\theta_{\text{int}}$ not too small). The steepness or smoothness of the decreases observed are a non-trivial combination of the different structural parameters and relative importance of the substructures (that dominate the signal at large radii). The ordering of the curves (from larger to smaller contrast) follows the ranking established in the main text: at $\theta_{\text{int}} = 0.5^\circ$, for most configurations shown in the different panels, dSphs (solid lines) reach signal contrasts
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Generalised J-factor contrast between the targets (color-coded lines) and the MW, $J_S^{\text{target}}/J_S^{\text{MW}}$, as a function of the integration angle $\theta_{\text{int}}$ around the target position. The three line styles correspond to the three family of targets (solid lines for dSphs, dash-dotted lines for dIrrs, and dotted lines for galaxy clusters). The left (resp. right) panels show the results for s-wave (resp. p-wave) annihilations, with or without Sommerfeld enhancement (top and bottom panels respectively); the value of the Sommerfeld parameter is fixed to $\epsilon_{\phi} \approx 10^{-2}$ and is chosen off resonance. The black dash-dot-dotted line separates the region between a favorable (> 1) and unfavorable (< 1) contrast. See text for discussion.</figcaption></figure>

as high as ~10–50, but only ~1–5 for galaxy clusters (dotted lines), whereas dIrrs (dash-dotted lines) are below the MW signal, except for IC10 (dash-dotted grey line). The main difference is for p-wave annihilation without Sommerfeld enhancement (bottom right panel), where only galaxy clusters outshine the MW.

To know whether these conclusions are generic or not, we show again $J_S^{\text{target}}/J_S^{\text{MW}}$ in Fig. 10, but now as a function of $\epsilon_{\phi}$ (for $\theta_{\text{int}} = 0.5^{\circ}$). We see a strong and complicated dependence that depends on the targets considered and whether s-wave (left panel) or p-wave (right panel) annihilations are considered. The most dramatic dependence is observed for dSphs (solid lines), for which we can go from favorable (> 1) to unfavorable (< 1) contrasts; for the p-wave case in particular (right panel), the contrast for Sculptor goes from ~10 (for $\epsilon_{\phi} \sim 10^{-3}$) to ~10$^{-2}$ (for $\epsilon_{\phi} > 1$).

These results illustrate the fact that devising an optimal data analysis for the various targets considered is not simple. Beside the usual considerations about the instrument characteristics and astrophysical backgrounds, the above figures stress that the optimal signal
<figure><img src="image_11.png" /></figure>

**Figure 10.** Generalised J-factor contrast between the targets (color-coded lines) and the MW, $J_S^{\text{target}}/J_S^{\text{MW}}$, as a function of the Sommerfeld parameter $\epsilon_\phi$. The three line styles correspond to the three family of targets (solid lines for dSphs, dash-dotted lines for dIrrs, and dotted lines for galaxy clusters). We show results for s-wave (left panel) and p-wave (right panel) at a fixed integration angle $\theta_{\text{int}} = 0.5^\circ$. The black dash-dot-dotted line separates the region between a favorable (> 1) and unfavorable (< 1) contrast. See text for discussion.

region is both target- and ‘particle physics model’-dependent—at least if we wish to consistently analyse the DM signal of some targets, which can possibly lay well below the MW own DM signal [73]. This should be kept in mind when searching for DM signals or setting limits on DM candidates from non-detection.

## References

[1] J. Lavalle and P. Salati, *Dark matter indirect signatures*, *Comptes Rendus Physique* **13** (July, 2012) 740–782, [1205.1004].

[2] T. Bringmann and C. Weniger, *Gamma ray signals from dark matter: Concepts, status and prospects*, *Physics of the Dark Universe* **1** (Nov., 2012) 194–217, [1208.5481].

[3] J. M. Gaskins, *A review of indirect searches for particle dark matter*, *Contemporary Physics* **57** (Oct., 2016) 496–525, [1604.00014].

[4] M. Di Mauro and M. W. Winkler, *Multimessenger constraints on the dark matter interpretation of the Fermi-LAT Galactic center excess*, *Phys. Rev. D* **103** (2021) 123005, [2101.11027].

[5] E. Charles, M. Sánchez-Conde, B. Anderson, R. Caputo, A. Cuoco, M. D. Mauro et al., *Sensitivity projections for dark matter searches with the Fermi large area telescope*, *Physics Reports* **636** (2016) 1–46.

[6] The Fermi-LAT Collaboration and W. B. Atwood, *The Large Area Telescope on the Fermi Gamma-ray Space Telescope Mission*, *APJ* (2009), [0902.1089v1].

[7] R. K. Leane, T. R. Slatyer, J. F. Beacom and K. C. Y. Ng, *GeV-scale thermal WIMPs: Not even slightly ruled out*, *Phys. Rev. D* **98** (2018) 023016, [1805.10305].

[8] M. Doro, M. A. Sánchez-Conde and M. Hütten, *Fundamental Physics Searches with IACTs*, *arXiv e-prints* (Nov., 2021) arXiv:2111.01198, [2111.01198].

[9] CTA Consortium collaboration, *Introducing the CTA concept*, *Astropart. Phys.* **43** (Mar., 2013) 3–18.
[10] The CTA Consortium, *Sensitivity of the Cherenkov Telescope Array to a dark matter signal from the Galactic centre*, Journal of Cosmology and Astroparticle Physics **2021** (Jan, 2021) 057-057.

[11] J. Abdallah, H. Araujo, A. Arbey, A. Ashkenazi, A. Belyaev, J. Berger et al., *Simplified Models for Dark Matter Searches at the LHC*, ArXiv e-prints (June, 2015), [1506.03116].

[12] Y. Zhao, X.-J. Bi, H.-Y. Jia, P.-F. Yin and F.-R. Zhu, *Constraint on the velocity dependent dark matter annihilation cross section from fermi-lat observations of dwarf galaxies*, Phys. Rev. D **93** (Apr., 2016) 083513, [1601.02181].

[13] Y. Zhao, X.-J. Bi, P.-F. Yin and X. Zhang, *Constraint on the velocity dependent dark matter annihilation cross section from gamma-ray and kinematic observations of ultrafaint dwarf galaxies*, Phys. Rev. D **97** (Mar., 2018) 063013, [1711.04696].

[14] K. K. Boddy, J. Kumar, A. B. Pace, J. Runburg and L. E. Strigari, *Effective J-factors for Milky Way dwarf spheroidal galaxies with velocity-dependent annihilation*, Phys. Rev. D **102** (July, 2020) 023029, [1909.13197].

[15] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer and N. Weiner, *A theory of dark matter*, Phys. Rev. D **79** (Jan., 2009) 015014, [0810.0713].

[16] M. Pospelov and A. Ritz, *Astrophysical Signatures of Secluded Dark Matter*, Phys. Lett. B **671** (2009) 391-397, [0810.1502].

[17] A. Sommerfeld, *Über die Beugung und Bremsung der Elektronen*, Annalen der Physik **403** (Jan., 1931) 257-330.

[18] J. Hisano, S. Matsumoto and M. M. Nojiri, *Explosive dark matter annihilation*, Phys. Rev. Lett. **92** (2004) 031303, [hep-ph/0307216].

[19] J. Hisano, S. Matsumoto, M. M. Nojiri and O. Saito, *Nonperturbative effect on dark matter annihilation and gamma ray signature from the galactic center*, Phys. Rev. D **71** (Mar., 2005) 063528, [hep-ph/0412403].

[20] S. Profumo, *TeV gamma-rays and the largest masses and annihilation cross sections of neutralino dark matter*, Phys. Rev. D **72** (2005) 103521, [astro-ph/0508628].

[21] M. Cirelli, A. Strumia and M. Tamburini, *Cosmology and Astrophysics of Minimal Dark Matter*, Nucl. Phys. B **787** (2007) 152-175, [0706.4071].

[22] J. March-Russell, S. M. West, D. Cumberbatch and D. Hooper, *Heavy Dark Matter Through the Higgs Portal*, JHEP **07** (2008) 058, [0801.3440].

[23] J. Silk and A. Stebbins, *Clumpy cold dark matter*, Astrophys. J. **411** (July, 1993) 439-449.

[24] L. Bergström, J. Edsjö, P. Gondolo and P. Ullio, *Clumpy neutralino dark matter*, Phys. Rev. D **59** (Feb., 1999) 043506, [astro-ph/9806072].

[25] C. Calcaneo-Roldan and B. Moore, *The Surface brightness of dark matter: Unique signatures of neutralino annihilation in the galactic halo*, Phys. Rev. D **62** (2000) 123005, [astro-ph/0010056].

[26] V. Berezinsky, V. Dokuchaev and Y. Eroshenko, *Small - scale clumps in the galactic halo and dark matter annihilation*, Phys. Rev. D **68** (2003) 103003, [astro-ph/0301551].

[27] F. Stoehr, S. D. M. White, V. Springel, G. Tormen and N. Yoshida, *Dark matter annihilation in the halo of the Milky Way*, Mon. Not. Roy. Astron. Soc. **345** (2003) 1313, [astro-ph/0307026].

[28] J. Lavalle, J. Pochon, P. Salati and R. Taillet, *Clumpiness of dark matter and the positron annihilation signal*, Astron. Astroph. **462** (Feb., 2007) 827-840, [arXiv:astro-ph/0603796].

[29] M. Kuhlen, J. Diemand and P. Madau, *The Dark Matter Annihilation Signal from Galactic Substructure: Predictions for GLAST*, Astrophys. J. **686** (Oct., 2008) 262-278, [0805.4416].
[30] L. Pieri, J. Lavalle, G. Bertone and E. Branchini, *Implications of High-Resolution Simulations on Indirect Dark Matter Searches*, Phys. Rev. D 83 (Jan., 2011) 023518, [0908.0195].

[31] E. Nezri, R. White, C. Combet, J. A. Hinton, D. Maurin and E. Pointecouteau, $\gamma$-rays from annihilating dark matter in galaxy clusters: stacking versus single source analysis, MNRAS 425 (Sept., 2012) 477-489, [1203.1165].

[32] M. A. Sánchez-Conde and F. Prada, *The flattening of the concentration-mass relation towards low halo masses and its implications for the annihilation signal boost*, MNRAS 442 (Aug., 2014) 2271-2277, [1312.1729].

[33] V. Bonnivard, M. Hütten, E. Nezri, A. Charbonnier, C. Combet and D. Maurin, CLUMPY: Jeans analysis, $\gamma$-ray and $\nu$ fluxes from dark matter (sub-)structures, Computer Physics Communications 200 (Mar., 2016) 336-349, [1506.07628].

[34] Á. Moliné, M. A. Sánchez-Conde, S. Palomares-Ruiz and F. Prada, Characterization of subhalo structural properties and implications for dark matter annihilation signals, MNRAS 466 (Apr., 2017) 4974-4990, [1603.04057].

[35] M. Stref and J. Lavalle, *Modeling dark matter subhalos in a constrained galaxy: Global mass and boosted annihilation profiles*, Phys. Rev. D 95 (Mar., 2017) 063003, [1610.02233].

[36] M. Lattanzi and J. Silk, *Can the WIMP annihilation boost factor be boosted by the Sommerfeld enhancement?*, Phys. Rev. D 79 (Apr., 2009) 083523, [0812.0360].

[37] J. Bovy, *Substructure boosts to dark matter annihilation from sommerfeld enhancement*, Phys. Rev. D 79 (Apr., 2009) 083539, [0903.0413].

[38] M. Kuhlen, P. Madau and J. Silk, *Exploring Dark Matter with Milky Way Substructure*, Science 325 (Aug., 2009) 970, [0907.0005].

[39] M. Kamionkowski, S. M. Koushiappas and M. Kuhlen, Galactic substructure and dark-matter annihilation in the Milky Way halo, Phys. Rev. D 81 (Feb., 2010) 043532, [1001.3144].

[40] B. E. Robertson and A. R. Zentner, Dark matter annihilation rates with velocity-dependent annihilation cross sections, Phys. Rev. D 79 (Apr., 2009) 083525, [0902.0362].

[41] R. Essig, N. Sehgal, L. E. Strigari, M. Geha and J. D. Simon, *Indirect dark matter detection limits from the ultrafaint Milky Way satellite Segue 1*, Phys. Rev. D 82 (Dec., 2010) 123503, [1007.4199].

[42] F. Ferrer and D. R. Hunter, *The impact of the phase-space density on the indirect detection of dark matter*, JCAP 9 (Sept., 2013) 05, [1306.6586].

[43] K. K. Boddy, J. Kumar, L. E. Strigari and M.-Y. Wang, *Sommerfeld-enhanced $j$-factors for dwarf spheroidal galaxies*, Phys. Rev. D 95 (June, 2017) 123008, [1702.00408].

[44] S. Bergström, R. Catena, A. Chiappo, J. Conrad, B. Eurenius, M. Eriksson et al., *$J$-factors for self-interacting dark matter in 20 dwarf spheroidal galaxies*, Phys. Rev. D 98 (Aug., 2018) 043017, [1712.03188].

[45] M. Petač, P. Ullio and M. Valli, *On velocity-dependent dark matter annihilations in dwarf satellites*, JCAP 2018 (Dec., 2018) 039, [1804.05052].

[46] S. Ando and K. Ishiwata, *Sommerfeld-enhanced dark matter searches with dwarf spheroidal galaxies*, Phys. Rev. D 104 (2021) 023016, [2103.01446].

[47] E. J. Baxter, J. Kumar, A. B. Pace and J. Runburg, *Prospects for measuring dark matter microphysics with observations of dwarf spheroidal galaxies*, JCAP 07 (2021) 030, [2103.11646].

[48] K. K. Boddy, J. Kumar and L. E. Strigari, *The effective $j$-factor of the galactic center for velocity-dependent dark matter annihilation*, Phys. Rev. D 98 (May, 2018) 063012, [1805.08379].
[49] C. Johnson, R. Caputo, C. Karwin, S. Murgia, S. Ritz, J. Shelton et al., Search for gamma-ray emission from p-wave dark matter annihilation in the Galactic Center, Phys. Rev. D 99 (May, 2019) 103007, [1904.06261].

[50] M. Petač, Equilibrium axisymmetric halo model for the Milky Way and its implications for direct and indirect dark matter searches, Phys. Rev. D 102 (Dec., 2020) 123028, [2008.11172].

[51] E. Board, N. Bozorgnia, L. E. Strigari, R. J. J. Grand, A. Fattahi, C. S. Frenk et al., Velocity-dependent J-factors for annihilation radiation from cosmological simulations, JCAP 2021 (Apr., 2021) 070, [2101.06284].

[52] B.-Q. Lu, Y.-L. Wu, W.-H. Zhang and Y.-F. Zhou, Constraints on the Sommerfeld-enhanced dark matter annihilation from the gamma rays of subhalos and dwarf galaxies, JCAP 2018 (Apr., 2018) 035, [1711.00749].

[53] J. Runburg, E. J. Baxter and J. Kumar, Constraining Dark Matter Microphysics with the Annihilation Signal from Subhalos, arXiv e-prints (June, 2021) arXiv:2106.10399, [2106.10399].

[54] A. Abramowski, F. Acero, F. Aharonian, A. G. Akhperjanian, G. Anton, A. Balzer et al., Search for Dark Matter Annihilation Signals from the Fornax Galaxy Cluster with H.E.S.S., Astrophys. J. 750 (May, 2012) 123, [1202.5494].

[55] S. Campbell, B. Dutta and E. Komatsu, Effects of velocity-dependent dark matter annihilation on the energy spectrum of the extragalactic gamma-ray background, Phys. Rev. D 82 (Nov., 2010) 095007, [1009.3530].

[56] S. Campbell and B. Dutta, Effects of p-wave annihilation on the angular power spectrum of extragalactic gamma-rays from dark matter annihilation, Phys. Rev. D 84 (Oct., 2011) 075004, [1106.4621].

[57] J. Zavala, M. Vogelsberger and S. D. M. White, Relic density and CMB constraints on dark matter annihilation with Sommerfeld enhancement, Phys. Rev. D 81 (2010) 083502, [0910.5221].

[58] S. Hannestad and T. Tram, Sommerfeld Enhancement of DM Annihilation: Resonance Structure, Freeze-Out and CMB Spectral Bound, JCAP 01 (2011) 016, [1008.1511].

[59] J. Hisano, M. Kawasaki, K. Kohri, T. Moroi, K. Nakayama and T. Sekiguchi, Cosmological constraints on dark matter models with velocity-dependent annihilation cross section, Phys. Rev. D 83 (2011) 123511, [1102.4658].

[60] G. Facchinetti, M. Stref, T. Lacroix, J. Lavalle, J. Pérez-Romero, D. Maurin et al., Analytical insight into dark matter subhalo boost factors for Sommerfeld-enhanced s- and p-wave γ-ray signals, 2203.16491.

[61] S. Cassel, Sommerfeld factor for arbitrary partial wave processes, Journal of Physics G Nuclear Physics 37 (Oct., 2010) 105009, [0903.5307].

[62] R. Iengo, Sommerfeld enhancement: general results from field theory diagrams, Journal of High Energy Physics 2009 (May, 2009) 024, [0902.0688].

[63] T. R. Slatyer, The Sommerfeld enhancement for dark matter with an excited state, JCAP 2010 (Feb., 2010) 028, [0910.5713].

[64] J. L. Feng, M. Kaplinghat and H.-B. Yu, Sommerfeld enhancements for thermal relic dark matter, Phys. Rev. D 82 (Oct., 2010) 083525, [1005.4678].

[65] K. Blum, R. Sato and T. R. Slatyer, Self-consistent calculation of the Sommerfeld enhancement, JCAP 2016 (June, 2016) 021, [1603.01383].

[66] K. Petraki, M. Postma and M. Wiechers, Dark-matter bound states from Feynman diagrams, JHEP 06 (2015) 128, [1505.00109].
[67] K. Petraki, M. Postma and J. de Vries, *Radiative bound-state-formation cross-sections for dark matter interacting via a Yukawa potential*, JHEP 04 (2017) 077, [1611.01394].

[68] A. S. Eddington, *The distribution of stars in globular clusters*, MNRAS 76 (May, 1916) 572-585.

[69] J. Binney and S. Tremaine, *Galactic Dynamics: Second Edition*. Princeton University Press, 2008.

[70] T. Lacroix, M. Stref and J. Lavalle, *Anatomy of Eddington-like inversion methods in the context of dark matter searches*, JCAP 1809 (2018) 040, [1805.02403].

[71] M. Hütten, M. Stref, C. Combet, J. Lavalle and D. Maurin, *γ-ray and ν searches for dark-matter subhalos in the milky way with a baryonic potential*, Galaxies 7 (May, 2019) 60, [1904.10935].

[72] M. Stref, T. Lacroix and J. Lavalle, *Remnants of galactic subhalos and their impact on indirect dark-matter searches*, Galaxies 7 (June, 2019) 65, [1905.02008].

[73] G. Facchinetti, J. Lavalle and M. Stref, *Statistics for dark matter subhalo searches in gamma rays from a kinematically constrained population model. I: Fermi-LAT-like telescopes*, arXiv e-prints (July, 2020) arXiv:2007.10392, [2007.10392].

[74] G. Facchinetti, M. Stref and J. Lavalle, *Tidal stripping of dark matter subhalos by baryons from analytical perspectives: disk shocking and encounters with stars*, arXiv e-prints (Jan., 2022) arXiv:2201.09788, [2201.09788].

[75] T. Bringmann, *Particle Models and the Small-Scale Structure of Dark Matter*, New J. Phys. 11 (2009) 105027, [0903.0189].

[76] G. Lake, *Detectability of gamma-rays from clumps of dark matter*, Nature 346 (July, 1990) 39-40.

[77] N. W. Evans, F. Ferrer and S. Sarkar, *A travel guide to the dark matter annihilation signal*, Phys. Rev. D 69 (June, 2004) 123501, [astro-ph/0311145].

[78] MAGIC, FERMI-LAT collaboration, M. L. Ahnen et al., *Limits to Dark Matter Annihilation Cross-Section from a Combined Analysis of MAGIC and Fermi-LAT Observations of Dwarf Satellite Galaxies*, JCAP 02 (2016) 039, [1601.06590].

[79] A. Albert, B. Anderson, K. Bechtol, A. Drlica-Wagner, M. Meyer, M. Sánchez-Conde et al., *Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-Lat*, Astrophys. J. 834 (Jan., 2017) 110, [1611.03184].

[80] S. Hoof, A. Geringer-Sameth and R. Trotta, *A global analysis of dark matter signals from 27 dwarf spheroidal galaxies using 11 years of Fermi-LAT observations*, JCAP 2020 (Feb., 2020) 012, [1812.06986].

[81] L. E. Strigari, *Dark matter in dwarf spheroidal galaxies and indirect detection: a review*, Reports on Progress in Physics 81 (May, 2018) 056901, [1805.05883].

[82] K. Bechtol, A. Drlica-Wagner, E. Balbinot, A. Pieres, J. D. Simon, B. Yanny et al., *Eight New Milky Way Companions Discovered in First-year Dark Energy Survey Data*, Astrophys. J. 807 (July, 2015) 50, [1503.02584].

[83] S. E. Koposov, V. Belokurov, G. Torrealba and N. W. Evans, *Beasts of the Southern Wild: Discovery of Nine Ultra Faint Satellites in the Vicinity of the Magellanic Clouds.*, Astrophys. J. 805 (June, 2015) 130, [1503.02079].

[84] A. Drlica-Wagner, K. Bechtol, E. S. Rykoff, E. Luque, A. Queiroz, Y. Y. Mao et al., *Eight Ultra-faint Galaxy Candidates Discovered in Year Two of the Dark Energy Survey*, Astrophys. J. 813 (Nov., 2015) 109, [1508.03622].
[85] B. P. M. Laevens, N. F. Martin, R. A. Ibata, H.-W. Rix, E. J. Bernard, E. F. Bell et al., A New Faint Milky Way Satellite Discovered in the Pan-STARRS1 $3\pi$ Survey, Astrophys. J. Lett. 802 (Apr., 2015) L18, [1503.05554].

[86] B. P. M. Laevens, N. F. Martin, E. J. Bernard, E. F. Schlafly, B. Sesar, H.-W. Rix et al., Sagittarius II, Draco II and Laevens 3: Three New Milky Way Satellites Discovered in the Pan-STARRS 1 $3\pi$ Survey, Astrophys. J. 813 (Nov., 2015) 44, [1507.07564].

[87] D. Homma, M. Chiba, S. Okamoto, Y. Komiyama, M. Tanaka, M. Tanaka et al., A New Milky Way Satellite Discovered in the Subaru/Hyper Suprime-Cam Survey, Astrophys. J. 832 (Nov., 2016) 21, [1609.04346].

[88] D. Homma, M. Chiba, S. Okamoto, Y. Komiyama, M. Tanaka, M. Tanaka et al., Searches for new Milky Way satellites from the first two years of data of the Subaru/Hyper Suprime-Cam survey: Discovery of Cetus III, Pub. Astron. Soc. Jap. 70 (Jan., 2018) S18, [1704.05977].

[89] D. Homma, M. Chiba, Y. Komiyama, M. Tanaka, S. Okamoto, M. Tanaka et al., Boötes. IV. A new Milky Way satellite discovered in the Subaru Hyper Suprime-Cam Survey and implications for the missing satellite problem, Pub. Astron. Soc. Jap. 71 (Oct., 2019) 94, [1906.07332].

[90] G. Torrealba, V. Belokurov, S. E. Koposov, T. S. Li, M. G. Walker, J. L. Sanders et al., The hidden giant: discovery of an enormous Galactic dwarf satellite in Gaia DR2, MNRAS 488 (Sept., 2019) 2743-2766, [1811.04082].

[91] A. Drlica-Wagner, K. Bechtol, S. Mau, M. McNanna, E. O. Nadler, A. B. Pace et al., Milky Way Satellite Census. I. The Observational Selection Function for Milky Way Satellites in DES Y3 and Pan-STARRS DR1, Astrophys. J. 893 (Apr., 2020) 47, [1912.03302].

[92] B. Mutlu-Pakdil, D. J. Sand, D. Crnojević, A. Drlica-Wagner, N. Caldwell, P. Guhathakurta et al., Resolved Dwarf Galaxy Searches within 5 Mpc with the Vera Rubin Observatory and Subaru Hyper Suprime-Cam, Astrophys. J. 918 (Sept., 2021) 88, [2105.01658].

[93] L. E. Strigari, S. M. Koushiappas, J. S. Bullock and M. Kaplinghat, Precise constraints on the dark matter content of Milky Way dwarf galaxies for gamma-ray experiments, Phys. Rev. D 75 (Apr., 2007) 083526, [astro-ph/0611925].

[94] G. D. Martinez, J. S. Bullock, M. Kaplinghat, L. E. Strigari and R. Trotta, Indirect Dark Matter detection from Dwarf satellites: joint expectations from astrophysics and supersymmetry, JCAP 2009 (June, 2009) 014, [0902.4715].

[95] A. Charbonnier, C. Combet, M. Daniel, S. Funk, J. A. Hinton, D. Maurin et al., Dark matter profiles and annihilation in dwarf spheroidal galaxies: prospects for present and future $\gamma$-ray observatories - I. The classical dwarf spheroidal galaxies, MNRAS 418 (Dec., 2011) 1526-1556, [1104.0412].

[96] A. Geringer-Sameth, S. M. Koushiappas and M. Walker, Dwarf Galaxy Annihilation and Decay Emission Profiles for Dark Matter Experiments, Astrophys. J. 801 (Mar., 2015) 74, [1408.0002].

[97] V. Bonnivard, C. Combet, M. Daniel, S. Funk, A. Geringer-Sameth, J. A. Hinton et al., Dark matter annihilation and decay in dwarf spheroidal galaxies: the classical and ultrafaint dSphs, MNRAS 453 (Oct., 2015) 849-867, [1504.02048].

[98] N. W. Evans, J. L. Sanders and A. Geringer-Sameth, Simple J-factors and D-factors for indirect dark matter detection, Phys. Rev. D 93 (May, 2016) 103512, [1604.05599].

[99] J. L. Sanders and J. Binney, A review of action estimation methods for galactic dynamics, MNRAS 457 (Apr., 2016) 2107-2121, [1511.08213].

[100] A. B. Pace and L. E. Strigari, Scaling relations for dark matter annihilation and decay profiles in dwarf spheroidal galaxies, MNRAS 482 (Jan., 2019) 3480-3496, [1802.06811].
[101] A. Chiappo, J. Cohen-Tanugi, J. Conrad and L. E. Strigari, *Dwarf spheroidal J-factor likelihoods for generalized NFW profiles*, MNRAS 488 (Sept., 2019) 2616–2628, [1810.09917].

[102] A. Alvarez, F. Calore, A. Genina, J. Read, P. D. Serpico and B. Zaldivar, *Dark matter constraints from dwarf galaxies with data-driven J-factors*, JCAP 2020 (Sept., 2020) 004, [2002.01229].

[103] R. R. Muñoz, P. Côté, F. A. Santana, M. Geha, J. D. Simon, G. A. Oyarzún et al., *A MegaCam Survey of Outer Halo Satellites. III. Photometric and Structural Parameters*, Astrophys. J. 860 (June, 2018) 66, [1806.06891].

[104] J. D. Simon, *The Faintest Dwarf Galaxies*, Ann. Rev. Astron. Astrophys. 57 (Aug., 2019) 375–415, [1901.05465].

[105] J. D. Simon, T. S. Li, D. Erkal, A. B. Pace, A. Drlica-Wagner, D. J. James et al., *Birds of a Feather? Magellan/IMACS Spectroscopy of the Ultra-faint Satellites Grus II, Tucana IV, and Tucana V*, Astrophys. J. 892 (Apr., 2020) 137, [1911.08493].

[106] S. Jenkins, T. S. Li, A. B. Pace, A. P. Ji, S. E. Koposov and B. Mutlu-Pakdil, *VLT Spectroscopy of Ultra-Faint Dwarf Galaxies. I. Bo{ö}tes 1, Leo IV, Leo V*, arXiv e-prints (Dec., 2020) arXiv:2101.00013, [2101.00013].

[107] V. Bonnivard, C. Combet, D. Maurin, A. Geringer-Sameth, S. M. Koushiappas, M. G. Walker et al., *Dark Matter Annihilation and Decay Profiles for the Reticulum II Dwarf Spheroidal Galaxy*, Astrophys. J. Lett. 808 (Aug., 2015) L36, [1504.03309].

[108] M. Hütten, C. Combet and D. Maurin, *CLUMPY v3: γ-ray and ν signals from dark matter at all scales*, Computer Physics Communications 235 (Feb., 2019) 336–345, [1806.08639].

[109] V. Bonnivard, C. Combet, D. Maurin and M. G. Walker, *Spherical Jeans analysis for dark matter indirect detection in dwarf spheroidal galaxies - impact of physical parameters and triaxiality*, MNRAS 446 (Jan., 2015) 3002–3021, [1407.7822].

[110] M. L. Mateo, *Dwarf Galaxies of the Local Group*, Ann. Rev. Astron. Astrophys. 36 (Jan., 1998) 435–506, [astro-ph/9810070].

[111] A. W. McConnachie, *The Observed Properties of Dwarf Galaxies in and around the Local Group*, Astron. J. 144 (July, 2012) 4, [1204.1562].

[112] S.-H. Oh et al., *High-resolution mass models of dwarf galaxies from LITTLE THINGS*, Astron. J. 149 (2015) 180, [1502.01281].

[113] S.-H. Oh, W. J. G. de Blok, E. Brinks, F. Walter and R. C. Kennicutt, Jr, *Dark and luminous matter in THINGS dwarf galaxies*, Astron. J. 141 (2011) 193, [1011.0899].

[114] G. Gentile, P. Salucci, U. Klein and G. L. Granato, *NGC 3741: Dark halo profile from the most extended rotation curve*, Mon. Not. Roy. Astron. Soc. 375 (2007) 199–212, [astro-ph/0611355].

[115] M. Winter, G. Zaharijas, K. Bechtol and J. Vandenbroucke, *Estimating the GeV Emission of Millisecond Pulsars in Dwarf Spheroidal Galaxies*, Astrophys. J. Lett. 832 (2016) L6, [1607.06390].

[116] V. Gammaldi, E. Karukes and P. Salucci, *Theoretical predictions for dark matter detection in dwarf irregular galaxies with gamma rays*, Phys. Rev. D 98 (2018) 083008, [1706.01843].

[117] V. Gammaldi, J. Pérez-Romero, J. Coronado-Blázquez, M. Di Mauro, E. V. Karukes, M. A. Sánchez-Conde et al., *Dark matter search in dwarf irregular galaxies with the Fermi Large Area Telescope*, Phys. Rev. D 104 (Oct., 2021) 083026, [2109.11291].

[118] P. Salucci, A. Lapi, C. Tonini, G. Gentile, I. Yegorova and U. Klein, *The Universal Rotation Curve of Spiral Galaxies. 2. The Dark Matter Distribution out to the Virial Radius*, Mon. Not. Roy. Astron. Soc. 378 (2007) 41–47, [astro-ph/0703115].
[119] J. F. Navarro, C. S. Frenk and S. D. M. White, *The Structure of cold dark matter halos*, Astrophys. J. **462** (1996) 563-575, [astro-ph/9508025].

[120] J. F. Navarro, C. S. Frenk and S. D. M. White, *A Universal Density Profile from Hierarchical Clustering*, Astrophys. J. **490** (Dec., 1997) 493-508, [astro-ph/9611107].

[121] J. Einasto, *On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters*, Trudy Astrofizicheskogo Instituta Alma-Ata **5** (1965) 87-100.

[122] G. A. Gómez-Vargas, M. A. Sánchez-Conde, J.-H. Huh, M. Peiró, F. Prada, A. Morselli et al., *Constraints on WIMP annihilation for contracted dark matter in the inner Galaxy with the Fermi-LAT*, JCAP **10** (2013) 029, [1308.3515].

[123] M. Schaller, C. S. Frenk, R. G. Bower, T. Theuns, A. Jenkins, J. Schaye et al., *Baryon effects on the internal structure of ΛCDM haloes in the EAGLE simulations*, Mon. Not. Roy. Astron. Soc. **451** (2015) 1247-1267, [1409.8617].

[124] S. Peirani, Y. Dubois, M. Volonteri, J. Devriendt, K. Bundy, J. Silk et al., *Density profile of dark matter haloes and galaxies in the HORIZON-AGN simulation: the impact of AGN feedback*, MNRAS **472** (Dec., 2017) 2153-2169, [1611.09922].

[125] S. Bose et al., *No cores in dark matter-dominated dwarf galaxies with bursty star formation histories*, Mon. Not. Roy. Astron. Soc. **486** (2019) 4790-4804, [1810.03635].

[126] A. Benítez-Llambay, C. S. Frenk, A. D. Ludlow and J. F. Navarro, *Baryon-induced dark matter cores in the EAGLE simulations*, MNRAS **488** (Sept., 2019) 2387-2404, [1810.04186].

[127] A. Burkert, *The Structure of dark matter halos in dwarf galaxies*, IAU Symp. **171** (1996) 175, [astro-ph/9504041].

[128] G. M. Voit, *Tracing cosmic evolution with clusters of galaxies*, Rev. Mod. Phys. **77** (Apr, 2005) 207-258.

[129] R. J. van Weeren, F. de Gasperin, H. Akamatsu, M. Brüggen, L. Feretti, H. Kang et al., *Diffuse Radio Emission from Galaxy Clusters*, Space Sci. Rev. **215** (2019) 16, [1901.04496].

[130] FERMI-LAT collaboration, M. Ackermann et al., *Search for cosmic-ray induced gamma-ray emission in Galaxy Clusters*, Astrophys. J. **787** (2014) 18, [1308.5654].

[131] M. Colavincenzo, X. Tan, S. Ammazzalorso, S. Camera, M. Regis, J.-Q. Xia et al., *Searching for gamma-ray emission from galaxy clusters at low redshift*, Mon. Not. Roy. Astron. Soc. **491** (2020) 3225-3244, [1907.05264].

[132] M. Ackermann, M. Ajello, A. Albert, W. B. Atwood, L. Baldini, J. Ballet et al., *Search for Gamma-Ray Emission from the Coma Cluster with Six Years of Fermi-LAT Data*, Astrophys. J. **819** (Mar., 2016) 149, [1507.08995].

[133] S.-Q. Xi, X.-Y. Wang, Y.-F. Liang, F.-K. Peng, R.-Z. Yang and R.-Y. Liu, *Detection of gamma-ray emission from the Coma cluster with Fermi Large Area Telescope and tentative evidence for an extended spatial structure*, Phys. Rev. D **98** (Sept., 2018) 063006, [1709.08319].

[134] R. Adam, H. Goksu, S. Brown, L. Rudnick and C. Ferrari, *γ-ray detection toward the Coma cluster with Fermi-LAT: Implications for the cosmic ray content in the hadronic scenario*, Astron. Astroph. **648** (Apr., 2021) A60, [2102.02251].

[135] V. Baghmanyan, D. Zargaryan, F. Aharonian, R. Yang, S. Casanova and J. Mackey, *Detailed study of extended gamma-ray morphology in the vicinity of the Coma cluster with Fermi-LAT*, arXiv e-prints (Oct., 2021) arXiv:2110.00309, [2110.00309].

[136] P. Blasi, S. Gabici and G. Brunetti, *Gamma rays from clusters of galaxies*, Int. J. Mod. Phys. A **22** (2007) 681-706, [astro-ph/0701545].
[137] A. Pinzke and C. Pfrommer, *Simulating the γ-ray emission from galaxy clusters: a universal cosmic ray spectrum and spatial distribution*, MNRAS 409 (Dec., 2010) 449–480, [1001.5023].

[138] T. E. Jeltema, J. Kehayias and S. Profumo, *Gamma rays from clusters and groups of galaxies: Cosmic rays versus dark matter*, Phys. Rev. D 80 (July, 2009) 023005, [0812.0597].

[139] M. Ackermann, M. Ajello, A. Allafort, L. Baldini, J. Ballet, G. Barbiellini et al., *Constraints on dark matter annihilation in clusters of galaxies with the Fermi large area telescope*, JCAP 2010 (May, 2010) 025, [1002.2239].

[140] A. Pinzke, C. Pfrommer and L. Bergström, *Prospects of detecting gamma-ray emission from galaxy clusters: Cosmic rays and dark matter annihilations*, Phys. Rev. D 84 (Dec., 2011) 123509, [1105.3240].

[141] C. Combet, D. Maurin, E. Nezri, E. Pointecouteau, J. A. Hinton and R. White, *Decaying dark matter: Stacking analysis of galaxy clusters to improve on current limits*, Phys. Rev. D 85 (Mar., 2012) 063517, [1203.1164].

[142] S. Ando and D. Nagai, *Fermi-LAT constraints on dark matter annihilation cross section from observations of the Fornax cluster*, JCAP 2012 (July, 2012) 017, [1201.0753].

[143] FERMI-LAT collaboration, M. Ackermann et al., *Search for extended gamma-ray emission from the Virgo galaxy cluster with Fermi-LAT*, Astrophys. J. 812 (2015) 159, [1510.00004].

[144] MAGIC collaboration, V. A. Acciari et al., *Constraining Dark Matter lifetime with a deep gamma-ray survey of the Perseus Galaxy Cluster with MAGIC*, Phys. Dark Univ. 22 (2018) 38–47, [1806.11063].

[145] D. Maurin, C. Combet, E. Nezri and E. Pointecouteau, *Disentangling cosmic-ray and dark-matter induced γ-rays in galaxy clusters*, Astron. Astroph. 547 (Nov., 2012) A16, [1203.1166].

[146] M. A. Sánchez-Conde, M. Cannoni, F. Zandanel, M. E. Gómez and F. Prada, *Dark matter searches with Cherenkov telescopes: nearby dwarf galaxies or local galaxy clusters?*, JCAP 2011 (Dec., 2011) 011, [1104.3530].

[147] T. H. Reiprich and H. Böhringer, *The Mass Function of an X-Ray Flux-limited Sample of Galaxy Clusters*, Astrophys. J. 567 (Mar., 2002) 716–740, [astro-ph/0111285].

[148] R. Piffaretti, M. Arnaud, G. W. Pratt, E. Pointecouteau and J. B. Melin, *The MCXC: a meta-catalogue of x-ray detected clusters of galaxies*, Astron. Astroph. 534 (Oct., 2011) A109, [1007.1916].

[149] S. L. Snowden, R. M. Mushotzky, K. D. Kuntz and D. S. Davis, *A Catalog of Galaxy Clusters Observed by XMM-Newton*, Astron. Astrophys. 478 (2008) 615–658, [0710.2241].

[150] A. Vikhlinin et al., *Chandra Cluster Cosmology Project II: Samples and X-ray Data Reduction*, Astrophys. J. 692 (2009) 1033–1059, [0805.2207].

[151] G. Schellenberger and T. H. Reiprich, *HICOSMO – cosmology with a complete sample of galaxy clusters – I. Data analysis, sample selection and luminosity-mass scaling relation*, Mon. Not. Roy. Astron. Soc. 469 (2017) 3738–3761, [1705.05842].

[152] G. W. Pratt, M. Arnaud, A. Biviano, D. Eckert, S. Ettori, D. Nagai et al., *The galaxy cluster mass scale and its impact on cosmological constraints from the cluster population*, Space Sci. Rev. 215 (2019) 25, [1902.10837].

[153] Y.-Y. Zhang, T. H. Reiprich, P. Schneider, N. Clerc, A. Merloni, A. Schwope et al., *HIFLUGCS: X-ray luminosity-dynamical mass relation and its implications for mass calibrations with the SPIDERS and 4MOST surveys*, Astron. Astroph. 599 (Mar., 2017) A138, [1608.06585].
[154] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont et al., Planck 2015 results. XXVII. The second Planck catalogue of Sunyaev-Zeldovich sources, Astron. Astroph. 594 (Sept., 2016) A27, [1502.01598].

[155] Y. Chen, T. H. Reiprich, H. Bohringer, Y. Ikebe and Y. Y. Zhang, Statistics of X-ray observables for the cooling-core and non-cooling core galaxy clusters, Astron. Astrophys. 466 (2007) 805, [astro-ph/0702482].

[156] R. Adam, H. Goksu, A. Leingartner-Goth, S. Ettori, R. Gnatyk, B. Hnatyk et al., MINOT: Modeling the intracluster medium (non-)thermal content and observable prediction tools, Astron. Astrophys. 644 (2020) A70, [2009.05373].

[157] M. Kuhlen and D. Malyshev, Atic, pamela, hess, and fermi data and nearby dark matter subhalos, Phys. Rev. D 79 (June, 2009) 123517, [0904.3378].

[158] L. Pieri, M. Lattanzi and J. Silk, Constraining the dark matter annihilation cross-section with Cherenkov telescope observations of dwarf galaxies, MNRAS 399 (Nov., 2009) 2033–2040, [0902.4330].

[159] T. R. Slatyer, N. Toro and N. Weiner, Sommerfeld-enhanced annihilation in dark matter substructure: Consequences for constraints on cosmic-ray excesses, Phys. Rev. D 86 (Oct., 2012) 083534, [1107.3546].

[160] J. Zavala and N. Afshordi, Clustering in the phase space of dark matter haloes - II. Stable clustering and dark matter annihilation, MNRAS 441 (June, 2014) 1329–1339, [1311.3296].

[161] K. K. Boddy, J. Kumar, J. Runburg and L. E. Strigari, Angular distribution of gamma-ray emission from velocity-dependent dark matter annihilation in subhalos, Phys. Rev. D 100 (Sept., 2019) 063019, [1905.03431].

[162] E. Piccirillo, K. Blanchette, N. Bozorgnia, L. E. Strigari, C. S. Frenk, R. J. J. Grand et al., Velocity-dependent annihilation radiation from dark matter subhalos in cosmological simulations, 2203.08853.

[163] S. Palomares-Ruiz and J. M. Siegal-Gaskins, Annihilation vs. decay: constraining dark matter properties from a gamma-ray detection, JCAP 2010 (July, 2010) 023, [1003.1142].

[164] B. Boucher, J. Kumar, V. B. Le and J. Runburg, J-factors for Velocity-dependent Dark Matter, arXiv e-prints (Oct., 2021) arXiv:2110.09653, [2110.09653].

[165] Fermi LAT Collaboration, Limits on dark matter annihilation signals from the Fermi LAT 4-year measurement of the isotropic gamma-ray background, JCAP 2015 (Sept., 2015) 008, [1501.05464].

[166] P. J. McMillan, The mass distribution and gravitational potential of the Milky Way, MNRAS 465 (Feb., 2017) 76–94, [1608.00971].

[167] M. Cirelli, N. Fornengo and A. Strumia, Minimal dark matter, Nucl. Phys. B 753 (2006) 178–194, [hep-ph/0512090].

[168] C. Garcia-Cely, A. Ibarra, A. S. Lamperstorfer and M. H. G. Tytgat, Gamma-rays from Heavy Minimal Dark Matter, JCAP 2015 (Oct., 2015) 058, [1507.05536].

[169] J. Choquette, J. M. Cline and J. M. Cornell, p-wave annihilating dark matter from a decaying predecessor and the Galactic Center excess, Phys. Rev. D 94 (July, 2016) 015018, [1604.01039].

[170] D. P. Finkbeiner, L. Goodenough, T. R. Slatyer, M. Vogelsberger and N. Weiner, Consistent Scenarios for Cosmic-Ray Excesses from Sommerfeld-Enhanced Dark Matter Annihilation, JCAP 05 (2011) 002, [1011.3082].

[171] Z.-P. Liu, Y.-L. Wu and Y.-F. Zhou, Sommerfeld enhancements with vector, scalar and pseudoscalar force-carriers, Phys. Rev. D 88 (2013) 096008, [1305.5438].
[172] Y.-C. Ding, Y.-L. Ku, C.-C. Wei and Y.-F. Zhou, *Consistent explanation for the cosmic-ray positron excess in p-wave Sommerfeld-enhanced dark matter annihilation*, JCAP 09 (2021) 005, [2104.14881].

[173] S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo et al., *Closing the window on wimp dark matter*, Eur. Phys. J. C 82 (2022) 31, [2107.09688].

[174] T. Hambye and L. Vanderheyden, *Minimal self-interacting dark matter models with light mediator*, JCAP 2020 (May, 2020) 001, [1912.11708].

[175] A. Das and B. Dasgupta, *Selection rule for enhanced dark matter annihilation*, Phys. Rev. Lett. 118 (June, 2017) 251101, [1611.04606].

[176] J. Binney, *The radius-dependence of velocity dispersion in elliptical galaxies*, MNRAS 190 (Mar., 1980) 873-880.

[177] R. Wojtak, E. L. Łokas, G. A. Mamon, S. Gottlöber, A. Klypin and Y. Hoffman, *The distribution function of dark matter in massive haloes*, MNRAS 388 (Aug., 2008) 815-828, [0802.0429].

[178] M. Henon, *Numerical Experiments on the Stability of Spherical Stellar Systems*, Astron. Astroph. 24 (Apr., 1973) 229.

[179] S. M. Kent and J. E. Gunn, *The dynamics of rich clusters of galaxies. i - the coma cluster*, Astron. J. 87 (July, 1982) 945-971.

[180] P. Cuddeford, *An analytic inversion for anisotropic spherical galaxies*, MNRAS 253 (Dec., 1991) 414-426.

[181] N. W. Evans and J. H. An, *Distribution function of the dark matter*, Phys. Rev. D73 (2006) 023524, [astro-ph/0511687].

[182] J. Diemand, B. Moore and J. Stadel, *Velocity and spatial biases in cold dark matter subhalo distributions*, MNRAS 352 (Aug., 2004) 535-546, [astro-ph/0402160].

[183] J. Diemand, M. Kuhlen and P. Madau, *Formation and Evolution of Galaxy Dark Matter Halos and Their Substructure*, Astrophys. J. 667 (Oct., 2007) 859-877, [astro-ph/0703337].

[184] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, A. Jenkins, A. Helmi et al., *The Aquarius Project: the subhaloes of galactic haloes*, MNRAS 391 (Dec., 2008) 1685-1711, [0809.0898].

[185] A. J. Benson, *The mass function of unprocessed dark matter haloes and merger tree branching rates*, MNRAS 467 (May, 2017) 3454-3466, [1610.01057].

[186] C. Giocoli, L. Pieri and G. Tormen, *Analytical approach to subhalo population in dark matter haloes*, MNRAS 387 (June, 2008) 689-697, [0712.1476].

[187] F. Jiang and F. C. van den Bosch, *Statistics of Dark Matter Substructure: I. Model and Universal Fitting Functions*, MNRAS 458 (May, 2016) 2848-2869, [1403.6827].

[188] S. Ando, T. Ishiyama and N. Hiroshima, *Halo substructure boosts to the signatures of dark matter annihilation*, Galaxies 7 (July, 2019) 68, [1903.11427].

[189] T. Ishiyama and S. Ando, *The abundance and structure of subhaloes near the free streaming scale and their impact on indirect dark matter searches*, MNRAS 492 (Mar., 2020) 3662-3671, [1907.03642].

[190] A. Charbonnier, C. Combet and D. Maurin, *CLUMPY: A code for γ-ray signals from dark matter structures*, Computer Physics Communications 183 (Mar., 2012) 656-668, [1201.4728].

[191] A. J. Benson, *Galacticus: A semi-analytic model of galaxy formation*, New Astronomy 17 (Feb., 2012) 175-197, [1008.1786].
[192] M. D. Weinberg, *Adiabatic invariants in stellar dynamics. 1: Basic concepts*, Astron. J. **108** (Oct., 1994) 1398-1402, [astro-ph/9404015].

[193] O. Y. Gnedin, L. Hernquist and J. P. Ostriker, *Tidal Shocking by Extended Mass Distributions*, Astrophys. J. **514** (Mar., 1999) 109-118.

[194] E. Hayashi, J. F. Navarro, J. E. Taylor, J. Stadel and T. Quinn, *The Structural Evolution of Substructure*, Astrophys. J. **584** (Feb., 2003) 541-558, [astro-ph/0203004].

[195] F. C. van den Bosch, G. Ogiya, O. Hahn and A. Burkert, *Disruption of dark matter substructure: fact or fiction?*, MNRAS **474** (Mar., 2018) 3043-3066, [1711.05276].

[196] S. B. Green, F. C. van den Bosch and F. Jiang, *The tidal evolution of dark matter substructure - ii. the impact of artificial disruption on subhalo mass functions and radial profiles*, MNRAS (Mar., 2021), [2103.01227].

[197] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, A. A. Klypin et al., *Profiles of dark haloes: evolution, scatter and environment*, MNRAS **321** (Mar., 2001) 559-575, [astro-ph/9908159].

[198] A. V. Macciò, A. A. Dutton and F. C. van den Bosch, *Concentration, spin and shape of dark matter haloes as a function of the cosmological model: WMAP1, WMAP3 and WMAP5 results*, MNRAS **391** (Dec., 2008) 1940-1954, [0805.1926].

[199] F. Prada, A. A. Klypin, A. J. Cuesta, J. E. Betancort-Rijo and J. Primack, *Halo concentrations in the standard Λ cold dark matter cosmology*, MNRAS **423** (July, 2012) 3018-3030, [1104.5130].

[200] A. A. Dutton and A. V. Macciò, *Cold dark matter haloes in the Planck era: evolution of structural parameters for Einasto and NFW profiles*, MNRAS **441** (July, 2014) 3359-3374, [1402.7073].

[201] F. C. van den Bosch and F. Jiang, *Statistics of Dark Matter Substructure: II. Comparison of Model with Simulation Results*, MNRAS **458** (May, 2016) 2870-2884, [1403.6835].

[202] J. Han, S. Cole, C. S. Frenk and Y. Jing, *A unified model for the spatial and mass distribution of subhaloes*, MNRAS **457** (Apr., 2016) 1208-1223, [1509.02175].

[203] R. S. Somerville and T. S. Kolatt, *How to plant a merger tree*, MNRAS **305** (May, 1999) 1-14, [astro-ph/9711080].

[204] H. Parkinson, S. Cole and J. Helly, *Generating dark matter halo merger trees*, MNRAS **383** (Jan., 2008) 557-564, [0708.1382].

[205] J. R. Bond, S. Cole, G. Efstathiou and N. Kaiser, *Excursion set mass functions for hierarchical Gaussian fluctuations*, Astrophys. J. **379** (Oct., 1991) 440-460.

[206] C. Lacey and S. Cole, *Merger rates in hierarchical models of galaxy formation*, MNRAS **262** (June, 1993) 627-649.

[207] S. Cole, C. G. Lacey, C. M. Baugh and C. S. Frenk, *Hierarchical galaxy formation*, MNRAS **319** (Nov., 2000) 168-204, [astro-ph/0007281].

[208] Planck Collaboration, *Planck 2018 results. vi. cosmological parameters*, Astron. Astroph. **641** (Sept., 2020) A6, [1807.06209].

[209] G. Facchinetti, *Analytical study of particle dark matter structuring on small scales and implications for dark matter searches*, Ph.D. thesis, University of Montpellier, 2021.