Analytical derivation of the radial distribution function in spherical dark matter halos

TL;DR

This paper presents an analytical method to derive the velocity distribution function of dark matter in spherical halos, complementing numerical simulations and enhancing understanding of dark matter behavior near Earth.

Contribution

It introduces a differential equation approach to analytically determine the full velocity distribution function in spherical, equilibrated dark matter halos, extending previous work on tangential components.

Findings

Analytical derivation matches numerical simulations below 0.7 v_esc.

Complete velocity distribution can be obtained analytically for spherical halos.

Method enhances understanding of dark matter velocity profiles.

Abstract

The velocity distribution of dark matter near the Earth is important for an accurate analysis of the signals in terrestrial detectors. This distribution is typically extracted from numerical simulations. Here we address the possibility of deriving the velocity distribution function analytically. We derive a differential equation which is a function of radius and the radial component of the velocity. Under various assumptions this can be solved, and we compare the solution with the results from controlled numerical simulations. Our findings complement the previously derived tangential velocity distribution. We hereby demonstrate that the entire distribution function, below 0.7 v_esc, can be derived analytically for spherical and equilibrated dark matter structures.

Figures

Equations

$$\frac{df}{dt}=\frac{\partial f}{\partial t}+\vec{v}\cdot \nabla f-\nabla \Phi \cdot \frac{\partial f}{\partial \vec{v}}=0,$$

$$\int v_r\frac{df}{dt}{d}^{3}v=0,$$

$$\frac{GM(r)}{r}=-\overline{v_r^2}\left(\frac{d\ln (\rho )}{d\ln (r)}+\frac{d\ln (\overline{v_r^2})}{d\ln (r)}+2\beta \right),$$

$$\int \frac{df}{dt}{d}^2{v}_{\theta ,\varphi }=0.$$

$$v_r\frac{\partial f}{\partial r}+\left(\frac{v_{\theta }^2+v_{\varphi }^2}{r}-\frac{\partial \Phi }{\partial r}\right)\frac{\partial f}{\partial v_r}-\frac{v_r{v}_{\theta }}{r}\frac{\partial f}{\partial v_{\theta }}-\frac{v_r{v}_{\varphi }}{r}\frac{\partial f}{\partial v_{\varphi }}=0$$

$$v_r\frac{\partial }{\partial r}\int fd{v}_{\theta ,\varphi }+\frac{1}{r}\frac{\partial }{\partial v_r}\int (v_{\theta }^2+v_{\varphi }^2)fd{v}_{\theta ,\varphi }$$

$$-\frac{\partial \Phi }{\partial r}\frac{\partial }{\partial v_r}\int fd{v}_{\theta ,\varphi }+2\frac{v_r}{r}\int fd{v}_{\theta ,\varphi }$$

$$v_r\frac{\partial F_R}{\partial r}+\frac{1}{r}\frac{\partial }{\partial v_r}\left(⟨v_{\theta }^2{⟩}_r+⟨v_{\varphi }^2{⟩}_r\right)-\frac{d\Phi }{dr}\frac{\partial F_R}{\partial v_r}+2\frac{v_r{F}_R}{r}=0.$$

$$f(r,v_r,v_t)\propto f_t(v_t,k(r))\times R(r)\times f(v_r),$$

$$f(E)=\frac{1}{\sqrt{8}\pi }{\int }_0^E\frac{d^2\rho }{d{\Psi }^2}\frac{d\Psi }{\sqrt{E-\Psi }},$$

$$f(v_t)\propto {\left(1+\frac{v_{\theta }^2+v_{\varphi }^2}{3k^2}\right)}^{-5/2}.$$

$$f(\vec{v},r)\propto {\left(1+\frac{v_{\theta }^2+v_{\varphi }^2}{3k^2}\right)}^{-5/2}R(r)f(v_r).$$

$$F_R\propto {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\left(1+\frac{v_{\theta }^2+v_{\varphi }^2}{3k^2}\right)}^{-5/2}R(r)f(v_r)d{v}_{\theta }d{v}_{\varphi },$$

$$⟨v_{\theta }^2{⟩}_r\propto {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\left(1+\frac{v_{\theta }^2+v_{\varphi }^2}{3k^2}\right)}^{-5/2}{v}_{\theta }^{2}f(v_r)R(r)d{v}_{\theta }d{v}_{\varphi }.$$

