6.6. (Sub-)haloes (clumps.h & stat.h)

CLUMPY needs to manipulate and calculate many quantities related to generic DM (sub-)haloes. Their most general description takes into account that the total DM distribution is the sum of a smooth contribution and a distribution of sub-halos. The sub-halos can then be seen themselves as scale-down versions of their parent (“host”) halo. Hence, multi-level of substructures should be accounted for. For the specific case of DM decay, if only the mean contribution is calculated, all these details are irrelevant (the signal is basically proportional to the mass, which remains the same whether we put some of it in sub-halos or not). For DM annihilation, they are important as sub-halos can boost the signal (note however that after the second level, the contributions are generally negligible).

The libraries clumps.h and stat.h contain functions to calculate,

  • mass of DM halo \(M_\Delta\), \([M_\Delta]={\rm M}_\odot\);
  • the properties of subhaloes (distributions, mass-concentration, etc.);
  • the intrinsic luminosity (see below), and J-factor for any host halo, spherical and triaxial, with or without sub-structures;
  • the average and variance of galactic subhaloes or field haloes (not yet implemented for the extragalactic calculation).

6.6.1. Overdensity \(\Delta\) (\(R_\Delta\) and \(M_\Delta\))

All integrations for (sub-)haloes must be carried out to their maximum radius. However, the radius of a halo is an ill-defined quantity. In general, the mass \(M_{\Delta}\) of a halo is connected to its size, \(R_{\Delta}\), via the relation

(6.8)\[R_{\Delta}( M_{\Delta}, z) =\left( \frac{3\, M_{\Delta}}{ 4\pi \times \Delta(z) \times \varrho_{\rm c}(z)}\right)^{1/ 3} \times (1 + z)\,,\]

where the subscript \(\Delta\) denotes a characteristic collapse overdensity above the critical density of the Universe, \(\varrho_{\rm c}= 3H^2(z)/8\pi G\). All calculations can be chosen to be performed with respect to any of the following overdensity definitions:

\[\begin{split}\Delta(z) & = {\rm const.},\\ \Delta(z) &= {{\rm const.} \times \Omega_{\rm m}(z)} =: {\Delta_{\rm m} \times \Omega_{\rm m}(z)},\\ \Delta(z) &= 18\pi^2 + 82\,[\Omega_{\rm m}(z) - 1)]- 39\,[\Omega_{\rm m}(z) - 1]^2 \,.\end{split}\]

If a mass-concentration \(c_{\Delta}(M_{\Delta},z)\) or halo mass function \({\mathrm{d}n}/{\mathrm{d}M}(M,\,z)\) are natively provided for a \(\Delta\) different from the user’s choice, \(c_{\Delta}(M_{\Delta},z)\) or \({\mathrm{d}n}/{\mathrm{d}M}(M,\,z)\) are internally rescaled to the user-chosen Delta using the algorithm described in appendix A of Hütten et al. (2018). Note that this rescaling presumes a given halo profile (see Section 6.5).

6.6.2. Concentration: \(c_\Delta-M_\Delta\)


The concentration \(c_\Delta\), at a given characteristic overdensity \(\Delta\), is defined to be

\[c_\Delta\equiv \frac{R_\Delta}{r_{-2}},\]

where \(R_\Delta\) is the radius of the DM halo for which the density equals this overdensity (see above), and \(r_{-2}\) is the position where the slope of the DM halo density profile reaches \(-2\) (see Section 6.5). A list of all implemented \(c_\Delta-M_\Delta\) descriptions is given in Mass-concentration relations.

Implementation in CLUMPY

We chose to code the generic \(c_\Delta-M_\Delta\) relationship (with \(\Delta=200,\, {\rm vir},\,...)\) and then rely on a conversion function to translate between different choices of \(\Delta\). The redshift dependence of the concentration is taken into account; haloes formed at earlier times are less concentrated. There is generally a dependence on the environment and an intrinsic scatter of the \(c_\Delta-M_\Delta\) relationship (see below). However, effects of halo selection (relaxed or not…) are not implemented. They can lead to \(\sim 10\%\) differences in the \(c_\Delta-M_\Delta\) relationship (see, e.g. Duffy et al. 2008 and App. of Prada et al. 2011).


