Halo Substructure Boosts to the Signatures of Dark Matter Annihilation

: 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 (cid:38) 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.


Introduction
One of the most popular candidates for dark matter is weakly interacting massive particles (WIMPs) [1,2].They are motivated by beyond-the-standard-model physics such as supersymmetry [3] or universal extra-dimensions [4], although the non-discovery of new physics at the TeV scale with the Large Hadron Collider puts these models to serious test [5].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 [6].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 but with very large range of possible values of ∼10 −12 -10 −3 M [7][8][9][10][11][12][13][14].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 ρ 2 ≥ ρ 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 .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 ) [37,51,52], one can resolve only down to 10 4 -10 5 M , 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 [53] 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 [54] 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 [55]), 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 [35,[46][47][48][49][50].This relatively mild amount of the annihilation boost makes the prospect of indirect dark matter searches less promising compared with earlier more optimistic predictions [36,40,41,56].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 [47] 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.

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).

Subhalo Boost Factor
The rate of dark matter annihilation is proportional to dark matter density squared, ρ 2 χ , where χ represents the dark matter particle.In the presence of substructure, ρ χ is divided into two terms: where ρ sm and ρ sh represent smooth and subhalo components, respectively.(In the following, we omit the subscript χ.)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 We assume that the smooth component ρ sm is characterized by the following Navarro-Frenk-White (NFW) profile [57,58]: where ρ s is a characteristic density and r s is a scale radius.These parameters, ρ s and r s , are evaluated such that the volume integral of ρ NFW yields the total halo mass M, and thus we have , where f sh is defined as the mass fraction in the subhalos.The first term is then simply where V = 4πr 3 vir /3 is the volume of the host out to its virial radius r vir , c vir ≡ r vir /r s is the concentration parameter.The parameters characterizing the host profile-ρ s , r s , and c vir -are all functions of M and z. 1  Next, we evaluate the second term of Equation ( 2), ρ 2 sh (x) .We characterize each subhalo i with the location of its center x i and mass m i .Density due to all the subhalos at a coordinate x is written as a sum of the density profile around the seed of each subhalo, that is, where δ N D is the N-dimensional Dirac delta function, and u sh (r|m) defines the density profile of the subhalo with mass m and is normalized to one after the volume integral. 2We define the ensemble average of the product of these delta functions as We note, however, that the concentration c vir has a scatter, which is often characterized by a log-normal distribution, whose mean cvir is the function of M and z.We will include this in the latter sections. 2 For 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.
its volume integral over the host halo as and call both dn sh /dm and dN sh /dm the subhalo mass function.We also obtain the mass fraction in the subhalos as By multiplying Equation ( 5) by itself and taking both the ensemble and the volume averages, we have where at the last equality we adopted the NFW function for the subhalo density profile mu sh (r|m) with the scale radius r s,sh , characteristic density ρ s,sh , and tidal truncation radius r t,sh ≡ c t,sh r s,sh beyond which the subhalo density abruptly decreases to zero.At the third equality of Equation ( 9), we ignored the term arising from j = i as we evaluate the quantity at one point 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 [27,38] 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 where the subscript 0 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 [47,[59][60][61][62]. Following a similar procedure as in Equation ( 9), we have 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 sh /dm out of x integration adopting x ≈ x as its spatial variable.Second, we performed volume integral for u(x − x |m) over x 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: where P msh (x)d 3 x represents the probability of finding a subhalo in a volume element d 3 x around x.With this and Equation (8), we have For simplicity, we assume that the subhalos are distributed following the smooth NFW component.In this case, we have ρ sm (x) = (1 − f sh )MP sh (x), and The luminosity from the smooth component in the presence of the subhalos (L sm ) is related to the host luminosity in the subhalos' absence via L sm = (1 − f sh ) 2 L host,0 , because the density in the smooth component gets depleted by a factor of 1 − f sh , if there are subhalos.Thus, the total luminosity from both the smooth component and the subhalos are given by Often, L total /L 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 sh .Thus, we replace L sh with (1 − f 2 ssh + B ssh )L sh , where the subscript "ssh" represents the contribution from the sub-subhalos.If the subhalo mass fraction f sh and the boost factor B sh depend only on the host mass, then one can assume f ssh (m) = f sh (m) and B ssh (m) = B sh (m), and repeat the calculations in an iterative manner.See, however, Section 4.1 for a more realistic treatment.

Characterization of Dark Matter Halos
We shall discuss the density profile of dark matter halos that are characterized by the virial radius r vir , the scale radius r s , and the characteristic density ρ 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: is the present critical density, H 0 = 100h km s −1 Mpc −1 is the Hubble constant, Ω m and Ω Λ are the density parameters for matter and the cosmological constant, respectively.In CDM cosmology with the cosmological constant, this critical value is given as [63] ∆ vir (z) = 18π 2 + 82d(z) − 39d 2 (z), (16) where Given the virial mass M and the redshift z of the halo of interest, r vir is therefore obtained by solving Alternatively, one can define M 200 and r 200 via M 200 is often adopted to define halo masses in N-body simulations.
The concentration parameter c vir ≡ r vir /r s (or c 200 ≡ 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 cvir (M, z) (e.g., Reference [64]) and the standard deviation of σ log c ≈ 0.13 [65].The mean cvir 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 vir is drawn from the distribution, it is used to obtain r s = r vir /c vir .Finally, ρ s is obtained through the condition of having mass M within r vir : 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 [29,32,44,66], 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 [46] and will be discussed in Section 4 (see also References [48,49]).

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 ρ(r), the concentration-mass relation, 3 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.
The concentration-mass relation is defined as the average concentration parameter as a function of halo mass.

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. [68] and Moore et al. [69] 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 [70][71][72]).Even though it triggered many studies attempting to reduce small-scale structures by imposing other non-CDM candidates such as warm dark matter [73] and self-interacting dark matter [74], it is also possible to solve it with standard baryonic physics including early reionization [75] and self-regulation of star formation in low-mass halos [76][77][78][79].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 [29,37,51,80]), indicating that it obeys a power law where the slope −α 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 [80,81].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 [65,[80][81][82][83][84][85][86][87]).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 ∼10 5 M [37,51], 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 [46,50,85,[88][89][90]), which can shed light on the resolution issue.Hiroshima et al. [50] 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 [44,45,48].Assuming that ∼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 −α = −2 and −1.9 [44].More extensive simulations are needed to obtain the subhalo mass function in wide mass ranges and also to compare with analytic models [50].

