6.8. Jeans analysis (jeans_analysis.h)

The Jeans equation relates the dark matter content of an object (potential well) to its light content (dynamics in the potential). The main steps and ingredients to derive the Jeans equation are recalled below. We recall that the user has to provide the data for the velocity dispersion data (\(\sigma_p^2(R)\)) and the parameters for the light profile, \(I(r)\), Table 6.7, and the anisotropy profile, \(\beta_{\rm ani}(r)\), Table 6.8. Jeans analysis calculations are implemented in the library jeans_analysis.h in CLUMPY.

A reference book for stellar dynamics is Binney and Tremaine's Galactic Dynamics, where the Jeans equation is described — see also Battaglia et al. (2013) and Courteau et al. (2014) for reviews on dynamics in the context of dSphs and the Galaxy respectively.


All the expressions below are valid for the spherical-halo case only.

6.8.1. Jeans equation

Assuming virial equilibrium and a non-rotating system, the Jeans equation is obtained by

  1. Using the collisionless Boltzmann equation for the 6D space-phase of stellar density \(f(t,\vec{x},\vec{v})\);
  2. Assuming steady-state and spherical symmetry;
  3. Extracting and combining the first two moments \(\int\!\!\int\!\!\int f \,v_i^n\, d^3\vec{v}\) for \(n=0\) (density of stars) and \(n=1\) (mean stellar velocity along the \(i\)-th component).

We can then link the enclosed mass \(M(r)\) in radius \(r\) and moments of the stellar distribution:

\[\frac{1}{\nu}\frac{\mathrm{d}}{\mathrm{d}r}(\nu \bar{v_r^2})+2\frac{\beta_{\rm ani}(r)\bar{v_r^2}}{r}=-\frac{GM(r)}{r^2} \quad\quad {\rm [Jeans~Equation]},\]


  • \(\nu(r)\), also denoted \(\rho(r)\): spatial density (spherical symmetry) of the light profile;
  • \(\bar{v_r^2}(r)\): radial velocity dispersion of baryons (e.g. stars in dSphs, or galaxies in galaxy clusters);
  • \(\beta_{\rm ani}(r)\equiv 1-\bar{v_{\theta}^2}(r)/\bar{v_r^2}(r)\): orbital anisotropy of the stellar component (depending on the radial \(v_r^2\) and tangential \(v_{\theta}^2\) velocity dispersion).

Note that higher orders of the Jeans equation can be derived (integrating on higher moments of the distribution function \(f\)), but it is not considered in this version of CLUMPY.

6.8.2. (De-)projection & Abel

In practice, we only have access to projected quantities on the sky, as illustrated in Fig. 6.15.


Fig. 6.15 Framework for Abel transform (projection).

For spherically symmetric systems, projection and de-projection are given by the Abel and inverse Abel transform (the quantity \(R\) is the projected radius, and \(y=\sqrt{r^2-R^2}\), see Fig. 6.15):

  • Projection: \(f(r)\rightarrow F(R)\)

    \[F(R) = 2 \int_R^\infty \frac{f(r)\,r\,\mathrm{d}r}{\sqrt{r^2-R^2}} = 2\int_0^\infty f(y) \, \mathrm{d}y\]


    \[y=\sqrt{r^2-R^2}\rightarrow \mathrm{d}y=\frac{r\,\mathrm{d}r}{\sqrt{r^2-R^2}}.\]
  • De-projection: \(F(R)\rightarrow f(r)~~~~~~~~[\) valid if \(f(r)\) drops to zero more quickly than \(1/r]\)

    \[f(r) = -\frac{1}{\pi} \int_r^\infty \frac{\mathrm{d}F}{\mathrm{d}R}(R) \times \frac{\mathrm{d}R}{\sqrt{R^2-r^2}} = -\frac{1}{\pi} \int_0^\infty \frac{\mathrm{d}F}{\mathrm{d}R}(R=\sqrt{(Y^2+r^2)}) \times \frac{\mathrm{d}Y}{R}\]


    \[Y=\sqrt{R^2-r^2}\rightarrow \mathrm{d}Y=\frac{R\,\mathrm{d}R}{\sqrt{R^2-r^2}}.\]

In practice, because of the singularity of the integrand for \(r=R\), the integration better behaves numerically (faster convergence) if it is performed over the variable \(y\) (for projection, see Fig. 6.15) or \(Y\) (for de-projection): this is what we use in CLUMPY.

6.8.3. Formal solutions

We can write the formal solution of the Jeans equation thanks to \(\nu(r)\bar{v_r^2}(r)\) and \(I(R)\sigma_p^2(R)\).

Solving for deprojected (3D): \(\nu(r)\bar{v_r^2}(r)\)

The above Jeans equation is a first order linear non-homogeneous differential equation with non constant coefficients. Its general solution is given for instance here, and for the Jeans equation this gives:

