6.5. Halo profiles (profiles.h)

6.5.1. Formal description

A DM profile is characterised by a normalisation \(\rho_{\rm s}\), a scale radius \(r_{\rm s}\), and several shape parameters \(\vec{\alpha}\). For a spherically symmetric DM halo, we have formally

\[\rho(r) = f(\rho_{\rm s},r_{\rm s}, \vec{\alpha})\,.\]

The physical size of the DM halo is usually taken to be \(R_{\rm \Delta}\), the radius at overdensity \(\Delta\). All profiles and their corresponding input parameter keyword names are gathered in Section 10.2. Their calculation is implemented in the library profiles.h.

6.5.2. Saturation radius \(r_{\rm sat}\)

Some of these DM profiles are steep enough in the inner regions to lead to a singularity of the \(\gamma\)-ray luminosity at the centre of the halo. However, a cut-off radius \(r_{\rm sat}\) naturally appears, within which the DM density saturates due to the balance between the annihilation rate and the gravitational infalling rate of DM particles. The saturation density reads (Berezinsky et al., 1992):

\[\rho_{\rm sat}= 3.10^{18} \left(\frac{m_\chi}{100~\rm GeV}\right) \times \left( \frac{10^{-26} {\rm cm}^3~{\rm s}^{-1}} {\langle \sigma v\rangle}\right) M_\odot~{\rm kpc}^{-3}.\]

The saturation radius is defined to be the radius for which \(\rho(r_{\rm sat})=\rho_{\rm sat}\), so that \(\rho(r<r_{\rm sat})=\rho_{\rm sat}\). This is a very crude description, but this cut-off only arises for the steepest profiles which are generally disfavoured. For the various family of profiles we use in CLUMPY, we have:

  • kEINASTO/kEINASTO_N family: for any model in Einasto, the profile is not cuspy, so \(\rho_{\rm sat}\) is assumed to never be reached.
  • kZHAO family: \(r_{\rm sat} = r_s \times \left(\frac{\rho_{\rm s}}{\rho_{\rm sat}}\right)^{1/\gamma}\).
  • kDPDV_GAO04: if \(\beta-3<0\): \(r_{\rm sat} = r_{200} \times \left\{\left(\frac{\rho_{\rm sat}}{\rho_{200}}\right) \frac{\left[\beta-\alpha ac/(1+ac)\right]}{\beta(1+ac)}\right\}^{1/(\beta-3)}\;\).


In CLUMPY, the user cannot select \(r_{\rm sat}\) for a particular halo, but must globally set a physically motivated value for \(\rho_{\rm sat}\), \([\rho_{\rm sat}]={\rm M}_\odot\,{\rm kpc}^{-3}\), for all (sub)halos involved in the simulation. This is done by the gDM_RHOSAT input parameter.

6.5.3. Radius \(r_{-2}\)

The concentration of a DM halo, defined in Section 6.6, depends on \(r_{-2}\), the radius for which the slope is \(-2\):

\[\left. \frac{{\rm d}~[r^2 \rho(r)~]}{{\rm d}r}\right|_{r=r_{-2}} = 0\quad \Leftrightarrow\quad \left.\frac{{\rm d}\log \rho}{{\rm d}\log r}\right|_{r=r_{-2}} = -2.\]

For any specific profile, it can be related to the scale radius (\(r_{\rm s},\,r_{0},\,r_{\rm e},...\)) defining a given halo’s dimension, as long as the outer slope of the profile is steeper than \(-2\). Values for \(r_{-2}\) are gathered in Section 10.2. To highlight the position of \(r_{-2}\) in various profiles, Fig. 6.8 shows the quantity \(r^2\rho(r)\) for Milky-Way like halos.


Fig. 6.8 Radius \(r_{-2}\) for various DM profiles (see Section 10.2) at the maximum of the curves.

