Halo Substructure Boosts to the Signatures of Dark Matter Annihilation

TL;DR

Dark matter substructure significantly enhances annihilation signals, but accurately estimating this boost factor across all scales remains challenging due to uncertainties in extrapolations from simulations and the need for well-calibrated models.

Contribution

This paper reviews recent advances in estimating the dark matter subhalo boost factor using simulations and analytic models, highlighting the importance of accurate assessments for dark matter detection.

Findings

Numerical simulations suggest a large boost factor (>100) for galaxy-sized halos.

Analytic models predict a more modest boost of order unity for galaxy halos.

Extrapolations to small-scale subhalos involve significant uncertainties.

Abstract

The presence of dark matter substructure will boost the signatures of dark matter annihilation. We review recent progress on estimates of this subhalo boost factor---a ratio of the luminosity from annihilation in the subhalos to that originating the smooth component---based on both numerical $N$-body simulations and semi-analytic modelings. Since subhalos of all the scales, ranging from the Earth mass (as expected, e.g., the supersymmetric neutralino, a prime candidate for cold dark matter) to galaxies or larger, give substantial contribution to the annihilation rate, it is essential to understand subhalo properties over a large dynamic range of more than twenty orders of magnitude in masses. Even though numerical simulations give the most accurate assessment in resolved regimes, extrapolating the subhalo properties down in sub-grid scales comes with great uncertainties---a straightforward extrapolation yields a very large amount of the subhalo boost factor of $\gtrsim$100 for galaxy-size halos. Physically motivated theoretical models based on analytic prescriptions such as the extended Press-Schechter formalism and tidal stripping modeling, which are well tested against the simulation results, predict a more modest boost of order unity for the galaxy-size halos. Giving an accurate assessment of the boost factor is essential for indirect dark matter searches and thus, having models calibrated at large ranges of host masses and redshifts, is strongly urged upon.

Equations

$${\rho }_{\chi }(\mathbf{x})={\rho }_{\chi ,\mathrm{sm}}(\mathbf{x})+{\rho }_{\chi ,\mathrm{sh}}(\mathbf{x}),$$

$$⟨{\rho }^2(\mathbf{x}){⟩}_{M,z}=⟨{\rho }_{\mathrm{sm}}^2(\mathbf{x}){⟩}_{M,z}+⟨{\rho }_{\mathrm{sh}}^2(\mathbf{x}){⟩}_{M,z}+2⟨{\rho }_{\mathrm{sm}}(\mathbf{x}){\rho }_{\mathrm{sh}}(\mathbf{x}){⟩}_{M,z}.$$

$${\rho }_{\mathrm{NFW}}(r)=\frac{{\rho }_s}{(r/r_s)(1+r/r_s{)}^2},$$

$$⟨{\rho }_{\mathrm{sm}}^2(\mathbf{x}){⟩}_{M,z}=\frac{1}{V}\int d^{3}x{\rho }_{\mathrm{sm}}^2(r)=\frac{4\pi (1-f_{\mathrm{sh}}{)}^2}{3V}{\rho }_s^2{r}_s^3\left[1-\frac{1}{(1+c_{\mathrm{vir}}{)}^3}\right],$$

$${\rho }_{\mathrm{sh}}(\mathbf{x})=\int d{m}^{\prime }\int d^3{x}^{\prime }\sum _i{\delta }_D(m^{\prime }-m_i){\delta }_D^3({\mathbf{x}}^{\prime }-{\mathbf{x}}_i)m^{\prime }{u}_{\mathrm{sh}}(\mathbf{x}-{\mathbf{x}}^{\prime }|m^{\prime }),$$

$$\frac{d{n}_{\mathrm{sh}}(\mathbf{x},m)}{dm}=⟨\sum _i{\delta }_D(m-m_i){\delta }_D^3(\mathbf{x}-{\mathbf{x}}_i)⟩,$$

$$\frac{d{N}_{\mathrm{sh}}}{dm}=\int d^{3}x\frac{d{n}_{\mathrm{sh}}(\mathbf{x},m)}{dm},$$

$$f_{\mathrm{sh}}(M,z)=\frac{1}{M}\int dmm\frac{d{N}_{\mathrm{sh}}}{dm}.$$

$$⟨{\rho }_{\mathrm{sh}}^2(\mathbf{x}){⟩}_{M,z}$$

$$B_{\mathrm{sh}}(M,z)=\frac{1}{L_{\mathrm{host},0}(M,z)}\int dm\frac{d{N}_{\mathrm{sh}}(m|M,z)}{dm}{L}_{\mathrm{sh}}(m),$$

$$2⟨{\rho }_{\mathrm{sm}}(\mathbf{x}){\rho }_{\mathrm{sh}}(\mathbf{x}){⟩}_{M,z}$$

$$\frac{d{n}_{\mathrm{sh}}(\mathbf{x},m)}{dm}=P_{\mathrm{sh}}(\mathbf{x})\frac{d{N}_{\mathrm{sh}}(m)}{dm},$$

$$2⟨{\rho }_{\mathrm{sm}}(\mathbf{x}){\rho }_{\mathrm{sh}}(\mathbf{x}){⟩}_{M,z}=\frac{2f_{\mathrm{sh}}M}{V}\int d^{3}x{\rho }_{\mathrm{sm}}(\mathbf{x})P_{\mathrm{sh}}(\mathbf{x}).$$

$$2⟨{\rho }_{\mathrm{sm}}(\mathbf{x}){\rho }_{\mathrm{sh}}(\mathbf{x}){⟩}_{M,z}=\frac{2f_{\mathrm{sh}}}{1-f_{\mathrm{sh}}}⟨{\rho }_{\mathrm{sm}}^2(\mathbf{x})⟩=\frac{8\pi f_{\mathrm{sh}}(1-f_{\mathrm{sh}})}{3V}{\rho }_s^2{r}_s^3\left[1-\frac{1}{(1+c_{\mathrm{vir}}{)}^3}\right].$$

$$L_{\mathrm{total}}$$

$${\Delta }_{\mathrm{vir}}(z)=18{\pi }^2+82d(z)-39d^2(z),$$

$$M=\frac{4\pi }{3}{\Delta }_{\mathrm{vir}}(z){\rho }_c(z)r_{\mathrm{vir}}^3.$$

$$M_{200}=\frac{4\pi }{3}200{\rho }_c(z)r_{200}^3.$$

$$M$$

$$f(x)$$

$$\frac{d{N}_{\mathrm{sh}}}{dm}\propto m^{-\alpha },$$

$$\rho (r)=\frac{{\rho }_s}{(r/r_s{)}^{\gamma }{\left[1+{\left(r/r_s\right)}^{\eta }\right]}^{(\beta -\gamma )/\eta }}.$$

$$\rho (r)={\rho }_s\exp \left\{-\frac{2}{{\alpha }_E}\left[{\left(\frac{r}{r_s}\right)}^{{\alpha }_E}-1\right]\right\},$$

$$c_V=\frac{\bar{\rho }(<r_{\max })}{{\rho }_{\mathrm{c}}}=2{\left(\frac{V_{\max }}{H_0{r}_{\max }}\right)}^2,$$

$${\tau }_{\mathrm{dyn}}(z)=1.628h^{-1}\mathrm{Gyr}{\left[\frac{{\Delta }_{\mathrm{vir}}(z)}{178}\right]}^{1/2}{\left[\frac{H(z)}{H_0}\right]}^{-1},$$

$$\dot{m}(z)=-A\frac{m(z)}{{\tau }_{\mathrm{dyn}}(z)}{\left[\frac{m(z)}{M(z)}\right]}^{\zeta },$$

$$\log A$$

$$\zeta$$

