Public Data Access Overview / Discussion Forum

Hello.
I'm working on the construction of density radial profiles for gas, dark matter and star particles (and also their combination, i.e., dm + gas + stars). I followed the instructions here: dark matter and stellar binned masses.. My goal is to reproduce Figure 2. of Y. Wang et al.:

but, instead of obtaining a dm and total profiles density of order ~1e7 (it's clearly superior, but I'm using this for comparison), I got a result of order ~1e10 (c.f. fig. bellow)

I wish to know if I'm doing something wrong or missing some step. The script I constructed to generate this plot is shown below.

#for all particles
with h5py.File(saved_filename) as f:
  dx_dm = f['PartType1']['Coordinates'][:,0] - sub['pos_x']
  dy_dm = f['PartType1']['Coordinates'][:,1] - sub['pos_y']
  dz_dm = f['PartType1']['Coordinates'][:,2] - sub['pos_z']

  dx_stars = f['PartType4']['Coordinates'][:,0] - sub['pos_x']
  dy_stars = f['PartType4']['Coordinates'][:,1] - sub['pos_y']
  dz_stars = f['PartType4']['Coordinates'][:,2] - sub['pos_z']

  dx_gas = f['PartType0']['Coordinates'][:,0] - sub['pos_x']
  dy_gas = f['PartType0']['Coordinates'][:,1] - sub['pos_y']
  dz_gas = f['PartType0']['Coordinates'][:,2] - sub['pos_z']


  mass_stars = f['PartType4']['Masses'][:]* 1e10 / 0.6774


  mass_gas = f['PartType0']['Masses'][:]* 1e10 / 0.6774


  mass_dm = (sub['mass_dm']/sub['len_dm'])* 1e10 / 0.6774


  rr_dm = np.sqrt(dx_dm**2 + dy_dm**2 + dz_dm**2)
  rr_dm *= scale_factor/little_h # ckpc/h -> physical kpc

  rr_stars = np.sqrt(dx_stars**2 + dy_stars**2 + dz_stars**2)
  rr_stars *= scale_factor/little_h

  rr_gas = np.sqrt(dx_gas**2 + dy_gas**2 + dz_gas**2)
  rr_gas *= scale_factor/little_h

bin_m_stars = []
bin_m_gas = []
bin_m_dm = []

bin_r_stars = []
bin_r_gas = []
bin_r_dm = []

rho_stars = []
rho_gas = []
rho_dm = []

phys_hmr = sub['halfmassrad_stars']*scale_factor/little_h

#the following loop divides the interval (0.4R*, 4R*) into 100 equally sized bins. Here R* is the stellar half-mass radius.
for n in range(100):
  w_pt4 = np.where( (rr_stars >= (n*0.036+0.4)*phys_hmr ) & (rr_stars <= ((n+1)*0.036+0.4)*phys_hmr ) )
  w_pt0 = np.where( (rr_gas >= (n*0.036+0.4)*phys_hmr ) & (rr_gas <= ((n+1)*0.036+0.4)*phys_hmr ) )
  w_pt1 = np.where( (rr_dm >= (n*0.036+0.4)*phys_hmr ) & (rr_dm <= ((n+1)*0.036+0.4)*phys_hmr ) )


  bin_m_stars.append(np.sum(mass_stars[w_pt4])) #for stars/gas
  bin_m_gas.append(np.sum(mass_gas[w_pt0]))
  bin_m_dm.append(len(w_pt1[0])*mass_dm) # for dark matter
  bin_r_stars.append(np.mean(rr_stars[w_pt4])/phys_hmr )
  bin_r_gas.append(np.mean(rr_gas[w_pt0])/phys_hmr )
  bin_r_dm.append(np.mean(rr_dm[w_pt1])/phys_hmr )



for i in range(100):
  rho_pt4 = (3*bin_m_stars[i])/(4*pi*(bin_r_stars[i]**3))
  rho_pt0 = (3*bin_m_gas[i])/(4*pi*(bin_r_gas[i]**3))
  rho_pt1 = (3*bin_m_dm[i])/(4*pi*(bin_r_dm[i]**3))
  rho_stars.append(rho_pt4)
  rho_gas.append(rho_pt0)
  rho_dm.append(rho_pt1)

rho_all = np.array(rho_stars)+np.array(rho_gas)+np.array(rho_dm)


plt.figure(figsize=(7,7))
plt.scatter(bin_r_stars,rho_stars,label='stars')
plt.scatter(bin_r_gas,rho_gas,label='gas')
plt.scatter(bin_r_dm,rho_dm,label='dm')
plt.scatter(bin_r_dm,rho_all,label='all')
plt.xscale('log');plt.yscale('log')
plt.axvline(x=3)
plt.legend()

At last, this is the subhalo used: (http://www.tng-project.org/api/TNG100-1/snapshots/z=0/subhalos/214452)

Thanks in advance.

The units look ok, I suspect this is just a different volume normalization definition. Do I understand correctly, you are normalizing by the volume of a sphere (i.e. cumulative) instead of a shell (i.e. differential)? Usually you would make a radial density profile by diving by the volume of the spherical shell between a radius r0 and r1 for a given bin.

Yes, not using the volume of the spherical shell was my mistake. After correcting this part, I got the image below.

the difference in some points is due to the fact that the authors binned the radial interval in a logarithmic interval, while I binned it in a linear interval.