I am wondering is there a way to access the mass loss rate due to stellar evolution between every two snapshot for one subhalo. I implement the following method but I don't know if there are any data field in IllustrisTNG which records the mass loss.
You're right, this isn't recorded any way at the group catalog level, so you'll have to derive it from the star particles themselves.
What I would probably do is look at the subhalo at snapshot N and load all Masses of all stars. Comparison to GFM_InitialMass will give the integrated mass loss since birth. You could compute the average using the ages, or else for the instantaneous, you should load also the IDs of those stars, then the IDs of all stars in snapshot N-1, find the same set of stars (all must exist at the previous snapshot, since star particles do not disappear). The difference of their Masses sum would be the mass loss between these two times.
One more question: as you said, since the stellar ID never disappeared in the simulation. If some IDs is missing in the catalog of a subhalo, it will mean that these IDs probably are stripped away. Is my assumption correct here? I am curious about if a stellar particle is totally destroyed by supernova explosion, is IllustrisTNG always keep that the mass of that ID 0 and not remove it from the simulation?
You're right that a star with a given ID might not be in the same subhalo (i.e. its progenitor or descendant) at a different snapshot, since it could have been stripped for instance. But it will be somewhere in the simulation volume, thus why you should in general search over all stars (not such a problem, as there aren't too many stars, as compared to e.g. gas or DM).
The physical model for SN explosions in Illustris[TNG] does not ever completely destroy star particles, because star particles represent entire populations, not just single stars.