(6.12)\[\nu(r)\bar{v_r^2}(r) = \frac{1}{f(r)}\times \int_r^\infty f(s) \nu(s) \frac{GM(s)}{s^2}\, \mathrm{d}s\]


(6.13)\[f(r) = f_{r_1} \exp\left[\int_{r_1}^r \frac{2}{t}\beta_{\rm ani}(t)\, \mathrm{d}t \right].\]

The integration on the lower boundary \(r_1\) leads to a normalisation factor that disappears in the solution \(\nu(r)\bar{v_r^2}(r)\), because the function \(f\) appears both in the numerator and denominator.

Solving for projected (2D): \(I(R)\sigma_p^2(R)\)

To compare the velocity dispersion to data, one needs to project the above equation (only line-of-sight projected velocity dispersion are available to the observer). Using the Abel transform, we have:

(6.14)\[\begin{split}I(R)\sigma_p^2(R)&=2\displaystyle \int_{R}^{\infty}\biggl (1-\beta_{\rm ani}(r)\frac{R^2}{r^2}\biggr ) \frac{\nu(r) \,\bar{v_r^2}(r)\,r}{\sqrt{r^2-R^2}}\mathrm{d}r\\ &= 2\displaystyle \int_{0}^{\infty}\biggl (1-\beta_{\rm ani}(r)\frac{R^2}{r^2}\biggr ) \nu(r) \,\bar{v_r^2}(r)\mathrm{d}y .\end{split}\]

using the change of variable \(y=\sqrt{r^2-R^2}\), which gives \(dy=r\,dr/\sqrt{r^2-R^2}\).

Depending on the spatial distribution for the anisotropy \(\beta_{\rm ani}(r)\), the above equations tell us that three integrations are required (which is quite time consumming): one to get \(f(r)\), one to get \(\nu(r)\bar{v_r^2}(r)\), and one last to get \(I(R)\sigma_p^2(R)\). In practice, there are several special cases for which:

  • \(f(r)\) has an analytical form, reducing the number of integration to two — one for \(\nu(r)\bar{v_r^2}(r)\) and one for \(I(R)\sigma_p^2(R)\);

  • relying on the so-called Kernel \({\cal K}\) (Mamon and Lokas, 2005), the solution can be rewritten as a single integration:

    (6.15)\[I(R)\sigma_p^2(R)=2 G\times \int_{R}^{\infty} {\cal K}\left(\frac{r}{R},\,\frac{r_a}{R}\right) \nu(r) M(r) \frac{\mathrm{d}r}{r}.\]

We present the anisotropy profiles used and the corresponding solutions for \(f(r)\) and/or \(I(R)\sigma_p^2(R)\) below.

6.8.4. Inputs in CLUMPY

From the above equations, four ingredients are required to perform the Jeans analysis. They are described below along with the parametrisations available in CLUMPY.

1. Dark matter density profile \(\rho_{\rm DM}\) (and total mass)

Any of the profiles discussed in Section 6.5 can be used. Note that the baryonic mass should enter the total mass in the calculation in principle, but it is neglected in dSphs, where the mass-to-light ratio is of the order of hundreds.

2. Surface brightness profiles \(I(R)\equiv\Sigma(R)\) and its deprojection \(\nu(r)\equiv\rho(r)\)

Several light profiles are routinely used in the literature (Plummer, King, Sérsic…). They are implemented in CLUMPY. These profiles have from 2 to 5 free parameters (normalisation, scale radius, plus some structural parameters). A difficulty is to decide whether the analytical form used for the profile should be assumed for the surface brightness (projected 2D profile) or for the density distribution (un-projected 3D profile). In the first case, de-projection is required (inverse Abel transform), whereas in the second case, projection is (Abel transform). The extension 2D or 3D is added at the end of the profile keyword to explicitly indicate which case it is, as listed in Table 6.7.