Most of the parametrisations deal with galaxy cluster or down to galaxy size objects, \(M_{\Delta}\gtrsim [10^{10}-10^{15}]\,M_\odot\). Only kB01_VIR, kENS01_VIR, kGIOCOLI12_VIR, and kSANCHEZ14_200 can be extrapolated down to the lowest halo masses, \(M_{\Delta}\ll 10^{10}\,M_\odot\).


In Fig. 6.12, we plot the \(c_\Delta-M_\Delta\) relationships for selected prescriptions available in CLUMPY, and in Fig. 6.13 the intrinsic luminosity \({\cal L}(M_{vir})\) for all profiles:


Fig. 6.12 Concentration-mass relations.


Fig. 6.13 Intrinsic halo luminosities.

Fig. 6.12 and Fig. 6.13 can be reproduced via

$ clumpy -e2 -D

which creates a default parameter file. Then the gEXTRAGAL_FLAG_CDELTAMDELTA_LIST can be expanded within the parameter file or simply overwritten in the command line:

$ clumpy -e2 -i clumpy_params_e2.txt --gEXTRAGAL_FLAG_CDELTAMDELTA_LIST=kB01_VIR,kENS01_VIR,kNET007_200,kDUFFY08F_VIR,kDUFFY08F_200,kDUFFY08F_MEAN,kETTORI10_200,kPRADA12_200,kGIOCOLI12_VIR,kSANCHEZ14_200

As calculating ten times the substructure contribution for each description may be rather time-consuming, the signal boost from substructures can be disregarded by overwriting the default value in the parameter file via the command line:

$ clumpy -e2 -i clumpy_params_e2.txt --gDM_SUBS_NUMBEROFLEVELS=0 --gEXTRAGAL_FLAG_CDELTAMDELTA_LIST=kB01_VIR,kENS01_VIR,kNET007_200,kDUFFY08F_VIR,kDUFFY08F_200,kDUFFY08F_MEAN,kETTORI10_200,kPRADA12_200,kGIOCOLI12_VIR,kSANCHEZ14_200


ROOT needs to be installed and be linked against CLUMPY to take advantage of CLUMPY’s built-in pop-up graphics. See the Download and Installation section for details.

6.6.3. Subhaloes: distributions

The distribution of subhaloes (or substructures) in a host halo is a function of their distance from the host centre and of their mass. There are some environment effects (e.g. subhaloes are expected to be disrupted by tidal forces if baryons are present) dependent on the position, that correlate the spatial and mass distribution. However, in this release of CLUMPY, we assume that these distributions are not correlated. In that case, we can express the total distribution as a product of three probability distribution functions (PDF) and a normalisation

\[\frac{{\rm d} N}{{\rm d}^{3}\vec{r}\ {\rm d}M} = N_{tot}\, \frac{{\rm d}{\cal P}_V(r)}{{\rm d}V}\, \times \frac{{\rm d}{\cal P}_M(M)}{{\rm d}M} \times \frac{{\rm d}{\cal P}_c}{{\rm d}c}(M,c)\ \textrm{,}\]

where \(N_{tot}\) is the total number of subhaloes within the virial radius of the parent halo \(R_{\Delta}\), and the other quantities are distributions for the mass, position, and concentration of the DM haloes, i.e.

\[\int_{M_{\rm min}}^{M_{\rm max}} \frac{{\rm d}{\cal P}_M(M)}{{\rm d}M}\, {\rm d} M = 1, \quad \int_{0}^{R_{vir}} \frac{{\rm d}{\cal P}_V(r)}{{\rm d}V}\, {\rm d}^{3}\vec{r} = 1, \quad {\rm and} \int_{c_{\rm min}(M)}^{c_{\rm max}(M)} \frac{{\rm d}{\cal P}_c(c,M)}{{\rm d}c}\, {\rm d}c= 1.\]

Mass distribution \({\rm d}{\cal P}_M/{\rm d}M\)

