Public Data Access Overview / Discussion Forum

Dear Dylan,

I am trying to compare the maximum velocity computed from the dynamical mass with Vmax from the Subfind Subhalo (galaxies) catalog. However, they do not match.
Could you elaborate on the method of computing Vmax in the Subhalo catalog? I have followed the reference of Lovell et al. (2019), and this is the approach we want to use for our rotation curves (RCs) or dynamical mass profiles.

I compute dynamical mass of an aligned galaxy (x-y plane, face-on) as following:

###Magnitude of radius or projected radius
Gas_Radius = np.linalg.norm(gas_pos_align , axis=1)
Star_Radius = np.linalg.norm(stars_pos_align , axis=1)
DM_Radius = np.linalg.norm(DM_pos_align , axis=1)
Gas_Mass = gas_masses = gas['Masses'][:] #In solar masses; 
Star_Mass = star_masses = stars['Masses'][:] #In solar masses; 
###count the number of DM particles using DM positions and multipy them by  $3.1 \times 10^5 M_\odot/h$  or  $0.45 \times 10^6 M_\odot$  given in TNG5-1 web : https://www.tng-project.org/data/docs/background/
DM_Mass = np.full(len(DM_Radius), 0.45*1e6) # in solar masses
def component_Mdyn_optimized(Xarray, Yarray, bin_length=None):
    ###radial bins
    if bin_length is None:
        bin_length = np.max(Xarray) / 2
    bin_edges = np.geomspace(0.1, 100, num=int(bin_length) + 1)
    ###Determine the bin index for each point
    indices = np.digitize(Xarray, bin_edges) - 1   #Subtract 1 to align with Python's 0-based indexing
    ###Initialize arrays to store the sum of Yarray values in each bin
    binned_M = np.zeros(len(bin_edges))
    ###Calculate the sum of Yarray values for each bin
    for i in range(len(bin_edges)):
        mask = indices <= i
        binned_M[i] = np.sum(Yarray[mask])
    return bin_edges, binned_M

Gas_R, Gas_M = component_Mdyn_optimized(Gas_Radius, Gas_Mass, bin_length=100)
Star_R, Star_M = component_Mdyn_optimized(Star_Radius, Star_Mass, bin_length=100)
DM_R, DM_M = component_Mdyn_optimized(DM_Radius, DM_Mass, bin_length=100)
Total_mass = Gas_M + Star_M + DM_M

Then I assume the spherical symmetry and compute the circular velocity curve for each component and spherically average them

def spherical_sym_RCs(R, M):
    G = 4.302*10**(-6) # units (kpc/Msun)*(km2/s2)
    Vc = np.sqrt(G * M / R)
    return Vc 
Gas_Vc = spherical_sym_RCs(Gas_R, Gas_M)
Star_Vc = spherical_sym_RCs(Star_R, Star_M)
DM_Vc = spherical_sym_RCs(DM_R, DM_M)
BH_Vc = spherical_sym_RCs(DM_R, BH_masses)
#Speharically averaged Rotation curve
Vc_tot_mass = spherical_sym_RCs(DM_R, Total_mass)
Vc_tot = np.sqrt(Gas_Vc**2 + Star_Vc**2 + DM_Vc**2)

I am not sure whats wrong in my calculations

Even After adding black hole, discrepency is about 7km/s
For your reference: Snapshot_number 89 and subhalo_id 411594

The calculation of SubhaloVmax and SubhaloVmaxRad are going to be similar (identical) to the Subfind code you find in the public version of AREPO.