$$\rho(r)=\left\{\begin{array}[]{lcc}\rho_{s}r_{s}^{3}/[r(r+r_{s})]^{2},&{\rm for}&r<r_{t},\\ 0,&{\rm for}&r\geq r_{t}.\end{array}\right.$$

$$\frac{V_{\max }(z)}{V_{\max ,a}}$$

$$\frac{r_{\max }(z)}{r_{\max ,0}}$$

$$r_s$$

Full text

\articlenumber

x \doinum10.3390/—— \pubvolumexx

\historyReceived: 27 March 2019; Accepted: 26 June 2019; Published: 1 July 2019 \TitleHalo Substructure Boosts to the Signatures of Dark Matter Annihilation \AuthorShin’ichiro Ando *1,2,**\orcidA, Tomoaki Ishiyama 3 and Nagisa Hiroshima 4,5,6 \AuthorNamesShin’ichiro Ando, Tomoaki Ishiyama and Nagisa Hiroshima

\corresCorrespondence: s.ando@uva.nl

1 Introduction

One of the most popular candidates for dark matter is weakly interacting massive particles (WIMPs) Bertone et al. (2005); Bringmann and Weniger (2012). They are motivated by beyond-the-standard-model physics such as supersymmetry Jungman et al. (1996) or universal extra-dimensions Hooper and Profumo (2007), although the non-discovery of new physics at the TeV scale with the Large Hadron Collider puts these models to serious test Boveia and Doglioni (2018). In addition, WIMPs can naturally explain the relic dark matter density with thermal freezeout mechanisms, where the WIMPs following the weak-scale physics were in chemical equilibrium until freezeout—when the expansion of the Universe became faster than the annihilation rate Steigman et al. (2012). Since dark matter is often the lightest particle in an extended sector, it can self-annihilate only into the standard-model particles, which end up producing gamma rays, charged cosmic rays, and neutrinos. Indirect detection of dark matter annihilation is therefore a direct test of the thermal freezeout of WIMPs.

WIMPs are also a subcategory of cold dark matter (CDM), where they were nonrelativistic when structure formation started. In the CDM framework, it is known that the structures form hierarchically, from smaller to larger ones. These virialized structures are referred to as halos and they are nearly spherically symmetric. Typical size of the smallest structure is highly model dependent. In the case of the supersymmetric neutralino that is one of the most popular WIMP candidates, the smallest halos tend to be of the Earth mass, 10${}^{-6}M_{\odot}$ but with very large range of possible values of $\sim$10*-12*–10${}^{-3}M_{\odot}$ Hofmann et al. (2001); Green et al. (2004); Loeb and Zaldarriaga (2005); Bertschinger (2006); Profumo et al. (2006); Bringmann (2009); Cornell and Profumo (2012); Diamanti et al. (2015). Smaller halos collapse at higher redshifts when the Universe was denser, and hence they are of higher density. A larger dark matter halo today contains lots of substructures (or subhalos) of all mass scales, which can go down to the Earth masses or even smaller and hence denser.

Since the annihilation rate depends on the dark matter density squared (and $\langle\rho^{2}\rangle\geq\langle\rho\rangle^{2}$), the presence of the subhalos will boost the gamma-ray signatures from dark matter annihilation. This subhalo boost of dark matter annihilation, in relation with the smallest-scale subhalos, has been a topic of interest for very many years Silk and Stebbins (1993); Bergstrom et al. (1999, 1998); Calcaneo-Roldan and Moore (2000); Tasitsiomi and Olinto (2002); Stoehr et al. (2003); Koushiappas et al. (2004); Baltz et al. (2007); Ando (2005); Oda et al. (2005); Pieri et al. (2005); Koushiappas (2006); Ando et al. (2007); Pieri et al. (2008); Diemand et al. (2007); Berezinsky et al. (2007); Lavalle et al. (2008); Siegal-Gaskins (2008); Ando et al. (2008); Lee et al. (2009); Kamionkowski and Koushiappas (2008); Springel et al. (2008, 2008); Ando (2009); Kamionkowski et al. (2010); Pinzke et al. (2011); Gao et al. (2012); Fornasa et al. (2013); Ng et al. (2014); Sánchez-Conde and Prada (2014); Ishiyama (2014); Bartels and Ando (2015); Stref and Lavalle (2017); Moliné et al. (2017); Zavala and Afshordi (2016); Hiroshima et al. (2018). The main difficulty is the fact that subhalos of all the scales ranging from the Earth mass (or even smaller) to larger masses (a significant fraction of their host’s mass) give a substantial contribution to the annihilation rate. Covering this very large dynamic range is challenging even with the state-of-the-art numerical simulations. In simulations of Milky-Way-size halos ($10^{12}M_{\odot}$) Diemand et al. (2008); Springel et al. (2008); Stadel et al. (2009), one can resolve only down to $10^{4}$–$10^{5}M_{\odot}$, and there still remains more than ten orders of magnitude to reach.

We will review recent progress on the subhalo contribution to dark matter annihilation. (See also Reference Berezinsky et al. (2014) for a review on generic processes that subhalos undergo.) We first discuss approaches using the numerical $N$-body simulations and estimate of the annihilation boost factor by adoping the results and extrapolation down to very-small-mass ranges. To complement the approach based on simulations, we then review an analytical approach. In the CDM framework, fraction of halos that collapse is described with the Press-Schechter formalism Press and Schechter (1974) based on spherical or ellipsoidal collapse models. This has been further extended to accommodate collapsed regions within larger halos (excursion set or extended Press-Schechter formalism Bond et al. (1991)), which can be applied to address statistics of halo substructure. More recent literature suggests that the annihilation boost factor, defined as the luminosity due to subhalos divided by the host luminosity, is modest, ranging from order of unity to a few tens for galaxy-size halos Kamionkowski and Koushiappas (2008); Bartels and Ando (2015); Zavala and Afshordi (2016); Moliné et al. (2017); Stref and Lavalle (2017); Hiroshima et al. (2018). This relatively mild amount of the annihilation boost makes the prospect of indirect dark matter searches less promising compared with earlier more optimistic predictions Springel et al. (2008); Pinzke et al. (2011); Gao et al. (2012); Anderson et al. (2016). We note that our focus is mainly on subhalo boost factors in extragalactic halos. For the subhalo boosts in the Galactic halo, on the other hand, we need to assess the spatial distribution of the subhalos too. The $N$-body simulations described in Section 3 can address this issue but again are subject to resolution issues as well as the baryonic effect. See, for example, Reference Stref and Lavalle (2017) for an alternative approach adopting analytical prescription.

This review is organized as follows. In Section 2, we introduce basic concepts of density profiles, mass functions, and the annihilation boost factors of the subhalos, starting with simple formulations. Here we make some simplifying assumptions, which are to be addressed in later sections. In Section 3, we summarize the progress from the numerical simulations for the subhalos and the annihilation boost factors. Section 4 presents more recent approaches based on realistic formulation than Section 2. In Section 4.1, we first show new analytic models that predict the subhalo mass functions well in agreement with the results from the numerical simulations for various ranges of the host masses and redshifts, and that the annihilation boost factors are on the order of unity even for cluster-size halos. Then, we summarize other semi-analytic approaches for computing the annihilation boost factors, based on self-similarity (Section 4.2) and universal phase-space clustering (Section 4.3) of the subhalos. We conclude the review in Section 5. Finally, for convenience, we summarize fitting functions for the subhalo mass functions, and annihilation boost factors that can be applicable to nearly arbitrary masses and redshifts in Appendix A.

2 Formulation

In this section, we introduce several important quantities such as density profiles, subhalo mass function, and the annihilation boost factors. This section is based on a simplified analytic model, which in several aspects are unrealistic but sets the basis for the latter discussions according to numerical simulations (Section 3) and more sophisticated semi-analytical models (Section 4).

2.1 Subhalo Boost Factor

The rate of dark matter annihilation is proportional to dark matter density squared, $\rho_{\chi}^{2}$, where $\chi$ represents the dark matter particle. In the presence of substructure, $\rho_{\chi}$ is divided into two terms:

[TABLE]

where $\rho_{{\rm sm}}$ and $\rho_{{\rm sh}}$ represent smooth and subhalo components, respectively. (In the following, we omit the subscript $\chi$.) The volume average of the density squared in a host halo characterized by its virial mass $M$ and redshift $z$, which is the relevant quantity for the indirect dark matter researches, is therefore written as

[TABLE]

We assume that the smooth component $\rho_{\rm sm}$ is characterized by the following Navarro-Frenk-White (NFW) profile Navarro et al. (1996, 1997):

[TABLE]

where $\rho_{s}$ is a characteristic density and $r_{s}$ is a scale radius. These parameters, $\rho_{s}$ and $r_{s}$, are evaluated such that the volume integral of $\rho_{\rm NFW}$ yields the total halo mass $M$, and thus we have $\rho_{\rm sm}(r)=(1-f_{\rm sh})\rho_{\rm NFW}(r)$, where $f_{\rm sh}$ is defined as the mass fraction in the subhalos. The first term is then simply

[TABLE]

where $V=4\pi r_{\rm vir}^{3}/3$ is the volume of the host out to its virial radius $r_{\rm vir}$, $c_{\rm vir}\equiv r_{\rm vir}/r_{s}$ is the concentration parameter. The parameters characterizing the host profile – $\rho_{s}$, $r_{s}$, and $c_{\rm vir}$ – are all functions of $M$ and $z$.111We note, however, that the concentration $c_{\rm vir}$ has a scatter, which is often characterized by a log-normal distribution, whose mean $\bar{c}_{\rm vir}$ is the function of $M$ and $z$. We will include this in the latter sections.

Next, we evaluate the second term of Equation (2), $\langle\rho_{\rm sh}^{2}(\bm{x})\rangle$. We characterize each subhalo $i$ with the location of its center $\bm{x}_{i}$ and mass $m_{i}$. Density due to all the subhalos at a coordinate $\bm{x}$ is written as a sum of the density profile around the seed of each subhalo, that is,

[TABLE]

where $\delta_{D}^{N}$ is the $N$-dimensional Dirac delta function, and $u_{\rm sh}(r|m)$ defines the density profile of the subhalo with mass $m$ and is normalized to one after the volume integral.222For the sake of simplicity for analytic expressions, we assume that the suhbalo mass is the only parameter characterizing its density profile. One can introduce many more parameters to make the model more realistic. We define the ensemble average of the product of these delta functions as

[TABLE]

its volume integral over the host halo as

[TABLE]

and call both $dn_{\rm sh}/dm$ and $dN_{\rm sh}/dm$ the subhalo mass function. We also obtain the mass fraction in the subhalos as

[TABLE]

By multiplying Equation (5) by itself and taking both the ensemble and the volume averages, we have

[TABLE]

where at the last equality we adopted the NFW function for the subhalo density profile $mu_{\rm sh}(r|m)$ with the scale radius $r_{s,{\rm sh}}$, characteristic density $\rho_{s,{\rm sh}}$, and tidal truncation radius $r_{t,{\rm sh}}\equiv c_{t,{\rm sh}}r_{s,{\rm sh}}$ beyond which the subhalo density abruptly decreases to zero. At the third equality of Equation (9), we ignored the term arising from $j\neq i$ as we evaluate the quantity at one point $\bm{x}$ and assume that subhalos do not overlap. We note, however, that such a term becomes relevant for obtaining the two-point correlation function, or the power spectrum; see References Ando et al. (2007); Ando (2009) for more details.

We define the subhalo boost factor as the ratio of the total luminosity from dark matter annihilation in the subhalos and that from the smooth component in the case that there is no substructure. By comparing Equations (4) and (9), and remembering that the luminosity is proportional to the volume integral of the density squared, the boost factor is simply written as

[TABLE]

where the subscript [math] shows that this is a quantity in the case of no subhalo contributions. Equation (10) is also valid for any other spherically symmetric density profiles than the NFW.

Finally, we evaluate the last cross-correlation term in Equation (2). See also References Charbonnier et al. (2012); Nezri et al. (2012); Ando and Komatsu (2013); Stref and Lavalle (2017); Hütten et al. (2016). Following a similar procedure as in Equation (9), we have

[TABLE]

where in the last equality, we first used the fact that the subhalo density profile is much more sharply peaked than their spatial distribution, and take $dn_{\rm sh}/dm$ out of $\bm{x}^{\prime}$ integration adopting $\bm{x}^{\prime}\approx\bm{x}$ as its spatial variable. Second, we performed volume integral for $u(\bm{x}-\bm{x}^{\prime}|m)$ over $\bm{x}^{\prime}$ variable, which simply returns one, to reach the last expression of Equation (11). Then we assume that the spatial distribution of the subhalos is independent of their masses:

[TABLE]

where $P_{\rm msh}(\bm{x})d^{3}x$ represents the probability of finding a subhalo in a volume element $d^{3}x$ around $\bm{x}$. With this and Equation (8), we have

[TABLE]

For simplicity, we assume that the subhalos are distributed following the smooth NFW component. In this case, we have $\rho_{\rm sm}(\bm{x})=(1-f_{\rm sh})MP_{\rm sh}(\bm{x})$, and

[TABLE]

The luminosity from the smooth component in the presence of the subhalos ($L_{\rm sm}$) is related to the host luminosity in the subhalos’ absence via $L_{\rm sm}=(1-f_{\rm sh})^{2}L_{\rm host,0}$, because the density in the smooth component gets depleted by a factor of $1-f_{\rm sh}$, if there are subhalos. Thus, the total luminosity from both the smooth component and the subhalos are given by

[TABLE]

Often, $L_{\rm total}/L_{\rm host,0}$ is also referred to as the subhalo boost factor in the literature. Note, however, that we have not included the effect of sub-subhalos (and beyond) yet in this formalism. In order to accommodate it, in the right-hand side of Equation (10), we need to include the sub-subhalo boost to the subhalo luminosity $L_{\rm sh}$. Thus, we replace $L_{\rm sh}$ with $(1-f_{\rm ssh}^{2}+B_{\rm ssh})L_{\rm sh}$, where the subscript “ssh” represents the contribution from the sub-subhalos. If the subhalo mass fraction $f_{\rm sh}$ and the boost factor $B_{\rm sh}$ depend only on the host mass, then one can assume $f_{\rm ssh}(m)=f_{\rm sh}(m)$ and $B_{\rm ssh}(m)=B_{\rm sh}(m)$, and repeat the calculations in an iterative manner. See, however, Section 4.1 for a more realistic treatment.

2.2 Characterization of Dark Matter Halos

We shall discuss the density profile of dark matter halos that are characterized by the virial radius $r_{\rm vir}$, the scale radius $r_{s}$, and the characteristic density $\rho_{s}$. The halo is virialized when a mean density within a region reaches some critical value times the critical density of the Universe at that time: $\Delta_{\rm vir}(z)\rho_{c}(z)$, where $\rho_{c}(z)=\rho_{c,0}[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}]$, $\rho_{c,0}=3H_{0}^{2}/(8\pi G)$ is the present critical density, $H_{0}=100h$ km s*-1* Mpc*-1* is the Hubble constant, $\Omega_{m}$ and $\Omega_{\Lambda}$ are the density parameters for matter and the cosmological constant, respectively. In CDM cosmology with the cosmological constant, this critical value is given as Bryan and Norman (1998)