Table 6.7 Global enumerators available in gENUM_LIGHTPROFILE (for exponential profiles, \(K_0\) and \(K_1\) are respectively the modified Bessel function of the second kind of order 0 and 1).
Keyword Surface brightness \(\Sigma(R)\equiv I(R)\)   Density profile \(\rho(r) \equiv\nu(r)\) # free param. References
kEXP2D \(\mathbf{\Sigma_0 \times \exp\left(-\frac{R}{r_c}\right)}\) \(\rightarrow\) \(\frac{\Sigma_0}{\pi r_c}\times K_0\left(\frac{r}{r_c}\right)\) 2 Evans, An, and Walker (2009)
kEXP3D \(2\rho_0 \,R\times K_1\left(\frac{R}{r_c}\right)\) \(\leftarrow\) \(\mathbf{\rho_0 \times \exp\left(-\frac{r}{r_c}\right)}\) 2 Evans, An, and Walker (2009)
kKING2D \(\mathbf{\Sigma_0\times\bigg[\left(1+\frac{R^2}{r_c^2}\right)^{-1/2}}\) \(\mathbf{-\left(1+\frac{r_{\rm lim}^2}{r_c^2}\right)^{-1/2}\bigg]^2}\) \(\rightarrow\) \(\frac{\Sigma_0}{\pi r_c}\times\frac{\cos^{-1}(z)/z-\sqrt{1-z^2}}{z^2(1+r_{\rm lim}^2/r_c^2)^{3/2}}\) with \(z^2 = \frac{1+r^2/r_c^2}{1+r_{\rm lim}^2/r_c^2}\) 3 King (1962), Strigari et al. (2008)
kPLUMMER2D \(\mathbf{\frac{\Sigma_0}{\pi r_c^2}\times \left(1+\frac{R^2}{r_c^2}\right)^{-2}}\) \(\rightarrow\) \(\frac{3\Sigma_0}{4\pi r_c^3}\times \left(1+\frac{r^2}{r_c^2}\right)^{-5/2}\) 2 Plummer (1911), Evans et al. (2009)
kSERSIC2D \(\mathbf{\Sigma_0\!\times\!\exp\left\{\!-b_n\!\left[\left(\frac{R}{r_c}\right)^{\frac{1}{n}}\!-\!\!1\right]\right\}}\) with \(b_n=2n-1/3+0.009876/n\) \(\rightarrow\) \(-\frac{1}{\pi} \int_r^\infty \frac{{\rm d}\Sigma(R)}{{\rm d}R}\times \frac{{\rm d}R}{\sqrt{R^2-r^2}}\) 3 Sérsic (1968), Prugniel & Simien (1997), Graham and Driver (2005), Merritt et al. (2006)
kZHAO3D \(2 \int_R^\infty \rho(r)\,r\times\frac{{\rm d}r}{\sqrt{r^2-R^2}}\) \(\leftarrow\) \(\mathbf{\rho_0\times\frac{(r/r_s)^{-\gamma}}{\left[1+\left(\frac{r}{r_s}\right)^\alpha\right]^{(\beta-\gamma)/\alpha}}}\) 5 Hernquist (1990), Zhao (1996)


3. Anisotropy profile \(\beta_{\rm ani}(r)\)

We recall that the anisotropy profile is given by a combination of the radial and tangential velocity dispersion (and that \(\beta_{\rm ani}(r)<1\) for any \(r\)):

\[\beta_{\rm ani}(r)\equiv 1-\frac{\bar{v_{\theta}^2}(r)}{\bar{v_r^2}(r)}\,.\]

Due to the lack of information on the anisotropy profile, the first profiles discussed in the literature were based on analytical studies to build dynamical models (in spherical symmetry) having realistic anisotropy profile. Very few of these models are known: the ones considered are either completely isotropic (\(\beta_{\rm ani}=0\)), or have a constant anisotropy or an anisotropy profile of Osipkov-Merritt type (isotropic in the centre and completely radially anisotropic at large radii). More recently, indications of radial anisotropy in the outer regions of dark matter haloes have been obtained from numerical simulations. In the inner region, dynamical formation and evolution processes can lead to a strong anisotropy. To better describe these profiles, Baes and van Hese (2007) introduced a technique to construct dynamical models with an arbitrary density and an arbitrary anisotropy profile.

We gather in Table 6.8 these various profile and the function \(f(r)\) and kernel \({\cal K}(u,u_a)\) required to solve Jeans Eqs (6.12) and (6.15), i.e.

\[\nu(r)\bar{v_r^2}(r) = \frac{1}{f(r)}\times \int_r^\infty f(s) \nu(s) \frac{GM(s)}{s^2}\, \mathrm{d}s\]


\[I(R)\sigma_p^2(R)=2 G\times \int_{R}^{\infty} {\cal K}\left(\frac{r}{R},\,\frac{r_a}{R}\right) \nu(r) M(r) \frac{\mathrm{d}r}{r},\]

