I am wondering that how Subfind identifies subhalo in the simulation. Is there any specific paper talk about the algorithm. Thanks!
You can find the details in Springel+ 2001.
Maybe also: https://arxiv.org/abs/1203.3695 is interesting.
I have a question on Subfind algorithm. Reading Springel+2001, it appears that the Subfind algorithm is applied to dark matter overdensities within FoF subhalos. It is less clear how the algorithm changes if baryons are introduduced. Should I rely on Dolag+2009? Thanks
Yes, Subfind runs on a "primary" mass component which is DM. Then baryons are attached to their nearest DM particles. Furthermore, gas temperature is accounted for during the gravitational unbinding procedure, as described in Dolag+09.