2025MNRAS.538..553L

We study accretion remnants around Milky Way analogues in the IllustrisTNG simulations to determine how well commonly used distribution functions (DFs) describe their phase-space distributions. We identify 30 Milky Way analogues and 116 remnants from mergers with stellar mass ratios greater than 1:20. Two power density profiles, as well as rotating constant-anisotropy and Osipkov-Merritt DFs, are fit to the remnants. We determine that the remnants are suitable for equilibrium modelling by assessing them in the context of the Jeans equation. Each of the models we consider is reasonably able to fit the stellar remnant energy and angular momentum distribution, as well as the magnitude and shape of velocity dispersion profiles. Case studies matched to two well-known merger remnants in the stellar halo - Gaia-Sausage/Enceladus (GS/E) and Sequoia - are explored in more depth. We find good evidence that remnants with high anisotropy $\beta$, such as GS/E, are better modelled with a superposition of two Osipkov-Merritt DFs than either a constant-anisotropy model or a single Osipkov-Merritt DF. We estimate an Osipkov-Merritt profile with a scale radius between 2 and 4 kpc would be a good first-order representation of GS/E, and comment on existing observational evidence for this as well as studies that could demonstrate it. Overall, we find that DF-based models work well for describing the kinematics of large merger remnants. Our results will be an important reference for future studies that seek to constrain both the spatial and kinematic properties of merger remnants in the Milky Way stellar halo.