Density Profile of Dark Matter Halos
By the end of the 1980's, it was already known in both analytic [91] and numerical [92] studies that the density profiles were described by power-law functions.Reference Dubinski and Carlberg [93] 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 [94].
Navarro et al. [57,58] 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 ρ(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 [95,96]) performed simulations with better mass resolutions.Whereas previous studies [57,58] used only ∼10,000 particles, they used ∼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 [97][98][99][100][101][102][103][104].In most cases, the slope of density became shallower as the radius went inward.A different approach was adopted by Jing and Suto [105], who used the triaxial model for describing the central structures.Moore et al. [96] and Diemand et al. [106] considered a more general profile, If β = 3, γ = 1, and η = 1, the profile is the same as NFW.More recent studies [37,51,52] archived one of the highest resolution dark matter only simulations for galaxy-size halos with mass resolution better than 10 4 M .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 ∼0.001r 200 .Besides, Springel et al. [37] and Stadel et al. [52] fitted the density using the Einasto profile [107] where α E is a free parameter.Note that r s and ρ s are not the same parameters as those in Equation ( 3).
Although we can obtain the density profile down to the radius ∼0.001r 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 [44,45,64,65,[108][109][110][111][112][113][114][115][116][117].
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 [108,118].
Traditionally, the concentration-mass relation has been calibrated with cosmological N-body simulations for relatively massive halos (10 10 M M 200 10 15 M ).Because the mass dependence of the concentration is weak for these mass halos, it is found that a single power-law function, c ∝ M −α c 200 , with slope α c in the range of 0.08 to 0.13, gives reasonable fits [108,[110][111][112]114].However, the dependence gradually becomes weaker toward less massive halos, and a clear flattening emerges [44,45,64,65,117,119,120], 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., [37,48,121]) where r max is the radius at which the circular velocity reaches its maximum value V max .This definition is also used to estimate the annihilation boost factor (e.g., [48]).
Even with the highest resolution simulations for galaxy-scale halos, the smallest resolved subhalo mass is around ∼10 5 M [37,51].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 [36], and ∼1000 for cluster-scale halos [40,41,56].
Another approach is adopting analytic models or fitting functions that can reproduce flattening of the concentration-mass relation (e.g., [44,64,115,116]).
In contrast to using the power-law extrapolation, the resulting boost factor is rather modest, three to a few tens [44,45,47] for Milky-Way halo, and less than ∼100 for cluster-scale halos [44,45,56].
The density profile at fixed halo mass shows a significant halo-to-halo scatter [122], possibly making a big impact on the annihilation signal.Inferring from the cosmological Millennium simulation [123], the effect of this non-universality on the annihilation flux is a factor of ∼3 [122], 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.

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., [37,48,97,108,121]).For example, Bullock et al. [108] 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., [124]).Diemand et al. [121] 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. [125] derived the concentration-mass relation of subhalos in Milky-Way-size halos by analyzing high resolution cosmological simulations [37,51] 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. [48] 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 [128] 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 [129,130]).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., [50]).

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 [11,12,14].Such halos are sometimes referred to as "microhalos." The density profiles of the microhalos have been investigated using cosmological N-body simulations [45, 66,[131][132][133].Diemand et al. [131] 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, ρ(r) ∝ r −γ , with a slope γ 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. [132] 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 ρ(r) ∝ r −1.5 , which is supported by follow-up cosmological simulations [45, 66,133] and cold-collapse simulations [134].Ishiyama [45] 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 [133,135]).Similar density structures are observed in recent simulations of ultracompact minihalos [136][137][138] and warm dark matter [139].The self-similar gravitational collapse models (e.g., [91,[140][141][142][143]) 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. [132] 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 ∼0.2 deg yr −1 , which are, however, stringently constrained with the diffuse gamma-ray backgrond [33].Gravitational perturbations to the millisecond pulsars might be detectable with future observations by pulsar timing arrays [132,[144][145][146]. Anderhalden and Diemand [66] have assumed a transition from the NFW to steeper cusps at scales corresponding to ∼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. [108].
Ishiyama [45] 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 [44] (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., [44,45,64,65,117,119]), 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.Ishiyama [45].Two subhalo mass functions, dn/dm = A/M vir (m/M vir ) −ξ , are used (A = 0.012, ξ = 2.0, and A = 0.030, ξ = 1.9 [44,48]).Thick dotted curves are for the NFW profile, where the empirical concentration-mass relation of field halos [44] are assumed for the full hierarchy of subhalos.Including the effect of steeper cusp of halos near the free streaming scale gives thick dashed curves.Besides, thick solid curves are results of incorporating the concentration of these halos derived from cosmological simulations [45].For comparison, boost factors obtained in other studies are shown with thin dashed curves [48] (two subhalo mass functions are used), thin solid curves [41], and crosses [36,51].

Models Based on Structure Formation and Tidal Evolution
In an analytical approach in Section 2, the subhalo luminosity L 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 [46] and Hiroshima et al. [50] 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 ca (m a , z a ) [64] and σ log c = 0.13 [65].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 ρ 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, as follows [90]: where , m(z) and M(z) are the subhalo and host-halo masses at z, respectively.Following Jiang and van den Bosch [90], Hiroshima et al. [50] adopted simple Monte Carlo simulations to estimate ṁ 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 [90] found A = 0.81 and ζ = 0.04, Hiroshima et al. [50] extended the mass and redshift ranges of applicability and found that these parameters are weakly dependent on both M and z: 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. [50] adopted a fitting formula given by Correa et al. [147].
The subhalo density profile after accretion is also well described with the NFW profile with a sharp truncation at r t : This is indeed a good approximation found in the simulations [37].In addition to r t , Peñarrubia et al. [129] found that the internal structure changes.If the inner profile is ∝ r −1 just like NFW, the maximum circular velocty V max and its corresponding radius r max evolve as respectively.After computing V max and r max at z, one can convert them to ρ s and r s through which are valid for the NFW profile.Finally by solving the condition the truncation radius r t is obtained.Hiroshima et al. [50] omitted subhalos with r t < 0.77r s from the subsequent calculations assuming that they were tidally disrupted [148].This criterion, however, might be a numerical artifact [130].Either case, Hiroshima et al. [50] 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 , ρ 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 sh /(dm a dz a ), is given by the excursion set or the extended Press-Schechter formalism [55].Especially Yang et al. [89] 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. [50] adopted their model III.
The subhalo mass function is obtained as 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. [50] adopted the log-normal distribution with the mean ca (m a , z a ) [64] and the standard deviation σ log c a = 0.13 [65].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 space halo finder [149].The bound mass is used as the subhalo mass, which nearly corresponds to the tidal mass [48].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 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 sh is larger.Again, a good agreement in f sh is found between the analytic model and the simulation results by Giocoli et al. [86].
The annihilation boost factor is then 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 sh (m a , z a , c a |M, z), is proportional to the volume integral of density squared ρ 2 sh (r) out to the truncation radius, where ρ s,sh , r s,sh and r t,sh are functions of (m a , z a , c a ) as well as M and z.
Then, the effect of sub n -subhalos (for n ≥ 1) can be estimated iteratively.At nth iteraction, when a subhalo accreted onto its host at z a with m a , it is assigned a sub-subhalo boost factor B (n−1) sh (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. [50] assumed that the sub-subhalos were distributed within the subhalo following n ssh (r) ∝ [1 + (r/r s ) 2 ] −3/2 .The luminosity due to sub-subhalos within a radius r is therefore proportional to their enclosed number and it gets suppressed by a factor of N ssh (< r t |r s )/N ssh (< r vir |r s,a ) due to the tidal stripping. 4 The luminosity due to the smooth component also decreases as L Thus the sub-subhalo boost after the nth iteration, B ssh is obtained by 4 We note that in estimating the effect of sub n -subhalos in the boost factors, Reference Hiroshima et al. [50] ignored the changes of ρ s and r s and hence did not include the factor of r 3 s in Equation (39) and ρ s r 3 s in Equation (40).In addition, in Equation (43), they multiplied L (0) sh by a factor of 1 + B We correct for all these effects in this review.
Similarly, the sub-subhalo mass fraction f ssh is obtained by f ssh (m a , z a , c a |M, z) = f sh (m a , z a ) N ssh (< r t |r s )/N ssh (< r vir |r s,a ) m sm (< r t |ρ s , r s )/m sm (< r vir |ρ s,a , r s,a ) , (42) where f sh (m a , z a ) is obtained with Equations ( 36) and m sm (< r|ρ s , r s ) ∝ ρ s r 3 s f (r/r s ) is the enclosed mass within r of the smooth component of the subhalo.The subhalo boost factor after nth iteration is obtained with Equation (37) by replacing sh [see discussions below Equation (15)]: dm a dz a d 2 N sh dm a dz a dc a P(c a |m a , z a ) The host luminosity in the absence of the subhalos L host,0 (M, z) is defined by marginalizing over the concentration parameter c vir : with the log-normal distribution P(c vir |M, z).
Figure 5 shows the subhalo boost factors B 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 ∼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 sub n -subhalos, which is saturated after the second iteration.The contribution to the boost factors due to sub-subhalos and beyond is 10% for the hosts with M host ≥ 10 13 M .The bottom left panel of Figure 5 shows the luminosity ratio L total /L host,0 = 1 − f 2 sh + B 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 [41,44,48].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 × 10 12 M halo by Hiroshima et al. [50] and the result of Springel et al. [37] 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.  [50].The effect of sub n -subhalos, up to n = 3, is shown in the right panel in the case of z = 0. Note that the three curves except for n = 0 overlap with each other.The bottom left panel shows the ratio between the total luminosity including the subhalo boost and the luminosity in absence of subhalos, L total /L host,0 = 1 − f 2 sh + B sh .The bottom right panel shows comparison of B sh between several models at z = 0: G12 [41], SC14 [44] and M17 [48] are based on N-body calculations while H18 [50] is on analytic calculations.The subhalo mass function for the N-body results is assumed to be dN sh /dm ∝ m −α .

Models for Self-Similar Subhalos
Assuming a self-similarity of the subhalos, Kamionkowski and Koushiappas [35] developed a fully analytic formulation for the probability distribution function of the dark matter density, P(ρ).Then Kamionkowski et al. [39] 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: where f sm (r) is the volume fraction occupied by the smooth component and ρ max is the highest dark matter subhalo density.Through the calibration with the numerical simulations, Kamionkowski et al. [39] found δ f = 0.2, α K = 0 and that the subhalo fraction was given by where κ = 0.007.Fornasa et al. [42] then suggested a larger value of κ = 0.15-0.2 to obtain a larger boost consistent with earlier work [36,40,41].The maximum subhalo density ρ max is estimated as , where c and z f are the concentration parameter and collapse redshift of the smallest halos.Kamionkowski et al. [39] adopted c = 3.5 and z f = 40 and ρ max = 80 GeV cm −3 .On the other hand, Ng et al. [43] obtained a smaller ρ max ≈ 20 GeV cm −3 even with a very small cutoff masses of M min = 10 −12 M Within the virial radius of the Milky-Way halo, the subhalo boost factor for dark matter annihilation is found to be no greater than ∼10 [39].

Universal Clustering of Dark Matter in Phase Space
Zavala and Afshordi [49] 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 [37], 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 [151] down to very small free-streaming cutoff scales.Based on the power spectrum, they obtained the subhalo boost factor greater than ∼30-100 for the Milky-Way size halos, which is significantly larger than the values obtained with other analytic work [39,46,50].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.

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 [37,150] 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 ∼100 (∼1000) for galaxy (cluster) size halos [36,41].
As a complementary approach, analytic models have been investigated.They are based on self similar propertiese of the subhalos [35,39], universal phase-space distribution [49] and extended Press-Schechter formalism combined with tidal stripping modeling [46,50].(More recent numerical approach also adopts the concentration-mass relation calibrated for the subhalos in order to take the tidal effects into account [48].)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. [50] predict the subhalo mass functions for various host masses and redshifts, which are found to be in good agreement with the simulation results [37,86,126,150].The annihilation boost factors based on these analytic models tend to be more modest, O(1) for galaxy-size halos and 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 .Simulations of microhalos with 10 −6 M suggest cuspier profiles towards halo centers such as r −1. 5 [45] 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 [152][153][154][155]) 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.

Figure 1 .
Figure 1.Dark matter distribution of a Milky-Way-size halo taken from a high-resolution cosmological N-body simulation [67].

Figure 3 .
Figure 3. Boost factor as a function of halo virial mass.The data used in six thick curves are taken from Ishiyama [45].Two subhalo mass functions, dn/dm = A/M vir (m/M vir ) −ξ , are used (A = 0.012,

Figure 4 .
Figure 4. Subhalo mass function for galaxy (M 200 = 1.8 × 10 12 M ) and cluster (M 200 = 5.9 × 10 14 M ) halos at z = 0 (top left), halos with 2.3 × 10 12 M at z = 2 and 4.7 × 10 11 M at z = 4 both of which would evolve to M 200 = 10 13 M at z = 0 (top right) and smaller halos of M 200 = 10 6 M and 10 7 M at z = 5 (bottom left).Results of the analytic models by Hiroshima et al. [50] are compared with those from the numerical simulations of similar halos and other fitting functions: Springel et al. [37], Diemand et al. [150] (top left), ν 2 GC H2 [126,127], Giocoli et al. [86] (top right) and Phi-2 (Ishiyama et al., in preparation; bottom left).The bottom right panel shows the subhalo mass fraction f sh as a function of the host mass M 200 for various values of redshift z.The thin solid curve is for z = 0 but with lower mass threshold of 1.73 × 10 10 h −1 M to be compared with Giocoli et al.[86] results shown as the squares.

Figure 5 .
Figure 5.The subhalo boost factor B sh as a function of the host mass M 200 for various values of redshift z (top left) based on the analytic models by Hiroshima et al.[50].The effect of sub n -subhalos, up to n = 3, is shown in the right panel in the case of z = 0. Note that the three curves except for n = 0 overlap with each other.The bottom left panel shows the ratio between the total luminosity including the subhalo boost and the luminosity in absence of subhalos, L total /L host,0 = 1 − f 2 sh + B sh .The bottom right panel shows comparison of B sh between several models at z = 0: G12[41], SC14[44] and M17[48] are based on N-body calculations while H18[50] is on analytic calculations.The subhalo mass function for the N-body results is assumed to be dN sh /dm ∝ m −α .

Figure A1 .
Figure A1.mass function m 2 dN/dm at z = 0.Each line corresponds to a different host halo mass.The fitting formula is applicable for the host mass from M host 10 −4 to 10 14 M and redshifts up to ∼6.