The distribution \({\rm d}{\cal P}_M/{\rm d}M\) is assumed to be a power law

\[\frac{{\rm d}{\cal P}_M}{{\rm d}M} = B M^{-\alpha_M}\,,\]

where \(\alpha_M\) is the slope \(\sim 1.9\). All functions using \({\rm d}{\cal P}_M/{\rm d}M\) assume that the normalisation \(B\) has been correctly set (so as to have a probability).

Spatial distribution \({\rm d}{\cal P}_V/{\rm d}V\)

\({\rm d}{\cal P}_V/{\rm d}V\) provides the spatial distribution of substructures in a host halo, the latter being the Galactic DM halo or any other DM halo (a dSph galaxy, another galaxy, a galaxy cluster). The distribution can be chosen from any DM profiles as used for the DM density profile, and is provided by the input parameters ending with SUBS_DPDV_FLAG_PROFILE, see Section and Section 10.2.

Concentration distribution \({\rm d}{\cal P}_c/{\rm d}c\)

\({\rm d}{\cal P}_c/{\rm d}c\) provides the concentration distribution of the substructures in the host halo, the latter being Galactic DM halo or any other halo (e.g. a dSph galaxy). It is triggered by the value set for gDM_LOGCDELTA_STDDEV:

  • gDM_LOGCDELTA_STDDEV = 0: use the mean concentration \(\bar{c}(M)\), as in CLUMPY version 1:

    \[\frac{{\rm d}{\cal P}_c}{{\rm d}c}(M,c) = \delta\left(\bar{c}(M)\right).\]
  • gDM_LOGCDELTA_STDDEV > 0: log-normal distribution

    \[\frac{{\rm d}{\cal P}_c}{{\rm d}c}(M,c) = \frac{\exp^{\displaystyle -\left [\frac{\ln{c} - \ln({\bar{c}(M)})}{\sqrt{2}\sigma_c(M)} \right ]^2}}{\sqrt{2 \pi}\,\,c\,\,\sigma_c(M)},\]

    i.e., the concentration of a halo of mass \(M\) is randomly drawn from the above distribution around the mean concentration \(\bar{c}(M)\) and using a typical scatter \(\sigma_c(M)\sim 0.2\) (Bullock et al. 2001). This distribution is normalised, so it is directly a PDF. Their works give a \(c_{200}-M_{200}\) relationship instead, this time working at the radius within which the average density is 200 times the critical density \(\varrho_{\rm crit}(z)\) of the universe (or 200 times the mean background density \(\varrho_{\rm mean}(z)\)). Subtle variations of some definitions exist in the literature and care is needed when trying to compare these parametrisations (see appendix B of Giocoli et al. (2010) for an example).

6.6.4. Luminosity \(\cal L\) and boost

The intrinsic luminosity \(\cal L\) of a halo is defined to be the total signal over the halo volume, i.e.

\[{\cal L} \equiv \int_{V_{\rm cl}} \rho^i(r)\; {\rm d}V,\]

with \(i=1\) for decay (in which case \({\cal L}\) is the mass, \({\cal L} = M_\Delta\)) and \(i=2\) for annihilation, with, e.g., \([{\cal L}] = \rm M_\odot^2\,\,kpc^{-3}\) (internal units of CLUMPY).


For annihilation, the luminosity of a halo is boosted for high redshifts, as the Universe was denser (w.r.t. a comoving description) at earlier times. Therefore, one has to carefully distinguish between a proper and comoving luminosity. In CLUMPY’s \(\tt -e2\) module (see Section 7.3), we output the proper (“physical”) luminosity. The comoving luminosity is defined as

\[{\cal L}_{\rm comoving} = \frac{{\cal L}_{\rm proper}}{(1+z)^3}\,,\]

obtained when calculating the luminosity in comoving coordinates, i.e., using \(R_\Delta\) as in (6.8). However, the observable quantity is the proper luminosity.

