Clustering in the Phase Space of Dark Matter Haloes. II. Stable Clustering and Dark Matter Annihilation
Abstract
We present a model for the structure of the particle phase space average density () in galactic haloes, introduced recently as a novel measure of the clustering of dark matter. Our model is based on the stable clustering hypothesis in phase space, the spherical collapse model, and tidal disruption of substructures, which is calibrated against the Aquarius simulations. Using this model, we can predict the behaviour of in the numerically unresolved regime, down to the decoupling mass limit of generic WIMP models. This prediction can be used to estimate signals sensitive to the small scale structure of dark matter. For example, the dark matter annihilation rate can be estimated for arbitrary velocitydependent cross sections in a convenient way using a limit of to zero separation in physical space. We illustrate our method by computing the global and local subhalo annihilation boost to that of the smooth dark matter distribution in a MilkyWaysize halo. Two cases are considered, one where the cross section is velocity independent and one that approximates Sommerfeldenhanced models. We find that the global boost is , which is at the low end of current estimates (weakening expectations of large extragalactic signals), while the boost at the solar radius is below the percent level. We make our code to compute publicly available, which can be used to estimate various observables that probe the nanostructure of dark matter haloes.
keywords:
cosmology: dark matter methods: analytical, numerical1 Introduction
Despite being the most dominant type of matter in the Universe, dark matter is evident so far only through its gravitational effects on luminous matter. The nature of dark matter will thus continue to be elusive unless it has detectable nongravitational interactions. Weakly Interacting Massive Particles (WIMPs) are among the favourite dark matter candidates and are expected to give promising signals that can be detected experimentally either directly through scattering off nuclei in laboratory detectors, or indirectly through the byproducts of their selfannihilation into standard model particles (e.g. photons and electron/positrons pairs). Several experiments are pursuing such a discovery and their sensitivity is reaching the level of interaction predicted by popular WIMP models (e.g. Abazajian et al., 2012; Aprile et al., 2012).
Interestingly, there are tantalizing anomalies that might be caused by new dark matter physics. For instance, although the excess of GeV cosmic ray positrons established firmly by the PAMELA and Fermi Satellites (Adriani et al., 2009; Abdo et al., 2009), and recently confirmed with high precision by the AMS collaboration (Aguilar et al., 2013), can be explained by “ordinary” astrophysical sources (e.g. nearby pulsars, Linden & Profumo 2013), it can also be interpreted as dark matter annihilation. Large annihilation rates, of above the expected rate of a TeV mass thermal relic, and primarily leptonic final states are needed, however, to explain the data (Bergström et al., 2009). These requirements can be satisfied by leptophilic dark matter models coupled to light force carriers that enhance the annihilation cross section via a Sommerfeld mechanism (e.g. ArkaniHamed et al., 2009; Pospelov & Ritz, 2009). An anomalous extended gammaray emission at intermediate galactic latitudes (peaking at GeV) has also been interpreted as a signal of dark matter annihilation (e.g. Hooper & Slatyer, 2013).
Predictions for the hypothetical nongravitational signatures of dark matter are highly dependent on the clustering of dark matter at small scales. Although the steady progress of numerical body simulations over the past few decades has given us a detailed picture of the spatial dark matter distribution from large ( Gpc) to subgalactic scales ( pc), the regime most relevant for certain dark matter detection efforts (those based on extragalactic signals) remains below the resolution of current simulations (M) at ; Springel et al. 2008; Diemand et al. 2008; Stadel et al. 2009). This is because the current cold dark matter (CDM) paradigm of structure formation predicts a hierarchical scenario with a very large mass range of gravitationally bound dark matter structures, from M clustersize haloes down to a decoupling mass limit of M (e.g. Bringmann, 2009). The characteristic sizes of these small haloes (commonly called microhaloes) varies from times the size of a Milky Way (MW) type halo (we call this the nanostructure of dark matter haloes).
Despite their limited resolution, current simulations clearly suggest that the contribution of lowmass (sub)haloes is dominant for dark matter annihilation signals in the case of extragalactic sources, such as gammarays from galaxy clusters (e.g. Gao et al., 2012) and from integrated backgrounds (e.g. Zavala et al., 2010; Fornasa et al., 2013). Within the solar radius, it seems that the dominant signal should come from the diffuse distribution of dark matter rather than by small scale subclumps (Springel et al., 2008). Substructures within the MW halo can be detected in gammarays either individually as the hosts of satellite galaxies (Ackermann et al., 2011), dark satellites devoid of stars (Ackermann et al., 2012) or as a dominant contribution to the angular power spectrum of a fullsky gammaray signal (e.g. SiegalGaskins, 2008).
Extrapolations below resolved scales are therefore needed to obtain a prediction of the expected signals from dark matter annihilation. For a given host halo, these extrapolations ultimately depend on the survivability of the smallest substructures as they are tidally disrupted by the host. Low mass subhaloes form earlier in the hierarchical scenario and, being the densest, are expected to survive tidal disruption contributing heavily to the annihilation signal. Current estimates on this contribution mostly rely on assumptions about the abundance, spatial distribution, and internal properties of unresolved subhaloes that could lead to significant over/underestimations. Typical estimates calibrated to simulation data vary by up to an order of magnitude (e.g. Springel et al., 2008; Kuhlen et al., 2008; Kamionkowski et al., 2010). The role of substructure, and therefore the uncertainty in the extrapolations to the unresolved regime, is magnified in Sommerfeldenhanced models where the annihilation cross section scales as the inverse velocity making the cold and dense subhaloes dominant not only in extragalactic structures (Zavala et al., 2011) but plausibly also locally (Slatyer et al., 2012).
In this paper, we present a novel method to study dark matter clustering that can be used to estimate the dark matter annihilation rate calibrated with simulation data. This method is based on a physicallymotivated model inspired by an extension into phase space of the stable clustering hypothesis (Afshordi et al., 2010), originally introduced in position space by Davis & Peebles (1977). The novelty of the method relies on using a more complete picture of dark matter clustering based on a new quantity, the particle phase space average density (), a coarsegrained phase space density that can be used straightforwardly to compute the annihilation rate. was introduced in a companion paper (Zavala & Afshordi, 2013, hereafter Paper I) where we studied its main features in galacticsize haloes using the simulation suite from the Aquarius project (Springel et al., 2008). Remarkably, we found it to be nearly universal at small scales (i.e. small separations in phase space) across time and in regions of substantially different ambient densities. This near universality is also preserved across the different Aquarius haloes, having similar masses but diverse mass accretion histories: they possess only slightly different s at small scales without any rescaling. We argue in this paper that this behaviour can be roughly described by a refinement of the stable clustering hypothesis through a simple model that incorporates tidal disruption of substructures.
In cases other than swave selfannihilation, the dependence of the interaction on the relative velocity of the annihilating particles has to be averaged over the velocity distribution of the dark matter particles. Although it is common to assume a Maxwellian velocity distribution, current simulations have shown that there are significant deviations from this assumption related to the dark matter assembly history (e.g. Vogelsberger et al., 2009). Our method is particularly useful for models where the cross section is velocity dependent (such as in Sommerfeldenhanced models) because the particle phase space average density () that we introduce does not rely on assumptions about the velocity distribution and can deal with these cases in a natural way.
This paper is organised as follows: In Section 2, we summarize the main results of Paper I, introducing the clustering of dark matter through and its nearly universal structure at small scales in galactic haloes. In Section 3, we describe how the annihilation rate is directly connected to the limit of zero physical separation of . In Section 4, we describe our model of the small scale structure of based on the spherical collapse model and the stable clustering hypothesis in phase space, refined by a subhalo tidal disruption prescription. In Section 5, we fit our model to the simulation data while in Section 6 we illustrate how this can be used to compute the global and local subhalo boost to the smooth annihilation rate in a MWsize halo. Finally, a summary and our main conclusions are given in Section 7.
2 The Particle Phase Space Average Density () on small scales
We follow the same definitions as in Paper I and study the clustering in phase space through , defined as the massweighted average (over a volume in phase space) of the coarsegrained phase space density of dark matter, on spheres of radius and , in position and velocity spaces, respectively:
(1) 
where is the phase space distribution function at the phase space coordinates and .
In Paper I of this series we implemented and tested an estimator of in body simulations based on pair counts and applied to the set of Aquarius haloes (Springel et al., 2008). In this paper, we refer exclusively to the results found for halo AqA2 that has a virial mass and radius, defined with a mean overdensity of 200 times the critical value, of M and kpc, respectively.
One of the main results that we obtained is that the structure of averaged within the virialized halo is clearly separated by two regimes: (i) a region at large scales (i.e. large separations in phase space) dominated by the smooth dark matter distribution where varies in time due to the insideout growth of the dark matter halo, and (ii) a region at small scales (i.e. small separations in phase space) dominated by gravitationally bound substructures where is nearly universal across time and ambient density. In the reminder of this paper we consider only the small scale regime, which is the one of relevance in estimating the impact of substructure on dark matter annihilation signals.
Although the small scale regime can be roughly described by a subhalo model with abundance and properties given by scaling laws fitted to the simulation data (see Fig. 5 of Paper I), we found that a better description is given by a family of superellipse contours of constant , with parameters that can be adjusted to fit the variations of in redshift and halocentric distance (see Tables 2 and 3 of Paper I):
(2) 
where is a shape parameter and and are the axes of the superellipse:
(3)  
Fig. 1 shows contours of log() in the small scale regime averaged within the virialized region of AqA2 at (solid lines). The dotted lines are the fitting contours through Eqs. (23) with the parameters as given in Table 1. In section 4, we present the model that motivates this fitting formula. We note that the shaded area on the left corner encompasses the region where our results could be affected by resolution conservatively by at the most (see Appendix of Paper I); we exclude this region from the fitting procedures and from subsequent analyses.
3 Dark matter annihilation rate: substructure boost
The number of dark matter selfannihilation events per unit time is given in terms of the phase space distribution at the phase space coordinates and :
(4)  
where is the total dark matter mass within the phase space volume , is the product of the annihilation cross section and the relative velocity between pairs^{1}^{1}1Note that in our notation, the relative velocity is , thus, . Throughout our paper we choose to use the former since this is the common choice in the literature., and we have used Eq.(1) to introduce .
The annihilation rate in a region of spatial volume is typically written as:
(5) 
where is the dark matter particle mass, is the local dark matter density (that includes contribution from the smooth dark matter distribution and from substructure), and is an average of over the velocity distribution of the dark matter particles (typically assumed to be Maxwellian). Instead of separating the dark matter spatial and velocity distributions, Eq. (4) gives directly, without any assumptions, the annihilation rate as an integral over the relative velocity at the limit of zero spatial separation of , henceforth abbreviated (). Since in general is an arbitrary function of , we can use Eq. (4) to easily accommodate any velocity dependent annihilation cross section. Notice also that Eq. (4) adapts to the region of interest for the annihilation rate by simply using averaged over that region.
The near universality of at small scales (i.e. small separations in phase space) makes the functional shape in Eqs. (23) valid across regions of substantially different ambient density with only slight changes to the fitting parameters (see sections 3.3 and 3.4 of Paper I) and can therefore be straightforwardly applied in Eq. (4) to calculate the annihilation rate in gravitationally bound substructure either globally in the whole halo, or locally in a certain region of the halo. Since we also found that at small scales is nearly insensitive to the assembly history of a particular halo (see Fig. 10 of Paper I), the results we find for the particular initial conditions of AqA2 are also a very good approximation for different initial conditions.
The integral in Eq. (4) can also be easily transformed into an integral over using Eq. (2) and a change of variables:
(6) 
Notice that because of the limit , the annihilation rate is at the end only sensitive to and . The limits of the integral over () correspond to the maximum and minimum values of the separation in velocity where substructures contribute. Note that we can only use Eqs. (23) down to:
(7) 
below which falls off more rapidly than the power law in Eq. (3), and more importantly, the smooth distribution of dark matter dominates in general (Paper I). We will then take this value as the transition between the smooth and subhaloes dominance of the annihilation rate.
To compute the contribution from the smooth component we substitute in Eq. (5) by , the spherically averaged Einasto density profile:
(8) 
where and are the density and radius at the point where the logarithmic density slope is 2, and is the Einasto shape parameter. We take the parameter values from the fit to the AqA2 halo given in Navarro et al. (2010): , , . We also assume that the velocity distribution of dark matter particles in the smooth component is Maxwellian. The substructure boost is then simply defined as:
(9) 
As an example, let us take the case where and compute the total annihilation rate in resolved substructures at within a MWsize halo. Since (see Table 1), according to Eq. (6) we have:
(10) 
If we take the value of given in Eq. (7) and the maximum value of we can resolve in AqA2, , then we can estimate the “resolved” substructure boost to the total annihilation rate in the AqA2 halo:
(11) 
The simulation particle mass of AqA2 is M and the minimum subhalo mass that we can trust in terms of global subhalo properties is M ( particles; see Figs. 26 and 27 of Springel et al. 2008). By using a subhalo model (see Section 3.2 of Paper I), we have checked that this minimum “resolved” subhalo mass also corresponds roughly to the value of we can resolve. As a consistency check, we can therefore compare the value in Eq. (11) to the subhalo boost computed in Springel et al. (2008) (using the same simulation data) by summing up the contribution of all resolved subhaloes above M within the virial radius of the AqA1 halo, and assuming a NFW density profile for each of them (see second green line from top to bottom in Fig. 3 of Springel et al. 2008). The resolved subhalo boost they found is , quite similar to ours.
Before continuing with the analysis of the applications of to compute the dark matter annihilation rate, we describe in the next section a physicallymotivated model that we will use to extrapolate the behaviour of to the unresolved regime.
4 Improved stable clustering hypothesis and spherical collapse model
The stable clustering hypothesis was originally introduced by Davis & Peebles (1977) to study the galaxy twopoint correlation function in the strongly nonlinear regime. The hypothesis proposes that the number of neighbouring dark matter particles within a fixed physical separation becomes constant (i.e. there is no net streaming motion between particles in physical coordinates) on sufficiently small scales. Jing (2001) found that the hypothesis is valid when averaging the mean pair velocity between particles in many simulated haloes but found that this is not generally satisfied within one single virialized halo.
The hypothesis can be extended to phase space (Afshordi et al., 2010) by stating that, on average, the number of particles within the physical distance and physical velocity of a given particle does not change with time for sufficiently small and . In Paper I we found evidence for the validity of the stable clustering hypothesis in phase space through the analysis of at small scales finding that it typically varies by a factor of a few in regions of substantially different ambient densities (by nearly four orders of magnitude).
The small scale structure of today within the virialized region of the halo is given by a collection of gravitationally bound merging substructures that collapsed earlier than the host halo^{2}^{2}2From here onwards, we will follow closely the notation given in Afshordi et al. (2010).. Each of these substructures has a characteristic phase space density imprinted at the collapse time, when the Hubble’s constant has a value , and each has a collapse mass . In the absence of tidal disruption, would be preserved from the time of collapse until today. The stable clustering hypothesis can then be used and applied to write a simplified solution to the collisionless Boltzmann equation at the phase space coordinates of the following form (for a spherically symmetric gravitational potential, Afshordi et al. 2010):
(12)  
where is the characteristic density of the collapsed subhalo, which is roughly times the critical density at the collapse time. Since we know that subhaloes will be subjected to tidal disruption as they merge and move through denser environments, Eq. (12) accounts for tidal stripping modifying the stable clustering hypothesis prediction (i.e. ) by introducing , which can be interpreted as the mean fraction of particles that remain bound. To make the connection with the simulation results, we equate Eq. (1) (in the smallscale regime) with Eq. (12). We can then associate a constant value of within the simulated MWsize halo to a subhalo that formed with a typical mass in the past.
The characteristic phase space density of a given (sub)halo at formation time can be estimated within the spherical collapse model using the characteristic densities and velocities of the collapsed object: and (e.g. Afshordi & Cen, 2002):
(13) 
In the spherical collapse model, the subhalo collapses when the r.m.s tophat linear overdensity (mass variance) crosses the linear overdensity threshold at an epoch given by:
(14) 
The mass variance is defined by:
(15) 
where is the tophat window function and is the linear power spectrum. We use the fitting function given in Taruya & Suto (2000) which is accurate to a few percent in the mass range we use here:
(16) 
where and M, where with and being the contributions from matter and baryons to the mass energy density of the Universe. The mass variance is normalized to the value at Mpc spheres at redshift zero. We assume the same cosmological parameters as in the Aquarius simulations (those of a WMAP1 flat cosmology): , , , and (the spectral index of the primordial power spectrum).
Using Eq. (14) we can then calculate the epoch of collapse of a given halo of mass and estimate its phase space density using Eq. (13). Since the area enclosed by the ellipse: is simply then we can equate the phase space volume encompassed by the ellipsoid to the volume of the collapsed halo :
(17) 
Thus, we can finally give the prediction of the improved stable clustering hypothesis for the curves of contours of constant :
(18) 
where and .
Comparing this prediction with the simulation data, we find that ellipses are not a good description, instead generalizing the ellipses to superellipses (Lamé curves) provides a good fit to the simulated MW halo at small (i.e., to ):
(19) 
Also, we find that a better fit to the simulated data is obtained from a mass dependent tidal disruption parameter instead of the constant value taken by Afshordi et al. (2010). In what follows, we introduce a model for tidal disruption that captures this mass dependence.
4.1 Tidal stripping
We introduce a simplified model of tidal stripping in which we assume that the characteristic density in a given subhalo changes with time according to:
(20) 
where is the characteristic freefall time of the host as the subhalo starts being stripped. If we assume that , then the solution to Eq. (20) as a function of redshift is:
(21)  
where
(22) 
with being the relevant redshift of first infall and we have assumed that and . If we take as an ansatz that evolves in a similar way as , then we can write:
(23) 
The value of at the time of first infall is uncertain, but considering that structures that formed earlier would be more resilient to tidal stripping, and that tidal disruption begins as the subhalo infalls into a larger structure (not necessary the final host halo) with a mass (), we model the initial condition as:
(24) 
where and are free parameters.
The infall redshift can be estimated using the Extended PressSchechter formalism. We are interested in the probability that a halo of mass , had progenitors of mass :
(25) 
where (at early times) is the overdensity barrier required for spherical collapse. Using the method of steepest descent, we can approximate the average infall time:
(26) 
Since , where is the linear growth factor (for an approximation formula see e.g. Carroll et al., 1992), we can then obtain the collapse redshift from the mass variance and finally obtain the value of for a given .
We introduce a halocentric distance dependence in by considering that the average density of the host within a given distance was established at an epoch when its characteristic density had that value: (this defines the characteristic free fall time). We also use the radial diagonal part of the Hessian of the potential as the quantity that drives tidal disruption, rather than simply .
After , the density does not dilute anymore as the Universe expands. If we then have an additional solution to Eq. (20) for :
(27)  
where
(28) 
and the subscript means that the quantity is renormalized to the radius where the average solution applies, i.e., to . If , then only this “second mode” (Eq.27) of stripping occurs.
In the end, we have a model of tidal disruption with free parameters: , , , and , in addition to the improved stable clustering prediction which includes additional free parameters (which in principle could depend on redshift and in halocentric distance): , and accounting for the stretching in phase space due to tidal shocking (Eq.19). This model is able to describe the smallscale behaviour of at different redshifts and at different regions within the halo as we show in the following section.
5 Simulation data and model comparison
5.1 Average behaviour within the virialized halo and redshift dependence
We first compare the model developed in Section 4 with the value of averaged over all particles that are gravitationally bound to the AqA2 halo, i.e., those in the smooth and substructure components (for a more detailed definition see the first paragraphs of section 3 of Paper I). Fig. 2 shows the small scale behaviour of in the simulation at different redshifts (solid lines) and fits by our model in dotted lines. Although the fits are poorer for (particularly at ), it is clearly a reasonable description at smaller scales (i.e., separations in phase space), which are the ones that matter the most to extrapolate the behaviour to the unresolved scales.
The variations across redshift can be accommodated by introducing a slight redshift dependence of the parameters , and . The former is given by a simple linear relation as in Paper I:
(29) 
while and are given in Table 2.
Redshift  