[TABLE]

where $d(z)=\Omega_{m}(1+z)^{3}/[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}]-1$. Given the virial mass $M$ and the redshift $z$ of the halo of interest, $r_{\rm vir}$ is therefore obtained by solving

[TABLE]

Alternatively, one can define $M_{200}$ and $r_{200}$ via

[TABLE]

$M_{200}$ is often adopted to define halo masses in $N$-body simulations.

The concentration parameter $c_{\rm vir}\equiv r_{\rm vir}/r_{s}$ (or $c_{200}\equiv r_{200}/r_{s}$) has been studied with numerical simulations and found to be a function of $M$ and $z$. It follows a log-normal distribution with the mean of $\bar{c}_{\rm vir}(M,z)$ (e.g., Reference Correa et al. (2015)) and the standard deviation of $\sigma_{\log c}\approx 0.13$ Ishiyama et al. (2013). The mean $\bar{c}_{\rm vir}$ has been calibrated at both large (galaxies, clusters) and very small (of Earth-mass size) halos, and found to decreases as a function of $M$ and $z$. Once $c_{\rm vir}$ is drawn from the distribution, it is used to obtain $r_{s}=r_{\rm vir}/c_{\rm vir}$. Finally, $\rho_{s}$ is obtained through the condition of having mass $M$ within $r_{\rm vir}$:

[TABLE]

where the second equality of Equation (19) holds in the case of the NFW profile.

In the case of the subhalos, the procedures above cannot be adopted. This is because they are subject to tidal effects from the host, which strip masses away from the subhalos. However, the regions well inside the scale radius $r_{s}$—because of strong self-gravity—is resilient against the tidal force and hence the annihilation rate hardly changes. These tidal processes, therefore, make the subhalos more concentrated and hence effectively brighter compared with the field halos of the same mass. In many analytical studies in the literature Diemand et al. (2007); Siegal-Gaskins (2008); Anderhalden and Diemand (2013); Sánchez-Conde and Prada (2014), however, the effect of tidal stripping was ignored and the concentration-mass relation of the field halos was adopted, which resulted in underestimate of the annihilation boost factor. This has been pointed out by Reference Bartels and Ando (2015) and will be discussed in Section 4 (see also References Zavala and Afshordi (2016); Moliné et al. (2017)).

3 Estimates of Annihilation Boost with Numerical Simulations