We also give the annihilation luminosity in the \(\tt -e2\) module in units of \(M_\Delta \times \varrho_{\rm m,0}\). For a NFW profile (see Section 6.5.5), the (comoving) annihilation luminosity, \({\cal L}(M_\Delta,z) = {\cal L}[c_\Delta(M_\Delta,z), z]\), can be given in a concise analytic form:

\[\frac{{\cal L}_{\rm comoving}(M_\Delta,z)}{M_\Delta \times \varrho_{\rm m,0}} = \frac{\Delta(z)}{9\,\Omega_{\rm m}(z)}\times f(c_\Delta)\,,\]

with the term

(6.9)\[f(c_\Delta) = c_\Delta^3 \times \frac{1-(1+c_\Delta)^{-3}}{\left[\ln(1+c_\Delta)-\frac{c_\Delta}{1+c_\Delta}\right]^2}\,.\]


The numerator \(1-(1+c_\Delta)^{-3}\) in (6.9) is often discarded in the literature, as it close to \(1\) for any \(c_\Delta\gtrsim 1\). However, it becomes significant for large redshifts, for which concentrations are expected to decrease with approximately \(c(z)= c_0/(1+z)\). Therefore, \(c \rightarrow 0\) for \(z\rightarrow \infty\), and interestingly,

\[\lim_{c \to 0}\, f(c) = 12 \,,\]

such that in the limit, a luminosity of \(4/3\) times the one of constant-density spheres of \(\Delta\) times the critical density is reached for the NFW case. In turn, for density profiles regular at their centers (\({\rm d} \rho/{\rm d}r(r\rightarrow 0) < \infty\)), constant-density spheres of exactly \(\Delta\) times the critical density are expected.

Also note that the proper luminosity consequently always diverges with at least \(\propto (1+z)^{3}\).

For a point source, the J-factor is simply given by \(J={\cal L}/d^2\) (\(d\): distance of the halo). Subhaloes in their host halo can generally be assumed to be point-like (size encompassed by angular resolution of instrument), and this property is taken advantage of to calculate statistical properties for a population of subhaloes. The intrinsic luminosity \({\cal L}_n\) from including \(n>1\) levels of substructures is recursively obtained from \({\cal L}_{n-1}\),