$$F_R\propto 2\pi k^2(r)f(v_r)R(r),$$

$$⟨v_{\theta }^2{⟩}_r=⟨v_{\varphi }^2{⟩}_r\propto 6\pi k^4(r)f(v_r)R(r).$$

$$\frac{R^{\prime }(r)}{R(r)}+\frac{1}{k^2(r)}\frac{d{k}^2(r)}{dr}+\frac{2}{r}=\left(\frac{GM(r)}{r^2}-\frac{6k^2(r)}{r}\right)\frac{f^{\prime }(v_r)}{v_{r}f(v_r)}.$$

$$\frac{\frac{R^{\prime }(r)}{R(r)}+\frac{1}{k^2(r)}\frac{d{k}^2(r)}{dr}+\frac{2}{r}}{\frac{GM(r)}{r^2}-\frac{6k^2(r)}{r}}=b_c,$$

$$\frac{f^{\prime }(v_r)}{v_{r}f(v_r)}=b_c,$$

$$R(r)=c_2\exp \left[-\int \left(b_c\left(6k^2(r)-v_c^2\right)+2+\frac{d\ln k^2(r)}{d\ln r}\right)d\ln r\right],$$

$$f(v_r)=c_3{e}^{\frac{1}{2}{b}_c{v}_r^2},$$

Full text

††pagerange: Analytical derivation of the radial distribution function in spherical dark matter halos–References††pubyear: 2012

1 Introduction

The many gravitational observations of dark matter (Clowe et al., 2006; Planck Collaboration et al., 2015) makes us eager to measure the dark matter either in underground observatories or indirectly from the decay or annihilation of dark matter. The analysis of a potential signal in underground observatories requires an assumption of the velocity distribution function of the dark matter near the Earth (Bernabei et al., 2013; Akerib et al., 2016), however, this velocity distribution remains unclear. This problem has in particular been addressed via numerical simulations (for a list of recent references, see Beraldo e Silva et al. (2015); Cerdeno et al. (2016); Kelso et al. (2016); Pillepich et al. (2014)).

A wide range of theoretical ideas have been suggested to understand the dynamical origin of the distribution function e.g. (Le Delliou & Henriksen, 2003; Williams et al., 2004; Lu et al., 2006; Salvador-Solé et al., 2007; Ascasibar et al., 2007; Dalal et al., 2010). Also several papers have attempted statistical mechanical approaches, e.g. Ogorodnikov (1957); Lynden-Bell (1967); Hansen et al. (2005); Hjorth & Williams (2010).

The most dominating guiding principle for theoretical investigations has been the Jeans theorem, which states that any function of the integrals of motion yields a solution to the steady-state Boltzmann equation (Binney & Tremaine, 2008). In other words, any function of energy and angular momentum, $f(E,L)$, will be a solution to the collisionless Boltzmann equation. This has lead to a large number of authors studying exactly such functions (Eddington, 1916; Evans & An, 2005, 2006; An & Evans, 2006; Williams et al., 2014; Wojtak et al., 2008). An alternative approach is to separate the cosmological structure in radial bins, and then consider the distribution of velocities in each bin. This way it was shown that the shape of the distribution function found in numerically simulated strutures is highly non-trivial and differs from a Gaussian (Hansen et al., 2006). In Beraldo e Silva et al. (2015) it was found that a non-separable distribution function was linearly related to the slope of the density profile.

In this paper we will take a new approach, in order to analytically derive the distribution function. We will attempt to derive the shape of the radial components of the distribution function. The idea is to integrate the Boltzmann equation over all the variables we are not interested in, leaving us with a differential equation, Eq. (7), over the variables $r$ and $v_{r}$. We will solve this equation under various assumptions, and we can finally compare the solution with the results of numerically simulated halos of collisionless particles.

Our conclusions are, that when we are given the density profile of a dark matter halo, then we can analytically derive both the radial and tangential distribution function, and hence we know the full distribution function of dark matter particles. The comparisons with numerically simulated structures show that this holds true for velocities below $\sim 0.7$ times the escape velocity.

2 Deriving the radial equation