In order to assess the annihilation boosts, one has to have reasonably good ideas on the density profiles $\rho(r)$, the concentration-mass relation,333The concentration-mass relation is defined as the average concentration parameter as a function of halo mass. and the subhalo mass function. Cosmological $N$-body simulations have been a powerful tool for probing all of them because once a halo collapses from initial density fluctuations, it evolves under a strongly nonlinear environment. They have indeed demonstrated that there are a large amount of surviving subhalos (see Figure 1) in halos and halos have cuspy density profiles.

3.1 Subhalo Abundance

Cosmological $N$-body simulations predict that there are many surviving subhalos in host halos as a consequence of hierarchical structure formation. Klypin et al. (1999) and Moore et al. (1999) performed high-resolution cosmological $N$-body simulations for the formation and evolution of galaxy-scale halos. They demonstrated that too many subhalos existed in simulated halos in comparison with the number of observed dwarf galaxies in the Local Group. This discrepancy is known as the “missing satellite problem” that has been investigated by a number of follow-up simulation studies (e.g., References Diemand et al. (2004); Reed et al. (2005); Kase et al. (2007)). Even though it triggered many studies attempting to reduce small-scale structures by imposing other non-CDM candidates such as warm dark matter Kamionkowski and Liddle (2000) and self-interacting dark matter Spergel and Steinhardt (2000), it is also possible to solve it with standard baryonic physics including early reionization Susa and Umemura (2004) and self-regulation of star formation in low-mass halos Benson et al. (2002); Stoehr et al. (2002); Kravtsov et al. (2004); Okamoto et al. (2010). Hence, it is no longer regarded as a serious problem of the CDM model.

These studies suggest a large number of “dark satellites” exist in halos, which do not contain optically visible components such as gases and stars. The population of dark satellites is more abundant in host halos than visible satellite galaxies and could enhance annihilation boosts significantly. To estimate subhalo boosts to annihilation signals accurately, understanding abundance of subhalos as well as their density structure is crucial.

A number of studies have calculated the subhalo mass function in halos using cosmological $N$-body simulations (e.g., References Diemand et al. (2007); Ishiyama et al. (2009); Diemand et al. (2008); Springel et al. (2008)), indicating that it obeys a power law

[TABLE]

where the slope $-\alpha$ ranges from $-2$ to $-1.8$, although no consensus has yet emerged. There is also a large halo-to-halo scatter for the subhalo abundances Ishiyama et al. (2009); Mao et al. (2015). The subhalo abundance at a fixed mass halo depends on their accretion history. Namely, it increases with the mass of halo and decreases with the halo concentration (e.g., References Kravtsov et al. (2004); Gao et al. (2004); Zentner et al. (2005); van den Bosch et al. (2005); Ishiyama et al. (2009); Giocoli et al. (2010); Gao et al. (2011); Ishiyama et al. (2013); Mao et al. (2015)).

Due to the limitation of currently available computational resources, simulations cannot resolve the full hierarchy of subhalos from the smallest to the most massive scales, which ranges more than twenty orders of magnitude in the mass. Even in the highest resolution simulations for galaxy-scale halos, the smallest resolved subhalo mass is around $\sim$10${}^{5}M_{\odot}$ Diemand et al. (2008); Springel et al. (2008), which is still more than ten orders of magnitude more massive than that of the cutoff scale. To study subhalo boosts to annihilation signals, a single-power-law subhalo mass function, Equation (21), is traditionally extrapolated beyond the resolution.

Another approach is to use some analytical models (e.g., References van den Bosch et al. (2005); Giocoli et al. (2008); Yang et al. (2011); Jiang and van den Bosch (2016); Bartels and Ando (2015); Hiroshima et al. (2018)), which can shed light on the resolution issue. Hiroshima et al. (2018) developed a model of the subhalo evolution calibrated with cosmological $N$-body simulations and found that the power-law index of the subhalo mass function is in a rather narrow range between $-2$ and $-1.8$ with a vast range of subhalo mass from $z=0$ to 5. This picture is more or less consistent with the assumption of the subhalo mass function of the single power law. More details on the analytic approach are discussed in the following section.

Note that the annihilation boost factors strongly depend on the underlying subhalo mass function Sánchez-Conde and Prada (2014); Ishiyama (2014); Moliné et al. (2017). Assuming that $\sim$10% of the halo mass is within subhalos, the difference of the boost factors in the Milky-Way-size halo could be as large as a factor of ten between the slope of $-\alpha=-2$ and $-1.9$ Sánchez-Conde and Prada (2014). More extensive simulations are needed to obtain the subhalo mass function in wide mass ranges and also to compare with analytic models Hiroshima et al. (2018).

3.2 Density Profile of Dark Matter Halos

By the end of the 1980’s, it was already known in both analytic Gunn and Gott (1972) and numerical Quinn et al. (1986) studies that the density profiles were described by power-law functions. Reference Dubinski and Carlberg (1991) studied the density profiles of dark matter halos using cosmological $N$-body simulations and argued that the profiles were well described by a Hernquist model Hernquist (1990).

Navarro et al. Navarro et al. (1996, 1997) simulated the structures of CDM halos systematically with masses in the range of galaxy to rich cluster size. They claimed that the radial density profile $\rho(r)$ could be described by a simple universal profile, Equation (3), the so-called NFW profile. They also claimed that the shape of the profile was universal, independent of cosmological parameters, the primordial power spectrum, and the halo mass. Today, the NFW profile has been extensively used to model halos analytically for various purposes.

After the work of NFW, A number of subsequent studies (e.g., References Fukushige and Makino (1997); Moore et al. (1999)) performed simulations with better mass resolutions. Whereas previous studies Navarro et al. (1996, 1997) used only $\sim$10,000 particles, they used $\sim$1,000,000 particles for a halo, and found that the slope was steeper than $-1$. In the original results of the NFW, the numerical two-body relaxation effects due to the small number of particles affected the structures of central regions and led to form a shallower cusp. Higher resolution simulations could resolve more inner structures of halos Ghigna et al. (2000); Jing (2000); Fukushige and Makino (2001); Klypin et al. (2001); Power et al. (2003); Fukushige et al. (2004); Hayashi et al. (2004); Kazantzidis et al. (2006). In most cases, the slope of density became shallower as the radius went inward. A different approach was adopted by Jing and Suto (2000), who used the triaxial model for describing the central structures. Moore et al. (1999) and Diemand et al. (2004) considered a more general profile,

[TABLE]

If $\beta=3$, $\gamma=1$, and $\eta=1$, the profile is the same as NFW.

More recent studies Diemand et al. (2008); Springel et al. (2008); Stadel et al. (2009) archived one of the highest resolution dark matter only simulations for galaxy-size halos with mass resolution better than $10^{4}M_{\odot}$. Their results are in agreement in that the density slope cannot be described by a single power law and the slope is around $-1$ at the radius $\sim$0.001$r_{200}$. Besides, Springel et al. (2008) and Stadel et al. (2009) fitted the density using the Einasto profile Einasto (1965)

[TABLE]

where $\alpha_{E}$ is a free parameter. Note that $r_{s}$ and $\rho_{s}$ are not the same parameters as those in Equation (3).

Although we can obtain the density profile down to the radius $\sim$0.001$r_{200}$, the result does not converge to a single power law. In addition, the physical origin of this flattening towards the center is not understood at all. However, the importance of understanding the central structures is increasing. In particular, if we would like to detect signals from dark matter annihilation, the central structure of the dark matter halo is essential.

The most important parameter to describe the halo profile is the concentration parameter. Assuming the universal NFW profile regardless of the halo mass, the concentration-mass relation gives the annihilation rate as a function of the halo mass. Combined with assumed subhalo mass functions, they enable to estimate the annihilation boost factor. The concentration-mass relation of halos has been widely investigated and a number of fitting functions has been suggested Bullock et al. (2001); Zhao et al. (2003); Maccio’ et al. (2007); Neto et al. (2007); Maccio’ et al. (2008); Zhao et al. (2009); Klypin et al. (2011); Prada et al. (2012); Ishiyama et al. (2013); Sánchez-Conde and Prada (2014); Klypin et al. (2016); Ishiyama (2014); Correa et al. (2015); Ludlow et al. (2016). The concentration shows a weak dependence on the halo mass. The average concentration at fixed halo mass becomes smaller with increasing halo mass because the central density is tightly correlated with the cosmic density at the halo formation epoch, reflecting the hierarchical structure formation Bullock et al. (2001); Wechsler et al. (2002).