6.5.4. Triaxial profiles

  1. Definition:

    Triaxial profiles are described for instance in Jing & Suto, ApJ 574, 538 (2002). We refer the reader to the Wikipedia page for a general description of an ellipsoid: the points \((a,0,0)\), \((0,b,0)\) and \((0,0,c)\) lie on the surface and the line segments from the origin to these points are called the semi-principal axes of length \(a,b,c\). They correspond to the semi-major axis and semi-minor axis of the appropriate ellipses.

  2. Convention in CLUMPY:

    In CLUMPY, we allow for a general triaxial model \((a\neq b\neq c)\), where the major axis is \(a\) the intermediate axis \(b\), and the minor axis is \(c\). Cosmological simulations have shown that these axes can vary as a function of the iso-density radius \(R_{\rm iso}\) (Jing & Suto, ApJ 574, 538 2002), but it is not considered here. In any case, the position \((X,Y,Z)\) in the coordinate framework attached to a given DM halo centre (see Fig. 6.2) corresponds to the iso-density radius

    \[R_{\rm iso} =\sqrt{\frac{X^2}{a^2} + \frac{Y^2}{b^2} + \frac{Z^2}{c^2}}.\]

    The density is given by \(\rho(R)\), where \(\rho\) is any spherical profile discussed above.

  3. Initialisation parameters:

    Triaxial parameters consist of 7 parameters (three axes for the shape, three Halo rotation (Euler angles), and one boolean to decide whether to take into account or not the triaxiality). Depending of the type of halo considered, triaxiality is implemented in different fashion in CLUMPY:

    • Galactic halo:

      Triaxial parameters for the Galactic halo gMW_XXX are set directly in the parameter file or over the command line as given in Table 6.5:

      Table 6.5 Setting triaxiality for the Milky Way halo
      Quantities (triaxial profile) Galactic halo: parameter names
      Is triaxiality enabled? (0 or 1) gGAL_TRIAXIAL_IS
      \(a\;\;[-]\) gGAL_TRIAXIAL_AXES_0
      \(b\;\;[-]\) gGAL_TRIAXIAL_AXES_1
      \(c\;\;[-]\) gGAL_TRIAXIAL_AXES_2
      \(\alpha_{\rm Euler}\;\;[\rm deg]\) gGAL_TRIAXIAL_ROTANGLES_0
      \(\beta_{\rm Euler}\;\;[\rm deg]\) gGAL_TRIAXIAL_ROTANGLES_1
      \(\gamma_{\rm Euler}\;\;[\rm deg]\) gGAL_TRIAXIAL_ROTANGLES_2

    • User-defined list of halos:

      Triaxial parameters are set directly in the file of user-defined haloes as extra-columns in the file (regardless of the type of the DM halo, be it a galaxy, a dwarf spheroidal galaxy, or a cluster of galaxy), see Format of halo definition files for the correct format.

    • Substructures:

      Triaxiality is not enabled/implemented in CLUMPY, since it would require a statitical distribution of axis ratios and orientations. First, the orientation is not important if the integration angle encompasses the halo volume (and this is the case for all subhaloes but a tiny fraction). Second, the exact distribution of axis ratios is not known, but the impact in the calculation is certainly negligible compared to other uncertainties (concentration, profile, etc.). In the end, it is not worth including this option in the calculation.

6.5.5. Family of profiles Zhao (kZHAO)

Many DM profiles used in the literature can be written in terms of a three parameter \((\alpha,\beta,\gamma)\) family (Hernquist 1990 and Zhao 1996):

\[\rho\left(r\,|\,r_{\rm s},\rho_{\rm s};\alpha,\beta,\gamma\right)=\frac{2^{\frac{\beta-\gamma}{\alpha}}\,\times\,\rho_{\rm s}}{\left(\frac{r}{r_{\rm s}}\right)^\gamma \times \left[1+\left(\frac{r}{r_{\rm s}}\right)^\alpha\right]^{\frac{\beta-\gamma}{\alpha}}}\;.\]