The Boltzmann equation is the differential equation describing the flow of particles in phase space (Binney & Tremaine, 2008):

[TABLE]

where the function $f(\vec{r},\vec{v},t)$ is the distribution function which gives the density of particles in phase space. In this paper we will discuss the dependence of this function on the radial velocity component $v_{r}$.

Solving this equation is impossible in most cases, and one therefore has to make approximations. The most well-known approach is to get rid of virtually all the velocity information, by integrating the Boltzmann equation over all velocities, for instance

[TABLE]

which leads to one of the Jeans equations

[TABLE]

where $\beta\equiv 1-\frac{\overline{v_{\theta}^{2}}}{\overline{v_{r}^{2}}}$ and we have assumed $\overline{v_{\phi}^{2}}=\overline{v_{\theta}^{2}}$. We have also assumed sphericity, such that $\frac{d\Phi}{dr}=\frac{GM(r)}{r^{2}}$. This is the Jeans equation most commonly used for describing collisionless structures like galaxy clusters and dwarf galaxies.

2.1 Deriving an equation for the radial VDF

The problem with the Jeans equation is that, by averaging over the entire velocity space in calculating the radial velocity dispersion $\sigma_{r}^{2}\equiv\overline{v_{r}^{2}}$, we lose information about the detailed shape of the velocity distribution function (VDF).

We will therefore integrate the Boltzmann equation only over the two angular velocity components instead of all three $v_{i}$’s

[TABLE]

Thereby the radial velocity will remain a free variable. The result, as we shall see, is a differential equation in the variables $v_{r}$ and $r$

We again start with the Boltzmann equation under the assumption of spherical symmetry and staticity:

[TABLE]

Integrating over the two angular velocity components in phase space gives us

[TABLE]

To simplify this equation we define $F_{R}\equiv\int fdv_{\theta,\phi}$, which is the density of particles in a given volume of space with a given radial velocity $v_{r}$. We also define $\langle v_{i}^{2}\rangle_{r}\equiv\int fv_{i}^{2}dv_{\theta,\phi}$, which is the weighted sum of squared i-velocity components in the phase plane ($v_{\theta}$ , $v_{\phi}$), with f acting as weight. If divided by $F_{R}$, this would be completely analogous to the averages $\overline{v_{i}^{2}}$ in the Jeans equation. Using this, our equation becomes

[TABLE]

This is a differential equation in two variables $v_{r}$ and $r$, and with four unknown functions, $F_{R}$, $\langle v_{\theta}^{2}\rangle_{r}$, $\langle v_{\phi}^{2}\rangle_{r}$ and $\Phi(r)$. As it is, the equation has far too many unknowns to be solved directly.

In order to proceed to solve this differential equation, we will assume that the distribution function can be separated into the components

[TABLE]

where $k(r)$ is an $r$-dependent normalisation parameter. Next we need to calculate approximate expressions for $F_{R}$ and $\langle v_{i}^{2}\rangle_{r}$ as functionals of $R(r)f(v_{r})$. This will be done using the shape of the tangential velocity distribution derived in Hansen et al. (2005), which was demonstrated to be in good agreement with numerical simulations (Hansen & Sparre, 2012). In this way, the number of unknown functions in the above equation can be reduced to two, namely $f_{r}(v_{r})$ and $R(r)$, and the equation can then be solved using regular separation of variables.

2.2 The tangential VDF

Historically, attempting to derive a general expression for VDFs in static, spherical systems, Eddington realised that in order for the system to be static, the inflow from an outer spherical shell must equal the flow out of the smaller inner radial bin. The flux of particles as a function of average velocity clearly also depends on the density in the various bins, and therefore the density profile of the entire structure. Using the requirement of staticity, Eddington was able to relate the velocity distribution function to the density profile of the spherically symmetric structure (Eddington, 1916)

[TABLE]

where $-\Psi(r)=\Phi(r)$ is the potential, $\rho(r)$ is the density profile, and $E=\Psi-\frac{v^{2}}{2}$ is the relative energy. We have also assumed $\beta=0$. This method for finding the velocity distribution function is know as the Eddington inversion method (Binney & Tremaine, 2008). One unfortunate problem with this method is, that the resulting distribution functions do not agree with the ones observed in numerical simulations.

