# Not All Subhaloes Are Created Equal: Modelling the Diversity of Subhalo
Density Profiles in TNG50

Felix M. Heinze,1 Giulia Despali2,3 and Ralf S. Klessen4,5
1Friedrich-Schiller-Universität Jena, Theoretisch-Physikalisches Institut,
Fröbelstieg 1, D-07743 Jena, Germany
2Alma Mater Studiorum - Università di Bologna, Dipartimento di Fisica e
Astronomia "Augusto Righi", Via Gobetti 93/2, I-40129 Bologna, Italy
3INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Gobetti 93/3, I-40129 Bologna, Italy
4Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Straße 2, D-69120 Heidelberg, Germany
5Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches
Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this work, we analyse the density profiles of subhaloes with masses
$M_{\mathrm{sh}}\geq 1.4\times 10^{8}$ M⊙ in the TNG50 simulation, with the
aim of including baryonic effects. We evaluate the performance of frequently
used models, such as the standard 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_{\mathrm{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^{8}$ M⊙ $\leq M_{\mathrm{sh}}\leq$ $10^{13}$
M⊙ that can be influenced by the presence of baryons.

###### keywords:

galaxies: haloes – cosmology: dark matter – methods: numerical

††pubyear: 2023††pagerange: Not All Subhaloes Are Created Equal: Modelling the
Diversity of Subhalo Density Profiles in TNG50–A

##  1 Introduction

The $\Lambda$ Cold Dark Matter ($\Lambda$CDM) model has so far been widely
successful in predicting and explaining numerous observations on cosmological
scales, such as the existence and properties of the cosmic microwave
background (CMB), the expansion history of the Universe, the abundances of the
chemical elements, as well as the formation and distribution of cosmic large-
scale structure (Planck Collaboration et al., 2020). Another key prediction of
the $\Lambda$CDM model is the presence of a large number of small dark matter
subhaloes within larger host haloes (Zavala & Frenk, 2019), which have first
been studied in $N$-body simulations by Moore et al. (1999b); Klypin et al.
(1999a). The increase in available computational resources over the last few
decades allowed for ever more precise theoretical predictions on galactic and
subgalactic scales by using high-resolution simulations. On these small
scales, a number of discrepancies between the theoretical model and
observations have been identified in the past (Bullock & Boylan-Kolchin,
2017). A few common examples are the cusp-core problem (Flores & Primack,
1994; Moore, 1994), the missing satellites problem (Klypin et al., 1999b;
Moore et al., 1999b) and the too big to fail problem (Boylan-Kolchin et al.,
2011), which point out mismatches in the abundances, central densities and
inner density profile shapes for subhaloes in simulations and dwarf galaxies
in the Local Group. However, there are still a number of uncertainties in both
the model predictions and the observations, which led to various attempts of
resolving these small-scale issues, both within the $\Lambda$CDM framework
(for example by including baryonic feedback, see e.g. Mashchenko et al., 2008;
Pontzen & Governato, 2012; Brooks & Zolotov, 2014; Read et al., 2016; Sawala
et al., 2016; Wetzel et al., 2016) and also by invoking other dark matter
models, such as warm dark matter (WDM, see e.g. Bode et al., 2001; Lovell et
al., 2014; Despali et al., 2020; Lovell & Zavala, 2023) or self-interacting
dark matter (SIDM, see e.g. Vogelsberger et al., 2012; Rocha et al., 2013;
Despali et al., 2019, 2022; Mastromarino et al., 2023; Robles et al., 2019;
Lovell & Zavala, 2023). Yet another idea to solve these problems is to modify
the theory of gravity, which is done for example in Modified Newtonian
Dynamics (MOND, see e.g. Milgrom, 1983; Sanders & McGaugh, 2002; Milgrom,
2002; Famaey & McGaugh, 2012; Kroupa, 2012). Up to the present day, these
controversies and debates are still ongoing, and in order to find the most
suitable model it is necessary to probe these subgalactic scales in ever more
detail.

Investigating the properties of dwarf galaxies in the Local Group has for a
long time been the primary method for probing these small scales
observationally (Mateo, 1998; Tolstoy et al., 2009). More recently, other
methods have been used, which even allow for studying the properties of small
(sub-)haloes that do not contain enough luminous matter to be detectable via
emission of radiation. In the Milky Way, these dark subhaloes can be detected
by the perturbations they cause in stellar tidal streams (Carlberg, 2012;
Erkal & Belokurov, 2015), which led to the detection of a subhalo with a mass
of $10^{6}$-$10^{8}$ M⊙ in the GD-1 stellar stream by Bonaca et al. (2019).
Extragalactic subhaloes can be detected by the effect they have on the flux
ratios of multiply-imaged quasars (Gilman et al., 2019; Hsueh et al., 2020) or
by the perturbations they cause in strong galaxy-galaxy lens systems. Using
this method, four subhalo detections have been claimed up to this day (Vegetti
et al., 2010a, b, 2012; Hezaveh et al., 2016), with masses between $1.9\times
10^{8}$ and $2.7\times 10^{10}$ M⊙, but many more are expected to follow from
the upcoming Euclid and LSST surveys (Collett, 2015) and high-resolution
follow-up images. Minor et al. (2021) and Şengül & Dvorkin (2022) reanalysed
some of the previously detected subhaloes and estimated that their
concentrations are exceptionally high when modelled with NFW profiles, making
them 2$\sigma$ outliers of the CDM model. The inferred concentration of the
GD-1 stream perturber in Bonaca et al. (2019) is also higher than the average
NFW concentration, pointing in the same direction.

In order to reliably test the $\Lambda$CDM model with these observations,
accurate theoretical predictions for the properties of these subhaloes are
needed, such as their mass function, their internal structure, their dynamical
properties, as well as their structural evolution as they experience tidal
stripping due to the presence of their host halo and tidal shocks due to close
encounters with other subhaloes. Unfortunately, state-of-the-art simulations
are limited in resolution and are therefore not able to resolve the entire
possible mass range of subhaloes down to earth masses. The finite resolution
also leads to several unwanted effects, such as artificial disruption (van den
Bosch et al., 2018; van den Bosch & Ogiya, 2018; Green et al., 2021), which
can severely limit the ability to obtain reliable data for the subhalo
properties.

Moliné et al. (2017, 2023) recently investigated concentrations and other
subhalo properties in great detail in $N$-body simulations for a wide range of
different subhalo and host halo masses and how they evolve over cosmic time.
Furthermore, $N$-body simulations have shown that subhaloes generally exhibit
much higher central densities and are on average more concentrated than
isolated haloes of the same mass (Ghigna et al., 2000; Bullock et al., 2001).

The density profiles of dark matter subhaloes have also been studied
extensively before. However, this has mostly been done in zoom-in simulations,
which provide a high resolution but fewer objects for a proper statistical
analysis. Green & van den Bosch (2019) studied the evolution of the subhalo
density profiles in idealised simulations as they experience tidal stripping
and showed that the structural evolution solely depends on the initial subhalo
concentration and the fraction of mass being stripped. Di Cintio et al. (2013)
used an approach similar to the one presented in this paper in order to
investigate the density profiles of subhaloes in zoom-in simulations of the
Local Group. They showed that the Einasto profile (Einasto, 1965) provides a
better model than the commonly used Navarro-Frenk-White (NFW) profile (Navarro
et al., 1996, 1997) and some modifications of it. They further demonstrated
that the shape parameter of the Einasto profile strongly depends on the total
subhalo mass, and is also being reduced by tidal stripping. Similar results
have also been obtained before for isolated field haloes (Merritt et al.,
2006; Navarro et al., 2010; Dutton & Macciò, 2014). Di Cintio et al. (2011)
already pointed out that the NFW profile might not be an accurate model for
the density profiles of subhaloes. Springel et al. (2008) investigated the
density profiles, concentrations, mass function and radial distribution of
subhaloes in the dark-matter-only simulation Aquarius. They found that the
density profiles of subhaloes show a similar behaviour to those of main haloes
and that the Einasto profile provides a much better fit than the NFW or Moore
profile (Moore et al., 1999a), even with fixed shape parameter.

In this paper, we investigate the density profiles of subhaloes in the state-
of-the-art cosmological simulation TNG50 of the IllustrisTNG project, which
offers a good compromise between size and resolution. Recent studies have
analysed the impact of baryons on the structural properties of (sub-)haloes as
well as their mass function and found that they can have a significant impact,
making haloes more spherical, reducing the abundance of lower mass subhaloes
and leading to higher concentrations and different slopes in the inner regions
of (sub-)halo density profiles (Mollitor et al., 2015; Zhu et al., 2017;
Despali & Vegetti, 2017; Garrison-Kimmel et al., 2017; Nadler et al., 2018).
TNG50 includes the effect of baryons, together with many other physical model
components. Our main goal is to provide a simple three-parameter analytical
fit for the density profiles that can be used for a wide range of subhaloes
with different properties.

The paper is structured as follows: In Section 2, we give a brief overview of
the TNG50 simulation and its relevant technical details. We also discuss the
selection of suitable subhaloes and the computation of their density profiles.
In Section 3, we present commonly used analytical models for describing the
density profiles of dark matter haloes and subhaloes (the NFW profile, the
Einasto profile and the truncated NFW profile) and further introduce new
improved models. We compare the performances of all these models in Section 4
for subhaloes with various different kinds of properties. In Section 5, we
investigate how both the individual parameters of the best performing model
and the average density profiles vary with different subhalo properties,
including mass, Vmax, distance from the host halo centre, baryon fraction and
infall time. Furthermore, we present simple scaling relations for how the
model parameter values depend on the subhalo mass and Vmax, and discuss how
baryons and tidal stripping affect the subhalo density profiles. Finally, our
summary and conclusions are presented in Section 6.

##  2 Methods

###  2.1 The TNG50 Simulation

The subhaloes analysed in this paper are drawn from IllustrisTNG, a suite of
large-volume cosmological gravo-magnetohydrodynamical simulations (Nelson et
al., 2019a). All of these simulations were performed using the moving-mesh
code Arepo (Springel, 2010), which calculates the gravitational forces using a
Particle-Mesh-Tree method and solves the ideal magneto-hydrodynamic equations
using a finite volume method on an adaptive mesh. Each simulation starts at a
redshift of $z=127$ and finishes at the present day $z=0$, with cosmologically
motivated initial conditions and a cosmology in agreement with Planck 2015,
given by $\Omega_{\Lambda,0}=0.6911$, $\Omega_{\mathrm{m,0}}=0.3089$,
$\Omega_{\mathrm{b,0}}=0.0486$, $\sigma_{8}=0.8159$, $n_{\mathrm{s}}=0.9667$
and $h=0.6774$ (Planck Collaboration et al., 2016). They also include a large
number of physical processes such as primordial and metal-line cooling,
heating by the extragalactic UV background, stochastic star formation,
evolution of stellar populations, feedback from supernovae and AGB stars, as
well as supermassive black hole formation and feedback. It has been shown that
the TNG simulations produce results which are largely consistent with
observations, beyond the regimes which were used to calibrate the model. For
more information about the underlying model and its numerical details, see
Weinberger et al. (2017) or Pillepich et al. (2018).

In order to be able to reliably study the small subgalactic scales, we make
use of the TNG50 simulation (Nelson et al., 2019b; Pillepich et al., 2019),
since it provides the highest resolution. TNG50-1 has a cubic simulation
volume with a comoving side length of 51.7 Mpc, a dark matter mass resolution
of $3.1\times 10^{5}$ M${}_{\odot}h^{-1}$, a baryonic mass resolution of
$5.7\times 10^{4}$ M${}_{\odot}h^{-1}$ and a gravitational softening length of
0.288 kpc at $z=0$.

Haloes are identified using the friends-of-friends (FoF) group finder
algorithm (Huchra & Geller, 1982) with a linking length of 0.2. Subhaloes are
subsequently identified using the Subfind algorithm (Springel et al., 2001;
Dolag et al., 2009). Additionally, merger trees have been created using
SubLink (Rodriguez-Gomez et al., 2015) and LHaloTree (Springel et al., 2005).
The halo and subhalo finding algorithm has an impact on the on the total mass
and radial extent of the structures (Onions et al., 2012; Behroozi et al.,
2015). However, we will show in Section 4 and 5 that our newly proposed
density profile model deviates significantly from the NFW profile in the inner
regions, where the density is well constrained and where the effects on
observations are the strongest.

###  2.2 Subhalo Selection

The TNG50-1 group catalogue of the snapshot at $z=0$ contains 5,688,113
Subfind groups in total. However, due to the limitations in resolution, not
all of these objects are suitable for the analysis of the density profiles.
Furthermore, some of the objects identified by the Subfind algorithm do not
have a cosmological origin, in the sense that they have not formed due to the
process of hierarchical structure formation. Instead, some of these objects
might be fragments of already formed galaxies, produced by baryonic processes
such as disk instabilities. These have been flagged in the simulation data and
we exclude them from the analysis of the subhalo density profiles. In this
work, we also exclude subhaloes with less than 300 dark matter particles.
Below this threshold, the scatter in the parameter values of the analytical
fit functions increases drastically, which we attribute to the limitations in
resolution. Other authors use a higher number of particles as a threshold
(e.g. 1000, see Di Cintio et al. 2013), but this would only affect the first
two mass bins of our sample. We also eliminate the main haloes as well as
subhaloes which belong to a Subfind group that has a total mass of less than
1010 M⊙. After excluding all the objects which are not suitable for further
analysis, 112,048 subhaloes remain in total. Their total masses lie in the
range between 1.4 $\times$ 108 M⊙ and 8.5 $\times$ 1012 M⊙, which we divide up
into 17 mass bins that cover a large part of the subhalo mass range that can
currently be probed observationally.

###  2.3 Computation of the Density Profiles

We computed the density profiles in spherical shells centered around the
minimum gravitational potential. The corresponding 40 logarithmically spaced
radial bins lie in the range between 10-1 and 103 ckpc/$h$. We chose the lower
limit such that the radial bins reach a bit below the gravitational softening
length, which marks the limit below which the density profile can no longer be
well resolved. The upper limit was determined by the largest subhalo. The
fixed number of radial bins in this range allows for an easy computation of
average density profiles, and it is large enough to guarantee that the
smallest subhaloes with an extent of a few kpc are still well resolved. For
computing the density value of each spherical shell we used all matter
particles, including dark matter, gas, stars, and black holes. We further used
the Poissonian error to estimate the uncertainties for the density values due
to the finite number of particles that sample the density distribution. Figure
1 shows the density profiles of four representative subhaloes with different
masses, together with some of the analytical models described in Section 3.

##  3 Analytical Density Profile Models

###  3.1 Traditional Models

####  3.1.1 The NFW Profile

Navarro et al. (1996, 1997) demonstrated in their $N$-body simulations that
haloes consisting of dark matter only have a universal density profile shape,
which is independent of their mass and can be well described by the Navarro-
Frenk-White (NFW) profile:

$\rho(r)=\frac{\rho_{0}}{\frac{r}{r_{\mathrm{s}}}\left(1+\frac{r}{r_{\mathrm{s}}}\right)^{2}}.$
(1)

The profile has an inner logarithmic slope of -1 and an outer log-slope of -3,
with the transition happening at the scale radius $r_{\mathrm{s}}$, which is
the only free parameter besides the overall normalisation $\rho_{0}$.
Therefore, the profile is completely characterised by the virial mass
$M_{200}$ and the concentration $c_{200}=r_{200}/r_{\mathrm{s}}$ of the halo.
The mass and concentration are not completely independent but correlated, with
more massive haloes being less concentrated in general. Simulations have also
shown that at fixed mass, the halo concentration is correlated with assembly
time, which introduces a large halo-to-halo scatter due to the different mass
assembly histories. This concentration-mass relation has been studied both for
haloes (Wechsler et al., 2002; Zhao et al., 2003, 2009; Dutton & Macciò, 2014;
Diemer & Kravtsov, 2015) and for subhaloes (Bullock et al., 2001; Moliné et
al., 2017, 2023).

####  3.1.2 The Truncated NFW Profile

Since subhaloes can experience severe tidal truncation (Kazantzidis et al.,
2004; Diemand et al., 2007), it has been pointed out in the past that subhalo
density profiles might not be well described by the NFW profile. One way to
take the truncation into account is by using a smoothly truncated version of
the NFW profile (tNFW) with a logarithmic slope of -5 well beyond the
truncation radius $r_{\mathrm{t}}$:

$\rho(r)=\frac{\rho_{0}}{\frac{r}{r_{\mathrm{s}}}\left(1+\frac{r}{r_{\mathrm{s}}}\right)^{2}}\cdot\left(\frac{r_{\mathrm{t}}^{2}}{r^{2}+r_{\mathrm{t}}^{2}}\right).$
(2)

This model has been used e.g. by Minor et al. (2021) for the analysis of
subhalo properties inferred from strong gravitational galaxy-galaxy lensing.

####  3.1.3 The Einasto Profile

The Einasto profile (Einasto, 1965) is another commonly used three-parameter
model for describing the density profiles of dark matter haloes. Its
functional form is given by:

$\rho(r)=\rho_{-2}\exp\left\\{\frac{-2}{\alpha}\left[\left(\frac{r}{r_{-2}}\right)^{\alpha}-1\right]\right\\}.$
(3)

$r_{-2}$ is the radius at which the log-slope is equal to -2, $\rho_{-2}$ is
the corresponding density and $\alpha$ controls the shape of the profile. The
Einasto profile is no longer a double power law but instead the logarithmic
slope is a power law function of radius:

$\frac{d\log\rho}{d\log r}=-2\left(\frac{r}{r_{-2}}\right)^{\alpha}.$ (4)

It has been shown in the past that the Einasto profile generally provides a
much better model for the density profile of both haloes and subhaloes
(Merritt et al., 2006; Navarro et al., 2010; Di Cintio et al., 2011; Dutton &
Macciò, 2014).

Figure 1: Four example density profiles of representative subhaloes with
different masses, including the best fits for the NFW, the Einasto and the
modified NFW profile. The softening length is highlighted in purple. Subhaloes
with masses well above $10^{11}$ M☉ (upper left) tend to have a pronounced
bump feature in the central regions which flattens towards the centre and
forms a core. Subhaloes with masses in the range between $10^{10}$ and
$10^{11}$ M☉ tend to have an inner log-slope of approximately -2 all the way
to the centre and either a soft transition to a slightly steeper outer slope
(upper right) or a sharper truncation (lower left). Subhaloes with masses
below $10^{9}$ M☉ generally have a more regular shape, with a continuously
changing logarithmic slope and can therefore often be better described by the
Einasto profile.

###  3.2 Additional Models

In Figure 1, one can see that the previously mentioned density profile models
do not generally provide a good fit for most of the subhaloes, especially the
more massive ones. Furthermore, the best-fit parameter values of the NFW
profile take extreme and unrealistic values. Because of that, we introduce
some other fit functions that better capture the inner and outer logarithmic
slopes as well as the truncation.

####  3.2.1 The Modified NFW Profile

A generalisation of the NFW profile is the so-called $(\alpha,\beta,\gamma)$
model (see e.g. Zhao, 1996), which has the following functional form:

$\rho(r)=\frac{\rho_{0}}{\left(\frac{r}{r_{\mathrm{s}}}\right)^{\gamma}\left[1+\left(\frac{r}{r_{\mathrm{s}}}\right)^{\alpha}\right]^{(\beta-\gamma)/\alpha}}.$
(5)

Here $\gamma$ and $\beta$ set the inner and outer log-slope and $\alpha$
controls the sharpness of the transition from one slope to the other. At the
scale radius $r_{\mathrm{s}}$, the log-slope is the average of the inner and
outer slopes, which is $-(\gamma+\beta)/2$. For $(\alpha,\beta,\gamma)$ = (1,
3, 1) we get the NFW profile, but also other commonly used profiles can be
obtained from it.

The $(\alpha,\beta,\gamma)$ model has five parameters in total. When looking
at many of the computed density profiles, one can observe that profiles with a
steeper outer log-slope generally also have a sharper transition from the
inner to the outer slope. This motivates the following model with the
functional form:

$\rho(r)=\frac{\rho_{0}}{\left(\frac{r}{r_{\mathrm{s}}}\right)^{2}\left[1+\left(\frac{r}{r_{\mathrm{s}}}\right)^{\alpha}\right]^{6}}.$
(6)

Figure 2: Subhalo mass dependence of the density profile model goodness of fit
values for the individual (left) and the average density profiles (right). The
light blue histograms in the background indicate the number of subhaloes in
each mass bin. For the individual subhaloes, the mean (solid line) and median
(dashed line) values of the goodness of fit per unit mass have been computed
for each mass bin. The scatter in $\chi^{2}$ has been omitted, since it is
dominated by failing fits.

Here, $\alpha$ controls both the outer log-slope and the sharpness of the
slope transition. It therefore determines the shape of the density profile. We
set the inner log-slope to -2, and we determined the value of the other
exponent by optimising the goodness of fit for the individual and average
profiles. A value of 6 turns out to be the best choice for the subhaloes in
our TNG50 sample. The scale radius $r_{\mathrm{s}}$ marks the radius at which
the log-slope changes. In our new parametrisation it corresponds, in practice,
to the truncation radius, and it could therefore be used to determine the
extent of the inner part of the profile. As the log-slope beyond the
truncation radius is very steep, most of the subhalo mass should be contained
within $r_{\mathrm{s}}$.

Throughout the paper, we refer to the model in Equation (6) as the _modified
NFW profile_.

####  3.2.2 The Modified Schechter Profiles

Instead of modelling the truncation with a double power law like we did in the
modified NFW profile, one can use a simple power law profile with an
exponential cutoff and an additional parameter $n$ that controls the sharpness
of the truncation:

$\rho(r)=\rho_{0}\left(\frac{r}{r_{\mathrm{t}}}\right)^{-\gamma}\cdot\exp\left[-\left(\frac{r}{r_{\mathrm{t}}}\right)^{n}\right].$
(7)

This looks very similar to the functional form of the Schechter luminosity
function (Schechter, 1976) with the additional parameter $n$. Because of that,
we will refer to it as the _modified Schechter (MS) profile_. Here $\rho_{0}$
is again the overall normalisation factor, $r_{\mathrm{t}}$ is the truncation
radius and $\gamma$ is the power law log-slope. From this four-parameter model
we construct two three-parameter models by fixing one parameter and leaving
the other parameters as free parameters.

In one case, we fix the inner log-slope to -2, similarly to the modified NFW
profile and leave $n$ and the other parameters free. From now on, we will
refer to this model as MS-(2, $n$).

For the other model we instead fix the sharpness of the truncation $n$ to 2
and leave $\gamma$ as a free parameter. From now on, we will refer to this
model as MS-($\gamma$, 2). Our analysis of the parameters of the MS-(2, $n$)
model in the appendix (see Figure 24, 27 and Table 1), suggests that the mean
values of $n$ indeed lie very close to 2 for most of the mass bins. Allowing
for more freedom to adjust the inner log-slope can lead to a better fit in the
case of an additional bump in the central regions of the density profile,
which is often present for subhaloes of higher mass (see Figure 1 in the upper
left corner). Similar parametrisations with an exponential truncation for
subhaloes have also been used in the past (see e.g. Kazantzidis et al., 2004;
Errani & Navarro, 2021).

Other frequently used radial profile models are: the Moore profile (Moore et
al., 1999a), the Hernquist profile (Hernquist, 1990), the Burkert profile
(Burkert, 1995) and the (Pseudo-)Jaffe profile (Jaffe, 1983; Muñoz et al.,
2001). While they have been used for modelling the density profiles of haloes
and subhaloes in previous works, we do not consider them since the new models
presented in this section provide a much better fit for our sample.

##  4 Comparison of the Model Performances

###  4.1 Fitting Process

We obtain the optimal fit for each model and each density profile by using the
non-linear least squares fitting algorithm of the Python library SciPy,
together with the Poissonian density uncertainties. For calculating the
goodness of fit, we use a reduced chi-squared statistic to include the
influence of the uncertainties as well as the different numbers of radial bins
and parameters:

$\chi^{2}=\frac{1}{N-p}\cdot\sum_{i=1}^{N}\frac{(\log\rho_{i}-\log
f_{i})^{2}}{(\sigma_{i}/\rho_{i})^{2}},$ (8)

with the total number of nonzero data points $N$, the number of fitted
parameters $p$, the difference between the log-values of the density
$\rho_{i}$ and the model fit value $f_{i}$, and the uncertainties of the log-
profile $\sigma_{i}/\rho_{i}$ at the data point $i$. We exclude the data
points below the softening length as well as data points where the density is
zero for both the fitting and the calculation of the goodness-of-fit.

We note that $\chi^{2}$ increases with subhalo mass, since lower-mass
subhaloes exhibit larger relative density uncertainties due to the limitations
in resolution. This makes it difficult to compare the goodness of fit for
subhaloes with different masses. However, the comparison of different model
performances within individual mass bins is not affected by this. We found
that $\chi^{2}$ is roughly proportional to $M_{\mathrm{sh}}$. Because of that,
we primarily look at $\chi^{2}$ per unit subhalo mass in order to have a
quantity that is relatively independent of mass. Otherwise, high values of
$\chi^{2}$ would only indicate where the most massive subhaloes can be found.

It should also be noted that it can make a significant difference whether one
computes the fit for a linear density profile or a log-density profile. The
difference in the obtained parameter values can be quite noticeable,
especially for the models described in Section 3.1, which do not perform well
in most cases. However, for the models that provide an accurate fit, the
difference is negligible. In order to be consistent, the further analysis will
be done with fits for the linear density profiles for all the fit functions as
they generally lead to more reliable results. Despite that, we still evaluate
the goodness of fit by testing how well these linear fits can reproduce the
overall shape of the log-density profile with the definition in Equation (8).

Figure 3: Median values of the goodness of fit per unit subhalo mass for
subhaloes with different distances from the host halo centre (on the left) and
different baryon fractions (on the right). Additionally, the radial
distribution of subhaloes and the distribution of baryon fractions is shown in
the histograms.

###  4.2 Performance for Different Subhalo Masses

In Figure 2, we show the performances of the fits for both the individual and
the average density profiles as a function of subhalo mass. The histogram
indicates the number of subhaloes in each mass bin.

The NFW profile performs the worst, followed by its smoothly truncated
version. The reason for this is that most of the subhalo density profiles in
TNG50 have an inner slope of -2 and not -1 (see Figure 1). A slope of -2
however can only be found at the scale radius $r_{\mathrm{s}}$ for the NFW
profile. Therefore, the parameters have to take extreme and unrealistic
values, stretching the profile to a large extent in order to produce a
reasonably good fit.

The Einasto profile has a better goodness of fit than the NFW profile. For
lower-mass subhaloes with $M_{\mathrm{sh}}\sim 10^{8}$ M☉ it even shows a
slightly better performance than the new models presented in Section 3.2, and
it also provides the best fit for the average density profiles of the lowest-
mass subhaloes. The reason for this can be seen in Figure 1: the lower-mass
subhaloes have a much more regular shape than the higher-mass ones, with no
sharp truncation or central bump and a gradually changing logarithmic slope.
Those are the properties reflected best by the Einasto profile.

The modified NFW and Schechter profiles fit the individual density profiles
equally well most of the time. However, the modified NFW profile is preferred
in general, since the mean goodness of fit values for the MS-(2, $n$) and the
MS-($\gamma$, 2) model are often a lot worse (solid lines in Figure 2), which
indicates that there are more outliers for which the fit fails. The modified
NFW profile also provides a very good fit for most of the average density
profiles. For subhaloes with masses above $10^{11}$ M⊙ the MS-(2, $n$) works
slightly better because the shape of the truncation can no longer be well-
described by the modified NFW. The average profiles of subhaloes at the
uppermost mass end are best described by the MS-($\gamma$, 2) profile, which
is mostly due to the fact that it allows for variations of the inner
logarithmic slope, and this can partially account for the central bumps that
often occur for these subhaloes.

Our results are compatible with those from Di Cintio et al. (2013), who looked
at the density profiles of subhaloes with masses above $2\times 10^{8}$
M${}_{\sun}h^{-1}$, with 56 in the hydro and 66 in the dark-matter-only run.
They computed the mean values of the goodness of fit for the entire sample of
subhaloes, where most of them had masses in the range $10^{8.5}-10^{9.5}$ M⊙,
and they found that among their tested models, the Einasto profile provides
the best fit for both the hydrodynamic and the dark-matter-only only run. They
have not analysed the goodness of fit for individual mass bins and they also
did not account for the density uncertainties.

###  4.3 Performance for Different Distances from the host halo centre

Figure 3 (on the left) shows the goodness of fit as a function of the subhalo
distance from the host centre (expressed as fraction of the host halo’s virial
radius R200). The goodness of fit values are relatively independent of
distance from the host halo centre. They only become better for subhaloes
closer than 0.1 R200 for most of the new models presented in Section 3.2 as
well as for the Einasto profile. The truncated NFW profile fit becomes worse
with decreasing distance. We found similar trends for the $\chi^{2}$-value
also when looking at individual mass bins. The radial distribution of
subhaloes also looks similar for the individual mass bins, with the exception
that more massive subhaloes are less likely to be found both closer to the
centre and further out than the virial radius R200.

###  4.4 Performance for Other Subhalo Properties

We also checked how the goodness of fit varies with baryon fraction,
triaxiality, group mass, concentration, infall time as well as
$V_{\mathrm{max}}$ and $R_{\mathrm{max}}$, i.e. the peak circular velocity and
the radius at which this velocity is reached. More details can be found in
Appendix A.

#### Baryon fraction

Figure 3 (on the right) presents the goodness of fit as a function of the
baryon fraction. It turns out that the fits become slightly better with higher
baryon fractions, both for our new models and for the Einasto profile. The
same is also true when looking at the majority of individual mass bins. For
the truncated NFW profile, the fits become slightly worse with a higher baryon
content. The typical baryon fractions that the most massive subhaloes have is
around 10-15 per cent (see also Figure 5).

Figure 4: Median values of the goodness of fit per unit subhalo mass for
subhaloes with different concentrations. Additionally, the distribution of
concentrations is shown in the histogram.

#### Triaxiality

We define the subhalo triaxiality as (Franx et al., 1991):

$T=\frac{a^{2}-b^{2}}{a^{2}-c^{2}},$ (9)

where $a$, $b$ and $c$ are the major, intermediate and minor axes of a
triaxial isodensity surface. We computed these values by diagonalising the
shape tensor, which is also sometimes called inertia tensor (Bailin &
Steinmetz, 2005; Zemp et al., 2011):

$S_{ij}=\frac{\sum_{k}m_{k}r_{k,i}r_{k,j}}{\sum_{k}m_{k}},$ (10)

with the mass $m_{k}$ and the the $i$-th and $j$-th coordinate $r_{k,i}$ and
$r_{k,j}$ of particle $k$. The axes $a$, $b$ and $c$ are the square roots of
the eigenvalues of $S_{ij}$. Using these values to bin the data, no change of
the goodness of fit with subhalo triaxiality has been found (see Figure 18 in
the Appendix A).

#### Concentration, Vmax and Rmax

Isolated haloes which are well described by the NFW profile have a
concentration defined as $c_{200}=r_{200}/r_{\mathrm{s}}$ with the virial
radius $r_{200}$ and the scale radius $r_{\mathrm{s}}$. This definition is
less suitable for subhaloes because the virial radius is not really well
defined due to tidal stripping (Moliné et al., 2017; Zavala & Frenk, 2019).
The concentration of a subhalo can be described independently of a density
profile by using the peak circular velocity $V_{\mathrm{max}}$, as:

$c_{V}=\frac{\bar{\rho}(<R_{\mathrm{max}})}{\rho_{\mathrm{c}}}=2\left(\frac{V_{\mathrm{max}}}{HR_{\mathrm{max}}}\right)^{2},$
(11)

where $\bar{\rho}$ is the mean density within the radius $R_{\mathrm{max}}$
corresponding to $V_{\mathrm{max}}$ in units of the critical density
$\rho_{\mathrm{c}}$ and $H$ is the Hubble parameter (Diemand et al., 2007;
Springel et al., 2008). We thus use this definition for the concentration of a
subhalo throughout this paper.

For our new models, the goodness of fit per unit subhalo mass slightly
decreases for $c_{V}>10^{5}$ (see Figure 4). The reason for this can be
inferred from Figure 2 and 5: the goodness of fit value per unit subhalo mass
is a bit smaller for higher-mass subhaloes and very high concentrations occur
primarily for more massive subhaloes. The performance of the new models also
stays relatively constant for different values of $R_{\mathrm{max}}$ and
$V_{\mathrm{max}}$, although it becomes slightly worse for values close to
$R_{\mathrm{max}}=1$ kpc. More details can be found in the appendix in Figure
21.

Figure 5: Subhalo mass vs. concentration with colour-coded baryon fraction for
all of the subhaloes. A small population of highly concentrated, massive and
baryon-rich subhaloes can be identified in the upper half of the diagram.

Figure 5 also shows that the subhaloes with $c_{V}>10^{5}$ generally also have
high baryon fractions $f_{\mathrm{b}}>0.10$ and one can additionally show that
they also have early infall times. This suggests that these high
concentrations are not only facilitated by tidal stripping but also by
baryonic processes, such as adiabatic contraction (Blumenthal et al., 1986;
Gnedin et al., 2004). These more concentrated subhaloes do not necessarily
have an additional bump in the central regions, as the one in Figure 1 (upper
left corner), but the majority of subhaloes with masses larger than $10^{11}$
M☉ has one, and it generally becomes a lot more pronounced for subhaloes with
even higher masses and baryon contents.

#### Infall time

In Figure 20 in the appendix, one can see how the goodness of fit value varies
with infall time. The infall time has been obtained by using the SubLink
merger tree in order to determine the snapshot at which both the subhalo and
its host halo start sharing the same group number. It should be noted, though,
that the subhalo could have also fallen into another halo at an even earlier
time, which later became a subhalo of the host halo at $z=0$. For these sub-
subhaloes the time at which severe tidal interactions happened is even
earlier. We did not consider these cases specifically and only used the infall
time into the main halo at $z=0$ as a first estimate of when tidal
interactions became relevant for each subhalo. One can see that the goodness
of fit does not depend much on the infall time for all models. We also found a
similar behaviour for the individual mass bins.

Figure 6: Subhalo mass dependence of the parameters of the modified NFW
profile from Equation (6) with colour-coded goodness of fit. The solid and
dashed orange lines show the mean and median values of the distribution for
each mass bin, together with the standard deviation. The red solid line
indicates the parameter values of the average density profiles for every mass
bin, together with their uncertainties. The average density profiles can be
found in Figure 10. Figure 7: Dependence of the modified NFW parameters on
$V_{\mathrm{max}}$ for five different mass bins with colour-coded subhalo
concentration. The solid and dashed lines indicate the mean and median of each
distribution for several $V_{\mathrm{max}}$ bins, together with the standard
deviations.

In summary, we were able to show that the new models proposed in Section 3.2
perform much better than the NFW profile and its smoothly truncated version,
both for the individual and for the average density profiles. The Einasto
profile performs similarly well or even better for subhaloes with masses
around $10^{8}$ M⊙. The new models can provide a good fit, largely independent
of concentration, baryon fraction, distance from the host halo centre, group
mass and triaxiality. The model which performs best for most of the subhaloes
is the modified NFW profile from Equation (6). In some regimes, other models
provide a better fit for the average density profiles, but the modified NFW
profile still has a comparable goodness of fit. We will therefore focus the
analysis in Section 5 on the modified NFW profile and only mention the other
models if they enable a better understanding of the density profiles in
general. A more detailed analysis of the other models can be found in the
Appendix A.

##  5 Scaling relations for the Modified NFW Profile and average density
profiles

In this section, we discuss in more detail the new modified NFW profile
introduced in Section 3.2.1, which we find to be the best fit to the simulated
subhaloes. We present the distribution of the fit parameters and derive
scaling relations as a function of the subhalo mass $M_{\mathrm{sh}}$ and
maximum circular velocity $V_{\mathrm{max}}$ that can be used to generate
analytical profiles for a population of subhaloes.

Finally, we discuss how the profile properties vary for individual subhaloes,
depending on their distances from the host halo centre, their concentrations,
their baryon fractions and their infall times, using both the distribution of
the parameter values of the modified NFW profile and the average density
profiles.

###  5.1 Scaling relations

Figure 6 shows how the parameters ($\rho_{0}$,$r_{\mathrm{s}}$,$\alpha$) of
the modified NFW profile (see Equation 6) depend on the subhalo mass. The
parameter $r_{\mathrm{s}}$ correlates strongly with mass, which is not
surprising since it describes the radial extent of the subhalo. $\alpha$ also
increases with mass, with a mean value of 2 for subhaloes with masses of $\sim
10^{8}$ M⊙ up to a value of 3.5 for subhaloes with masses above $10^{12}$ M⊙.

When looking at mass alone, the mean value of the normalisation $\rho_{0}$
stays relatively constant around $10^{4.5}$ M⊙ kpc-3, but with a very large
scatter. Fitting the mean values of the other distributions with a power law,
we obtain the following simple scaling relations:

$r_{\mathrm{s}}=(43.4\pm 3.2)\
\mathrm{kpc}\left(\frac{M_{\mathrm{sh}}}{10^{10}\
\mathrm{M_{\odot}}}\right)^{0.351\pm 0.016},$ (12) $\alpha=(2.647\pm 0.065)\
\left(\frac{M_{\mathrm{sh}}}{10^{10}\ \mathrm{M_{\odot}}}\right)^{0.0439\pm
0.0047}.$ (13)

Figure 8: $V_{\mathrm{max}}$ – subhalo mass relation with colour-coded subhalo
concentration. The fit to the mean values of each mass bin is represented by
the solid black line. The dashed line indicates the result from Springel et
al. (2008).

Figure 6 also shows the parameters obtained by fitting the average density
profile of each mass bin with the Modified NFW (see Figure 10). While these
values agree very well with the those of the individual profiles for the first
two parameters, one can see that $\alpha$ is close to 1.5 for the average
density profiles and therefore significantly differs from the mean values of
the distribution. This is because within each mass bin, subhaloes of various
different sizes go into the averaging process, which effectively creates a
log-slope transition that looks quite different from the ones in most of the
individual subhaloes.

We then obtain a more refined parametrisation by considering both the subhalo
mass and $V_{\mathrm{max}}$. We obtain relations with much less scatter for
both $\rho_{0}$ and $r_{\mathrm{s}}$, which can be seen in Figure 7. One can
see that for each mass bin $\rho_{0}$ increases with $V_{\mathrm{max}}$ while
$r_{\mathrm{s}}$ decreases with $V_{\mathrm{max}}$, which is the expected
behaviour for higher concentrations. The concentrations are colour-coded and
follow the increase in $V_{\mathrm{max}}$. Within a single mass bin, the
parameter $\alpha$ does not change much with $V_{\mathrm{max}}$. That is the
reason why we only describe it in terms of the subhalo mass. Analytic
relations that fit the other mean parameter values are given by:

$\rho_{0}=(4.9\pm
0.6)\mathrm{\frac{M_{\odot}}{kpc^{3}}}\left(\frac{M_{\mathrm{sh}}}{10^{13}\mathrm{M_{\odot}}}\right)^{-1.9\pm
0.1}\left(\frac{V_{\mathrm{max}}}{10^{2}\mathrm{km\ s^{-1}}}\right)^{5.2\pm
0.1}$ (14) $r_{\mathrm{s}}=(10.2\pm
0.1)\mathrm{kpc}\left(\frac{M_{\mathrm{sh}}}{10^{10}\mathrm{M_{\odot}}}\right)^{0.9\pm
0.1}\left(\frac{V_{\mathrm{max}}}{10^{2}\mathrm{km\ s^{-1}}}\right)^{-1.6\pm
0.1}.$ (15)

Besides that, there is also a general relation between $V_{\mathrm{max}}$ and
the subhalo mass. Springel et al. (2008) found:

$V_{\mathrm{max}}=10\,\mathrm{km\
s^{-1}}\left(\frac{M_{\mathrm{sh}}}{3.37\times
10^{7}\,\mathrm{M_{\odot}}}\right)^{0.29},$ (16)

which is close to the relation that we obtain:

$V_{\mathrm{max}}=10\ \mathrm{km\ s^{-1}}\left(\frac{M_{\mathrm{sh}}}{(9.0\pm
1.7)\times 10^{7}\ \mathrm{M_{\odot}}}\right)^{0.343\pm 0.007}.$ (17)

The relation is depicted in Figure 8. One can see that there is a noticeable
scatter due to different subhalo concentrations. $V_{\mathrm{max}}$ increases
with $M_{\mathrm{sh}}$ according to Equation (17) and within each mass bin
$V_{\mathrm{max}}$ increases with concentration (which can also be seen in
Figure 7). This increase in concentration and $V_{\mathrm{max}}$ leads to an
increase of $\rho_{0}$ and a decrease of $r_{\mathrm{s}}$ as described by the
Equations (14) and (15). If we plug Equation (17) into Equation (15) we
approximately obtain back Equation (12).

In the Appendix A, we also included similar plots and analytical expressions
for the scaling relations for the other density profile models. In Figure 9,
we want to additionally highlight how the inner log-slope parameter $\gamma$
of the MS-($\gamma$, 2) profile depends on the subhalo mass. One can see that
the mean values for the inner slope are close to -2 for all subhalo mass bins
and the same applies to the parameter values of the average profiles, although
with an offset to slightly higher values. The scatter becomes larger for
lower-mass subhaloes since the behaviour of the log-slopes changes as
discussed before.

Figure 9: Subhalo mass dependence of the $\gamma$-parameter of the
MS-($\gamma$, 2) profile from Equation (7) with colour-coded goodness of fit
as well as the mean, median and average density profile values for each mass
bin.

Figure 10: The average density profiles for all the mass bins can be seen in
the image on the left. On the right, the average density profiles for every
second mass bin are plotted together with the best modified NFW profile fit.
The residuals of the log-profiles are plotted below. The purple line indicates
the softening length.

###  5.2 Average Density Profiles

In Figure 10, we show the average density profiles for different mass bins,
both with (right image) and without the modified NFW profile fits (left
image). The average profiles have been computed by taking the average of the
density values of each radial bin for all subhaloes in the same mass bin
(including density values of zero). One can see that subhaloes with masses
below $10^{10}$ M⊙ have a gradually decreasing log-slope in the central
regions, while above $10^{11}$ M⊙ the central bump becomes more and more
prominent. In both regimes the density profile has a central core, while for
masses between $10^{10}$ M⊙ and $10^{11}$ M⊙ the log-slope stays constant all
the way to the centre, resulting in a cusp. The modified NFW profile performs
very well in almost all of the cases shown. The only regimes in which it
strongly deviates from the true values are in the central regions of very
massive subhaloes and the outer regions, where the standard deviations are
very large. In the Appendix A in Figure 22 one can see two examples of average
density profiles for given mass bins, together with the best modified NFW
profile fits as well as the combined density uncertainties and standard
deviations.

Figure 11: Left: dependence of the $r_{\mathrm{s}}$-parameter of the modified
NFW profile on host halo centre distance for different mass bins with colour-
coded concentration. The solid and dashed lines indicate the mean and median
of each distribution for different distance bins, together with the standard
deviations. A complete version of this plot can be found in Figure 29 in the
Appendix A. Right: average density profiles for different distances from the
host halo centre for subhaloes with masses between $3.5\times 10^{9}$ M⊙ and
$6.7\times 10^{9}$ M⊙. The purple region indicates the regime below the
softening length.

Within a single mass bin the profile shapes can be quite diverse. In the rest
of this section, we investigate the deviations from these average profiles due
to different subhalo distances from the host halo centre, different
concentrations, baryon fractions and infall times (see Figures 11, 13, 15 and
17).

###  5.3 Dependence of the Parameters on Host Halo Distance

The left panel in Figure 11 shows how the $r_{\mathrm{s}}$-parameter of the
modified NFW profile varies with distance to the host halo centre for
different mass bins. One can see that below the virial radius $R_{200}$ the
value of $r_{\mathrm{s}}$ decreases with distance in a very similar way for
all subhalo mass bins by up to a factor of 10. Beyond the virial radius (i.e.
for subhaloes associated with the FOF group, but located at $r>R_{200}$), the
mean values of $r_{\mathrm{s}}$ stay more or less constant, which is
consistent with the fact that these subhaloes are not yet subjected to tidal
stripping. The colour coding also shows that the subhalo concentration
increases towards the host halo centre. This behaviour is expected, since
subhaloes closer to the centre generally also have earlier infall times and
are therefore more tidally truncated. As $r_{\mathrm{s}}$ decreases,
$\rho_{0}$ increases and $\alpha$ slightly increases as well for all mass bins
in a similar way (see Figure 29 in the Appendix A).

The change in $r_{\mathrm{s}}$ can also be seen in the average profiles. The
right panel in Figure 11 shows the average density profiles for different
distances from the host halo centre in the mass bin between $3.5\times 10^{9}$
M⊙ and $6.7\times 10^{9}$ M⊙. Inside the virial radius, the effect of tidal
stripping clearly depends on distance. Besides that, the inner slope becomes
slightly steeper.

Figure 12: Dependence of the modified NFW parameters on the subhalo
concentration $c_{V}$ for five different mass bins with colour-coded infall
time $a_{\mathrm{in}}$. The solid and dashed lines indicate the mean and
median of each distribution for several $c_{V}$ bins, together with the
standard deviations.

Figure 13: The average density profiles for different subhalo concentration
bins are shown. On the left for subhaloes with masses between $2.6\times
10^{8}$ M⊙ and $5.0\times 10^{8}$ M⊙ and on the right for subhaloes with
masses between $1.7\times 10^{11}$ M⊙ and $3.3\times 10^{11}$ M⊙. The purple
line indicates the softening length.

###  5.4 Dependence of the Parameters on Subhalo Concentration

In Figure 12 the relationships between the modified NFW parameters and the
subhalo concentration are depicted for five different mass bins. The
distributions are not as clearly separated as for $V_{\mathrm{max}}$, they
also have a larger scatter and the mean values show more fluctuations.
Nonetheless, one can see that for higher-mass subhaloes the parameter values
change less with concentration than for lower-mass subhaloes. This can also be
seen in the average density profiles for different concentrations in Figure
13. For lower-mass subhaloes $r_{\mathrm{s}}$ becomes smaller and $\rho_{0}$
becomes larger for higher concentrations, as one would expect. However, for
higher-mass subhaloes $\rho_{0}$ and $r_{\mathrm{s}}$ do not significantly
vary with concentration and the increase in $c_{V}$ mainly comes from an
increase in the strength and radius of the central bump, which is not
parametrised by the model. Exceptions to this can be seen for the last two
bins of exceptionally high concentrations, where there is an additional
truncation besides the prominent central bump. However, concentrations above
$10^{9}$ are extremely rare and only occur for four subhaloes in the sample.

###  5.5 The Effect of Tidal Interactions

Once a subhalo gets sufficiently close to the centre of its host halo, the
tidal forces become strong enough to remove matter from the outer regions:
this starts to happen close to the virial radius $R_{200}$ (see Figure 11).
The subhalo then continually loses orbital energy by dynamical friction and
sinks closer towards the host halo centre. Matter gets removed beyond the
tidal radius $r_{\mathrm{tid}}$, which is defined as the radius where the
tidal force of the host halo is equal to the gravitational force of the
subhalo. Besides tidal stripping in a slowly varying external potential,
subhaloes can also experience tidal shocks if the external potential varies
very rapidly, which happens either close to the main halo centre or during
close encounters with other subhaloes (Zavala & Frenk, 2019). This transfers
orbital energy of the subhalo to internal energy of its particles, which can
alter the subhalo’s internal structure and make some of its particles unbound.
van den Bosch et al. (2018) found that tidal shock heating due to close
encounters with other subhaloes is negligible compared to the tidal effects of
the host halo and that tidal shocks which significantly exceed the subhalo’s
binding energy generally do not lead to its complete disruption.

Figure 14: Dependence of the modified NFW parameters on the infall time
$a_{\mathrm{in}}$ for five different mass bins with colour-coded concentration
$c_{V}$. The solid and dashed lines indicate the mean and median of each
distribution for several $a_{\mathrm{in}}$ bins, together with the standard
deviations.

Figure 15: The average density profiles for different subhalo infall time bins
are shown. On the left for subhaloes with infall masses between $6.9\times
10^{9}$ M⊙ and $1.6\times 10^{10}$ M⊙ and on the right for subhaloes with
masses between $1.9\times 10^{11}$ M⊙ and $4.3\times 10^{11}$ M⊙. The purple
line indicates the softening length.

Because of the expansion of the dark matter particle orbits during tidal
shocks, the inner density can be reduced, but the effect is usually not strong
enough to produce a central core (Hayashi et al., 2003). Besides that, tidal
effects mainly affect the outer regions of the subhalo, which can also be seen
in our average density profiles for different infall times in Figure 15.
However, the central density does not always decrease with earlier infall
times and the shape parameter $\alpha$ stays relatively constant and depends
primarily on the subhalo mass. Other studies, e.g. Kazantzidis et al. (2004)
and Green & van den Bosch (2019), find a decrease of the central density and
steeper truncations for individual subhalo density profiles as tidal stripping
progresses. Figure 14 shows that the parameters $\rho_{0}$ and
$r_{\mathrm{s}}$ on average change in a very similar way for subhaloes with
different masses and that the subhaloes become more concentrated for earlier
infall times. There is, however, also a large number of subhaloes with recent
infall times and a high concentration, which is not unexpected and reflects
the different assembly histories. Subhaloes that formed earlier, when the
Universe had a higher mean density, also have a higher concentration (Wechsler
et al., 2002).

A more sophisticated analysis of the effects of tidal stripping on subhaloes
with different properties at infall time would require following the evolution
of each individual subhalo, which goes beyond the scope of this study.

###  5.6 The Effect of Baryons

Besides tidal truncation due to early infall times or small distances to the
host halo centre, a high fraction of baryons is another indicator for a high
subhalo concentration. We already saw in Figure 5 that a high baryon fraction
is generally present for subhaloes with $c_{V}>10^{6}$. Baryons can cool, sink
to the subhalo centre and gravitationally attract more mass, including dark
matter. This leads to a significant density increase in the central regions,
which could also lead to an increase in the inner log-slope. If the subhalo is
massive and the gravitational potential well therefore deep enough, a large
amount of baryons can accumulate in the centre and form stars. This is also
the reason for the bump feature we see in subhaloes with masses larger than
$10^{11}$ M⊙ as these central bumps are mainly comprised of baryons. However,
the dark matter distribution also gets reshaped due to gravitational
interactions with the baryons, and thus the dark matter profile also shows an
inner log-slope of -2 and can have a small central bump. For a comparison of
the dark matter and baryon component of two subhalo density profiles, see
Figure 30. According to Mollitor et al. (2015), feedback is the reason why the
bump has a flat core in the centre. Tidal torques that the main galaxy disc
exerts on subhaloes can disrupt them and additionally affect the shape of
their density profiles (Garrison-Kimmel et al., 2017; Nadler et al., 2018).

Figure 16: Dependence of the modified NFW parameters on the subhalo baryon
fraction $f_{\mathrm{b}}$ for five different mass bins with colour-coded
concentration $c_{V}$. The solid and dashed lines indicate the mean and median
of each distribution for several $c_{V}$ bins, together with the standard
deviations.

Figure 17: The average density profiles for different subhalo baryon fraction
bins are shown. On the left for subhaloes with masses between $2.6\times
10^{8}$ M⊙ and $5.0\times 10^{8}$ M⊙ and on the right for subhaloes with
masses between $9.0\times 10^{10}$ M⊙ and $1.7\times 10^{11}$ M⊙. The purple
line indicates the softening length.

In Figure 16, the dependence of the modified NFW parameters on baryon fraction
is shown for five different mass bins. We see that $r_{\mathrm{s}}$ decreases
with increasing baryon fraction by up to a factor of 10 for very baryon-
dominated subhaloes. Most of the subhaloes with a baryon fraction larger than
0.4 have a concentration higher than $10^{6}$. Analogously, $\rho_{0}$
increases with baryon fraction in a similar way for all mass bins, leading to
a higher concentration.

Figure 17 shows the average density profiles for different baryon fractions
for subhaloes with masses in the range between $2.6\times 10^{8}$ M⊙ and
$5.0\times 10^{8}$ M⊙ (on the left) and for subhaloes in the mass bin between
$9.0\times 10^{10}$ M⊙ and $1.7\times 10^{11}$ M⊙ (on the right). For less
massive subhaloes $\rho_{0}$ increases and $r_{\mathrm{s}}$ decreases with
baryon fraction, making the subhalo more concentrated, and additionally the
inner log-slope becomes steeper. Here it should be noted, that tidal stripping
removes mass from dark matter dominated outer regions, which also increases
the ratio of baryon mass to total mass. For the more massive subhaloes,
$r_{\mathrm{s}}$ also decreases with baryon fraction and the central bump
becomes more pronounced.

Baryons can also have other effects, such as reducing the number of subhaloes
by facilitating their disruption, especially for small subhalo masses and
small distances from the host halo centre (Despali & Vegetti, 2017; Garrison-
Kimmel et al., 2017).

##  6 Summary & Conclusions

In this paper, we have presented a detailed study of the density profiles of
subhaloes with masses above $1.4\times 10^{8}$ M⊙ in the state-of-the-art
cosmological simulation TNG50. As an improvement to many other simulations
that have been used to study subhalo properties in previous works, TNG50
offers a good compromise between resolution and statistics, includes the
effect of baryons and makes use of a large variety of other physical model
components. The details of the outcomes presented in this paper could differ
for simulations that use different ways to model the hydrodynamics, but
overall, our main findings are mostly consistent with other studies.

We tested the performance of commonly used density profile models, such as the
standard NFW, the Einasto and the truncated NFW profiles and further
introduced new models, which show a much better performance on both the
individual and the average density profiles. Of the models analysed, the best
performing model is a modified version of the NFW profile, given by Equation
(6). The model parameter values follow simple scaling relations given by
Equations (12) - (15). For the lowest-mass subhaloes in the sample with masses
around $10^{8}$ M⊙, the Einasto profile shows a better performance compared to
the new models. For the average density profiles of more massive subhaloes,
the modified Schechter profiles, which model the truncation with an
exponential factor, provide a slightly better fit. The goodness of fit of the
new models is not significantly affected by other subhalo properties,
including baryon fraction, distance from the host halo centre, triaxiality,
concentration and infall time. This suggests that the new models are able to
describe a large variety of subhaloes in the simulation. We have not
considered models with more than three parameters in this study. It can be
expected that the generalised NFW profile given by Equation (5) provides an
even better fit than the three-parameter model we proposed, since it allows
for more variations in the inner log-slope and the exact shape of the
truncation. An investigation of a model that additionally parametrises the
strength and extent of the central bump that is present in the most massive
subhaloes could provide further insights on the effects of baryons on the
density profiles.

The other results of this study can be summarised as follows:

  1. 1.

The inner log-slope of most subhalo density profiles in TNG50 is close to -2,
which is significantly larger than the results obtained in previous studies
(e.g. Springel et al. (2008) and Di Cintio et al. (2013)). For subhaloes with
masses in the range between $10^{10}$ M⊙ and $10^{11}$ M⊙, the inner log-slope
stays relatively constant all the way to the centre, resulting in a central
cusp. Subhaloes with masses above $10^{11}$ M⊙ have an additional bump with a
flat core in the centre that becomes more pronounced for more massive
subhaloes with a higher baryon fraction. For subhalo masses smaller than
$10^{10}$ M⊙, the inner log-slopes continuously become smaller towards the
centre, also resulting in a core. This core can be of physical origin and
might be the result of the reduced amount of baryons, but limitations in the
resolution could also have an impact. The inner log-slope of -2 is also
present when looking only at the dark matter component of the density
profiles, whereas the baryonic component falls of more strongly. The majority
of matter in the central bump comes from the presence of baryons, but a small
central bump can also manifest in the dark matter component due to the
gravitational interactions with the baryons.

  2. 2.

The average density profiles for each mass bin have a nearly universal shape
with a modified NFW shape parameter of $\alpha\approx 1.6$, while the mean
shape parameters of the individual density profiles increase with mass
according to Equation (13). The outer subhalo shape and extent might be
affected by the subhalo finding algorithm. Nevertheless, the analysis of
observations is predominantly influenced by the behavior in the inner regions.

  3. 3.

The model parameters of the modified NFW profile start to change below the
virial radius of the host halo R200. For smaller distances, the parameter
$r_{\mathrm{s}}$ becomes smaller by up to a factor of 10 and the normalisation
$\rho_{0}$ becomes larger by up to a factor $10^{4}$, resulting in a higher
concentration for subhaloes in the same mass bin. This affects subhaloes of
all masses in the same way. For subhaloes closer to the centre, the inner log-
slope of the average density profiles can also increase.

  4. 4.

The parameters of the modified NFW profile vary more strongly with
concentration for lower-mass subhaloes. For subhaloes with masses above
$10^{11}$ M⊙, the parameter values do not change significantly with $c_{V}$
and an increase in concentration is mostly associated with a more pronounced
central bump feature. An additional tidal truncation can lead to exceptionally
high concentrations with $c_{V}>10^{9}$, which are only seen in 4 subhaloes in
the sample.

  5. 5.

Most subhalo concentrations lie in the range $10^{3}\lesssim c_{V}\lesssim
10^{5}$. Subhaloes with $c_{V}>10^{5}$ have in common that they generally have
a much higher baryon fraction ($f_{\mathrm{b}}>0.10$).

  6. 6.

A larger baryon fraction is generally accompanied by a decrease in
$r_{\mathrm{s}}$ by up to a factor of 10 and an increase in $\rho_{0}$,
resulting in a much higher concentration. For subhaloes above $10^{11}$ M⊙,
the baryon fraction also correlates with the size of the central bump.

  7. 7.

Tidal interactions mainly affect the outer regions of the subhalo density
profile, making the subhaloes more concentrated for earlier infall times. The
mean shape parameter $\alpha$ stays relatively constant for different infall
times and depends primarily on the subhalo mass. The average density profiles
for different infall times within different infall mass bins show no
consistent decrease in the central density with earlier infall times. This is
different from the results of other studies on the tidal stripping of
individual subhaloes (e.g. Kazantzidis et al. (2004) and Green & van den Bosch
(2019)), which find a decrease of the central density and steeper truncations
as tidal stripping progresses. However, a more sophisticated analysis of how
the density profiles change over time for different types of subhaloes
requires following their evolution individually.

The NFW and Einasto profiles, together with the concentration-mass relation,
are routinely used to predict the internal structure of dark matter subhaloes
and estimate their observable effects. In turn, the comparisons with
observational data are fundamental to constrain the properties of dark matter.
The new analytical expression for the best-fitting subhalo density profile
model presented in this paper (Equation 6), as well as the scaling relations
between the model parameters and subhalo properties (Equations 12 \- 15)
provide a more realistic description of subhaloes in the presence of baryonic
physics and can thus improve our understanding of the detectability of
substructures and their distribution, leading to more precise constraints on
the dark matter distribution in galaxies. In particular, our results suggest
that NFW profiles are too shallow to correctly represent the inner slope of
hydrodynamical subhaloes and this can lead to a mismodelling of the subhalo
concentrations and central densities.

## Acknowledgements

The authors acknowledge computing resources provided by the Ministry of
Science, Research and the Arts (MWK) of the State of Baden-Württemberg through
bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1
FUGG and 35/1597-1 FUGG. They also thank for data storage at SDS@hd funded
through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG.

The IllustrisTNG simulations were undertaken with compute time awarded by the
Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU
and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High
Performance Computing centre Stuttgart (HLRS), as well as on the machines of
the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.

GD was supported by a Gliese Fellowship.

RSK acknowledges financial support from the European Research Council via the
ERC Synergy Grant “ECOGAL” (project ID 855130), from the German Excellence
Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948)
“STRUCTURES”, and from the German Ministry for Economic Affairs and Climate
Action in project “MAINN” (funding ID 50OO2206).

## Data Availability

The data analysed in this work are available upon request to the corresponding
author. The entire data of the IllustrisTNG simulations, including TNG50, are
publicly available and accessible at www.tng-project.org/data (Nelson et al.,
2019a). This research made use of the public Python packages matplotlib
(Hunter, 2007), NumPy (Harris et al., 2020) and SciPy (Jones et al., 2001).

## References

  * Bailin & Steinmetz (2005) Bailin J., Steinmetz M., 2005, ApJ, 627, 647
  * Behroozi et al. (2015) Behroozi P., et al., 2015, MNRAS, 454, 3020
  * Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
  * Bode et al. (2001) Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93
  * Bonaca et al. (2019) Bonaca A., Hogg D. W., Price-Whelan A. M., Conroy C., 2019, ApJ, 880, 38
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40
  * Brooks & Zolotov (2014) Brooks A. M., Zolotov A., 2014, The Astrophysical Journal, 786, 87 pp.
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Bullock et al. (2001) Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V., Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559
  * Burkert (1995) Burkert A., 1995, ApJ, 447, L25
  * Carlberg (2012) Carlberg R. G., 2012, ApJ, 748, 20
  * Collett (2015) Collett T. E., 2015, ApJ, 811, 20
  * Despali & Vegetti (2017) Despali G., Vegetti S., 2017, MNRAS, 469, 1997
  * Despali et al. (2019) Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F., 2019, MNRAS, 484, 4563
  * Despali et al. (2020) Despali G., Lovell M., Vegetti S., Crain R. A., Oppenheimer B. D., 2020, MNRAS, 491, 1295
  * Despali et al. (2022) Despali G., Walls L. G., Vegetti S., Sparre M., Vogelsberger M., Zavala J., 2022, MNRAS, 516, 4543
  * Di Cintio et al. (2011) Di Cintio A., Knebe A., Libeskind N. I., Yepes G., Gottlöber S., Hoffman Y., 2011, MNRAS, 417, L74
  * Di Cintio et al. (2013) Di Cintio A., Knebe A., Libeskind N. I., Brook C., Yepes G., Gottlöber S., Hoffman Y., 2013, MNRAS, 431, 1220
  * Diemand et al. (2007) Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859
  * Diemer & Kravtsov (2015) Diemer B., Kravtsov A. V., 2015, ApJ, 799, 108
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Dutton & Macciò (2014) Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
  * Einasto (1965) Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
  * Erkal & Belokurov (2015) Erkal D., Belokurov V., 2015, MNRAS, 454, 3542
  * Errani & Navarro (2021) Errani R., Navarro J. F., 2021, MNRAS, 505, 18
  * Famaey & McGaugh (2012) Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10
  * Flores & Primack (1994) Flores R. A., Primack J. R., 1994, ApJ, 427, L1
  * Franx et al. (1991) Franx M., Illingworth G., de Zeeuw T., 1991, ApJ, 383, 112
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel S., et al., 2017, MNRAS, 471, 1709
  * Ghigna et al. (2000) Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 2000, ApJ, 544, 616
  * Gilman et al. (2019) Gilman D., Birrer S., Treu T., Nierenberg A., Benson A., 2019, MNRAS, 487, 5721
  * Gnedin et al. (2004) Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16
  * Green & van den Bosch (2019) Green S. B., van den Bosch F. C., 2019, MNRAS, 490, 2091
  * Green et al. (2021) Green S. B., van den Bosch F. C., Jiang F., 2021, MNRAS, 503, 4075
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hayashi et al. (2003) Hayashi E., Navarro J., Taylor J., Stadel J., Quinn T., 2003, ApJ, 584, 541
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hezaveh et al. (2016) Hezaveh Y. D., et al., 2016, ApJ, 823, 37
  * Hsueh et al. (2020) Hsueh J. W., Enzi W., Vegetti S., Auger M. W., Fassnacht C. D., Despali G., Koopmans L. V. E., McKean J. P., 2020, MNRAS, 492, 3047
  * Huchra & Geller (1982) Huchra J. P., Geller M. J., 1982, ApJ, 257, 423
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Jaffe (1983) Jaffe W., 1983, MNRAS, 202, 995
  * Jones et al. (2001) Jones E., Oliphant T., Peterson P., 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org
  * Kazantzidis et al. (2004) Kazantzidis S., Mayer L., Mastropietro C., Diemand J., Stadel J., Moore B., 2004, ApJ, 608, 663
  * Klypin et al. (1999a) Klypin A., Gottlöber S., Kravtsov A., Khokhlov A., 1999a, ApJ, 516, 530
  * Klypin et al. (1999b) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999b, ApJ, 522, 82
  * Kroupa (2012) Kroupa P., 2012, Publications of the Astronomical Society of Australia, 29, 395
  * Lovell & Zavala (2023) Lovell M. R., Zavala J., 2023, MNRAS, 520, 1567
  * Lovell et al. (2014) Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014, MNRAS, 439, 300
  * Mashchenko et al. (2008) Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174 pp.
  * Mastromarino et al. (2023) Mastromarino C., Despali G., Moscardini L., Robertson A., Meneghetti M., Maturi M., 2023, MNRAS, 524, 1515
  * Mateo (1998) Mateo M. L., 1998, ARA&A, 36, 435
  * Merritt et al. (2006) Merritt D., Graham A. W., Moore B., Diemand J., Terzić B., 2006, AJ, 132, 2685
  * Milgrom (1983) Milgrom M., 1983, ApJ, 270, 365
  * Milgrom (2002) Milgrom M., 2002, New Astron. Rev., 46, 741
  * Minor et al. (2021) Minor Q., Gad-Nasr S., Kaplinghat M., Vegetti S., 2021, MNRAS, 507, 1662
  * Moliné et al. (2017) Moliné Á., Sánchez-Conde M. A., Palomares-Ruiz S., Prada F., 2017, MNRAS, 466, 4974
  * Moliné et al. (2023) Moliné Á., et al., 2023, MNRAS, 518, 157
  * Mollitor et al. (2015) Mollitor P., Nezri E., Teyssier R., 2015, MNRAS, 447, 1353
  * Moore (1994) Moore B., 1994, Nature, 370, 629
  * Moore et al. (1999a) Moore B., Quinn T., Governato F., Stadel J., Lake G., 1999a, MNRAS, 310, 1147
  * Moore et al. (1999b) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999b, ApJ, 524, L19
  * Muñoz et al. (2001) Muñoz J. A., Kochanek C. S., Keeton C. R., 2001, ApJ, 558, 657
  * Nadler et al. (2018) Nadler E. O., Mao Y.-Y., Wechsler R. H., Garrison-Kimmel S., Wetzel A., 2018, ApJ, 859, 129
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Navarro et al. (2010) Navarro J. F., et al., 2010, MNRAS, 402, 21
  * Nelson et al. (2019a) Nelson D., et al., 2019a, Computational Astrophysics and Cosmology, 6, 2
  * Nelson et al. (2019b) Nelson D., et al., 2019b, MNRAS, 490, 3234
  * Onions et al. (2012) Onions J., et al., 2012, MNRAS, 423, 1200
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Pillepich et al. (2019) Pillepich A., et al., 2019, MNRAS, 490, 3196
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A1
  * Pontzen & Governato (2012) Pontzen A., Governato F., 2012, Monthly Notices of the Royal Astronomical Society, 421, 3464
  * Read et al. (2016) Read J. I., Agertz O., Collins M. L. M., 2016, Monthly Notices of the Royal Astronomical Society, 459, 2573
  * Robles et al. (2019) Robles V. H., Kelley T., Bullock J. S., Kaplinghat M., 2019, MNRAS, 490, 2117
  * Rocha et al. (2013) Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Sanders & McGaugh (2002) Sanders R. H., McGaugh S. S., 2002, ARA&A, 40, 263
  * Sawala et al. (2016) Sawala T., Frenk C. S., Fattahi et al. A., 2016, Monthly Notices of the Royal Astronomical Society, 457, 1931
  * Schechter (1976) Schechter P., 1976, ApJ, 203, 297
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Springel et al. (2008) Springel V., et al., 2008, MNRAS, 391, 1685
  * Tolstoy et al. (2009) Tolstoy E., Hill V., Tosi M., 2009, ARA&A, 47, 371
  * Vegetti et al. (2010a) Vegetti S., Czoske O., Koopmans L. V. E., 2010a, MNRAS, 407, 225
  * Vegetti et al. (2010b) Vegetti S., Koopmans L. V. E., Bolton A., Treu T., Gavazzi R., 2010b, MNRAS, 408, 1969
  * Vegetti et al. (2012) Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341
  * Vogelsberger et al. (2012) Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, The Astrophysical Journal Letters, 827, 6 pp.
  * Zavala & Frenk (2019) Zavala J., Frenk C. S., 2019, Galaxies, 7, 81
  * Zemp et al. (2011) Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2011, ApJS, 197, 30
  * Zhao (1996) Zhao H., 1996, MNRAS, 278, 488
  * Zhao et al. (2003) Zhao D. H., Jing Y. P., Mo H. J., Börner G., 2003, ApJ, 597, L9
  * Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Börner G., 2009, ApJ, 707, 354
  * Zhu et al. (2017) Zhu Q., Hernquist L., Marinacci F., Springel V., Li Y., 2017, MNRAS, 466, 3876
  * Şengül & Dvorkin (2022) Şengül A. Ç., Dvorkin C., 2022, MNRAS, 516, 336
  * van den Bosch & Ogiya (2018) van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066
  * van den Bosch et al. (2018) van den Bosch F. C., Ogiya G., Hahn O., Burkert A., 2018, MNRAS, 474, 3043

##  Appendix A Extra Material

Figures 18, 19, 20 and 21 depict how the density profile models perform for
different subhalo triaxialities, group masses, infall times,
$R_{\mathrm{max}}$ and $V_{\mathrm{max}}$.

Figure 22 shows two example average density profiles for given mass bins,
together with the best modified NFW profile fit.

Figures 23, 24 and 25 show the mass dependence of the parameter values for the
Einasto profile, the truncated power law and the modified Schechter profile.
Figures 26, 27 and 28 show the dependence of the parameter values on
$V_{\mathrm{max}}$ for different mass bins for the Einasto profile, the
truncated power law and the modified Schechter profile. One can see that for
all of the diagrams the overlap and scatter of the distributions is
significantly larger than for the modified NFW profile (Figures 6 and 7),
which is another thing that makes it the best model of the ones we have
analysed. We didn’t show the distributions for the NFW profile and its
smoothly truncated version, since the parameters take extreme values and do
not show nice scaling relations.

Figure 29 is the complete version of Figure 11 and shows how the parameter
values of the modified NFW profile depend on the distance from the host halo
centre.

Figure 30 shows the density profiles of the first two subhaloes in Figure 1
decomposed into their dark matter and baryon components.

Table 1 lists the equations for the scaling relations for the parameter values
of the Einasto profile, the truncated power law and the modified Schechter
profile.

Figure 18: Median values of the goodness of fit per unit subhalo mass for
subhaloes with different triaxialities. Additionally, the distribution of
triaxialities is shown in the histogram. Figure 19: Median values of the
goodness of fit per unit subhalo mass for subhaloes with different group
masses. Additionally, the distribution of group masses is shown in the
histogram. Figure 20: Median values of the goodness of fit per unit subhalo
mass for subhaloes with different infall times. Additionally, the distribution
of infall times is shown in the histogram.

Figure 21: Median values of the goodness of fit per unit subhalo mass for
subhaloes with different $R_{\mathrm{max}}$ (left) and $V_{\mathrm{max}}$
(right). Additionally, the distributions of $R_{\mathrm{max}}$ and
$V_{\mathrm{max}}$ are shown in the histograms.

Figure 22: Average density profiles with standard deviations for subhaloes with masses between $5.0\times 10^{8}$ M⊙ and $9.6\times 10^{8}$ M⊙ (left) as well as between $6.7\times 10^{9}$ M⊙ and $1.3\times 10^{10}$ M⊙ (right), together with the best modified NFW profile fits. The purple area indicates the regime below the softening length. Table 1: Scaling relations for the parameter values of the Einasto profile and the modified Schechter profiles. Density Profile Model | Normalisation Parameter | Radius Parameter | Shape Parameter
---|---|---|---
| [$\,\mathrm{M_{\sun}\ kpc^{-3}}$] | [kpc] |
Einasto Profile | $\rho_{-2}=930\ \left(\frac{M}{10^{13}\,\mathrm{M_{\odot}}}\right)^{-1.6}\left(\frac{V_{\mathrm{max}}}{10^{2}\,\mathrm{km\ s^{-1}}}\right)^{4.8}$ | $r_{-2}=1.6\ \left(\frac{M}{10^{10}\,\mathrm{M_{\odot}}}\right)^{0.9}\left(\frac{V_{\mathrm{max}}}{10^{2}\,\mathrm{km\ s^{-1}}}\right)^{-1.6}$ | $\alpha=0.64\ \left(\frac{M}{10^{10}\,\mathrm{M_{\odot}}}\right)^{0.047}$
MS-(2, $n$) | $\rho_{0}=1.2\ \left(\frac{M}{10^{13}\,\mathrm{M_{\odot}}}\right)^{-2.2}\left(\frac{V_{\mathrm{max}}}{10^{2}\,\mathrm{km\ s^{-1}}}\right)^{5.6}$ | $r_{\mathrm{t}}=4.4\ \left(\frac{M}{10^{10}\,\mathrm{M_{\odot}}}\right)^{0.9}\left(\frac{V_{\mathrm{max}}}{10^{2}\,\mathrm{km\ s^{-1}}}\right)^{-1.8}$ | $n=2.2\ \left(\frac{M}{10^{10}\,\mathrm{M_{\odot}}}\right)^{0.046}$
MS-($\gamma$, 2) | $\rho_{0}=13\ \left(\frac{M}{10^{13}\,\mathrm{M_{\odot}}}\right)^{-2.0}\left(\frac{V_{\mathrm{max}}}{10^{2}\,\mathrm{km\ s^{-1}}}\right)^{5.9}$ | $r_{\mathrm{t}}=5.5\ \left(\frac{M}{10^{10}\,\mathrm{M_{\odot}}}\right)^{0.9}\left(\frac{V_{\mathrm{max}}}{10^{2}\,\mathrm{km\ s^{-1}}}\right)^{-1.5}$ | $\gamma=1.94\ \left(\frac{M}{10^{10}\,\mathrm{M_{\odot}}}\right)^{0.001}$

Figure 23: Subhalo mass dependence of the parameters of the Einasto profile
with colour-coded goodness of fit. The solid and dashed orange lines show the
mean and median values of the distribution for each mass bin, together with
the standard deviation. The red solid line shows the parameter values of the
average density profiles for every mass bin, together with their
uncertainties. Figure 24: Subhalo mass dependence of the parameters of the
MS-(2, $n$) profile with colour-coded goodness of fit. The solid and dashed
orange lines show the mean and median values of the distribution for each mass
bin, together with the standard deviation. The red solid line shows the
parameter values of the average density profiles for every mass bin, together
with their uncertainties. Figure 25: Subhalo mass dependence of the parameters
of the MS-($\gamma$, 2) profile with colour-coded goodness of fit. The solid
and dashed orange lines show the mean and median values of the distribution
for each mass bin, together with the standard deviation. The red solid line
shows the parameter values of the average density profiles for every mass bin,
together with their uncertainties. Figure 26: Dependence of the Einasto
profile parameter values on $V_{\mathrm{max}}$ for five different mass bins
with colour-coded subhalo concentration. The solid and dashed lines indicate
the mean and median of each distribution for several $V_{\mathrm{max}}$ bins,
together with the standard deviations. Figure 27: Dependence of the MS-(2,
$n$) profile parameter values on $V_{\mathrm{max}}$ for five different mass
bins with colour-coded subhalo concentration. The solid and dashed lines
indicate the mean and median of each distribution for several
$V_{\mathrm{max}}$ bins, together with the standard deviations. Figure 28:
Dependence of the MS-($\gamma$, 2) parameter values on $V_{\mathrm{max}}$ for
five different mass bins with colour-coded subhalo concentration. The solid
and dashed lines indicate the mean and median of each distribution for several
$V_{\mathrm{max}}$ bins, together with the standard deviations. Figure 29:
Complete version of Figure 11, which shows the dependence of the modified NFW
parameter values on the subhalo distance from the host halo centre for five
different mass bins with colour-coded concentration $c_{V}$. The solid and
dashed lines indicate the mean and median of each distribution for different
distance bins together with the standard deviation.

Figure 30: The density profiles of the first two subhaloes from Figure 1
decomposed into their dark matter (orange) and baryon components (blue). The
purple area indicates the regime below the softening length. For these density
profiles we did not compute data points below the softening length.