Halo Catalog Analysis Example: calculating radial profiles

In this example, we’ll show how to start from a subhalo catalog and calculate how various halo properties vary as a function of 3d distance from some nearby halo, i.e. we’ll show how to use Halotools to calculate radial profiles. In particular, we’ll see how the radial_profile_3d function can be used to show that the mass accretion rate of halos decreases as those halos approach a nearby cluster. We’ll also show how to calculate the mean number density of some sample as a function of cluster-centric distance.

Before following this tutorial, make sure you understand Halo Catalog Analysis Example: halo properties as a function of host halo mass, as that tutorial covers basic material such as how halo data is organized into an Astropy Table object.

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

halotools/docs/notebooks/halocat_analysis/basic_examples/halo_catalog_analysis_tutorial2.ipynb

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.

Example 1: \(dM_{\rm vir}/dt\) vs. halo-centric distance

Let’s start out by selecting two samples of halos from the Bolshoi simulation: one sample of group-mass host halos, and a second sample of lower-mass halos.

Note that the lower-mass sample includes subhalos since we do not make a halo_upid cut.

from halotools.sim_manager import CachedHaloCatalog
halocat = CachedHaloCatalog(simname = 'bolshoi', redshift = 0)

host_mask = ((halocat.halo_table['halo_upid'] == -1) &
                (halocat.halo_table['halo_mvir'] > 1e13) &
                (halocat.halo_table['halo_mvir'] < 2e13))
group_mass_hosts = halocat.halo_table[host_mask]

low_mass_mask = ((halocat.halo_table['halo_mpeak'] > 5e11) &
                (halocat.halo_table['halo_mpeak'] < 6e11))
low_mass_halos = halocat.halo_table[low_mass_mask]

The first question we will ask of these is halos is the following: how does the mass accretion rate of (sub)halos vary as a function of the distance to the group-mass hosts?

As with all mock_observables functions that accept 3d positions for arguments, we first format our positions into an array of the expected shape.

from halotools.mock_observables import return_xyz_formatted_array

group_mass_hosts_pos = return_xyz_formatted_array(group_mass_hosts['halo_x'], group_mass_hosts['halo_y'], group_mass_hosts['halo_z'])
low_mass_halos_pos = return_xyz_formatted_array(low_mass_halos['halo_x'], low_mass_halos['halo_y'], low_mass_halos['halo_z'])

The radial_profile_3d function accepts two different kinds of inputs for the separation bins. If you pass in rbins_normalized and normalize_rbins_by, then this combination of arguments allows you to calculate how various quantities vary as a function of, for example, \(x = r / R_{\rm vir}.\) The way this works is that the normalize_rbins_by argument stores the value of \(R_{\rm vir}\) for each point in sample1, and the rbins_normalized argument will be interpreted as referring to the distance of points in sample2 scaled by this value.

In the following call to the radial_profile_3d function, we will calculate the mean mass accretion rate of the lower-mass halos in 20 bins linearly spaced between \(0.1 < r / R_{\rm vir} < 15.\)

from halotools.mock_observables import radial_profile_3d

rbins_normalized = np.linspace(0.1, 15, 20)
rbins_midpoints = (rbins_normalized[:-1] + rbins_normalized[1:])/2.

result = radial_profile_3d(group_mass_hosts_pos, low_mass_halos_pos, low_mass_halos['halo_dmvir_dt_tdyn'],
                           rbins_normalized = rbins_normalized,
                           normalize_rbins_by = group_mass_hosts['halo_rvir'],
                           period=halocat.Lbox)

plt.plot(rbins_midpoints, result, color='blue')

plt.xlabel(r'$r / R_{\rm vir}$', fontsize=20)
plt.ylabel(r'$\langle dM_{\rm vir} / dt \rangle$  $[M_{\odot}/{\rm yr}]$', fontsize=20)
plt.xticks(size=15); plt.yticks(size=15)

As we’ll see in the next example, if you instead want to calculate the profile of a quantity as a function of the absolute distance rather than some scaled distance, you can use the rbins_absolute argument instead.

Example 2: Number density vs host-centric distance

In this next example, we’ll calculate the answer to the following question: how does the abundance of (sub)halos vary as a function of host-centric distance to our sample of group-mass host halos?

The radial_profile_3d function has a return_counts argument that can be used to additionally return the number of objects as a function of the input distance. In the following call to the radial_profile_3d function, we will calculate the mean mass accretion rate of the lower-mass halos in 20 bins linearly spaced in r between \(0.5 {\\rm Mpc} < r < 10 {\\rm Mpc}.\)

rbins_absolute = np.linspace(0.5, 10, 20)
rbins_midpoints = (rbins_absolute[:-1] + rbins_absolute[1:])/2.

result, counts = radial_profile_3d(group_mass_hosts_pos, low_mass_halos_pos, low_mass_halos['halo_dmvir_dt_tdyn'],
                           rbins_absolute = rbins_absolute,
                           period=halocat.Lbox, return_counts = True)

shell_volumes = 4*np.pi*(rbins_midpoints**2)*np.diff(rbins_absolute)
mean_number_density = (counts/float(len(group_mass_hosts)))/shell_volumes

plt.plot(rbins_midpoints, mean_number_density, color='blue')

plt.xlabel(r'$r$  $[Mpc]$', fontsize=20)
plt.ylabel(r'$\langle dN/dV \rangle$  $[M_{\odot}/{\rm Mpc^{3}}]$', fontsize=20)
plt.xticks(size=15); plt.yticks(size=15)