It has been suggested that the tangential component of the distribution function should be particularly simple (Hansen et al., 2005). For sufficiently small velocities, the tangential velocity components only cause the particles to move around within the same radial bins. Given spherical symmetry, the potential and density at constant $r$ are themselves constant. In this vastly simplified case, the Eddington inversion method yields

[TABLE]

This function is known as a Tsallis distribution and has wide applications in non-extensive statistical mechanics, where functions of this type tend to maximise the generalised entropy. Clearly, the resulting analytical distribution function has a long tail of high energy particles, which cannot exist in finite structures, where the maximum velocity at any radius is the escape velocity, $v_{\rm esc}=\sqrt{2\left|\Phi\right|}$. This immediately tell us, that this function cannot provide a good fit at high velocities. Surprisingly enough, this function gives an excellent fit for all velocities smaller than $0.7\,v_{\rm esc}$ at all radii. This holds for structures formed in controlled numerical simulations (spherical collapse, head-on collisions, and various perturbed structures) and also in cosmological simulated structures inside a density slope around $\gamma=-2.5$.

3 Solving the radial equation

As clarified above, we need to make two assumptions to solve the equation for the radial distribution function, Eq. (7), namely that the full distribution function approximately is separable, $f(r,v_{t},v_{r})=R(r)\times f_{t}(v_{t})\times f_{r}(v_{r})$, and that the shape of the tangential distribution funtion is known. Substituting the tangential VDF in Eq. (10), we get

[TABLE]

For ease of notation we denote the $v_{r}$-component of the function as $f(v_{r})$ as opposed to the full function $f(v_{r},r)$. We then calculate $\langle v_{i}^{2}\rangle_{r}$ and $F_{R}$

[TABLE]

[TABLE]

These integrals give us

[TABLE]

[TABLE]

Plugging these into the radial velocity equation, we get

[TABLE]

Separation implies that we have the two ordinary differential equations

[TABLE]

[TABLE]

with the solutions

[TABLE]

[TABLE]

where we have that $f(v_{r},r)\propto R(r)f(v_{r})$.

4 Comparison with numerical simulations

In order to test the accuracy of our solution, we compare with the velocity distributions extracted from a range of controlled numerical simulations. The reasons for using controlled (instead of cosmological) simulations is that we thereby can assure that equilibrium is achieved. When a structure has reached equilibrium we divide it in radial bins, and then we plot the resulting radial and tangential velocity distribution in each radial bin.

We consider 3 different perturbation schemes. The first follows a scheme to resemble the effect of mergers by repeatedly changing the energies of individual particles by changing the depth of the potential. Practically we vary the gravitational potential, where the value of G is changed by a factor of 1.1 and 0.9 from its normal value, each for a few dynamical times. Thereby the particles are accelerated/deccelerated while the system is allowed to relax through phase-mixing (Sparre & Hansen, 2012).

Only the trustworthy region is considered, excluding the central resolution-limited region with radii smaller than 5 times the softening length. The outer non-fully-equilibrated region is identified by considering regions where $\gamma=d{\rm ln}\rho/d{\rm ln}r$ starts varying between individual perturbations. Effectively this removes regions outside radii where the density slope is -3. The issue of trustworthy regions was analysed in detail in Sparre & Hansen (2012).

In Figure 1 we show the distribution of velocities in 3 different radial bins, where the slope of density is about $\gamma=d{\rm ln}\rho/d{\rm ln}r=-1.8,-2.3,-3$ The bins are shifted vertically to improve readability. The blue stars are the radial velocities, and the red solid line is the Gaussian in Eq. (20). In comparison we also plot the tangential velocity distribution (green diamonds) which is well fitted with the shape in Eq. (10).

Figure 2 shows the same distribution in lin-log space in order to make the tails in the distribution much more visible.

The second perturbation scheme is a standard cold collapse: a spherical distribution of particles are placed at rest, and the violent relaxation during the collapse leads to a new equilibrated structure. The resulting velocity distributions at 3 different radii are shown in Figure 3. The cold collapse is a single major perturbation (as opposed to the repeated changing of G as described above), and it is therefore less likely that a perfectly beautiful distribution can be obtained. Infact, in Figure 3 we see that a Gaussian is not a perfect fit, since variations are even visible by eye on the lin-log figure.

