Galaxy Catalog Analysis Example: Galaxy-galaxy lensing

In this example, we’ll show how to calculate \(\Delta\Sigma(r),\) the galaxy-galaxy lensing signal of a mock catalog.

There is also a Jupyter Notebook in the following location that can be used as a companion to the material in this section of the tutorial:


By following this tutorial together with this notebook, you can play around with your own variations of the calculation as you learn the basic syntax.

Generate a mock galaxy catalog

Let’s start out by generating a mock galaxy catalog into an N-body simulation in the usual way. Here we’ll assume you have the z=0 rockstar halos for the bolshoi simulation, as this is the default halo catalog.

from halotools.empirical_models import PrebuiltSubhaloModelFactory
model = PrebuiltSubhaloModelFactory('behroozi10')
from halotools.sim_manager import CachedHaloCatalog
halocat = CachedHaloCatalog(simname = 'bolshoi', redshift = 0, halo_finder = 'rockstar')

Extract subsamples of galaxies and dark matter particles

Predictions for galaxy-galaxy lensing are calculated from the cross-correlation between the galaxy positions and the dark matter density field. For the latter, we will use the dark patter particle positions stored in the ptcl_table attribute of the mock.

px = model.mock.ptcl_table['x']
py = model.mock.ptcl_table['y']
pz = model.mock.ptcl_table['z']

As described in Formatting your xyz coordinates for Mock Observables calculations, functions in the mock_observables package such as delta_sigma take array inputs in a specific form: a (Npts, 3)-shape Numpy array. You can use the return_xyz_formatted_array convenience function for this purpose, which has a built-in mask feature that we’ll also demonstrate to select a random downsampling of \(10^{5}\) dark matter particles.

px = model.mock.ptcl_table['x']
py = model.mock.ptcl_table['y']
pz = model.mock.ptcl_table['z']

from halotools.mock_observables import return_xyz_formatted_array

import numpy as np
Nptcls_to_keep = int(1e5)
randomizer = np.random.random(len(model.mock.ptcl_table))
sorted_randoms = np.sort(randomizer)
ptcl_mask = np.where(sorted_randoms < sorted_randoms[Nptcls_to_keep])[0]
particle_positions = return_xyz_formatted_array(px, py, pz, mask = ptcl_mask)

Now we’ll extract the x, y, z positions of various subsamples of our galaxies.

x = model.mock.galaxy_table['x']
y = model.mock.galaxy_table['y']
z = model.mock.galaxy_table['z']

mstar11_mask = model.mock.galaxy_table['stellar_mass'] > 1e11
mstar11_positions = return_xyz_formatted_array(x, y, z, mask = mstar11_mask)

mstar105_mask = (model.mock.galaxy_table['stellar_mass'] > 10**10.3) & (model.mock.galaxy_table['stellar_mass'] < 10**10.7)
mstar105_positions = return_xyz_formatted_array(x, y, z, mask = mstar105_mask)

mstar105_central_mask = mstar105_mask * (model.mock.galaxy_table['halo_upid'] == -1)
mstar105_central_positions = return_xyz_formatted_array(x, y, z, mask = mstar105_central_mask)

mstar105_satellite_mask = mstar105_mask * (model.mock.galaxy_table['halo_upid'] != -1)
mstar105_satellite_positions = return_xyz_formatted_array(x, y, z, mask = mstar105_satellite_mask)

Calculate \(\Delta\Sigma(R_{\rm p})\)

As of Halotools 0.5, the delta_sigma function requires you to specify the particle masses of your simulation, as well as the factory by which you have randomly downsampled the particles in the snapshot to perform your calculation.

logrp_bins = np.linspace(-1,1,15)
rp_bins = 10**logrp_bins

particle_masses = halocat.particle_mass
downsampling_factor = (halocat.num_ptcl_per_dim**3)/float(len(particle_positions))
effective_particle_masses = particle_masses*downsampling_factor

result_mstar11_in_mpc = mean_delta_sigma(mstar11_positions, particle_positions,
                                    rp_bins, period)

result_mstar105_in_mpc = mean_delta_sigma(mstar105_positions, particle_positions,
                                    rp_bins, period)

result_mstar105_central_in_mpc = mean_delta_sigma(mstar105_central_positions, particle_positions,
                                    rp_bins, period)

result_mstar105_satellite_in_mpc = mean_delta_sigma(mstar105_satellite_positions, particle_positions,
                                    rp_bins, period)

Recall that all Halotools length units are comoving and in Mpc/h. However, the conventional units to plot \(\Delta\Sigma\) are \(h*M_{\odot}/pc^2\), since in those units the galaxy-galaxy lensing signal is roughly order unity for typical \(L_{\ast}\) galaxy samples. So now we convert units and plot the results.

length_unit_conversion = 1e6
result_mstar11_in_pc = result_mstar11_in_mpc/length_unit_conversion**2
result_mstar105_in_pc = result_mstar105_in_mpc/length_unit_conversion**2
result_mstar105_satellite_in_pc = result_mstar105_satellite_in_mpc/length_unit_conversion**2
result_mstar105_central_in_pc = result_mstar105_central_in_mpc/length_unit_conversion**2

Plot the results

logrp_mids = 0.5*(logrp_bins[:-1] + logrp_bins[1:])

fig, ax = plt.subplots(1, 1, figsize=(8, 6))


__=ax.plot(10**logrp_mids, result_mstar11_in_pc, label=r'All galaxies: $M_{\ast} > 10^{11}M_{\odot}$')
__=ax.plot(10**logrp_mids, result_mstar105_satellite_in_pc, label=r'Satellites: $M_{\ast} \approx 10^{10.5}M_{\odot}$')
__=ax.plot(10**logrp_mids, result_mstar105_in_pc, label=r'All galaxies: $M_{\ast} \approx 10^{10.5}M_{\odot}$')
__=ax.plot(10**logrp_mids, result_mstar105_central_in_pc, label=r'Centrals: $M_{\ast} \approx 10^{10.5}M_{\odot}$')

__=ax.set_xlim(0.1, 10)
__=ax.set_ylim(0.5, 200)

__=ax.set_xlabel(r'$R_{\rm p} $  $\rm{[Mpc / h]}$', fontsize=16)
__=ax.set_ylabel(r'$\Delta\Sigma(R_{\rm p})$  $[h M_{\odot} / {\rm pc}^2]$', fontsize=16)
__=ax.legend(loc='best', fontsize=13)

This tutorial continues with Galaxy Catalog Analysis Example: Projected galaxy clustering.