Public Data Access Overview / Discussion Forum

Hi,

I was trying to calculate comoving densities of star particles to evaluate the dynamical states of galaxies. So I started with the StellarHsml as the SPH kernel smoothing length to perform kernel estimation for star particles. But there are two concerns:

1. You said in the Data Specification that this radii are useful for visualization, so does it imply that it's not advisable to use them for scientific measurements of, say, densities? If so, what do you think is a better way to estimate stellar mass densities at each star particle's position?
2. I built KDTree's using the star particles and tried querying the nearest neighbors within StellarHsml. But for many star particles I get far fewer than 32\pm 1, e.g., 17 neighbors. So is there any ceiling or constraints for StellarHsml, when there are really not many neighbors around a star particle?

Thanks!

Dear Haowen,

StellarHsml is just as described in the documentation, you are free to use it if useful. It isn't intended as a density estimator, but it could be used as one.

(More often, people would just like to make e.g. a radial profile of stellar density or stellar surface density, in which case you can just sum up all the star particle masses and divide by the bin volume or area).

I'm surprised you found such a star particle with only 17 neighbors within StellarHsml. If you can put here the simulation, snapshot, and star particle ID or index I will take a look.

Hi Dylan,

Below is a code snipped I used to find the number of neighbor star particles, and a subset of the outputs:

kdtree_star = KDTree(star['Coordinates'], boxsize=box_size)
for i in range(len(star['Masses'])):
            if not i % 1000:
                print(i, '/', len(star['Masses']))
            ind = kdtree_star.query_ball_point(star['Coordinates'][i], r=star['StellarHsml'][i], workers=-1)
            if len(ind) < 10:
                print('Snapshot: %d, Subfind ID: %d, star particle index: %d, # of neighbors: ' % (n_snap, ID, i), len(ind))

outputs example:

Snapshot: 17, Subfind ID: 0, star particle index: 289691, # of neighbors:  5
Snapshot: 17, Subfind ID: 0, star particle index: 289715, # of neighbors:  6
Snapshot: 17, Subfind ID: 0, star particle index: 289731, # of neighbors:  6
Snapshot: 17, Subfind ID: 0, star particle index: 289785, # of neighbors:  5
Snapshot: 17, Subfind ID: 0, star particle index: 289802, # of neighbors:  9
Snapshot: 17, Subfind ID: 0, star particle index: 289820, # of neighbors:  7
Snapshot: 17, Subfind ID: 0, star particle index: 289842, # of neighbors:  6

So, I'm afraid the situation is much worse than you've thought, as there are particles with even fewer neighbors. The number of particles with fewer than 10 neighbors within StellarHsml in the subhalo mentioned above is 5508 (there are ~4e5 star particles in that subhalo). What makes sense, though, is that these particles are mostly in the outskirts of the subhalo, which makes it relatively hard to find neighbors. But If that's unexpected, could you please look into this issue? PS: This problem also applies to SubfindHsml, which often does not enclose the closest 64 +/- 1 DM particles.

This is snapshot 17 of what simulation?

Sorry, it was TNG50-1.

If I load all the stars of TNG50-1 at snapshot 17, then I see:

pos[289691,:]

Out[13]: array([ 5772.67829945, 25893.67652536, 25846.987316  ])

ids[289691]

Out[14]: 101456533949

StellarHsml[289691]

Out[15]: 1.0550892

In [16]: dists = periodicDists(pos[289691,:],pos)

In [17]: len( np.where(dists <= StellarHsml[289691])[0] )

Out[18]: 33

Perhaps it isn't the same star, if you had loaded just the stars of subhalo 0? But this is likely the problem - the StellarHsml value is computed "globally", considering all particles in the snapshot, not just the stars in one subhalo.

Hi Dylan,

Sorry for the ambiguity. The particle indices I listed were for subhalo 0 only. But Thank you very much for clarifying how StellarHsml is calculated! So to estimate the comoving stellar mass density, in principle I need to load all the star particles and use StellarHsml to query the tree, right?

In general that's true, to calculate any property from the particles or cells you need to be careful to load "beyond" a single subhalo or halo, if for example you are considering near the boundary of that object, where particles/cells exist in the volume, but are outside the object.