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_\mathrm{sh} \geq 1.4 \times 10^8$ M$_\odot$ 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$_\odot$ $\leq M_\mathrm{sh} \leq$ $10^{13}$ M$_\odot$ that can be influenced by the presence of baryons.


INTRODUCTION
The Λ Cold Dark Matter (Λ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 Λ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 -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 Λ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. 2019Despali et al. , 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(Vegetti et al. ,b, 2012;;Hezaveh et al. 2016), with masses between 1.9 × 10 8 and 2.7 × 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 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 Λ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. (2017Moliné et al. ( , 2023) ) recently investigated concentrations and other subhalo properties in great detail in -body simulations for a wide range of different subhalo and host halo masses and how they evolve over cosmic time.Furthermore, -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(Navarro et al. , 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, V max , 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 V max , and discuss how baryons and tidal stripping affect the subhalo density profiles.Finally, our summary and conclusions are presented in Section 6.

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  = 127 and finishes at the present day  = 0, with cosmologically motivated initial conditions and a cosmology in agreement with Planck 2015, given by Ω Λ,0 = 0.6911, Ω m,0 = 0.3089, Ω b,0 = 0.0486,  8 = 0.8159,  s = 0.9667 and ℎ = 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 × 10 5 M ⊙ ℎ −1 , a baryonic mass resolution of 5.7 × 10 4 M ⊙ ℎ −1 and a gravitational softening length of 0.288 kpc at  = 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.

Subhalo Selection
The TNG50-1 group catalogue of the snapshot at  = 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 10 10 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 × 10 8 M ⊙ and 8.5 × 10 12 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.

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 10 3 ckpc/ℎ.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. Navarro et al. (1996Navarro et al. ( , 1997) ) demonstrated in their -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:

The NFW Profile
(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  s , which is the only free parameter besides the overall normalisation  0 .Therefore, the profile is completely characterised by the virial mass  200 and the concentration  200 =  200 / 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. 2003Zhao et al. , 2009;;Dutton & Macciò 2014;Diemer & Kravtsov 2015) and for subhaloes (Bullock et al. 2001;Moliné et al. 2017Moliné et al. , 2023)).

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  t : (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.

The Einasto Profile
The Einasto profile (Einasto 1965) is another commonly used threeparameter model for describing the density profiles of dark matter haloes.Its functional form is given by: (3)  −2 is the radius at which the log-slope is equal to -2,  −2 is the corresponding density and  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: 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).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.

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.

The Modified NFW Profile
A generalisation of the NFW profile is the so-called (, , ) model (see e.g.Zhao 1996), which has the following functional form: Here  and  set the inner and outer log-slope and  controls the sharpness of the transition from one slope to the other.At the scale radius  s , the log-slope is the average of the inner and outer slopes, which is −( + )/2.For (, , ) = (1, 3, 1) we get the NFW profile, but also other commonly used profiles can be obtained from it.
The (, , ) 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: Here,  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  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  s .
Throughout the paper, we refer to the model in Equation ( 6) as the modified NFW profile.

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  that controls the sharpness of the truncation: This looks very similar to the functional form of the Schechter luminosity function (Schechter 1976) with the additional parameter .
Because of that, we will refer to it as the modified Schechter (MS) profile.Here  0 is again the overall normalisation factor,  t is the truncation radius and  is the power law log-slope.From this fourparameter 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  and the other parameters free.From now on, we will refer to this model as MS-(2, ).
For the other model we instead fix the sharpness of the truncation  to 2 and leave  as a free parameter.From now on, we will refer to this model as MS-(, 2).Our analysis of the parameters of the MS-(2, ) model in the appendix (see Figure A7, A10 and Table A1), suggests that the mean values of  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.

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: with the total number of nonzero data points , the number of fitted parameters , the difference between the log-values of the density   and the model fit value   , and the uncertainties of the log-profile   /  at the data point .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  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  2 is roughly proportional to  sh .Because of that, we primarily look at  2 per unit subhalo mass in order to have a quantity that is relatively independent of mass.Otherwise, high values of  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 logdensity 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).

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  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  sh ∼ 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, ) and the MS-(, 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, ) 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-(, 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 × 10 8 M ⊙ ℎ −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.

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 R 200 ).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 R 200 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  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 R 200 .

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  max and  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).

Triaxiality
We define the subhalo triaxiality as (Franx et al. 1991): where ,  and  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): with the mass   and the the -th and -th coordinate  , and  ,  of particle .The axes ,  and  are the square roots of the eigenvalues of    .Using these values to bin the data, no change of the goodness of fit with subhalo triaxiality has been found (see Figure A1 in the Appendix A).

Concentration, V max and R max
Isolated haloes which are well described by the NFW profile have a concentration defined as  200 =  200 / s with the virial radius  200 and the scale radius  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  max , as: where ρ is the mean density within the radius  max corresponding to  max in units of the critical density  c and  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   > 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  max and  max , although it becomes slightly worse for values close to  max = 1 kpc.More details can be found in the appendix in Figure A4.
Figure 5 also shows that the subhaloes with   > 10 5 generally also have high baryon fractions  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 A3 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  = 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  = 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.
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.

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  sh and maximum circular velocity  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.

