Analytic expressions for the Dekel-Zhao profile

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

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

Quantities expressed analytically:

• 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:

• 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-Zhao 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.

Fitting_functions.py provides simple functions to fit simulated density profiles with the Dekel-Zhao profile. The fits are constrained to avoid positive inner slopes.

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])