Traditionally, the concentration-mass relation has been calibrated with cosmological $N$-body simulations for relatively massive halos ($10^{10}M_{\odot}\lesssim M_{\rm 200}\lesssim 10^{15}M_{\odot}$). Because the mass dependence of the concentration is weak for these mass halos, it is found that a single power-law function, $c\propto M_{200}^{-\alpha_{c}}$, with slope $\alpha_{c}$ in the range of 0.08 to 0.13,  gives reasonable fits Bullock et al. (2001); Maccio’ et al. (2007); Neto et al. (2007); Maccio’ et al. (2008); Klypin et al. (2011). However, the dependence gradually becomes weaker toward less massive halos, and a clear flattening emerges Ishiyama et al. (2013); Sánchez-Conde and Prada (2014); Ishiyama (2014); Correa et al. (2015); Hellwing et al. (2016); Ludlow et al. (2016); Pilipenko et al. (2017), ruling out single power-law concentration-mass relation for the full hierarchy of subhalos.

These fitting functions are valid for the NFW density profile. More generally, the concentration can be defined independently of the density profile and the subhalo mass as (e.g.,  Diemand et al. (2007); Springel et al. (2008); Moliné et al. (2017))

[TABLE]

where $r_{\rm max}$ is the radius at which the circular velocity reaches its maximum value $V_{\rm max}$. This definition is also used to estimate the annihilation boost factor (e.g.,  Moliné et al. (2017)).

Even with the highest resolution simulations for galaxy-scale halos, the smallest resolved subhalo mass is around $\sim$10${}^{5}M_{\odot}$ Diemand et al. (2008); Springel et al. (2008). To estimate the annihilation boost from the full hierarchy of subhalos, we have to make some assumption of the concentration at unresolved scales, which has a significant impact on the result. One approach is extrapolating single power-law fittings to the smallest scale beyond the mass range calibrated with simulations, although the literature including the above cautions the risk of such extrapolations. With such extrapolations, the concentration of the smallest halo can reach more than 100, substantially enhancing the annihilation boost. A number of studies have computed the concentration in such a manner and the resulting boost factor is a few hundreds for Milky-Way halos Springel et al. (2008), and $\sim$1000 for cluster-scale halos Pinzke et al. (2011); Gao et al. (2012); Anderson et al. (2016).

Another approach is adopting analytic models or fitting functions that can reproduce flattening of the concentration-mass relation

(e.g., Prada et al. (2012); Sánchez-Conde and Prada (2014); Klypin et al. (2016); Correa et al. (2015)).

In contrast to using the power-law extrapolation, the resulting boost factor is rather modest, three to a few tens Sánchez-Conde and Prada (2014); Ishiyama (2014); Stref and Lavalle (2017) for Milky-Way halo, and less than $\sim$100 for cluster-scale halos Sánchez-Conde and Prada (2014); Ishiyama (2014); Anderson et al. (2016).

The density profile at fixed halo mass shows a significant halo-to-halo scatter Reed et al. (2011), possibly making a big impact on the annihilation signal. Inferring from the cosmological Millennium simulation Springel et al. (2005), the effect of this non-universality on the annihilation flux is a factor of $\sim$3 Reed et al. (2011), which indicates that the uncertainty of the concentration-mass relation for low-mass halos has a more significant effect.

These discussions are based on the universal density profile and the concentration-mass relation for field halos. There is a concern that whether or not we can apply the universal NFW profile for the full hierarchy of halos and subhalos beyond the range that cosmological simulations have been able to tackle. We discuss this issue in Section 3.4. More importantly, we have to use the concentration-mass relation for subhalos, not field halos. We also discuss this issue in the following section.

3.3 Density Profile of Dark Matter Subhalos

Density structures of subhalos are more challenging to be investigated than field halos because it requires much higher resolutions. Therefore, to evaluate the subhalo contribution to the annihilation signals, the universal NFW profile and the concentration-mass relation for field halos have been historically assumed to be the same for subhalos as a first approximation, although the underlying assumption is not well studied. Complex physical mechanisms relevant to subhalos could change their original density profiles, such as the tidal effect from host halos, the encounter with other subhalos, and denser environment than the field.

Cosmological simulations have been suggesting that the density profile of subhalos is cuspy in analogy with field halos. On the other hand, the average concentration of subhalos tend to be higher than those of field halos (e.g., Bullock et al. (2001); Ghigna et al. (2000); Diemand et al. (2007); Springel et al. (2008); Moliné et al. (2017)). For example, Bullock et al. (2001) showed that subhalos and halos tend to be more concentrated in dense environments than in the field, and the scatter of concentrations is larger. This result was taken into account to estimate the gamma-ray flux from dark matter annihilation (e.g., Ullio et al. (2002)). Diemand et al. (2007) showed that outer regions of subhalos tend to be tidally stripped by host halos, which gives higher concentrations. These results suggest that both earlier formation of halos/subhalos in dense environments and tidal effect are responsible for the increased concentration. Pieri et al. (2011) derived the concentration-mass relation of subhalos in Milky-Way-size halos by analyzing high resolution cosmological simulations Diemand et al. (2008); Springel et al. (2008) and showed that it depends on the location of subhalos relative to host halos. Subhalos have considerably large concentrations near the center than at the edge of host halos. Moliné et al. (2017) quantified the concentration of subhalos in Milky-Way-size halos as a function of not only subhalo mass but distance from host halo center, and found a factor 2–3 enhancement of the boost factor compared to the estimation that relied on the concentration-mass relation of field halos (see also Figure 2).

As shown in the literature in the above, higher concentrations of subhalos than field halos could have a big impact on the annihilation boost. However, van den Bosch and Ogiya (2018) argued that subhalos even in state-of-the-art cosmological simulations suffer from excessive mass loss and artificial tidal disruption due to inadequately large force softening (see also Peñarrubia et al. (2010); van den Bosch et al. (2018)). If that is the case, it might be possible that subhalos have larger concentrations than those ever considered. These issues should be addressed by extremely high resolution cosmological $N$-body simulations and analytic models (e.g., Hiroshima et al. (2018)).

3.4 Density Profile of Dark Matter Halos Near Cutoff Scales

In the CDM framework, smaller halos collapse first, and then they merge into more massive halos. Since the smallest halos contain no subhalos, their central structures might entirely differ from that observed in more massive halos. If the dark matter particle is the lightest supersymmetric particle such as the neutralino, the smallest halo mass is predicted to be around the Earth mass Profumo et al. (2006); Bringmann (2009); Diamanti et al. (2015). Such halos are sometimes referred to as “microhalos.”

The density profiles of the microhalos have been investigated using cosmological $N$-body simulations Diemand et al. (2005); Ishiyama et al. (2010); Anderhalden and Diemand (2013); Ishiyama (2014); Angulo et al. (2017). Diemand et al. (2005) simulated the formation of Earth-mass microhalos by means of cosmological $N$-body simulation. They claimed that a single power law could describe the density profiles of microhalos, $\rho(r)\propto r^{-\gamma}$, with a slope $\gamma$ in the range of 1.5 to 2. As a consequence of such steep slope, most microhalos could not be completely destructed by the Galactic tide and encounters with stars, even in the Galactic center.

Ishiyama et al. (2010) have performed $N$-body simulations with much higher resolution and showed that the density profile of microhalos had steeper cusps than the NFW profile. The central density scales as $\rho(r)\propto r^{-1.5}$, which is supported by follow-up cosmological simulations Anderhalden and Diemand (2013); Ishiyama (2014); Angulo et al. (2017) and cold-collapse simulations Ogiya and Hahn (2018). Ishiyama (2014) has also shown that the cusp slope gradually becomes shallower with increasing halo mass. Major merger of halos is responsible for the flattening, indicating that the process of violent relaxation plays a key role (see also Ogiya et al. (2016); Angulo et al. (2017)). Similar density structures are observed in recent simulations of ultracompact minihalos Gosenca et al. (2017); Delos et al. (2018a, b) and warm dark matter Polisensky and Ricotti (2015).

The self-similar gravitational collapse models (e.g., Gunn and Gott (1972); Gott (1975); Fillmore and Goldreich (1984); Bertschinger (1985); Gurevich and Zybin (1995)) can also give hints to understand the main physical origin of such steeper cusps, because the smallest halos do not contain smaller density fluctuations by definition and collapse from initially overdense patches.