For this profile, we have \(r_{-2} = r_{\rm s} \times \left( \frac{\beta-2}{2-\gamma} \right)^{-\frac{1}{\alpha}}\), with

  • \(r_{\rm s}\): scale radius, \([r_{\rm s}]={\rm kpc}\),
  • \(\rho_{\rm s}\): density normalisation \(\rho(r_{\rm s})\), \([\rho_{\rm s}]={\rm M}_\odot\,{\rm kpc}^{-3}\),
  • \(\alpha\): logarithmic transition slope (controls the sharpness of the transition between the inner and outer slopes),
  • \(\beta\): outer logarithmic slope,
  • \(\gamma\): inner logarithmic slope.

Widely-used special cases of the Hernquist-Zhao profile are: Einasto profiles

In the last decade, it has been realised that the inner slope of the DM halo should decrease as \(r\) decreases, which has given rise to several logarithmic inner-slope parametrisations, generally based on the Einasto profile which follows

(6.7)\[\rho(r)\propto \exp\left(-Ar^\alpha\right)\quad \Leftrightarrow \quad \frac{{\rm d}\log \rho}{{\rm d}\log r} = -\alpha A \,r^{\alpha},\]

and was first used by Einasto to describe the mass and light distributions in galaxies. Note that it is formally identical to Sersic’s law which is traditionally applied to projected (2D) density profiles rather that 3D density profiles. In the literature, two main expressions of the Einasto density profile of DM haloes have been formulated (see also Retana-Montenegro et al. (2012) for a nice and comprehensive discussion of the properties of Einasto profiles):

  1. Standard Einasto description (kEINASTO):

The standard Einasto is used for instance in Navarro et al. (2004) and Springel et al. (2008), where \(A = \frac{2}{\alpha}\, r_{-2}^{\;-\alpha}\)

\[\Rightarrow \quad \rho^{\rm EINASTO}\left(r\,|\,r_{-2},\rho_{-2};\alpha\right)=\rho_{-2}\, \exp\left\{-\frac{2}{\alpha}\left[\left(\frac{r}{r_{-2}}\right)^\alpha -1\right]\right\}\;,\]

and both teams find that the shape parameter \(\alpha=0.17\) as a good description of the density profile of DM halos. Springel et al. (2008) also find that \(\alpha=0.68\) describes well the spatial distribution of subhaloes \({\rm d}{\cal P}/{\rm d}V\) in the Milky-Way.

For the standard Einasto profile, we have \(\rho(r_{-2}) = \rho_{-2}\), with

  • \(r_{-2}\): radius, \([r_{-2}] ={\rm kpc}\), for which the slope is \(-2\),
  • \(\rho_{-2}\): normalisation, \([\rho_{-2}] = {\rm M}_\odot\,{\rm kpc}^{-3}\).
  • \(\alpha:\) parameter controlling the sharpness of the rolling logarithmic slope.
  1. Einasto ‘s \(r^{1/n}\) profile (kEINASTO_N):

Merritt et al. (2006) implemented a so-called Einasto’s “\(r^{1/n}\) model”,

\[\begin{split}\alpha = \frac{1}{n},\quad A&=d_n\, r_{\rm e }^{\;-1/n}\\[0.3cm] \Rightarrow \quad \rho^{\rm EINASTO_N}\left(r\,|\,r_{\rm e},\rho_{\rm e};n\right)& = \rho_{\rm e} \, \exp\left\{-d_n \times \left[ \left(\frac{r}{r_{\rm e}}\right)^{\frac{1}{n}} - 1 \right]\right\}\,,\end{split}\]

with \(d_n\) being given by the implicit equation

\[\Gamma(3n) = 2\,\gamma(3n,\,d_n)\,,\]

where \(\Gamma(a)\) is the complete gamma function and \(\gamma(a,x)\) the lower incomplete gamma function. For \(n\gtrsim 0.5\), \(d_n\) can be approximately written explicitly as a function of \(n\) (Merritt et al., 2006):

\[d_n \approx 3n-\frac{1}{3}+\frac{0.0079}{n}\,.\]