Table 6.8 Global enumerators available in gNAMES_ANISOTROPYPROFILE (kOSIPKOV is a special case of kBAES with \(\beta_0=0\), \(\beta_\infty=1\), and \(\eta=2\)).
Keyword Anisotropy \(\beta_{\rm anis}(r)\) \(f(r)\) to solve Eq. (6.12) \(=\nu(r)\bar{v_r^2}(r)\) Kernel \({\cal K}(u,u_a)\) to solve Eq. (6.14) \(=I(R)\sigma_p^2(R)\) Refs (for profile and/or solution)
kCONSTANT \(\beta_0\) \(r^{2\beta_0}\) \(\frac{\sqrt{1-u^{-2}}}{1-2\beta_0}+\frac{\sqrt{\pi}}{2}\frac{\Gamma(\beta_0-1/2)}{\Gamma(\beta_0)}u^{2\beta_0-1}\) \(\left(\frac{3}{2}-\beta_0\right)\left[1-I\left(\frac{1}{u^2},\beta_0+\frac{1}{2},\frac{1}{2}\right)\right]\) with \(I(x,a,b)=\) incomplete Beta function Mamon and Lokas (2005,2006)
kBAES \(\frac{\beta_0 + \beta_\infty (r/r_a)^\eta}{1+(r/r_a)^\eta}\) \(r^{2\beta_0}\left[1+\left(\frac{r}{r_a}\right)^\eta\right]^{2(\beta_\infty-\beta_0)/\eta}\) No analytical Kernel Baes & van Hese (2007)
kOSIPKOV \(\frac{r^2}{r^2+r_a^2}\) (special case of kBAES) \(\frac{r_a^2+r^2}{r_a^2}\) \(\displaystyle\frac{(u^2+u_a^2)(u_a+1/2)}{u(u_a^2+1)^{3/2}}\tan^{-1}\sqrt{\frac{u^2-1}{u_a^2+1}}\) \(-\displaystyle\frac{\sqrt{1-u^{-2}}}{2(u_a^2+1)}\) Osipkov (1979), Merritt (1985), Mamon and Lokas (2005, 2006)

6.8.5. Sketch of Jeans run

Fig. 6.16 gives a sketch on how to proceed with a Jeans analysis, with the inputs and outputs.


Fig. 6.16 Diagram of the Jeans analysis with CLUMPY. From a kinematic data file and a parameter file describing the free parameters, a statistical Jeans analysis can be run with \({\tt clumpy\_jeansMCMC}\) or \({\tt clumpy\_jeansChi2}\). A statistical output file is created, from which estimators of different quantities of interest (i.e., J-factors) can be obtained with \({\tt clumpy -s}\).

6.8.6. Log-likelihood (for fit)

To obtain the best-fit DM halo parameters solving the Jeans equation (see Section 8.1), we can rely on several likelihood functions. The details can be found in Bonnivard et al. (2015).

Binned and unbinned analysis

For the binned analysis, the velocity dispersion profiles \(\sigma_{\rm obs}(R)\) are built from the individual stellar velocities and the likelihood function is:

\[\mathcal{L}^{\rm bin}= \prod_{i=1}^{N_{\rm bins}} \frac{(2\pi)^{-1/2}}{\Delta{\sigma_i}(R_i)} \exp\biggl [-\frac{1}{2}\biggl (\frac{\sigma_{\rm obs}(R_{i}) -\sigma_{p}(R_i)}{\Delta {\sigma}_i(R_i)}\biggr )^{2}\biggr ].\]

The quantity \(\Delta {\sigma_{\rm i}}(R_i)\) is the error on the velocity dispersion at the radius \(R_i\).

For the unbinned analysis, we assume that the distribution of line-of-sight stellar velocities is Gaussian, centred on the mean stellar velocity \(\bar{v}\):

\[\mathcal{L}^{\rm unbin}= \prod_{i=1}^{N_{\rm stars}} \frac{(2\pi)^{-1/2}}{\sqrt{\sigma_{p}^2(R_i)+\Delta_{v_{ i}}^{2}}} \exp\biggl [-\frac{1}{2}\biggl (\frac{(v_{\rm i} -\bar{v})^{2}}{\sigma_{p}^2(R_i)+\Delta_{v_{i}}^{2}}\biggr )\biggr ],\]

where the dispersion of velocities at radius \(R_i\) of the \(i\)-th star comes from both the intrinsic dispersion \(\sigma_{p}(R_i)\) and the measurement uncertainty \(\Delta_{v_{i}}\)

Analysis with or without membership probabilities

Kinematic samples are often contaminated by interlopers from the Milky Way (MW) foreground stars. Different methods can be used in order to remove those contaminants, based e.g. on sigma-clipping, virial theorem, or expectation maximisation (EM) algorithms. The latter allows in particular the estimation of membership probabilities for each star of the object. These probabilities can be used as weights when building the velocity dispersion profile \(\sigma_{\rm obs}(R)\) in the binned case, or used directly in the unbinned likelihood, where the above equation becomes

\[\mathcal{L}^{\rm unbin}_{W}= \prod_{i=1}^{N_{\rm stars}} \left(\frac{(2\pi)^{-1/2}}{\sqrt{\sigma_{p}^2(R_i)+\Delta_{v_{i}}^{2}}} \exp\biggl [-\frac{1}{2}\biggl (\frac{(v_{\rm i} -\bar{v})^{2}}{\sigma_{p}^2(R_i)+\Delta_{v_{i}}^{2}} \biggr )\biggr ] \right)^{P_{i}},\]

with \(P_{i}\) the membership probability of the \(i\)-th star.