Such microhalos with steep cusps can cause a significant effect on indirect dark matter searches. Ishiyama et al. (2010) argued that the central parts of microhalos could survive against the encounters with stars except in the very Galactic center. The nearest microhalos could be observable via gamma rays from dark matter annihilation, with usually large proper motions of $\sim$0.2 deg yr*-1*, which are, however, stringently constrained with the diffuse gamma-ray backgrond Ando et al. (2008). Gravitational perturbations to the millisecond pulsars might be detectable with future observations by pulsar timing arrays Ishiyama et al. (2010); Baghram et al. (2011); Kashiyama and Oguri (2018); Dror et al. (2019). Anderhalden and Diemand (2013) have assumed a transition from the NFW to steeper cusps at scales corresponding to $\sim$100 times more massive than the cutoff and have found that such profiles can enhance moderately the annihilation boost of a Milky-Way-size halo by 5–12%. They also have found that concentrations of microhalos are consistent with a toy model proposed by Bullock et al. (2001).

Ishiyama (2014) showed that the steeper inner cusps of halos in the smallest scale and near the cutoff scale could increase the annihilation rate of a Milky-Way-size halo by 12–67%, compared with estimates adopting the universal NFW profile and an empirical concentration-mass relation Sánchez-Conde and Prada (2014) (see Figure 3). The value, however, depends strongly on the adopted subhalo mass function and concentration model. They have found that concentrations near the free-streaming scale show little dependence on the halo mass and corresponding conventional NFW concentrations are 60–70, consistent with the picture that the mass dependence is gradually becoming weaker toward less massive halos (e.g., Ishiyama et al. (2013); Sánchez-Conde and Prada (2014); Ishiyama (2014); Correa et al. (2015); Hellwing et al. (2016); Ludlow et al. (2016)), ruling out a single power-law concentration-mass relation.

As shown in the literature above, steep density cusps of halos near the free-streaming scale have an impact on the annihilation boost. However, these studies rely on the density structure seen in field halos, not subhalos. It is also important to quantify the structures of subhalos near the free-streaming scale by larger simulations. Another concern is that the cutoff in the matter power spectrum should suppress the number of subhalos near the free-streaming scale, which should weaken the annihilation signal. However, the shape of the mass function near the free-streaming scale is not understood well for the neutralino dark matter. The structure of subhalos and the subhalo mass function near the free-streaming scale should be explored by larger volume cosmological $N$-body simulations.

4 Semi-Analytic Approaches

4.1 Models Based on Structure Formation and Tidal Evolution

In an analytical approach in Section 2, the subhalo luminosity $L_{\rm sh}$ is characterized with the mass of the subhalo and the redshift of interest; see, for example, Equation (10). The mass and redshift, however, are not the only quantities that fully characterize the subhalo properties. Indeed, they depend on the accretion history and mass loss after they fall onto their host halo, that is, two subhalos that have the identical mass could have formed with different masses and accreted at different redshifts, evolved down to $z=0$ reaching the same mass. Bartels and Ando (2015) and Hiroshima et al. (2018) developed an analytical prescription to take these effects into account, which we follow in this section.

A subhalo is characterized with its mass and redshift when it accreted onto its host, ($m_{a}$, $z_{a}$). The concentration parameter $c_{a}$ is drawn from the log-normal distribution with mean $\bar{c}_{a}(m_{a},z_{a})$ Correa et al. (2015) and $\sigma_{\log c}=0.13$ Ishiyama et al. (2013). Since the subhalo was a field halo when it just accreted, one can use the relations in Section 2.2 to obtain $r_{s,a}$ and $\rho_{s,a}$ for the NFW profile.

After the accretion, the subhalos evolve by losing their mass through tidal forces. The mass-loss rate is typically characterized by a dynamical timascale at the redshift $z$,

[TABLE]

as follows Jiang and van den Bosch (2016):

[TABLE]

where $H(z)=H_{0}[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}]^{1/2}$, $m(z)$ and $M(z)$ are the subhalo and host-halo masses at $z$, respectively. Following Jiang and van den Bosch (2016), Hiroshima et al. (2018) adopted simple Monte Carlo simulations to estimate $\dot{m}$ based on the assumption that the subhalo loses all the masses beyond its tidal radius in one complete orbit at its peri-center passage. While Jiang and van den Bosch (2016) found $A=0.81$ and $\zeta=0.04$, Hiroshima et al. (2018) extended the mass and redshift ranges of applicability and found that these parameters are weakly dependent on both $M$ and $z$:

[TABLE]

One can solve Equation (26) to obtain the subhalo mass at a redshift of interest $z$, $m(z)$, with a boundary condition of $m(z_{a})=m_{a}$. For the evolution of the host, $M(z)$, Hiroshima et al. (2018) adopted a fitting formula given by Correa et al. (2015).

The subhalo density profile after accretion is also well described with the NFW profile with a sharp truncation at $r_{t}$:

[TABLE]

This is indeed a good approximation found in the simulations Springel et al. (2008). In addition to $r_{t}$, Peñarrubia et al. (2010) found that the internal structure changes. If the inner profile is $\propto r^{-1}$ just like NFW, the maximum circular velocty $V_{\rm max}$ and its corresponding radius $r_{\rm max}$ evolve as

[TABLE]

respectively. After computing $V_{\rm max}$ and $r_{\rm max}$ at $z$, one can convert them to $\rho_{s}$ and $r_{s}$ through

[TABLE]

which are valid for the NFW profile. Finally by solving the condition

[TABLE]

the truncation radius $r_{t}$ is obtained. Hiroshima et al. (2018) omitted subhalos with $r_{t}<0.77r_{s}$ from the subsequent calculations assuming that they were tidally disrupted Hayashi et al. (2003). This criterion, however, might be a numerical artifact van den Bosch et al. (2018). Either case, Hiroshima et al. (2018) checked that whether one implements this condition or not did not have impact on the results of, for example, subhalo mass functions.

Thus, given $(m_{a},z_{a},c_{a})$, one can obtain all the subhalo parameters after the evolution, $(m,r_{s},\rho_{s},r_{t})$, in a deterministic manner. The differential number of subhalos accreted onto a host with a mass $m_{a}$ and at redshift $z_{a}$, $d^{2}N_{\rm sh}/(dm_{a}dz_{a})$, is given by the excursion set or the extended Press-Schechter formalism Bond et al. (1991). Especially Yang et al. (2011) obtained analytical formulation for the distribution that provides good fit to the numerical simulation data over a large range of $m/M$ and $z$. Hiroshima et al. (2018) adopted their model III.

The subhalo mass function is obtained as

[TABLE]

where $P(c_{a}|m_{a},z_{a})$ is the probability distrbution for $c_{a}$ given $m_{a}$ and $z_{a}$, for which Hiroshima et al. (2018) adopted the log-normal distribution with the mean $\bar{c}_{a}(m_{a},z_{a})$ Correa et al. (2015) and the standard deviation $\sigma_{\log c_{a}}=0.13$ Ishiyama et al. (2013). We show the subhalo mass functions obtained with Equation (35) for various values of $M$ and $z$ in Figure 4, where comparison is made with simulation results of similar host halos.

Halos and subhalos formed in these simulations were identified with ROCKSTAR phase space halo finder Behroozi et al. (2013). The bound mass is used as the subhalo mass, which nearly corresponds to the tidal mass Moliné et al. (2017).

For all these halos, one can see remarkable agreement between the analytic model and the corresponding simulation results in resolved regimes. Successfully reproducing behaviors at resolved regimes, this analytic model is able to make reliable predictions of the subhalo mass functions below resolutions of the numerical simulations, without relying on extrapolating a single power-law functions, from which most of the previous studies in the literature had to suffer. The subhalo mass fraction is then obtained as

[TABLE]

and is shown in Figure 4 (bottom right) for various values of redshifts. The subhalo mass fraction is found to increase as a function of $M$ and $z$. At higher redshifts, since there is shorter time for the subhalos to experience tidal mass loss, $f_{\rm sh}$ is larger. Again, a good agreement in $f_{\rm sh}$ is found between the analytic model and the simulation results by Giocoli et al. (2010).

The annihilation boost factor is then

[TABLE]

which is to be compared with Equation (10) that was derived with a simpler (and unrealistic) discussion. The superscript (0) represents the quantity in the absense of sub-subhalos and beyond. The subhalo luminosity, $L_{\rm sh}^{(0)}(m_{a},z_{a},c_{a}|M,z)$, is proportional to the volume integral of density squared $\rho_{\rm sh}^{2}(r)$ out to the truncation radius,