\[{\cal L}_{n}(M,c) = {\cal L}_{\rm sm}(M) + {\cal L}_{\rm cross-prod}(M) + N_{\rm tot}(M)\int_{M_{\rm min}}^{M_{\rm max}(M)} {\cal L}_{n-1}(M')\frac{{\rm d}{\cal P}}{{\rm d}M'}(M') {\rm d}M\,,'\]


\[\begin{split}{\cal L}_{\rm sm}(M)&\equiv\int_{\rm V_{\rm cl}} \left[\rho_{\rm cl}^{\rm sm}(M)\right]^2 \,{\rm d} V\,;\\ {\cal L}_{\rm cross-prod}(M)&\equiv2\int_{\rm V_{\rm cl}} \rho_{\rm cl}^{\rm sm}(M)\,\langle\rho_{\rm subs}(M)\rangle \,{\rm d} V\,.\end{split}\]


\({\cal L}(M,c)\) depends on both the host mass and concentration (the latter being calculated from the mass. It can be calculated with or without contributions from substructures in the DM halo under scrutiny.

Below, we illustrate the effect of taking into account multiple levels of substructures and/or the distribution of \({\rm d}{\cal P}/{\rm d}c\) on the intrinsic luminosity (for annihilating DM). The three plots in Fig. 6.14 from left to right correspond to:

  • Calculation with (gDM_SUBS_NUMBEROFLEVELS = 1) and without (gDM_SUBS_NUMBEROFLEVELS = 0) substructures, that is the so-called boost-factor;
  • Impact of accounting for more than the first level of substructures (i.e. gDM_SUBS_NUMBEROFLEVELS > 1) w.r.t. the previous level, and also the impact of accounting for a distribution of concentrations (i.e. gDM_LOGCDELTA_STDDEV > 0).
  • Detailed view of the gDM_LOGCDELTA_STDDEV > 0 case, and how it shifts \(\langle {\cal L}\rangle_{c}\) w.r.t. the luminosity \({\cal L}\) in the gDM_LOGCDELTA_STDDEV > 0 case.

Fig. 6.14 Boost factor for kB01_VIR, kSANCHEZ14_200, and \({\rm d}P/{\rm d}M\) slope \(\alpha_{ M}=2\) (left); Impact of having \({\rm d}P/{\rm d}c\) distribution (blue), and adding extra-levels of subhaloes (middle); For two \(<c>\) values, \({\rm d}P/{\rm d}c\) log-normal distribution (blue) multiplied by \({\cal L}\) (magenta) shift the result (solid red line is the integrand of \(<{\cal L}>_c\)), hence the blue curve in the above plot (right). Click to enlarge the figures.

Fig. 6.14 can be reproduced via (provided ROOT is installed):

$ clumpy -e2 -D; clumpy -e2 -i clumpy_params_e2.txt

6.6.5. Mean and variance

A typical DM halo of \(10^{12}M_\odot\) (Milky-Way like) contains up to \(10^{14}\) substructures, which renders the explicit calculation of the signal summed over all haloes prohibitive. This huge number allows the use of the continuum limit as the subhalo positions and masses are random variables, drawn from distribution functions (described above) obtained by N-body numerical simulations and/or semi-analytical calculations

At first order, a random variable (e.g. the mass and position of substructures) is described by its average value and variance. Departure from the average can arise if a small number of objects contribute significantly to the total J-factor, which happens if a massive subhalo dominates, or if one of the smallest subhaloes (the smaller, the more numerous they are) is sitting almost at the observer location. The latter configuration only happens for subhaloes in the Galaxy, since substructures in dSphs or extragalactic objects are far away.

As presented in Section 6.1 for the galactic halo, J-factor skymaps rely on a combination of the calculation of the average signal and the calculation of individual drawn clumps above and below a critical distance \(l_{\rm crit}\): the critical distance is obtained by requiring the relative error of the signal integrated from \(l_{\rm crit}\) to remain lower than a user-defined precision requirement. This strategy ensures a controlled and extremely quick calculation of skymaps: the number of subhaloes to draw in the Galaxy is reduced from a few tens of thousands to a few hundreds depending on the configuration.

Average mass and the mean of (some power of) the distance:

\[\begin{split}\langle M\rangle \!&=\!\!\! \int_{M_{\rm min}}^{M_{\rm max}}\!\!\!\!\! M\frac{{\rm d} {\cal P}_M}{{\rm d} M} {{\rm d} M},\\ \langle l^n \rangle \!&=\!\!\! \int_0^{\Delta\Omega}\!\!\!\! \int_{l_{\rm min}}^{l_{\rm max}}\!\!\!\! l^{\,n} \frac{{\rm d} {\cal P}_V}{{\rm d} V} l^{\,2} {\rm d} l \, {\rm d}\Omega.\end{split}\]

Mean luminosity over \(M\) and \(c\):

\[\langle {\cal L} \rangle \!=\!\!\! \int_{M_{\rm min}}^{M_{\rm max}} \!\!\frac{{\rm d} {\cal P}_M}{{\rm d} M}(M) \!\!\! \int_{c_{\rm min}(M)}^{c_{\rm max}(M)} \!\!\frac{{\rm d} {\cal P}_c}{d c}(M,c) \,{\cal L}(M,c)\, {\rm d} c \,{\rm d} M,\]

Average and variance on J:

\[ \begin{align}\begin{aligned}\langle J_{\rm subs}\rangle &= N_{\rm tot}\!\!\int_0^{\Delta\Omega} \!\! \int_{l_{\rm min}}^{l_{\rm max}} \frac{{\rm d}{\cal P}_V}{{\rm d} V}(l,\Omega) \,{\rm d} l \,{\rm d}\Omega \langle {\cal L} \rangle\\\sigma^2_{\rm subs} &= \langle {\cal L}^2\rangle \left\langle\frac{1}{l^4}\right\rangle - \langle{\cal L}\rangle^2 \left<\frac{1}{ l^2}\right>^2.\end{aligned}\end{align} \]