For \(n\) an integer, which is mandatory in CLUMPY, \(M_{\rm tot} \sim 2\pi \,n\, r_{\rm e}^3\, \rho_{\rm e}\, e^{d_n}\, d_n^{-3n} (3n-1)!\) (Merritt et al. 2006). Like for the standard Einasto profile, the profile kEINASTO_N provides good fits to N-body simulations with \(n=6\) (\({\,n}^{-1}={1}/{6}\approx 0.17\)).

For this profile, we have \(r_{-2} = r_{\rm e} \times \left( \frac{2 n}{d_n}\right)^n\) (Graham et al. 2006), with which the profiles kEINASTO and kEINASTO_N can be translated one into the other, with

  • \(r_{\rm e}\): scale radius, \([r_{\rm e}] ={\rm kpc}\), which contains half the mass of the halo,
  • \(\rho_{\rm e}\): normalisation, \([\rho_{\rm e}] = {\rm M}_\odot\,{\rm kpc}^{-3}\). With the above definition, we have \(\rho(r_{\rm s}) = \rho_{\rm e}\),
  • \(n\): shape parameter (required to be an integer in CLUMPY). Burkert (kBURKERT)

The Burkert profile (Burkert, 1995) is a phenomenological density distribution that resembles an isothermal profile in the inner regions (\(r\ll r_0\)), and with an outer logarithmic slope of \(-3\) as in a NFW profile:

\[\rho\left(r\,|\,r_{0},\rho_{0}\right)= \frac{\rho_0}{\left(1+\frac{r}{r_0}\right) \times \left[1+\left(\frac{r}{r_0}\right)^2\right]}\;.\]

For this profile, \(\rho(r_0) = \rho_0/4\)., and \(r_{-2}\) is the root of \(x^3-x-2=0\) (with \(x=r_{-2}/r_0\), i.e. \(r_{-2} \approx 1.5213797068 \times r_0\)) with

  • \(r_0\): scale (or core) radius, \([r_{0}] ={\rm kpc}\),
  • \(\rho_0\): core density \(\rho(r=0)\), \([\rho_{0}] = {\rm M}_\odot\,{\rm kpc}^{-3}\). \({\rm d}\mathcal{P}/{\rm d}V\) distribution