[TABLE]

where $\rho_{s,{\rm sh}}$, $r_{s,{\rm sh}}$ and $r_{t,{\rm sh}}$ are functions of $(m_{a},z_{a},c_{a})$ as well as $M$ and $z$.

Then, the effect of subn-subhalos (for $n\geq 1$) can be estimated iteratively. At $n$th iteraction, when a subhalo accreted onto its host at $z_{a}$ with $m_{a}$, it is assigned a sub-subhalo boost factor $B_{\rm sh}^{(n-1)}(m_{a},z_{a})$. After the accretion, the outer region of the subhalo is stripped away by the tidal force and thus all the sub-subhalos within this stripped region will disappear, reducing the sub-subhalo boost accordingly. Hiroshima et al. (2018) assumed that the sub-subhalos were distributed within the subhalo following $n_{\rm ssh}(r)\propto[1+(r/r_{s})^{2}]^{-3/2}$. The luminosity due to sub-subhalos within a radius $r$ is therefore proportional to their enclosed number

[TABLE]

and it gets suppressed by a factor of $N_{\rm ssh}(<r_{t}|r_{s})/N_{\rm ssh}(<r_{\rm vir}|r_{s,a})$ due to the tidal stripping.444We note that in estimating the effect of subn-subhalos in the boost factors, Reference Hiroshima et al. (2018) ignored the changes of $\rho_{s}$ and $r_{s}$ and hence did not include the factor of $r_{s}^{3}$ in Equation (39) and $\rho_{s}r_{s}^{3}$ in Equation (40). In addition, in Equation (43), they multiplied $L_{\rm sh}^{(0)}$ by a factor of $1+B_{\rm ssh}^{(n)}$ instead of $1-f_{\rm ssh}^{2}+B_{\rm ssh}^{(n)}$. We correct for all these effects in this review. The luminosity due to the smooth component also decreases as $L_{\rm sh}^{(0)}(<r_{t}|\rho_{s},r_{s})/L_{\rm sh}^{(0)}(<r_{\rm vir}|\rho_{s,a},r_{s,a})$, where

[TABLE]

Thus the sub-subhalo boost after the $n$th iteration, $B_{\rm ssh}^{(n)}$ is obtained by

[TABLE]

Similarly, the sub-subhalo mass fraction $f_{\rm ssh}$ is obtained by

[TABLE]

where $f_{\rm sh}(m_{a},z_{a})$ is obtained with Equations (36) and $m_{\rm sm}(<r|\rho_{s},r_{s})\propto\rho_{s}r_{s}^{3}f(r/r_{s})$ is the enclosed mass within $r$ of the smooth component of the subhalo. The subhalo boost factor after $n$th iteration is obtained with Equation (37) by replacing $L_{\rm sh}^{(0)}$ with $[1-f_{\rm ssh}^{2}+B_{\rm ssh}^{(n)}]L_{\rm sh}^{(0)}$ [see discussions below Equation (15)]:

[TABLE]

The host luminosity in the absence of the subhalos $L_{\rm host,0}(M,z)$ is defined by marginalizing over the concentration parameter $c_{\rm vir}$:

[TABLE]

with the log-normal distribution $P(c_{\rm vir}|M,z)$.

Figure 5 shows the subhalo boost factors $B_{\rm sh}$ as a function of the host mass $M$ at various redshifts $z$ (top left). The boost factors are on the order of unity, while it can be as larger as $\sim$5 for cluster-size halos. It is also noted that they are larger at higher redshifts, because the subhalos have less time to be disrupted. The top right panel of Figure 5 shows the effect of subn-subhalos, which is saturated after the second iteration. The contribution to the boost factors due to sub-subhalos and beyond is $\lesssim$10% for the hosts with $M_{\rm host}\geq 10^{13}M_{\odot}$. The bottom left panel of Figure 5 shows the luminosity ratio $L_{\rm total}/L_{\rm host,0}=1-f_{\rm sh}^{2}+B_{\rm sh}$ (Equation 15) as a function of the host masses for various values of the redshifts. The bottom right panel of Figure 5 shows comparison with the results of the other work Gao et al. (2012); Sánchez-Conde and Prada (2014); Moliné et al. (2017). We note that the analytic models do not rely on the subhalo mass function prepared separately, as the models can provide them in a self-consistent manner. The resulting boost factors are, however, found to be more modest than the previous results. This is mainly because the subhalo mass function adopted in the literature is larger than the predictions of the analytic models. However, they might be larger because of halo-to-halo variance. See discrepancy between predictions of the subhalo mass function for the $1.8\times 10^{12}M_{\odot}$ halo by Hiroshima et al. (2018) and the result of Springel et al. (2008) shown in the top left panel of Figure 4.

Finally, for convenience of the reader who might be interested in using the results without going into details of the formalism, we provide fitting functions for both the subhalo mass functions and the annihilation boost factors. They are summarized in Appendix A.

4.2 Models for Self-Similar Subhalos

Assuming a self-similarity of the subhalos, Kamionkowski and Koushiappas (2008) developed a fully analytic formulation for the probability distribution function of the dark matter density, $P(\rho)$. Then Kamionkowski et al. (2010) applied the formulation to the result of cosmological $N$-body simulations, to obtain the fitting function for the Galactic local boost factor at Galactocentric radius $r$:

[TABLE]

where $f_{\rm sm}(r)$ is the volume fraction occupied by the smooth component and $\rho_{\rm max}$ is the highest dark matter subhalo density. Through the calibration with the numerical simulations, Kamionkowski et al. (2010) found $\delta_{f}=0.2$, $\alpha_{K}=0$ and that the subhalo fraction was given by

[TABLE]

where $\kappa=0.007$. Fornasa et al. (2013) then suggested a larger value of $\kappa=0.15$–0.2 to obtain a larger boost consistent with earlier work Springel et al. (2008); Pinzke et al. (2011); Gao et al. (2012). The maximum subhalo density $\rho_{\rm max}$ is estimated as $\rho_{\rm max}=[c^{3}/f(c)/12]200\rho_{c}(z_{f})$, where $c$ and $z_{f}$ are the concentration parameter and collapse redshift of the smallest halos. Kamionkowski et al. (2010) adopted $c=3.5$ and $z_{f}=40$ and $\rho_{\rm max}=80$ GeV cm*-3*. On the other hand, Ng et al. (2014) obtained a smaller $\rho_{\rm max}\approx 20$ GeV cm*-3* even with a very small cutoff masses of $M_{\rm min}=10^{-12}M_{\odot}$ Within the virial radius of the Milky-Way halo, the subhalo boost factor for dark matter annihilation is found to be no greater than $\sim$10 Kamionkowski et al. (2010).

4.3 Universal Clustering of Dark Matter in Phase Space

Zavala and Afshordi (2016) investigated the behavior of dark matter particles that belong to the halo substructure in the phase space of distance and velocity. Reconstructing the phase-space distribution using the Aquarius numerical simulations Springel et al. (2008), they found universality of the coarse-grained phase-space distribution ranging from dwarfs to clusters of galaxies. They developed physically motivated models based on the stable clustering hypothesis, spherical collapses and the tidal stripping of the subhalos and applied to the obtained phase-space distribution data from the simulations to find a good agreement. Then, they computed the nonlinear matter power spectrum based on the halo model Seljak (2000) down to very small free-streaming cutoff scales. Based on the power spectrum, they obtained the subhalo boost factor greater than $\sim$30–100 for the Milky-Way size halos, which is significantly larger than the values obtained with other analytic work Kamionkowski et al. (2010); Bartels and Ando (2015); Hiroshima et al. (2018). This discrepancy might come from the treatment at very small scales, where it is very hard to calibrate the analytic models against the results of the numerical simulations.

5 Conclusions

It is established that dark matter halos are made up with lots of substructures. Especially in the cold dark matter scenario, small structures form first and merge and accrete to create larger halos. If the dark matter is made of weakly interacting massive particles such as the supersymmetric neutralino, the smallest halos can be as light as or even lighter than the Earth. The rate of dark matter annihilation and hence its signatures such as gamma-ray fluxes are proportional to dark matter density squared and therefore, having small-scale “clumpy” subhalos will boost the signals.

