  We study accretion remnants around Milky Way analogs in the IllustrisTNG
simulations to determine how well commonly used distribution functions (DFs)
describe their phase-space distributions. We identify 30 Milky Way analogs 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 are 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 scale radius
between 2-4 kpc would be a good first-order representation of GS/E, and comment
on existing observational evidence for this as well as studies which 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 which seek to constrain both the spatial and
kinematic properties of merger remnants in the Milky Way stellar halo.
