Some example usage of how to build up a dataframe of galaxy cluster properties, including NFW halo profiles. Each cluster is treated as an individual, meaning we track its individual mass and redshift, and other properties. This is useful for fitting a stacked weak lensing profile, for example, where you want to avoid fitting a single average cluster mass. See the project documentation for more details.
In [1]:
from __future__ import print_function
import numpy as np
from astropy import units
from matplotlib import pyplot as plt
%matplotlib inline
from clusterlensing import ClusterEnsemble
In [2]:
z = [0.1,0.2,0.3]
c = ClusterEnsemble(z)
c.describe
Out[2]:
In [3]:
c.show()
In [4]:
n200 = np.ones(3)*20.
c.n200 = n200
c.show()
Notice that astropy units are present for the appropriate columns.
In [5]:
print('z: \t', c.z)
print('n200: \t', c.n200)
print('r200: \t', c.r200)
print('m200: \t', c.m200)
print('c200: \t', c.c200)
print('rs: \t', c.rs)
In [6]:
c.r200.value
Out[6]:
In [7]:
c.dataframe
Out[7]:
In [8]:
c.z = np.array([0.4,0.5,0.6])
c.dataframe
Out[8]:
In [9]:
c.m200 = [3e13,2e14,1e15]
c.show()
In [10]:
c.n200 = [20,30,40]
c.show()
In [11]:
c.massrich_slope = 1.5
c.show()
In [12]:
# Show basic table without Pandas formatting
c.show(notebook = False)
In [13]:
rmin, rmax = 0.1, 5. # Mpc
nbins = 20
rbins = np.logspace(np.log10(rmin), np.log10(rmax), num = nbins)
print('rbins range from', rbins.min(), 'to', rbins.max(), 'Mpc')
In [14]:
c.calc_nfw(rbins) # calculate the profiles
sigma = c.sigma_nfw # access the profiles
deltasigma = c.deltasigma_nfw
In [15]:
sigma[0]
Out[15]:
In [16]:
fig = plt.figure(figsize=(12,5))
fig.suptitle('Centered NFW Cluster Profiles', size=30)
first = fig.add_subplot(1,2,1)
second = fig.add_subplot(1,2,2)
for rich, profile in zip(c.n200,deltasigma):
first.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
first.set_xscale('log')
first.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
first.set_ylabel('$\Delta\Sigma(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$',
fontsize=20)
first.set_xlim(rbins.min(), rbins.max())
first.legend(fontsize=20)
for rich, profile in zip(c.n200,sigma):
second.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
second.set_xscale('log')
second.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
second.set_ylabel('$\Sigma(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$',
fontsize=20)
second.set_xlim(0.05, 1.)
second.set_xlim(rbins.min(), rbins.max())
second.legend(fontsize=20)
fig.tight_layout()
plt.subplots_adjust(top=0.88)
When the true underlying dark matter distribution is offset from the assumed cluster "center" (such as a BCG or some other center proxy), the weak lensing profiles measured around the assumed centers will be different than for the perfectly centered case. One way to account for this is to describe the cluster centroid offsets as a Gaussian distribution around the true centers. We say the probability of an offset is given by
$P(R_\mathrm{off}) = \frac{R_\mathrm{off}}{\sigma_\mathrm{off}^2}e^{-\frac{1}{2}\left(\frac{R_\mathrm{off}}{\sigma_\mathrm{off}}\right)^2}$,
which is parameterized by the width of the 2D offset distribution $\sigma_\mathrm{off} = \sqrt{\sigma_x^2 + \sigma_y^2}$. Then the measured surface mass density is given by
$\Sigma^\mathrm{sm}(R) = \int_0^\infty \Sigma(R | R_\mathrm{off})\ P(R_\mathrm{off})\ \mathrm{d}R_\mathrm{off}$,
where
$\Sigma(R | R_\mathrm{off}) = \frac{1}{2\pi} \int_0^{2\pi} \Sigma(r')\ \mathrm{d}\theta$,
and
$r' = \sqrt{R^2 + R_\mathrm{off}^2 - 2 R R_\mathrm{off} \cos{\theta}}$.
More details on the cluster miscentering problem can be found in Ford et al 2015, George et al 2012, and Johnston et al 2007.
To calculate the miscentered profiles, simply create an array of offsets in units of Mpc, and pass it to the calc_nfw method:
In [17]:
offsets = np.array([0.1,0.1,0.1]) #same length as number of clusters
c.calc_nfw(rbins, offsets=offsets)
deltasigma_off = c.deltasigma_nfw
sigma_off = c.sigma_nfw
In [18]:
fig = plt.figure(figsize=(12,5))
fig.suptitle('Miscentered NFW Cluster Profiles', size=30)
first = fig.add_subplot(1,2,1)
second = fig.add_subplot(1,2,2)
for rich, profile in zip(c.n200,deltasigma_off):
first.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
first.set_xscale('log')
first.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
ytitle = '$\Delta\Sigma^\mathrm{sm}(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$'
first.set_ylabel(ytitle, fontsize=20)
first.set_xlim(rbins.min(), rbins.max())
first.legend(fontsize=20)
for rich, profile in zip(c.n200,sigma_off):
second.plot(rbins, profile, label='$N_{200}=$ '+str(rich))
second.set_xscale('log')
second.set_xlabel('$r\ [\mathrm{Mpc}]$', fontsize=20)
ytitle = '$\Sigma^\mathrm{sm}(r)\ [\mathrm{M}_\mathrm{sun}/\mathrm{pc}^2]$'
second.set_ylabel(ytitle, fontsize=20)
second.set_xlim(rbins.min(), rbins.max())
second.legend(fontsize=20)
fig.tight_layout()
plt.subplots_adjust(top=0.88)
The centered profile calculations are straightforward, and this package uses the formulas given in Wright & Brainerd 2000 for this. However, as outlined above, the calculation of the miscentered profiles requires a double integration for $\Sigma^\mathrm{sm}(R)$, and there is a third integration for calculating $\Delta\Sigma^\mathrm{sm}(R) = \overline{\Sigma^\mathrm{sm}}(< R) - \Sigma^\mathrm{sm}(R)$.
For increased precision, you can adjust parameters specifying the number of bins to use in these integration (but note that this comes at the expense of increased computational time). See the documentation for details.