It is, however, of an extraordinary challenge to estimate this subhalo boost factor, that is, the ratio of luminosity from dark matter annihilation in the subhalos to that in the smoothly distributed main component. This is mainly because subhalos of all the mass scales ranging from Earth to galaxy masses can contribute to the boost factor nearly equally per decade in mass. In this review, we cover recent progress to overcome this issue to obtain realistic and unbiased estimates on the subhalo boost factor that will impact on interpretation of the measurements on particle physics parameters such as the annihilation cross section. While cosmological $N$-body simulations provide the most accurate avenue to study structures in highly nonlinear regime, it is inevitably limited by the numerical resolution. Even the state-of-the-art $N$-body simulations Springel et al. (2008); Diemand et al. (2006) can resolve subhalos ranging for only several decades, which is still more than ten orders of magnitude in short to resolve all the subhalos. Therefore, the boost estimates have to rely on extrapolation of the subhalo properties such as its mass function and concentration parameter, which are often well described with power-law functions. Danger of extrapolating trends found in resolved regime for other many orders of magnitude had been widely acknowledged but nevertheless, it was found that the estimates based on such extrapolations tended to give very large amount of boost factor of $\sim$100 ($\sim$1000) for galaxy (cluster) size halos Springel et al. (2008); Gao et al. (2012).

As a complementary approach, analytic models have been investigated. They are based on self similar propertiese of the subhalos Kamionkowski and Koushiappas (2008); Kamionkowski et al. (2010), universal phase-space distribution Zavala and Afshordi (2016) and extended Press-Schechter formalism combined with tidal stripping modeling Bartels and Ando (2015); Hiroshima et al. (2018). (More recent numerical approach also adopts the concentration-mass relation calibrated for the subhalos in order to take the tidal effects into account Moliné et al. (2017).) Most importantly, these are all calibrated with the cosmological $N$-body simulations at resolved regimes and proven to reproduce the simulation results such as the subhalo mass functions. For example, the most recent analytic models by Reference Hiroshima et al. (2018) predict the subhalo mass functions for various host masses and redshifts, which are found to be in good agreement with the simulation results Springel et al. (2008); Diemand et al. (2006); Giocoli et al. (2010); Ishiyama et al. (2015). The annihilation boost factors based on these analytic models tend to be more modest, $\mathcal{O}(1)$ for galaxy-size halos and $\lesssim\mathcal{O}(10)$ for cluster-size halos. However, none of these models have been tested against simulations at very small host halos that are less massive than $10^{6}M_{\odot}$. Simulations of microhalos with $10^{-6}M_{\odot}$ suggest cuspier profiles towards halo centers such as $r^{-1.5}$ Ishiyama (2014) and if this is the case for the subhalos too, it would boost the annihilation rate further.

It is known that including baryons in the simulations affects properties of subhalos such as spatial distribution and density profiles (e.g., References Zhu et al. (2016); Errani et al. (2017); Chua et al. (2017); Kelley et al. (2018)) and hence there might be some effect on the annihilation boost factors. This, however, remains largely unexplored and has to wait for future progress. However, since the subhalos of all masses ranging down to about the Earth mass contribute to the boost factors and the baryons will likely affect only halos of dwarf galaxies or larger, we anticipate that it is not a very important effect for the annihilation boost factors.

The subhalo boosts directly impact the obtained upper limits on the dark matter annihilation cross section from the extragalactic halo observations. Therefore, to obtain the most accurate estimates of the boost factor by reducing uncertainties on structure formation at small scales as well as the physics of tidal stripping is of extreme importance for the indirect searches for particle dark matter through self-annihilation with the current and near future observations of high-energy gamma-rays, neutrinos and charged cosmic rays.

\authorcontributions

SA wrote most of the review sections except for Section 3, which was contributed by TI including Figures 1–3. NH prepared Figures 4 and 5 and wrote Appendix A.

Acknowledgements.

This work was supported by JSPS KAKENHI Grant Numbers JP17H04836 (SA), JP18H04340 (SA and NH), JP18H04578 (SA), JP15H01030 (TI), JP17H04828 (TI), JP17H01101 (TI) and JP18H04337 (TI). TI has been supported by MEXT as “Priority Issue on Post-K computer” (Elucidation of the Fundamental Laws and Evolution of the Universe) and JICFuS. Numerical computations were partially carried out on the K computer at the RIKEN Advanced Institute for Computational Science (Proposal numbers hp150226, hp160212, hp170231, hp180180), and Aterui supercomputer at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. \conflictsofinterestThe authors declare no conflict of interest. \appendixtitlesyes \appendixsectionsmultiple

Appendix A Fitting Formulae

In this appendix, we provide fitting functions obtained with the analytical calculation in Hiroshima et al. (2018) that covers more than twenty orders of magnitude in the mass range and the redshift up to $\sim$10. The mass function of the subhalo $m^{2}dN_{\rm sh}/dm$, the luminosity of the subhalo separated from the particle physics factors $L_{\rm sh}$, and the boost factor $B_{\rm sh}=L_{\rm sh}/L_{\rm host,0}$ [see Equation (10)] as functions of the host mass and the redshift are provided here. We also summarize fitting functions for the boost factor $B_{\rm sh}$ at $z=0$ in the literature. Note that the host mass is always measured in units of the solar mass ($M_{\odot}$) in this appendix.

A.1 Subhalo Mass Function

The fitting formula is written in the follwing form:

[TABLE]

introducing a cutoff mass $m_{c}$; $a$, $b$, and $m_{c}$ in Equation (47) are functions of the host mass and the redshift:

[TABLE]

In Figure 6, we show the comparison between the mass function obtained in analytical calculations Hiroshima et al. (2018) and Equation (47). By integrating Equation (47), we obtain the mass fraction shown in Figure 7.

A.2 Subhalo Luminosity

The luminosity of a subhalo is written as

[TABLE]

assuming a subhalo of NFW profile. In this section, the characteristic density $\rho_{s}$ and the scale radius $r_{s}$ are measured in units of g/cm3 and cm, respectively. For simplicity, we do not show the integration over the distribution of the virial concentration parameter $P(c_{\rm vir}|M,z)$ in the above expression. The constant of proportionality is detemined by fixing the particle physics model. However, we do not include it in the following expression since it cancels in the calculation of the boost factor $B_{\rm sh}$ by taking the ratio of the host and subhalo luminosity $L_{\rm sh}/L_{\rm host}$.

The fitting formula for the luminosity takes the form of

[TABLE]

where

[TABLE]

A.3 Annihilation Boost Factor

The boost factor is sensitive to the models of the concentration-mass relation. We provide fitting functions for two different concentration-mass relation models here. One coresponds to the canonical model in Reference Hiroshima et al. (2018) which assumes the concentration-mass relation derived in Reference Correa et al. (2015). The other corresponds to the concentration-mass relation in Reference Okoli and Afshordi (2016). For both cases, the fitting function of the boost factor is written in a combination of two sigmoid functions, $f(x)=(1+e^{-ax})^{-1}$,

[TABLE]

Funcitons $X$, $Y$, $a$, $b$, $c$, $m_{1}$, and $m_{2}$ depend on the redshift but they do not on the host mass.

• •

For Correa’s concentration Correa et al. (2015)

[TABLE]

• •

For Okoli’s concentration Okoli and Afshordi (2016)

[TABLE]

All of these formulae are applicable for hosts at arbitrary redshifts up to $z\sim 7$.

A.4 Fitting Functions for the Boost Factor in the Literature

Several works provide the fitting function for the boost factor at $z=0$. We summarize functions provided in Gao et al. (2012); Sánchez-Conde and Prada (2014), and Moliné et al. (2017) here.

• •

Gao et al. (2012) have analyzed cluster scale halo of $M_{\rm host}=[5,20]\times 10^{14}h^{-1}M_{\odot}$ in the Phoenix project Gao et al. (2012). Subhalos down to $m\sim 10^{6}$ can be resolved in their calculations.

[TABLE]

• •

Sánchez-Conde and Prada (2014) derive the boost factor based on the concentration-mass relation in Reference Prada et al. (2012). The fitting function is provided for their fiducial model assuming the minimum halo mass to be $M_{\rm min}=10^{-6}M_{\odot}$ and the subhalo mass function $dN/dm\propto m^{-2}$. Each subhalo is assumed to be a field halo.

[TABLE]

with

[TABLE]

• •

Moliné et al. (2017) derive the boost factor taking the dependence of the survived halo properties on the distance from the host, that is, the host potential. The form of the function is similar to that of Sánchez-Conde and Prada (2014) but adopitng different log bases:

[TABLE]

and parameters $b_{i}$ are

[TABLE]

for $\alpha=2$ and

[TABLE]

for $\alpha=1.9$.

\reftitle

References