Besides the description of halo and subhalo profiles, we also categorise as profile the spatial distribution \({\rm d}\mathcal{P}/{\rm d}V\) of subhaloes in their host halo.

  • Anti-biased Einasto (kDPDV_SPRINGEL08_ANTIBIASED)

    Although not being a strict member of the family (6.7), a so-called anti-biased Einasto profile is sometimes used to described the population of substructures in its host halo (e.g., the Milky Way), where tidal effects tend to destroy haloes in the vicinity of the (Galactic) centre, hence the name anti-biased (Springel et al., 2008). In practice, this amounts to a simple multiplication by \(r/r_{-2}\), such that

    \[\begin{split}\rho^{\rm DPDV\_SPRINGEL08\_ANTIBIASED}\left(r\,|\,r_{\rm e},\rho_{\rm e};\alpha\right) &= \left(\frac{r}{r_{\rm e}}\right) \times \rho^{\rm EINASTO}(r) \\ &= \left(\frac{r}{r_{\rm e}}\right) \times \rho_{\rm e} \exp\left\{-\frac{2}{\alpha}\left[\left(\frac{r}{r_{\rm e}}\right)^\alpha -1\right]\right\}\;.\end{split}\]

    With this profile, we have we have \(r_{-2} = \left(\frac{3}{2}\right)^{1/\alpha} r_{\rm e}\), with

    • \(r_{\rm e}\): scale radius, \([r_{\rm e}] ={\rm kpc}\), for which the slope for the standard Einasto profile’s slope is \(-2\),
    • \(\rho_{\rm e}\): normalisation, \([\rho_{\rm e}] = {\rm M}_\odot\,{\rm kpc}^{-3}\). With the definition of the standard Einasto profile, we still have \(\rho(r_{\rm s}) = \rho_{\rm e}\),
    • \(\alpha\): shape parameter as for the standard Einasto profile.
  • Gao \((ac,\,\alpha,\beta)\) profile (kDPDV_GAO04)

    The fit function from Gao et al., (2004), see also Madau et al. (2008), with parameters \(a\), \(\alpha\), and \(\beta\) is used to model the distribution of \(N\) subclumps in a given host halo of concentration \(c(M)\) and scale radius \(r_{200}\). Note that these authors give

    \[N(<x) = \frac{(1+ac)\,x^\beta}{1+ac x^\alpha}, \quad {\rm with} \quad x=\frac{r}{r_{200}},\]

    so that using \(N(<x) = \int r^2 \rho(r)\,{\rm d}r\), with \(\rho(r)\) the profile we are looking for, we end up with

    \[\rho^{\rm GAO}\left(r\,|\,r_{200},\rho_{200};ac,\alpha,\beta\right)= \frac{\rho_{200}}{\beta-\alpha\frac{ac}{1+ac}}\times \frac{ (1+ac) \left(\frac{r}{r_{200}}\right)^{\beta -3} \left[\beta + ac (\beta - \alpha) \left(\frac{r}{r_{200}}\right)^\alpha\right]}{\left[1+ ac \left(\frac{r}{r_{200}}\right)^\alpha\right]^2}.\]

    Given this choice of normalisation, we have \(r_{-2}\) the root of a second order polynomial involving \((ac,\alpha,\beta)\), with

    • \(r_{200}\): scale radius, \([r_{200}] ={\rm kpc}\),
    • \(\rho_{200}\): normalisation \([\rho_{200}]={\rm M}_\odot\,{\rm kpc}^{-3}\), \(\rho_{200}=\rho(r_{200}) = \rho_{200}=200\times \rho_{\rm m}(z)\), 200 times the mean background matter density,
    • \(ac\): first ‘combined’ parameter of the Gao function: \(a\) is a parameter, \(c\) is the concentration of the host halo,
    • \(\alpha\):second parameter of the Gao fit (inner slope is \(\beta-3\)),
    • \(\beta\): third parameter of the Gao fit (outer slope is \(\beta-3-\alpha\)).

6.5.7. Plots \(\rho^2(r)\), \(r^2\rho(r)\),…

The following figures display a comparison of the various profiles in the central regions of the halos. The saturation density is reached for some of them \((1\,\rm GeV\,cm^{-3} = 2.63 \times 10^7~{\rm M}_\odot\,{\rm kpc}^{-3})\).


Fig. 6.9 DM smooth profile squared \(\rho^2(r)\) (normalised to \(\rho_0 = 0.3~\rm GeV\,cm^{-3}\) at \(r=8\,\rm kpc\)).


Fig. 6.10 Function \({\rm d}M/{\rm d}r \propto r^2\rho(r)\; \Rightarrow\;\) mass \(M =4\pi\int r^2 \rho(r)\, {\rm d}r\).


Fig. 6.11 Function \({\rm d}{\cal L}/{\rm d}r \propto r^2\rho^2(r) \; \Rightarrow\;\) intrinsic luminosity \({\cal L} = 4\pi\int r^2 \rho(r)^2 \,{\rm d}r\).


To calculate the mass \(M\) and luminosity \({\cal L}\) of a clump, the functions shown in the last two previous plots need to be integrated. If no analytical solution exists, we split the integration into several parts to speed up the calculation:

\[I \equiv \int\limits_0^{R_{\Delta}} f(r)\;{\rm d}r\; = \int\limits_0^{r=r_{\rm sat}} f(r)\;{\rm d}r \;\;+ \int\limits_{r_{\rm sat}}^{r=0.01\,r_{\rm s}} f(r)\;{\rm d}r\;\; + \int\limits_{r=0.01\,r_{\rm s}}^{R_{\Delta}} f(r)\;{\rm d}r\,.\]

The first part is integrated with a linear step integrator, and the second and third parts are integrated in adaptive logarithmic steps.