0.0  
0.95  
2.2  
3.5 
The full tidal stripping model described in Section 4.1 fits the simulated data with the following values for the five free parameters: , , , and . Although we did not explore exhaustively the large parameter space we found that large deviations from these values seem to give a poorer fit (this is especially the case for and ). We also note that the behaviour of as a function of mass and redshift can be roughly approximated by the following formula (up to ):
(30) 
5.2 Halocentric distance dependence
Fig. 3 shows the contours of log() averaged over particles that are located in concentric shells at different distances from the halo centre as described in the caption. The solid lines are the simulation data for the AqA2 halo while the dashed lines are the result of our model with the same value of the free parameters associated to as in the previous average case (i.e., , , , and ). We note that is a function of radius through Eq. (27) that can be approximated by:
(31) 
where the radialdependent parameters , and , as well as the shape parameters , and , that fit the simulation data are given in Table 3.
0.00.2  

0.20.4  
0.40.6  
0.60.8  
0.81.0 
6 Boost to Dark Matter Annihilation due to unresolved substructure
Once we have calibrated the model to the simulation data, we can use it to extrapolate the behaviour of to scales that are unresolved. We argue that this extrapolation method offer advantages over other commonly used methods in the computation of the boost to the annihilation rate due to gravitationally bound unresolved substructures. This is the case especially for methods that extrapolate the concentrationmass relation. The support for this argument is twofold:

The small scale modelling of is physically motivated by: the stable clustering hypothesis in phase space, the spherical collapse model, and a simplified tidal disruption description.

The structure of is directly connected to the annihilation rate (Eq.4), with the small scale behaviour reflecting the substructure contribution.
By calibrating the model to , we avoid the use of the subhalo model and instead of uncertainties on the abundance of subhaloes, their radial distribution, and their internal properties, our model is sensitive to the behaviour of and to the assumption that all other shape parameters remain unchanged at unresolved scales. However, as we show below, of all these, is the only one of relevance.
As a first example, let us compute the total annihilation rate due to substructures within a MWsize halo at in the case of :
(32)  
where we have used Eqs.(12),(14),(17), and (19). Notice that only and enter into the annihilation rate (i.e., and are irrelevant). The limits of the integral are the minimum and maximum subhalo masses at collapse that contribute to , which are nontrivially connected to the tidally disrupted subhalo masses today. With the values of () mentioned in Section 3 corresponding to resolved subhaloes we can directly estimate M (M) since . These masses are a factor of several above the minimum (maximum) subhalo masses at contributing to in the AqA2 halo.
The projection of along the direction can be seen in Fig. 4 for the whole virialized region of AqA2 at . This is the quantity of interest to compute the annihilation rate and we show with a black solid line the case when kpc (which is the minimum physical separation we can resolve). The dashed red line is a fit by our full model with the same cut in while the limit to in the model is shown with a black dashed line. The extrapolation using the fitting function that directly relates with (see Eq. 2) is shown with a blue dashed line. We have extended the horizontal axis to the corresponding velocity separations of WIMP models with M. From Fig. 4 we can see that the simple formula (see Eq. 2):
(33) 
provides a good approximation to our full modelling and it immediately suggests that the annihilation rate scales roughly logarithmically with (for ).
We use Eq. (32) to compute the subhalo boost (using the smooth Einasto distribution described in Section 3) as a function of the minimum collapsed mass down to the decoupling masses corresponding to WIMPs; this is shown in Fig. 5. The filled circle is the value in Eq. (11) corresponding to resolved substructures and found directly with the small scale fitting function of in Section 3 (i.e., Eq. 10). The star symbol on the right (left) is shown for reference, and corresponds to the subhalo boost reported by Springel et al. (2008) for a minimum subhalo mass () M. The extrapolation in this case was done under the assumption that the radial subhalo luminosity profile preserves its shape and that the rescaling of the normalisation with follows the resolved trend. Recall that these are not the masses at collapse so to put the circle and star symbols on the right of Fig. 5 we use MM as explained two paragraphs above, while the location of the star symbol on the left of the figure is somewhat uncertain. The triangle symbol is the extrapolation to lower masses made by Kamionkowski et al. (2010) using the probability distribution function of the density field (subhaloes imprint a powerlaw tail in this PDF) calibrated with the MWsized simulation Via Lactea II (Diemand et al., 2008). The results from our method are quite close to those found by Kamionkowski et al. (2010) and are an order of magnitude lower than the estimates by Springel et al. (2008).
Zavala et al. (2010) also estimated a subhalo boost with a statistical analysis of all haloes in the MillenniumII simulation (BoylanKolchin et al., 2009), implicitly extrapolating the subhalo mass function and concentrationmass relation to the unresolved regime. They found a large range of possible subhalo boosts, within , depending on the exact parameters of the extrapolation, for a M halo.
6.1 Sommerfeldenhanced models
If the annihilation crosssection is enhanced by a Sommerfeld mechanism (e.g. ArkaniHamed et al., 2009; Pospelov & Ritz, 2009), then the annihilation rate increases with lower relative velocities until saturating due to the finite range of the interaction between the particles prior to annihilation. Instead of using a specific Sommerfeldenhanced model, we generically approximate the cross section as:
(34) 
The value of is commonly near (the socalled “” boost), but it can reach near resonances; we only consider the former case. For Eq. (6.1) and using our modelling we obtain:
(35)  
where is the collapse mass corresponding to the characteristic velocity at saturation ; using Eq. (19):
(36) 
Fig. 6 shows the subhalo boost, relative to the smooth dark matter distribution of the host halo, for Sommerfeldlike models where the cross section scales as in Eq. (6.1) with . To compute the Sommerfeldenhanced smooth component, we simply took the annihilation rate corresponding to the Einasto profile, defined in Section 3, and scaled it up based on the characteristic velocity dispersion of the AqA2 halo: , with . The solid line is for ; below this value, all substructures are saturated and thus the boost is just a trivial scaled version of that shown in Fig. 5. The dotted line roughly approximates one of the benchmark point models defined in Finkbeiner et al. (2011) (BM1 in their Table 1) to satisfy a number of astrophysical constraints (relic abundance, CMB power spectrum, etc.) while at the same time providing a fit to the cosmic ray excesses observed by the PAMELA and Fermi satellites. Enhancements much larger than this value (e.g. dashed line, ) are likely ruled out by astrophysical constraints but illustrate the transition between the unsaturated and saturated regimes. We note that BM1 is still compatible with the new measurement of the positron excess reported by the AMS collaboration (see Cholis & Hooper, 2013). Since allowed Sommerfeldenhanced models in the parameterisation we have used have , then the enhanced substructure boost is at the end just a scaled version of the constant cross section boost: .
6.2 Halocentric distance boost
Finally, we can also apply our methodology to compute the substructure boost as a function of the distance to the halo centre. To do so, all that is needed is to replace the global values of the parameters and for the radial dependent values that we fit to the simulation data in Section 5, and to account for the renormalisation in mass of the phasespace volume where is averaged. A simple but good approximation for is given in Eq. (31). Using it, our results can be easily reproduced with the values given in Table 3.
Fig. 7 shows the substructure boost to the annihilation rate of the smooth halo as a function of the distance to the halo centre. Since in Section 5.2 we calibrated our fitting parameters in concentric shells with a thickness of , we plot the results with bars showing the extent of the shells in which was averaged. We remark that in this case is the total mass in a given shell, not as it has been so far in this section. In Fig. 7 we show the cases with (circles) and Sommerfeldlike enhancement corresponding to the benchmark point BM1 (stars) presented in Finkbeiner et al. (2011) (see Section 6.1). Our results are in good agreement with those of Kamionkowski et al. (2010) (see their Fig. 4), although we seem to predict lower boosts in the central regions.
Values specific to a certain radius can be approximated by interpolating and and considering that for small volumes one has and thus the subhalo boost can be computed independently of the specific value of . For instance, for kpc, we find that the boost is only for the constant cross section case, and for the Sommerfeldlike model BM1.
7 Summary and Conclusions
The clustering of dark matter at scales unresolved by current numerical simulations is a key ingredient in many predictions of nongravitational (and some gravitational) signatures of dark matter. The characteristic scale of the smallest haloes (commonly called microhaloes) contributing to these signals is times the size of the Milky Way halo, and we therefore refer to it as the nanostructure of dark matter haloes. The degree of uncertainty of this nanostructure clustering can be as much as two orders of magnitude for WIMP dark matter models, since the minimum bound haloes have masses orders of magnitude below the highest resolution simulations to date. In the case of dark matter annihilation for example, the small, cold and dense dark structures have the dominant contribution to the hypothetical extragalactic signals, and thus the bulk of the predicted annihilation rate comes from extrapolating, in diverse ways, the dark matter clustering from the resolved to the unresolved regime.
Most extrapolation methods rely on assumptions about the abundance, spatial distribution, and internal properties of (sub)haloes. For example, the total subhalo boost to the annihilation rate over the smooth dark matter distribution of a single halo, is typically computed by extrapolating, either explicitly or implicitly, the subhalo mass function and the concentrationmass relation. Extending these functions as powerlaws down to lower masses yields the higher boosts. In this paper, we present an alternative method based on a novel perspective on the clustering of dark matter introduced in Zavala & Afshordi (2013) (Paper I). The method rests upon writing the dark matter annihilation rate in a given volume as an integral over velocities of a new coarsegrained Particle Phase Space Average Density (; see Eqs. 1 and 4).
The structure of was analysed in detail in Paper I, where it was found to be nearly universal in time and across different ambient densities, in the regime dominated by substructures. Here we present a model of the structure of inspired by the stable clustering hypothesis and the spherical collapse model, and improved by incorporating a prescription for the tidal disruption of subhaloes.
Our modelling provides a physicallymotivated explanation of the twodimensional functional shape of , and calibrated to the MWsize Aquarius simulations, gives a firm basis to extrapolate into the unresolved substructure regime. In summary, the main advantages of our method are twofold:

The free parameters in our model are fitted to a single 2D function that is a very sensitive measure of cold small scale substructure and is directly connected to the annihilation rate.

The annihilation rate is written as an integral over relative velocities of the “zeroseparation” limit of multiplied by the annihilation cross section times the relative velocity (see Eq. 4). This allows us to accommodate any velocitydependence coming from particle physics models in a straightforward way.
Although our model has several free functions that are fitted to the MWsize Aquarius simulations, only two are of relevance for the prediction of the annihilation rate: (i) in Eq. (19) parametrizes the ratio of phase space volume in spherical collapse, to the characteristic phase space volume of subhaloes at the time of formation; and (ii) in Eq. (12) that can be interpreted as the dilution of the characteristic phase space density at collapse due to tidal disruption. For the simulated data we analised, the former is and we present tabulated fitted values across redshifts and distances to the MWhalo centre (Tables 2 and 3), while for the latter we present simple parameterisations (see Tables 2 and 3, Eqs. 2930).
As a sample application of our model, we computed the subhalo boost (over the smooth dark matter distribution) in a MWsize halo, both globally (i.e. over the whole virialized halo) and locally as a function of halocentric distance, for cases where the annihilation cross section is constant and for a generic Sommerfeldlike enhanced case where up to a maximum saturation value. We find that in the former, the global boost is for typical WIMP models (with decoupling masses in the range M).
Our estimate of the subhalo boost is in the low end of current estimates being in agreement with Kamionkowski et al. (2010) (based on a different simulation), and a factor of lower than the boost computed in Springel et al. (2008) (based on the same simulation suite than the one used here). The discrepancy with the latter is likely caused by their implicit extrapolation of the subhalo mass function and, perhaps more importantly, the concentrationmass relation. Evidence of this can be seen in the structure of in the resolved regime (without recurring to the modelling of the unresolved regime) where our analysis indicates that the annihilation rate (in the case) from substructures scales logarithmically with the maximum (coarsegrained) phase space density (Eq. 10) rather than as a power law. Since we argue that is a more direct measure of the annihilation signal, our results seem to disfavour large boost factors.
has also other potential applications in dark matter direct detection, pulsar timing and transient weak lensing (see Rahvar et al., 2013) that will be explored in the future.
We have made our code to compute (with our full model) publicly available online at http://spaces.perimeterinstitute.ca/p2sad/. Interested users should refer to that site for instructions on how to use the code.
Acknowledgments
We thank the members of the Virgo consortium for access to the Aquarius simulation suite and our special thanks to Volker Springel for providing access to SUBFIND and GADGET3 whose routines on the computation of the 2PCF in real space were used as a basis for the 2PCF code in phase space we developed. JZ and NA are supported by the University of Waterloo and the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation. JZ acknowledges financial support by a CITA National Fellowship. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca) and Compute/Calcul Canada. The Dark Cosmology Centre is funded by the DNRF.
References
 Abazajian et al. (2012) Abazajian K. N., Blanchet S., Harding J. P., 2012, Phys. Rev. D, 85, 043509
 Abdo et al. (2009) Abdo A. A., et al., 2009, Physical Review Letters, 102, 181101
 Ackermann et al. (2011) Ackermann M., et al., 2011, Physical Review Letters, 107, 241302
 Ackermann et al. (2012) Ackermann M., et al., 2012, ApJ, 747, 121
 Adriani et al. (2009) Adriani O., et al., 2009, Nature, 458, 607
 Afshordi & Cen (2002) Afshordi N., Cen R., 2002, ApJ, 564, 669
 Afshordi et al. (2010) Afshordi N., Mohayaee R., Bertschinger E., 2010, Phys. Rev. D, 81, 101301
 Aguilar et al. (2013) Aguilar M., et al., 2013, Physical Review Letters, 110, 141102
 Aprile et al. (2012) Aprile E., et al., 2012, Physical Review Letters, 109, 181301
 ArkaniHamed et al. (2009) ArkaniHamed N., Finkbeiner D. P., Slatyer T. R., Weiner N., 2009, Phys. Rev. D, 79, 015014
 Bergström et al. (2009) Bergström L., Edsjö J., Zaharijas G., 2009, Physical Review Letters, 103, 031103
 BoylanKolchin et al. (2009) BoylanKolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150
 Bringmann (2009) Bringmann T., 2009, New Journal of Physics, 11, 105027
 Carroll et al. (1992) Carroll S. M., Press W. H., Turner E. L., 1992, ARA&A, 30, 499
 Cholis & Hooper (2013) Cholis I., Hooper D., 2013, Phys. Rev. D, 88, 023013
 Davis & Peebles (1977) Davis M., Peebles P. J. E., 1977, ApJS, 34, 425
 Diemand et al. (2008) Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nature, 454, 735
 Finkbeiner et al. (2011) Finkbeiner D. P., Goodenough L., Slatyer T. R., Vogelsberger M., Weiner N., 2011, Journal of Cosmology and Astroparticle Physics, 5, 2
 Fornasa et al. (2013) Fornasa M., Zavala J., SánchezConde M. A., SiegalGaskins J. M., Delahaye T., Prada F., Vogelsberger M., Zandanel F., Frenk C. S., 2013, MNRAS, 429, 1529
 Gao et al. (2012) Gao L., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 419, 1721
 Hooper & Slatyer (2013) Hooper D., Slatyer T. R., 2013, arXiv:1302.6589
 Jing (2001) Jing Y. P., 2001, ApJ, 550, L125
 Kamionkowski et al. (2010) Kamionkowski M., Koushiappas S. M., Kuhlen M., 2010, Phys. Rev. D, 81, 043532
 Kuhlen et al. (2008) Kuhlen M., Diemand J., Madau P., 2008, ApJ, 686, 262
 Linden & Profumo (2013) Linden T., Profumo S., 2013, ApJ, 772, 18
 Navarro et al. (2010) Navarro J. F., Ludlow A., Springel V., Wang J., Vogelsberger M., White S. D. M., Jenkins A., Frenk C. S., Helmi A., 2010, MNRAS, 402, 21
 Pospelov & Ritz (2009) Pospelov M., Ritz A., 2009, Physics Letters B, 671, 391
 Rahvar et al. (2013) Rahvar S., Baghram S., Afshordi N., 2013, arXiv:1310.5412
 SiegalGaskins (2008) SiegalGaskins J. M., 2008, JCAP, 10, 40
 Slatyer et al. (2012) Slatyer T. R., Toro N., Weiner N., 2012, Phys. Rev. D, 86, 083534
 Springel et al. (2008) Springel V., Wang J., Vogelsberger M., Ludlow A., Jenkins A., Helmi A., Navarro J. F., Frenk C. S., White S. D. M., 2008, MNRAS, 391, 1685
 Springel et al. (2008) Springel V., White S. D. M., Frenk C. S., Navarro J. F., Jenkins A., Vogelsberger M., Wang J., Ludlow A., Helmi A., 2008, Nature, 456, 73
 Stadel et al. (2009) Stadel J., Potter D., Moore B., Diemand J., Madau P., Zemp M., Kuhlen M., Quilis V., 2009, MNRAS, 398, L21
 Taruya & Suto (2000) Taruya A., Suto Y., 2000, ApJ, 542, 559
 Vogelsberger et al. (2009) Vogelsberger M., Helmi A., Springel V., White S. D. M., Wang J., Frenk C. S., Jenkins A., Ludlow A., Navarro J. F., 2009, MNRAS, 395, 797
 Zavala & Afshordi (2013) Zavala J., Afshordi N., 2013, arXiv:1308.1098
 Zavala et al. (2010) Zavala J., Springel V., BoylanKolchin M., 2010, MNRAS, 405, 593
 Zavala et al. (2011) Zavala J., Vogelsberger M., Slatyer T. R., Loeb A., Springel V., 2011, Phys. Rev. D, 83, 123513