Scaling relations
Figure 6 shows how the parameters ( 0 , s ,) of the modified NFW profile (see Equation 6) depend on the subhalo mass.The parameter  s correlates strongly with mass, which is not surprising since it describes the radial extent of the subhalo. also increases with mass, with a mean value of 2 for subhaloes with masses of ∼ 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  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: = (2.647± 0.065)  sh 10 10 M ⊙ 0.0439±0.0047 .
(13) 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  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  max .We obtain relations with much less scatter for both  0 and  s , which can be seen in Figure 7.One can see that for each mass bin  0 increases with  max while  s decreases with  max , which is the expected behaviour for higher concentrations.The concentrations are colour-coded and follow the increase in  max .Within a single mass bin, the parameter  does not change much with  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: Besides that, there is also a general relation between  max and the subhalo mass.Springel et al. (2008) found: which is close to the relation that we obtain: The relation is depicted in Figure 8.One can see that there is a noticeable scatter due to different subhalo concentrations. max increases with  sh according to Equation ( 17) and within each mass bin  max increases with concentration (which can also be seen in Figure 7).This increase in concentration and  max leads to an increase of  0 and a decrease of  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  of the MS-(, 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.

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 A5 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.
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). 10 11 10 11

Dependence of the Parameters on Host Halo Distance
The left panel in Figure 11 shows how the  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  200 the value of  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  >  200 ), the mean values of  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  s decreases,  0 increases and  slightly increases as well for all mass bins in a similar way (see Figure A12 in the Appendix A).
The change in  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 × 10 9 M ⊙ and 6.7 × 10 9 M ⊙ .Inside the virial radius, the effect of tidal stripping clearly depends on distance.Besides that, the inner slope becomes slightly steeper.

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  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  s becomes smaller and  0 becomes larger for higher concentrations, as one would expect.However, for higher-mass subhaloes  0 and  s do not significantly vary with concentration and the increase in   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.

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  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  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.
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  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  0 and  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.

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   > 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 A13.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).
In Figure 16, the dependence of the modified NFW parameters on baryon fraction is shown for five different mass bins.We see that  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,  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 × 10 8 M ⊙ and 5.0 × 10 8 M ⊙ (on the left) and for subhaloes in the mass bin between 9.0 × 10 10 M ⊙ and 1.7 × 10 11 M ⊙ (on the right).For less massive subhaloes  0 increases and  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,  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).

SUMMARY & CONCLUSIONS
In this paper, we have presented a detailed study of the density profiles of subhaloes with masses above 1.4 × 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: (i) 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.
(ii) The average density profiles for each mass bin have a nearly universal shape with a modified NFW shape parameter of  ≈ 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.
(iii) The model parameters of the modified NFW profile start to change below the virial radius of the host halo R 200 .For smaller distances, the parameter  s becomes smaller by up to a factor of 10 and the normalisation  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.
(iv) 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   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   > 10 9 , which are only seen in 4 subhaloes in the sample.
(v) Most subhalo concentrations lie in the range 10 3 ≲   ≲ 10 5 .Subhaloes with   > 10 5 have in common that they generally have a much higher baryon fraction (  b > 0.10).
(vi) A larger baryon fraction is generally accompanied by a decrease in  s by up to a factor of 10 and an increase in  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.
(vii) 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  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 concentrationmass 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.show the dependence of the parameter values on  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.This paper has been typeset from a T E X/L A T E X file prepared by the author.
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.

Figure 2 .
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  2 has been omitted, since it is dominated by failing fits.

Figure 3 .
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.

Figure 4 .
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.

Figure 5 .
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 6 .
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 .
Figure 7. Dependence of the modified NFW parameters on  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  max bins, together with the standard deviations.

Figure 8 .
Figure 8.  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 9 .
Figure 9. Subhalo mass dependence of the -parameter of the MS-(, 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 .
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.
Figure 11.Left: dependence of the  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 A12 in the Appendix A. Right: average density profiles for different distances from the host halo centre for subhaloes with masses between 3.5 × 10 9 M ⊙ and 6.7 × 10 9 M ⊙ .The purple region indicates the regime below the softening length.

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

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

Figure 14 .Figure 15 .
Figure 14.Dependence of the modified NFW parameters on the infall time  in for five different mass bins with colour-coded concentration   .The solid and dashed lines indicate the mean and median of each distribution for several  in bins, together with the standard deviations.

Figure 16 .Figure 17 .
Figure 16.Dependence of the modified NFW parameters on the subhalo baryon fraction  b for five different mass bins with colour-coded concentration   .The solid and dashed lines indicate the mean and median of each distribution for several   bins, together with the standard deviations.

Figure A1 .
Figure A1.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
Figure A5shows two example average density profiles for given mass bins, together with the best modified NFW profile fit.

Figure
Figure A12 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.

FigureFigure A2 .Figure A3 .
FigureA13shows the density profiles of the first two subhaloes in Figure1decomposed into their dark matter and baryon components.

Figure A6 .
Figure A6.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 A7 .
Figure A7.Subhalo mass dependence of the parameters of the MS-(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 A8 .
Figure A8.Subhalo mass dependence of the parameters of the MS-(, 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 A9 .
Figure A9.Dependence of the Einasto profile parameter values on  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  max bins, together with the standard deviations.

Figure A10 .
Figure A10.Dependence of the MS-(2, ) profile parameter values on  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  max bins, together with the standard deviations.

Figure A11 .
Figure A11.Dependence of the MS-(, 2) parameter values on  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  max bins, together with the standard deviations.

Figure A12 .Figure A13 .
Figure A12.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   .The solid and dashed lines indicate the mean and median of each distribution for different distance bins together with the standard deviation.