It is important to keep in mind that in principle there are no free parameters in the fits. The width of the distribution function is given by the radial dispersion, and the height of the curve is given by the number of particles in the radial bin. However, since the curves are not fitted perfectly (the tail is wrong), we allow both the dispersion and normalization as free fitting parameters, and the resulting fits are a few percent different from the true values from the simulations.

In order to address the question if a Gaussian is really the best fit, we repeat the data in the lowest lines of Figure 3 in a lin-lin space, where we show two best fit lines: the (blue) Gaussian is visibly seen to provide a better fit than the (green) Tsallis shape: the Tsallis shape is too low around $v\sim\sigma$, and the Tsallis is too high at velocities above $v\approx 2\sigma$. For this radial bin the reduced chi-squares are about 2 and 6 respectively, when fitting velocities inside 1.8 times the dispersion (corresponding to 1.2 on the x-axis). Another way of quantizing the preference for the Gaussian is that the resulting error-bars on the fitting parameters (dispersion and magnitude) are approximately a factor of two smaller when fitted to a Gaussian than when fitted to a Tsallis shape. The statistical preference for the Gaussian is even bigger when fitting to higher velocities (because the Tsallis has a much longer tail), however, that may not be entirely reliable given the very few particles in the high-velocity tail. With a few percent error-bars (which requires about half a million particles in the structure) such Poisson error-bars are easily achieved. For the bins at largest radii there is such a small difference between the radial and tangential distributions, that with only half million particles in the entire structure, the statistical error-bars would allow a fit to both a Gaussian and a Tsallis shape, when one allows both dispersion and magnitude as free parameters, unless one includes the highest-velocity particles (in which case the Gaussian provides a statistically better fit).

The last perturbation scheme is a repeated pattern of kick-flow, in the sense that first the particle energies are perturbed (to resemble the effect of mergers) and secondly the system is allow to relax. Each energy kick is arranged to conserve energy in each radial bin (Hansen et al., 2010). Between each kick the system is allowed to relax through phase-mixing.

In Figure 5 we again plot the velocity distribution from 3 radial bins, and the lower bins correspond to a region near the centre, but well outside a region of 5 times the softening length.

These 3 very different perturbation schemes allow us to draw the same conclusion: all the radial distribution function are well fitted with the shape in Eq. (20) for velocities below approximately $0.7v_{\rm esc}$. This complements the known fact that the tangential VDF is well fit in the same region with the form in Eq. (10). One should keep in mind, that the majority of particles are having energies below $0.7v_{\rm esc}$, and that the departure from the Gaussian shape is not visible when plotting in linear-linear.

The reason why the theoretical predictions and the numerically simulated structures do not agree at high velocities is unknown to us. In the derivation we have assumed that the tangential velocities follow the shape in Eq. (10), which is known to be violated at high velocities. Assuming a cut-off in the tangential shape does not change the shape of the resulting radial distribution, but only a marginally different value for the radial dispersion. We therefore believe the origin should lie in the assumption of separability of the full distribution function, see Eq. (11). The full differential equation (7) could in principle be solved without assuming separability, however, the resulting solution is much less transparent and we will therefore not pursue this path here.

5 Conclusion

We derive the equation involving only the radius and the radial component of the velocity distribution function by integrating the Boltzmann equation over both the tangential velocities and the spatial angles. We solve this equation under the assumption that the distribution function is separable in the radial and tangential velocities. These solutions indicate that the radial distribution function is close to a Gaussian shape, Eq. (20), whereas the tangential is close to a Tsallis shape, Eq. (10).

We compare the solution of the radial distribution function with the results of numerical simulations, and we find that for velocities smaller than approximately 0.7 times the escape velocity, the solution fits rather well. In principle this is a prediction with zero free parameters, and hence the agreement with simulated data is not entirely trivial. Thus, for small velocities (which is the dominant component) the full velocity distribution function can be derived analytically.

Acknowledgement

It is a pleasure to thank Martin Sparre for providing data for the figures. We thank the anonymous referee for suggestions which improved the paper. We thanks Gary Mamon for constructive comments. This project is partially funded by the Danish council for independent research, under the project “Fundamentals of Dark Matter Structures”, DFF – 6108-00470.