Clustering in the Phase Space of Dark Matter Haloes. II. Stable Clustering and Dark Matter Annihilation

Clustering in the Phase Space of Dark Matter Haloes. II. Stable Clustering and Dark Matter Annihilation

Jesús Zavala and Niayesh Afshordi
Perimeter Institute for Theoretical Physics, 31 Caroline St. N., Waterloo, ON, N2L 2Y5, Canada
Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada
Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, Ontario M5S 3H8, Canada
Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, 2100 Copenhagen, Denmark

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 velocity-dependent 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 Milky-Way-size halo. Two cases are considered, one where the cross section is velocity independent and one that approximates Sommerfeld-enhanced 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.

cosmology: dark matter methods: analytical, numerical

1 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 non-gravitational 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 self-annihilation 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. Arkani-Hamed et al., 2009; Pospelov & Ritz, 2009). An anomalous extended gamma-ray 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 non-gravitational 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 sub-galactic 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 cluster-size 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 low-mass (sub)haloes is dominant for dark matter annihilation signals in the case of extragalactic sources, such as gamma-rays 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 gamma-rays 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 full-sky gamma-ray signal (e.g. Siegal-Gaskins, 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 Sommerfeld-enhanced 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 physically-motivated 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 coarse-grained 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 galactic-size 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 re-scaling. 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 s-wave self-annihilation, 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 Sommerfeld-enhanced 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 MW-size 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 mass-weighted average (over a volume in phase space) of the coarse-grained phase space density of dark matter, on spheres of radius and , in position and velocity spaces, respectively:


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 Aq-A-2 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 inside-out 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.

Figure 1: Contours of the logarithm of the particle phase space average density () of the Aq-A-2 halo (solid) at and for the fitting function defined in Eqs. (2-3) (dashed) with the parameters as given in Table 1. The dashed region in the left corner encompasses the region where resolution is a potential issue (see Appendix of Paper I).

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 halo-centric distance (see Tables 2 and 3 of Paper I):


where is a shape parameter and and are the axes of the superellipse:


Fig. 1 shows contours of log() in the small scale regime averaged within the virialized region of Aq-A-2 at (solid lines). The dotted lines are the fitting contours through Eqs. (2-3) 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.

Table 1: Values of the fitting parameters in Eq. (3) at over the virialized region of the Aq-A-2 halo. Together with Eqs. (2) and (3), these parameters provide a good description of  at small (see Fig. 1).

3 Dark matter annihilation rate: substructure boost

The number of dark matter self-annihilation events per unit time is given in terms of the phase space distribution at the phase space coordinates and :


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 pairs111Note 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:


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. (2-3) 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 Aq-A-2 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:


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. (2-3) down to:


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:


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 Aq-A-2 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:


As an example, let us take the case where and compute the total annihilation rate in resolved substructures at within a MW-size halo. Since (see Table 1), according to Eq. (6) we have:


If we take the value of given in Eq. (7) and the maximum value of  we can resolve in Aq-A-2, , then we can estimate the “resolved” substructure boost to the total annihilation rate in the Aq-A-2 halo:


The simulation particle mass of Aq-A-2 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 Aq-A-1 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 physically-motivated 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 two-point correlation function in the strongly non-linear 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 halo222From 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):


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 small-scale regime) with Eq. (12). We can then associate a constant value of within the simulated MW-size 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):


In the spherical collapse model, the subhalo collapses when the r.m.s top-hat linear overdensity (mass variance) crosses the linear overdensity threshold at an epoch given by:


The mass variance is defined by:


where is the top-hat 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:


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 :


Thus, we can finally give the prediction of the improved stable clustering hypothesis for the curves of contours of constant :


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


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:


where is the characteristic free-fall 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:




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:


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:


where and are free parameters.

The infall redshift can be estimated using the Extended Press-Schechter formalism. We are interested in the probability that a halo of mass , had progenitors of mass :


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:


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 halo-centric 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 .

Figure 2: Contours of the logarithm of the particle phase space average density () of the Aq-A-2 halo (solid) and for the model described in Section 4 (dotted) at different redshifts, clockwise from top left: (black), (red), (cyan) and (orange). The dashed regions in the left corners encompass the regions where resolution is a potential issue (see Appendix of Paper I).

After , the density does not dilute anymore as the Universe expands. If we then have an additional solution to Eq. (20) for :




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 halo-centric distance): , and accounting for the stretching in phase space due to tidal shocking (Eq.19). This model is able to describe the small-scale 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 Aq-A-2 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:


while and are given in Table 2.

Table 2: Values of the fitting parameters in Eq. (19) that together with the tidal disruption model explained in Section 4.1 provide a good description of  at small (see Fig. 2).

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


5.2 Halo-centric distance dependence

Figure 3: Contours of the logarithm of the particle phase space average density () for the Aq-A-2 halo at for samples of particles taken at different radial shells from the halo centre. Left panel: (red) and (blue). Right panel: (magenta), (cyan), and (violet). The dotted lines are fits given by our model described in section 4. The relevant parameters of the fits are given in Table 3. For reference, the dashed line is  averaged over the whole halo.

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 Aq-A-2 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:


where the radial-dependent parameters , and , as well as the shape parameters , and , that fit the simulation data are given in Table 3.

Table 3: The second, third and fourth columns are the values of the parameters in Eq. (3) for the fit to  averaged over different radial shells centered in the Aq-A-2 halo as indicated in the first column (see Fig.3). The last three columns are the parameters of a simple approximation (Eq. 31) of the parameter that controls tidal stripping (; see Section 4.1).

6 Boost to Dark Matter Annihilation due to unresolved sub-structure

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 concentration-mass relation. The support for this argument is two-fold:

  • 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 MW-size halo at in the case of :


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 non-trivially 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 Aq-A-2 halo.

Figure 4: Projection of the particle phase space average density () for the Aq-A-2 halo at in the direction fixing  kpc (black solid). The red dashed line is the fit with our model using the same cut in , while the black and blue dashed lines are extrapolations to using our full model and using the fitting function in Eq. (2), respectively. The horizontal axis is extended all the way to the value of corresponding to M, the typical mass of WIMP models.

The projection of  along the direction can be seen in Fig. 4 for the whole virialized region of Aq-A-2 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):


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 re-scaling 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 power-law tail in this PDF) calibrated with the MW-sized 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 Millennium-II simulation (Boylan-Kolchin et al., 2009), implicitly extrapolating the subhalo mass function and concentration-mass 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.

Figure 5: Total subhalo boost to the annihilation rate (assuming ) within the virialized region of a MW-size halo as a function of the minimum subhalo mass at collapse. The filled circle shows the value of the boost estimated directly from the fit to the simulation data (Eq. 11), while the star symbols are the values in the resolved (M) and extrapolated cases (M) taken from Springel et al. (2008). The filled triangle is the boost estimated by Kamionkowski et al. (2010) using a method calibrated with the Via Lactea II simulation (Diemand et al., 2008). The dotted region shows the theoretical expectation for the decoupling mass limit of WIMPs (Bringmann, 2009).

6.1 Sommerfeld-enhanced models

If the annihilation cross-section is enhanced by a Sommerfeld mechanism (e.g. Arkani-Hamed 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 Sommerfeld-enhanced model, we generically approximate the cross section as:


The value of is commonly near (the so-called “” boost), but it can reach near resonances; we only consider the former case. For Eq. (6.1) and using our modelling we obtain:


where is the collapse mass corresponding to the characteristic velocity at saturation ; using Eq. (19):

Figure 6: Same as Fig. 5 but for a Sommerfeld-like model where with and different saturation values as shown in the legend. The dotted line is an approximation to one of the benchmark points (BM1) presented in Finkbeiner et al. (2011) that avoid different astrophysical constraints and can account for high energy cosmic ray excesses.The dotted region shows the theoretical expectation for the decoupling mass limit of WIMPs (Bringmann, 2009).

Fig. 6 shows the subhalo boost, relative to the smooth dark matter distribution of the host halo, for Sommerfeld-like models where the cross section scales as in Eq. (6.1) with . To compute the Sommerfeld-enhanced 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 Aq-A-2 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 Sommerfeld-enhanced 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 Halo-centric 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 phase-space 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.

Figure 7: Subhalo boost to the annihilation rate as a function of halo-centric distance for a minimum subhalo mass at collapse of M. The boost was computed for concentric shells with a thickness of in the cases of constant cross section (circle symbols) and for an approximation of the Sommerfeld-enhanced model BM1 presented in Finkbeiner et al. (2011) (star symbols).

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 Sommerfeld-like 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 Sommerfeld-like 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 non-gravitational (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 concentration-mass relation. Extending these functions as power-laws 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 coarse-grained 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 physically-motivated explanation of the two-dimensional functional shape of , and calibrated to the MW-size Aquarius simulations, gives a firm basis to extrapolate into the unresolved substructure regime. In summary, the main advantages of our method are two-fold:

  • 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 “zero-separation” limit of  multiplied by the annihilation cross section times the relative velocity (see Eq. 4). This allows us to accommodate any velocity-dependence coming from particle physics models in a straightforward way.

Although our model has several free functions that are fitted to the MW-size 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 MW-halo centre (Tables 2 and 3), while for the latter we present simple parameterisations (see Tables 2 and 3, Eqs. 29-30).

As a sample application of our model, we computed the subhalo boost (over the smooth dark matter distribution) in a MW-size halo, both globally (i.e. over the whole virialized halo) and locally as a function of halo-centric distance, for cases where the annihilation cross section is constant and for a generic Sommerfeld-like 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 concentration-mass 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 (coarse-grained) 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 Interested users should refer to that site for instructions on how to use the code.


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 GADGET-3 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 ( and Compute/Calcul Canada. The Dark Cosmology Centre is funded by the DNRF.


  • 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
  • Arkani-Hamed et al. (2009) Arkani-Hamed 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
  • Boylan-Kolchin et al. (2009) Boylan-Kolchin 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ánchez-Conde M. A., Siegal-Gaskins 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
  • Siegal-Gaskins (2008) Siegal-Gaskins 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., Boylan-Kolchin 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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description