The radius of subhalo

Qin PENG
  • 17 Dec '20

Hi Dylan,

Could you tell me how to compute the r_200 and the r_500 of the subhalo? In general, are there formulas to compute those radii which the snapshot do not provide?

Dylan Nelson
  • 17 Dec '20

Hi Qin,

These are provided only for (FoF) halos, i.e. "groups", and not for subhalos.

In general, if you want these values for a central subhalo, then you just look up the value of its parent group (SubhaloGrNr).

If you want these values for a satellite subhalo, it is less clear that this makes physical sense, and these values are not in the catalog, but you could compute them yourself. Then I would suggest to make a cumulative enclosed mass radial profile, going outwards until the enclosed density equals the chosen overdensity, i.e. just follow the definition exactly.

Qin PENG
  • 19 May

Hi Dylan,

When you calculate the r200 of a halo, you consider all the particles within the sphere or only the particles belong to the FoF?

Thank you.

Dylan Nelson
  • 19 May

The values in the catalog (such as Group_R_Crit200) are the usual "spherical overdensity" definitions, i.e. they are computed considering all particles/cells in the simulation.

Qin PENG
  • 19 May

Thank you. I think I can have a good sleep.

Qin PENG
  • 2
  • 21 May

Hi Dylan,

I compared the cirtical density and the values calculated from Group_R_Crit200 and Group_M_Crit200. They are not consistent. So, what happen? The following is my code.

import illustris_python as il
import numpy as np
import astropy
import h5py
from astropy.cosmology import FlatLambdaCDM
from astropy import units as u
from astropy.constants import k_B, m_p, M_sun, G

def cal(snapNum,haloID):

    with h5py.File(il.snapshot.snapPath(basePath, snapNum), 'r') as f:
        header = dict(list(f['Header'].attrs.items()))

    h = header['HubbleParam']
    Om0 = header['Omega0']
    Ob0 = header['OmegaBaryon']
    OmL  = header['OmegaLambda']
    z = header['Redshift']
    BoxSize = header['BoxSize']

    H = h*100.

    H0 = 100.*h*u.km/(u.s*u.Mpc)

    cosmo = FlatLambdaCDM(H0=H,Om0=Om0,Ob0=Ob0)

    group_fields = il.groupcat.loadSingle(basePath=basePath,snapNum=snapNum,haloID=haloID)

    R_Crit200 = group_fields['Group_R_Crit200']
    Group_M_Crit200 = group_fields['Group_M_Crit200']

    R_200 = R_Crit200/h
    M_200 = Group_M_Crit200*1.e10/h

    rho_c = (3.*((cosmo.H(z))**2.)/(8.*np.pi*G)).to_value('g/cm3')

    print('rho_c: '+str(rho_c))

    Rho = M_200*M_sun.to_value('g')/(4/3*np.pi*((R_200*u.kpc.to('cm'))**3))/200

    print('Calculate: '+str(Rho))


cal(98,100)
cal(96,100)
cal(94,100)
cal(85,100)
cal(79,100)

The output is

rho_c: 8.695940038154481e-30
Calculate: 8.451114320552526e-30
rho_c: 8.897716120485159e-30
Calculate: 8.053119361338877e-30
rho_c: 9.114355308430861e-30
Calculate: 7.683802112321716e-30
rho_c: 1.0335492706884433e-29
Calculate: 6.28251581145008e-30
rho_c: 1.1453751646326055e-29
Calculate: 5.546595653884129e-30
Dylan Nelson
  • 21 May

Hello,

Because it seems correct at z=0 but not at high redshift, the issue is usually a missing scalefactor a to account for comoving units.

In this case, Group_R_Crit200 is comoving (as all all lengths).

  • Page 1 of 1