Analytic expressions for the Dekel+ profile

This python program enables to compute the analytical expressions of the Dekel+ dark matter density profile. It is based on the following article:

The Dekel+ profile: a mass-dependent dark-matter density profile with flexible inner slope and analytic potential, velocity dispersion, and lensing properties [PDF]

Quantities expressed analytically in this program:

• The density
• The average density
• The enclosed mass
• The orbital velocity
• The logarithmic density slope
• The gravitational potential
• The velocity dispersion

For the velocity dispersion, it provides both an expression in terms of incomplete beta functions (Freundlich et al. 2020b, Eq. 19) and two expressions in terms of finite sums (Zhao 1996, Eq. 19, and Freundlich et al. 2020a, Eqs. B8 and B10). Some of the quantities are also defined for other profiles (NFW, Einasto, Di Cintio+, generalised NFW with flexible inner slope).

Quantities integrated numerically in this program:

• The distribution function
• The projected surface density
• The projected mass
• The lensing convergence
• The lensing shear

Except the distribution function, these quantities can be expressed analytically in terms of Fox H functions and series expansions (Freundlich et al. 2020b, Section 2.3 and Appendix E).

The program also implements the mass-dependent prescriptions for the Dekel+ and the Di Cintio+ profiles, and provides functions to convert the shape parameters (a,c) into an inner slope s1 and a concentration c2 and vice-versa.

Example:

import numpy as np
import matplotlib.pyplot as plt
import Dekel_profile as prf

s1,c1=(1.,15.)
s2,c2=(1.,5.)
s3,c3=(0.,15.) 
s4,c4=(0.,5.)

rvir=1.
mvir=1.
rhovir=3.*mvir/(4.*np.pi*rvir**3)

p1=prf.params_from_s1c2(s1,c1,0.01*rvir,rvir,mvir)
p2=prf.params_from_s1c2(s2,c2,0.01*rvir,rvir,mvir)
p3=prf.params_from_s1c2(s3,c3,0.01*rvir,rvir,mvir)
p4=prf.params_from_s1c2(s4,c4,0.01*rvir,rvir,mvir)

color1='r'
color2='orangered'
color3='b'
color4='navy'
dashes2=(9,6)
dashes3=(5,4)
dashes4=(2,2)

x=np.linspace(-2,0,100)
y1=np.log10(prf.rho_function(10**x,p1,'dekel')/rhovir)
y2=np.log10(prf.rho_function(10**x,p2,'dekel')/rhovir)
y3=np.log10(prf.rho_function(10**x,p3,'dekel')/rhovir)
y4=np.log10(prf.rho_function(10**x,p4,'dekel')/rhovir)

plt.figure()
plt.plot(x,y1,color=color1,lw=2,label=r'$s_1=%.0f$, $c_{2}=%.0f$'%(s1,c1))
plt.plot(x,y2,color=color2,ls='--',dashes=dashes2,lw=2,label=r'$s_1=%.0f$, $c_{2}=%.0f$'%(s2,c2))
plt.plot(x,y3,color=color3,ls='--',dashes=dashes3,lw=2,label=r'$s_1=%.0f$, $c_{2}=%.0f$'%(s3,c3))
plt.plot(x,y4,color=color4,ls='--',dashes=dashes4,lw=2,label=r'$s_1=%.0f$, $c_{2}=%.0f$'%(s4,c4))

plt.legend(frameon=False,loc='lower left')
plt.xlabel(r'$\log(r/R_{\rm vir})$')
plt.ylabel(r'$\log(\rho/\overline{\rho_{\rm vir}})$')
plt.axis([-2,0,-1,4])