Not all subhaloes are created equal: modelling the diversity of subhalo density profiles in TNG50

In this work, we analyse the density profiles of subhaloes with masses M_sh ≥ 1.4 × 10⁸ M_⊙ in the TNG50 simulation, with the aim of including baryonic effects. We evaluate the performance of frequently used models, such as the standard Navarro-Frenk-White (NFW), the Einasto, and a smoothly truncated version of the NFW profile. We find that these models do not perform well for the majority of subhaloes, with the NFW profile giving the worst fit in most cases. This is primarily due to mismatches in the inner and outer logarithmic slopes, which are significantly steeper for a large number of subhaloes in the presence of baryons. To address this issue, we propose new three-parameter models and show that they significantly improve the goodness of fit independently of the subhalo's specific properties. Our best-performing model is a modified version of the NFW profile with an inner log-slope of -2 and a variable truncation that is sharper and steeper than the slope transition in the standard NFW profile. Additionally, we investigate how both the parameter values of the best density profile model and the average density profiles vary with subhalo mass, V_max, distance from the host halo centre, baryon content, and infall time, and we also present explicit scaling relations for the mean parameters of the individual profiles. The newly proposed fit and the scaling relations are useful to predict the properties of realistic subhaloes in the mass range 10⁸ M_⊙ ≤M_sh ≤ 10¹³ M_⊙ that can be influenced by the presence of baryons.