# An Accurate Comprehensive Approach to Substructure:
III. Masses and Formation Times of the Host Haloes

Eduard Salvador-Solé1, Alberto Manrique1, David Canales2 and Ignacio
Botella1,3
1Institut de Ciències del Cosmos, Universitat de Barcelona, Martí i Franquès
1, E-08028 Barcelona, Spain
2School of Aeronautics and Astronautics, Purdue University, 701 W. Stadium
Ave., West Lafayette, IN 47907-2045, USA
3Dept. of Astronomy, Graduate School of Science, Kyoto University,
Kitashirakawa, Oiwakecho, Sakyo-ku, Kyoto, 606-8502, Japan E-mail:
<EMAIL_ADDRESS>
###### Abstract

With this Paper we complete a comprehensive study of substructure in dark
matter haloes. In Paper I we derived the radial distribution and mass function
(MF) of accreted subhaloes (scaled to the radius and mass of the host halo)
and showed they are essentially universal. This is not the case, however, for
those of stripped subhaloes, which depend on halo mass and assembly history.
In Paper II we derived these latter properties in the simplest case of purely
accreting haloes. Here we extend the study to ordinary haloes having suffered
major mergers. After showing that all the properties of substructure are
encoded in the mean truncated-to-original subhalo mass ratio profile, we
demonstrate that the dependence of the subhalo MF on halo mass arises from
their mass-dependent concentration, while the shape of the subhalo radial
distribution depends on the time of the last major merger of the host halo. In
this sense, the latter property is a better probe of halo formation time than
the former. Unfortunately, this is not the case for the radial distribution of
satellites as this profile is essentially disconnected from subhalo stripping
and the properties of accreted subhaloes are independent of the halo formation
time.

###### keywords:

methods: analytic — galaxies: haloes, substructure — cosmology: theory, dark
matter — dark matter: haloes — haloes: substructure

##  1 INTRODUCTION

In the last decade, the subject of halo assembly history has attracted much
attention in connection with the so-called “missing satellite problem”, namely
that the abundance of satellite galaxies in the Milky Way (MW) (and Andromeda;
Tollerud, Boylan-Kolchin, & Bullock 2014) does not seem to conform with the
expected one in the favourite $\Lambda$CDM cosmology (Moore et al., 1999;
Klypin et al., 1999; Bullock & Boylan-Kolchin, 2017). Indeed, one possible
explanation for that problem is that the MW halo may not have had the typical
assembly history of haloes of its mass. Some aspects of the MW suggest,
indeed, that it has had a particularly quiescent history (e.g. Wyse 2001;
Deason et al. 2013; Ruchti et al. 2015; Lancaster, Belokurov, & Evans 2019).

Substructure in dark matter haloes is believed to harbour important
information on their assembly history as it is the direct consequence of the
way they have grown. Indeed, haloes undergo long periods of frequent minor
mergers (generically called accretion), separated by sporadic major mergers.
The difference between these two kinds of mergers is that the largest halo in
minor mergers (the accreting object) is much more massive than its partners
(the accreted objects), so it remains essentially in equilibrium during that
process and the accreted haloes survive within their host as subhaloes. On the
contrary, all haloes (usually two) partaking of a major merger are similarly
massive so that the event causes them to go out of equilibrium and to form a
new virialised halo after the system relaxes again. In this case the
progenitor haloes are thus destroyed, but their subhaloes are transferred to
the newborn halo. Subhaloes thus accumulate within haloes, where they are more
or less stripped and shock-heated depending on the characteristics of each
transient host. Consequently, their final properties are the result of their
past history.

Unfortunately, we do not know yet what are the typical properties of
substructure in haloes of different masses and formation times. None of the
analytic models of halo substructure so far developed (Taylor & Babul 2001;
Zentner & Bullock 2003; Sheth 2003; Lee 2004; Oguri & Lee 2004; Taylor & Babul
2004; Peñarrubia & Benson 2005; van den Bosch et al. 2005; Zentner et al.
2005; Kampakoglou & Benson 2007; Giocoli et al. 2008; Benson et al. 2013;
Pullen et al. 2014; Jiang & van den Bosch 2016; Griffen et al. 2016; van den
Bosch & Jiang 2016; van den Bosch et al. 2018; Green & van den Bosch 2019;
Font et al. 2020) have been able to provide definite answers to these
questions. Nor have high-resolution simulations, which have only been able to
draw the properties of substructure in a handful of haloes of the MW-mass
(Diemand et al., 2007; Springel et al., 2008; Boylan-Kolchin et al., 2010;
Wang et al., 2011) or of other masses (e.g. Angulo et al. 2009; Elahi et al.
2009; Giocoli et al. 2010; Klypin et al. 2011; Gao et al. 2011, 2012; Onions
et al. 2012; Lovell et al. 2014; Cautun et al. 2014; Ishiyama et al. 2020;
Lovell et al. 2021).

But things are rapidly changing. The incoming new data gathered by means of
the Gaia satellite (Gaia Collaboration et al., 2018) will allow to accurately
determine the MW’s substructure and it has now become feasible to observe
satellites in neighbouring MW analogues (Danieli et al., 2017; Geha et al.,
2017; Smercina et al., 2018; Bennet et al., 2019, 2020; Crnojević et al.,
2019; Carlsten et al., 2020; Mao et al., 2020; Carlsten et al., 2021). In
parallel, simulations have also greatly improved. The Copernicus Complexio
$N$-body simulations together with a semi-analytic galaxy formation model
(Hellwing et al. 2016; Bose et al. 2016, 2020) or the Apostle and Auriga
(Richings et al., 2020), FIRE-2 (Samuel et al., 2020), and Artemis (Font et
al., 2020; Font, McCarthy, & Belokurov, 2021; Engler et al., 2021)
hydrodynamic simulations have gathered a considerable number of simulated MW
(and Andromeda) analogs. On the other hand, it is now possible to reach very
high stellar mass resolutions (Grand et al., 2021) which significantly
improves the statistics of substructure at the level of faint and ultra-faint
satellites ($10^{2}$ M⊙ $<M_{\star}<10^{6}$ M⊙).

In addition, great progress has also been made on analytic grounds. A very
complete model of substructure formation has been built (Jiang et al., 2021)
that allows one to study the effects of different initial conditions in the
accretion and evolution of satellite galaxies. Likewise, the powerful
ConflUent System of Peak trajectories (CUSP) formalism (Manrique & Salvador-
Solé, 1995, 1996; Manrique et al., 1998; Salvador-Solé et al., 2012a, b; Juan
et al., 2014a, b), making the link between the properties of haloes and their
seeds (peaks) in the random Gaussian field of density perturbations (see
Salvador-Solé & Manrique 2021 for an overview) has been successfully applied
to the study of halo substructure (Salvador-Solé et al. 2021a and Salvador-
Solé et al. 2021b, hereafter Papers I and II, respectively).

In Paper I we determined the properties of accreted subhaloes, which act as
initial conditions in their evolution through stripping inside their host
haloes. To do this we took profit of the fact that, as shown in Salvador-Solé
et al. (2021c), all halo properties arising from their gravitational
clustering process do not depend on their particular assembly history so that
one has the right to focus on the simplest case of purely accreting haloes
evolving inside-out (see below for a brief explanation of this important
result). However, subhalo stripping as well as dynamical friction are two
(coupled) mechanisms acting on the dynamical evolution of subhaloes that are
not directly connected to gravitational clustering. Consequently, the
properties of stripped subhaloes do depend on the halo assembly history. In
Paper II we built a detailed model of tidal stripping and shock heating of
subhaloes as they orbit inside haloes which allowed us to derive their final
properties. However, this was done in the simplest case of purely accreting
haloes only and neglecting dynamical friction.

In the present Paper, we complete this study and extend the characterisation
of substructure to ordinary haloes, i.e. haloes having suffered major mergers,
paying special attention to the role of halo mass and formation time and
analysing the possible use of those properties as a probe for halo assembly
history. Our treatment does not include dynamical friction. However, the
predictions for low-mass subhaloes ($M_{\rm s}\la 10^{-4}M_{\rm h}$, where
$M_{\rm h}$ is the mass of the host halo) should not be affected by that
omission. On the other hand, it does not include baryons either. However, by
comparing our predictions to the results of simulations including them, it is
still possible to unravel to some extent the influence of baryons physics in
the properties of substructure.

The layout of the Paper is as follows. In Section 2 we remind the main results
of Papers I and II for purely accreting haloes of different masses. In Section
3 we extend those results to haloes having suffered major mergers. And in
Section 4 we analyse the properties of substructure in haloes of a fixed mass
and different formation times. Our results are summarised and discussed in
Section 5.

Some comments on the notation used in this Paper are in order. Unless
otherwise stated, when we refer to subhaloes without specifying their kind, we
mean stripped (or truncated) subhaloes. The halo formation time used is
defined as the time they suffered their last major merger. It thus differs
from the most usual definition: the time haloes reach 50% of their final mass.
We prefer the former not only because it is less arbitrary (why 50% and not,
say, 75%?), but also because it is physically better motivated. Indeed, as
mentioned, virialised haloes interrupt their identity in major mergers where
they are destroyed and a new virialised object appears (Salvador-Solé &
Manrique, 2021). Lastly, the notation we use for the cumulative or
differential abundances of subhaloes, dependent in general on subhalo mass,
$M_{\rm s}$, and radial location inside the host halo, $r$, is the same as in
Papers I and II. The cumulative number of stripped or accreted haloes out to
$r$ and down to $M_{\rm s}$ are denoted as ${\cal N}^{\rm stp}(<r,>M_{\rm s})$
and ${\cal N}^{\rm acc}(<r,>M_{\rm s})$, respectively. When one of the
arguments takes its maximum value, i.e. when the integrals over $r$ or $M_{\rm
s}$ are complete, we drop the corresponding argument. For instance, ${\cal
N}^{\rm stp}(>M_{\rm s})$ stands for ${\cal N}^{\rm stp}(<R_{\rm h},>M_{\rm
s})$, where $R_{\rm h}$ is the total halo radius, and ${\cal N}^{\rm stp}(<r)$
stands for ${\cal N}^{\rm stp}(<r,>0)$ (or for ${\cal N}^{\rm stp}(<r,<M_{\rm
h})$, where $M_{\rm h}$ is the total halo mass). Lastly, the differential form
with respect to any argument, $r$ or $M_{\rm s}$, of any of the previous
functions is denoted without the corresponding preceding inequality symbol.
For example, ${\cal N}^{\rm stp}(r,M_{\rm s})$ stands for the double
derivative with respect to $r$ and $M_{\rm s}$ and ${\cal N}^{\rm stp}(M_{\rm
s})$ stands for the differential subhalo mass function. This greatly
simplifies the notation as it avoids writing the dumb arguments $R_{\rm h}$ or
$M_{\rm h}$ (or 0) as well as the symbols of single and multiple derivatives
in most expressions.

Also like in Papers I and II, $M_{\rm h}$ is defined as the mass encompassed
by the virial radius $R_{\rm h}$ within which the inner mean density is equal
to the virial overdensity (Bryan & Norman, 1998; Henry, 2000) times the
current mean cosmic density. In particular, we assume the MW mass equal to
$M_{\rm h}=2.2\times 10^{12}$ M⊙. The cosmology adopted is that given by the
best WMAP7 parameters (Komatsu et al., 2011), with CDM spectrum according to
the prescription given by Bardeen et al. (1986) with the Sugiyama (1995) shape
parameter. The reader is referred to Papers I and II for the role of diffuse
dark matter (dDM) in the properties of substructure, just briefly referred to
in this Paper.

##  2 Purely Accreting haloes

During accretion haloes evolve inside-out because the later particles (and
subhaloes) fall onto them, the larger their initial turn-around radius as well
as their final apocentric radius due to the ordered virialisation process
taking place in this case (Salvador-Solé & Manrique, 2021). On the contrary,
the violent relaxation suffered by haloes in major mergers causes them to lose
the memory of their past history, so that their final properties are
indistinguishable from those of purely accreting haloes with the same mass
$M_{\rm h}$ at the same cosmic time $t_{\rm h}$. Salvador-Solé & Manrique
(2021) provide a formal proof for this important result, but the origin of it
is as follows. There is a one-to-one correspondence between haloes with mass
$M_{\rm h}$ at the cosmic time $t_{\rm h}$ and their seeds: peaks of density
contrast $\delta$ in the initial Gaussian random density field filtered with a
Gaussian window of scale $S$. That correspondence does not involve any other
halo or peak characteristic. In particular, it does not depend on how clumpy
the initial mass distribution is inside the initial patch encompassed by the
filter or, equivalently, on how lumpy the collapse of the peak is.

Therefore, all properties of haloes with $M_{\rm h}$ at $t_{\rm h}$ arising
from gravitational clustering, i.e. from their mass assembly, through
accretion and major mergers, are degenerate with respect to their formation
time. This is why to study them one has the right to assume pure accretion. In
fact, since gravitation is scale-free, such halo properties would be strictly
universal (i.e. independent of halo mass and formation time), except for the
length scale introduced by the power-spectrum of density fluctuations in the
CDM cosmology, at the base of the typical mass-concentration $M$–$c$ relation
(see Salvador-Solé et al. 2021c in progress).

However, subhalo stripping is not related to gravitational clustering and its
effects on subhaloes are not erased by violent relaxation. Consequently,
stripped subhaloes do retain the memory of the halo assembly history. This is
why substructure is expected to depend on halo formation time and, since
haloes with different masses have different typical formation times, on halo
mass as well. In this sense, the properties of substructure derived in Paper
II for purely accreting haloes might substantially differ from the properties
of ordinary haloes having suffered major mergers. Nevertheless, to understand
the latter we need first to comprehend the former from which they follow (see
Sec. 3). It is thus worthwhile reminding the results of Paper II.

As explained in Papers I and II, we distinguish between “the time of accretion
of a subhalo” onto the halo and “the time of its first crossing”. The latter
corresponds to the first time the subhalo orbits within the (virialised and
non-virialised parts of) the halo after reaching turnaround. During the first
few crossings of the system, subhalo orbits shrink due to their energy
exchange with the shells they cross, which causes the non-virialised part of
the subhalo to contract adiabatically. But, after these few initial crossings
and neglecting the effects of dynamical friction, subhalo orbits stabilise,
with the apocentre at the instantaneous virial radius of the (newly virialised
part of the) halo, which thus grows inside-out. The time at which subhalo
orbits become stable is what we adopt as the time of their accretion onto the
virialised halo.111Strictly speaking, the continuous arrival of new subhaloes
that cross the virialised halo causes it a slight abiabatic contraction. But
the characteristic time scale of this effect is very long and it can be safely
ignored. This is in fact the reason why we can see the inner halo as
virialised despite the continuous non-stabilised recent arrivals crossing it.
In what follows we concentrate in monitoring the stripping of subhaloes after
their accretion onto the host halo, i.e. once their orbits are fixed
(neglecting dynamical friction).

The stripped subhalo abundance per infinitesimal truncated mass and radius at
$M_{\rm s}^{\rm tr}$ and $r$ within a purely accreting halo with $M_{\rm h}$
at $t_{\rm h}$ is given by

$\displaystyle{\cal N}^{\rm stp}(r,M_{\rm s}^{\rm tr})={\cal N}^{\rm
tr}(r,M_{\rm s}^{\rm
tr})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle+\Bigg{\langle}\\!\\!\int_{M_{\rm s}}^{M(r)}\\!\\!{\rm
d}M\,{\cal N}^{\rm acc}(r,M)\\!\\!\int_{R^{\rm
tr}(v,r,M)}^{R(r,M)}\\!\\!\\!\\!{\rm d}r^{\prime}{\cal N}^{\rm
stp}_{\rm[M,t(r)]}(r^{\prime},M_{\rm s}^{\rm
tr})\\!\\!\Bigg{\rangle},\\!\\!\\!\\!$ (1)

where angular brackets indicate average over the tangential velocity $v$ of
subhaloes at their apocentre at $r$ where they spend most of the time. The
subindex [M,t] in the properties of subhaloes is to indicate that they refer
to host halos with $M$ at $t$. For simplicity in the notation, we have skipped
the subindex [Mh,th] for the host halo itself, but we will re-introduce it in
Sections 3 and 4 when dealing with haloes of different masses and times.

The first term on the right of equation (1), equal to

${\cal N}^{\rm tr}(r,M_{\rm s}^{\rm tr})=\mu(r,M_{\rm s}^{\rm tr})\,{\cal
N}^{\rm acc}(r,M_{\rm s}^{\rm tr}),$ (2)

gives the contribution directly arising from the stripping of accreted
subhaloes of suited mass $M_{\rm s}$, which have their apocentre at $r$ (see
Paper I) and whose abundance is

${\cal N}^{\rm acc}(r,M_{\rm s})\\!=\\!4\pi r^{2}\frac{\rho(r)}{M_{\rm
h}}\,{\cal N}^{\rm acc}(M_{\rm s}),$ (3)

where ${\cal N}^{\rm acc}(M_{\rm s})$ is their differential MF and $\rho(r)$
is the halo density profile. In equation (2), $\mu(r,M_{\rm s}^{\rm tr})$ is
the mean (averaged over $v$) truncated-to-original mass ratio of subhaloes
with $M_{\rm s}^{\rm tr}$ at $r$, calculated in Paper II by monitoring the
mass loss through repetitive stripping and shock heating of subhaloes accreted
at $t(r)$ when the host halo had radius $r$ and mass $M(r)$. And the second
term on the right of equation (1) gives the contribution arising from
subsubhaloes that were lying in accreted subhaloes with mass $M$ and
tangential velocities $v$ at $r$ and have been released into the intra-halo
medium when their hosts, with initial radius $R(r,M)$, have been truncated at
$R^{\rm tr}(v,r,M)$. See Paper II for the expression of the truncation radius
in objects with the NFW (Navarro et al., 1997) density profile.

Equation (1) can be rewritten in the simple form

${\cal N}^{\rm stp}(r,M_{\rm s}^{\rm tr})=[1+f_{\rm rel}(r,M_{\rm s}^{\rm
tr})]\,\mu(r,M_{\rm s}^{\rm tr})\,{\cal N}^{\rm acc}(r,M_{\rm s}^{\rm tr}),$
(4)

where $f_{\rm rel}(r,M_{\rm s}^{\rm tr})$ is the virtual fraction of accreted
subhaloes with mass $M_{\rm s}^{\rm tr}$ converted into stripped ones of that
mass at $r$ arising from subsubhaloes. This fraction is the solution of the
Fredholm integral equation of the second kind (with the boundary condition
$f_{\rm rel}(0,M_{\rm s}^{\rm tr})=0$) of the differential equation

$-\frac{{\rm d}}{{\rm d}r}\left[f_{\rm rel}(r,M_{\rm s}^{\rm
tr})\,\mu(r,M_{\rm s}^{\rm tr})\frac{R_{\rm
h}^{3}}{r^{3}}\right]\\!=\\!\Bigg{\langle}\frac{\partial{\cal R}}{\partial
r}\frac{{\cal N}^{\rm stp}({\cal R},M_{\rm s}^{\rm tr})}{{\cal N}^{\rm
acc}(M_{\rm s}^{\rm tr})}\Bigg{\rangle},$ (5)

with ${\cal R}(v,r)\equiv R^{\rm tr}(v,r,M_{\rm s})R_{\rm h}/R_{\rm
s}(r,M_{\rm s})$. The function $f_{\rm rel}$ is always less than a few
percent, so it can be safely neglected in front of unity, though, for the sake
of completeness, it is kept in all the expressions below. On the contrary, the
mean truncated-to-original subhalo mass ratio, $\mu(r,M_{\rm s}^{\rm tr})$,
depicted in Figure 1 plays a crucial role in the properties of substructure.
It is thus worth explaining its main features.

As the strength of stripping and shock heating depends on the concentration
$c$ of both subhaloes and the host halo (see Paper II) and $c$ depends on the
mass of the object through the well-known mass-concentration ($M$–$c$)
relation, $\mu$ depends on the masses of subhaloes and the host halo. Although
those dependencies are not gathered in Figure 1, which focuses on the $\mu$
profile for subhaloes with $10^{8}$ M⊙ in haloes with $2.2\times 10^{12}$ M⊙,
they are very important to understand the properties of substructure. Indeed,
the $\mu$ varies in a non-trivial way according to the mass of the host halo:
at $R_{\rm h}$ it is always close to unity, but its inwards decrease is less
steep in more massive haloes (see Fig. 5) because their concentration is
lower, which causes the pericentric radius reached by subhaloes with identical
$v$ at apocentre to be larger so that stripping and shock heating is less
marked (see Paper II). As a result, the mass integral of $\mu$ varies with
halo mass as $\propto M_{\rm h}^{0.08}$. Regarding the dependence of $\mu$ on
subhalo mass, it turns out that $\mu(r,M_{\rm s}^{\rm tr})$ is separable. The
reason for this is that subhaloes are truncated by tidal stripping at the
radius (dependent on their own concentration) where the inner mean density
equals that of the halo at the pericentre, which is the same for subhaloes of
all masses with identical $v$ at $r$ (see Paper II). The factor dependent on
subhalo mass is essentially proportional to $({M_{\rm s}^{\rm tr}})^{-0.03}$.
Both mass dependencies are weak, however, particularly that on subhalo mass
due to the fact that subhaloes accreted at any time $t(r)$ have similar
concentrations. In this sense the $\mu$ profile for any fixed halo mass can be
seen, in a first approximation, to depend only on $r$.

A more subtle issue is that, for the above mentioned reasons, the shape of the
$\mu$ profile will depend on the $M$–$c$ relation. Of course, the real $\mu$
profile predicted by CUSP implicitly follows from the $M$–$c$ relation that
can be derived within that framework (Salvador-Solé et al., 2021c). However,
if we want to reproduce the results of simulations, we must use the $M$–$c$
relation found in numerical studies by e.g. Gao et al. (2008) with a limited
mass resolution similar to that affecting those empirical results. The effect
of the limited mass resolution of simulations, which affects all (sub)haloes
at early times when they are little massive, is apparent in Figure 1. While
the $\mu$ profile derived from the CUSP $M$–$c$ relation with no limited mass
resolution is ever decreasing inwards, that found with the Gao et al. (2008)
$M$–$c$ relation stops decreasing at $r\sim 0.08R_{\rm h}$ and then begins to
increase again. The reason for this strange result is the following. In a
purely accreting halo evolving inside-out as considered here, its
concentration at earlier times decreases with decreasing $r$ as $r/r_{\rm s}$,
because of the fixed value of $r_{\rm s}$. Since the concentration of accreted
subhaloes also decreases with increasing $z$ (this is so in all empirical as
well as theoretical $M$–$c$ relations), stripping keeps on being effective at
small radii populated by subhaloes accreted at higher redshifts. However, due
to the limited mass resolution, the concentration of (sub)haloes in the Gao et
al. (2008) $M$–$c$ relation is bounded to a minimum value independent of mass
reached at$z\sim 3$ ($r\sim 0.3R_{\rm h}$). As a consequence, at radii smaller
than $0.3R_{\rm h}$, the concentration of the inside-out evolving host halo
continues to decrease, while that of subhaloes does not, and subhalo stripping
becomes (artificially) ineffective.

The previous discussion also illustrates that, in normal conditions, the $\mu$
profile decreases inwards despite the fact that the smaller the radius, the
lower the concentration of the host halo seen by subhaloes. Indeed, the main
cause shaping the $\mu$ profile is the time subhaloes at different radii have
been undergoing stripping rather than the different typical tidal forces they
see.

Figure 1: Mean truncated-to-original subhalo mass ratio profile predicted by
CUSP for subhaloes with $M_{\rm s}=10^{8}$ M⊙ in purely accreting MW-mass
haloes using the unbiased CUSP $M$–$c$ relation (solid red line) and the Gao
et al. (2008) empirical $M$–$c$ relation affected by the limited mass
resolution of simulations (long-dashed red line). The vertical dotted black
line marks the radius where the inwards decreasing behaviour of $\mu$ in the
latter case is inverted for the reason explained in the text.

(A colour version of this Figure is available in the online journal.)

Figure 2: Scaled number density profiles of subhaloes with different masses
(red lines) predicted for purely accreting haloes with current MW-mass using
the CUSP $M$–$c$ relation (solid red lines) and the Gao et al. (2008) $M$–$c$
relation affected by the limited mass resolution of simulations (long-dashed
red lines). For comparison we plot the fit by Han et al. (2016) to the profile
found in the halo A of the Level 1 Aquarius simulation (black dashed line)
affected by a similar resolution but having not evolved by accretion before
$z\sim 6$ corresponding to the radius marked with a vertical dotted black line
like in Fig. 1. The solid black line is the scaled halo density profile. To
better appreciate the effect of changing the subhalo masses, we plot the
predictions for subhaloes with $10^{-2}M_{\rm h}$ and $10^{-4}M_{\rm h}$ in
thick and thin red lines, respectively.

(A colour version of this Figure is available in the online journal.)

Figure 3: Cumulative MF of accreted subhaloes (upper curves) and of stripped
subhaloes (lower curves) (superindex “x” on ${\cal N}$ stands for “acc” or
“str”, respectively) predicted by CUSP for purely accreting haloes of the
different quoted masses $M_{\rm h}$. At the scale of the plot, the predictions
for the Gao et al. or CUSP $M$–$c$ relations coincide. The MFs of accreted
subhaloes for haloes with different masses overlap as their counterparts in
simulated ordinary haloes (solid black lines) and the same is true for the MFs
of stripped subhaloes provided they are multiplied by $[M_{\rm h}/(10^{12}$
M⊙$)]^{-0.08}$. These results agree with what is found in simulated ordinary
haloes (solid black lines) as derived by Han et al. 2018, here properly
normalised so as to include first-level subhaloes only (see Paper II).

(A colour version of this Figure is available in the online journal.)

Since both $\mu(r,M_{\rm s}^{\rm tr})$ and ${\cal N}^{\rm acc}(r,M_{\rm
s}^{\rm tr})$ are separable, so is also the radial abundance of stripped
subhaloes ${\cal N}^{\rm stp}(r,M_{\rm s}^{\rm tr})$ (eq. [4]). Thus, the
number density profile per infinitesimal mass of stripped subhaloes, $n^{\rm
stp}(r,M_{\rm s}^{\rm tr})\equiv{\cal N}^{\rm stp}(r,M_{\rm s}^{\rm tr})/(4\pi
r^{2})$, scaled to their total number in the halo,

$\frac{n^{\rm stp}(r,M_{\rm s}^{\rm tr})}{\bar{n}^{\rm stp}(R_{\rm h},M_{\rm
s}^{\rm tr})}\\!=\\!\frac{{\cal N}^{\rm stp}(r,M_{\rm s}^{\rm tr})/(4\pi
r^{2})}{3{\cal N}^{\rm stp}(M_{\rm s}^{\rm tr})/(4\pi R_{\rm h}^{3})},$ (6)

is independent of subhalo mass. From now on a bar on a quantity dependent on
$r$ means the radial average of that quantity inside $r$.

In Figure 2 we show the scaled number density profiles of subhaloes of two
masses ($10^{-2}M_{\rm h}$ and $10^{-4}M_{\rm h}$) found for the two above
mentioned $M$–$c$ relations. For comparison we plot the profile found by Han
et al. (2016) in the MW-mass halo A of the Level 1 Aquarius simulation
(Springel et al., 2008), also found to be roughly independent of subhalo mass
(except for the effects of dynamical friction; Han et al. 2018). This halo is
particularly well suited to the comparison with our predictions for purely
accreting haloes because it suffered the last major merger at $z\sim 6$
($r\sim 0.08R_{\rm h}$) and has been evolving by accretion (and growing
inside-out) since then. Note that, given its formation time, the concentration
$r/r_{\rm s}$ cannot be traced down to radii smaller than $r=0.08R_{\rm h}$.
This is the reason that its $\mu$ profile is not seen to increase at smaller
radii. In fact, as we will see in Section 3, it should be flat there, though
the $\mu$ profile of the simulated halo is not well determined at those radii
because it is dominated by orphan objects.222Accreted subhaloes are dubbed
“orphan” when their stripped subhaloes have masses below the mass resolution
of the simulation.

Like Figure 1, Figure 2 focuses on MW-mass haloes so it does not inform on the
dependence on halo mass of the scaled subhalo number density profile in purely
accreting haloes. However, our calculations show that the profiles for haloes
of different masses are quite similar. They do not overlap, however, because,
as mentioned, the corresponding $\mu$ profiles are not simply shifted with
respect to each other by a constant factor (see Fig. 6 dealing with ordinary
haloes).

Lastly, integrating over $r$ the subhalo abundance given in equation (4), we
obtain the subhalo differential MF, which takes the form

${\cal N}^{\rm stp}(M_{\rm s}^{\rm tr})=\overline{(1+f_{\rm rel})\mu}(R_{\rm
h})\,{\cal N}^{\rm acc}(M_{\rm s}^{\rm tr})\,.$ (7)

And integrating it over subhalo mass from $M_{\rm s}^{\rm tr}$, we arrive at
the cumulative MF, ${\cal N}^{\rm stp}(>M_{\rm s}^{\rm tr})$.

In Figure 3 we plot the cumulative MFs of accreted and stripped subhaloes
predicted for purely accreting haloes. As can be seen, the MF of accreted
subhaloes is universal, i.e. independent of halo mass, in agreement with the
results of simulations (Han et al. 2018 and references therein). This is
equivalent to say that the differential scaled subhalo abundance ${\cal
N}^{\rm acc}(M_{\rm s}/M_{\rm h})$ per infinitesimal $M_{\rm s}/M_{\rm h}$ is
also universal or that the subhalo abundance ${\cal N}^{\rm acc}(M_{\rm
s}/M_{\rm h})$ per infinitesimal $M_{\rm s}$ varies with halo mass as $M_{\rm
h}^{-1}$ (not to mix up with the differential subhalo abundance ${\cal N}^{\rm
acc}(M_{\rm s})$ per infinitesimal $M_{\rm s}$, which is roughly proportional
to $M_{\rm h}$; e.g. Fakhouri, Ma, & Boylan-Kolchin 2010). The subhalo
abundance ${\cal N}^{\rm acc}(M_{\rm s}/M_{\rm h})$ per infinitesinal $M_{\rm
s}$ which will have important consequences in Section 4. In Figure 3 we also
see that the scaled MF of stripped subhaloes is also universal provided it is
multiplied by $M_{\rm h}^{-0.08}$. Such a dependence on halo mass of the MF of
stripped subhaloes following from the above mentioned dependence of the
integral over $r$ of their $\mu$ profile as a consequence of the mass
dependence of halo concentration fully agrees with that observed in
simulations (the factor rendering the MF of simulated haloes universal is
$M_{\rm h}^{\eta}$ with $\eta=-0.1$; Han et al. 2018; Rodríguez-Puebla et al.
2016; see also Zentner et al. 2005; Giocoli et al. 2008; Gao et al. 2011).
Strictly speaking, the subhalo MF found in simulations is influenced at the
high-mass end by the effects of dynamical friction, ignored in our model. But
this effect is insignificant for subhaloes less massive than $10^{-4}M_{\rm
h}$ as represented here (see Fig. 10 in Paper II). There is one caveat,
however, in this agreement: while our predictions are for purely accreting
haloes, the results of simulations refer to ordinary haloes having suffered
major mergers. In other words, it is not clear whether the agreement is
accidental or it will persist when dealing with ordinary haloes.

The empirical result that less massive haloes are poorer (in the sense that
their subhalo MF are lower) than more massive ones is commonly interpreted as
due to the different formation times of haloes of different masses. Indeed,
the less massive a halo, the earlier it typically forms, so: 1) stripping has
more time to proceed and 2) it was more efficient because haloes are denser at
high redshifts. However, this explanation is at odds with the recent finding
by Bose et al. (2020) that the earlier MW-mass haloes form, the richer they
are. On the other hand, it is not supported either by the fact that purely
accreting haloes show the same dependence despite they all have the same
arbitrarily small formation time. Our results rather point to the fact that
such a dependence is due to the mass dependence of halo concentration, though
we must first confirm that the mass dependence of the MF in purely accreting
haloes is preserved in ordinary ones. In fact, our predictions also show that
the longer subhaloes have been stripped, the lower their final $\mu$ profile.
In other words, the different typical concentrations and typical formation
times of haloes of different masses go in the opposite direction, so we must
clarify which is the dominant effect and why in the richness of ordinary
haloes.

##  3 Ordinary Haloes of Different Masses

The density, accreted dDM mass fraction and mean abundance of accreted
subhaloes profiles, $\rho(r)$, $f_{\rm dDM}^{\rm acc}(r)$ and ${\cal N}^{\rm
acc}(r,M_{\rm s})$, respectively, are not related to tidal
stripping,333Stripping redistributes the dDM lost by subhaloes in a leading
arm and a trailing tail over their orbits. But, for the same reason that the
contribution from subhaloes to the halo density profile can be calculated
assuming they lie at their apocentre where they spend most of the time, the
contribution from the stripped dDM can also be calculated assuming it is
located at the apocentre of their orbit. Consequently, stripping does not
essentially alter the halo density profile. so they are the same in both
purely accreting and ordinary haloes having suffer major mergers (Salvador-
Solé & Manrique 2021). Thus, the only functions in the expression of the
radial distribution of stripped subhaloes (eq. [4]) that depend on stripping
are $\mu$ and $f_{\rm rel}$ defined in equations (2) and (5), respectively.
Therefore, to obtain the radial distribution of stripped subhaloes in ordinary
haloes we must first determine these two functions in such haloes.

Next we show how to obtain these properties in ordinary haloes from their
counterparts in purely accreting ones derived in Section 2, hereafter
distinguished with superindex PA. We will also use from now the scaled
arguments $x=r/R_{\rm h}$ and $m=M_{\rm s}/M_{\rm h}$ so that all functions of
those arguments should essentially coincide for purely accreting haloes of
different masses. Note that the properties of ordinary haloes of any given
mass derived next are their average over all formation times of such haloes.
Of course, those theoretical properties are hard to compare to the results of
simulations which at present provide the properties of substructure in a small
number of haloes of all masses. Nevertheless, these mean properties allow us
to elucidate the origin of their dependence on halo mass without being
disturbed by the statistical deviations of individual objects.

Be $F_{\rm[M_{\rm h},t_{\rm h}]}(x,m)$ the fraction of accreted subhaloes per
infinitesimal mass and radius that satisfy some condition in haloes of $M_{\rm
h}$ at $t_{\rm h}$ averaged over their formation time (i.e. the time of their
last major merger) and $F^{\rm PA}_{\rm[M_{\rm h},t_{\rm h}]}(x,m)$ its
counterpart in purely accreting haloes. Taking into account the inside-out
growth of haloes after their last major merger, we have the following relation
between the two quantities

$\displaystyle F_{\rm[M_{\rm h},t_{\rm h}]}(x,m)=\int_{0}^{t(x)}{\rm
d}t\,f_{\rm[M_{\rm h},t_{\rm h}]}(t)\,F^{\rm PA}_{\rm[M_{\rm h},t_{\rm
h}]}(x,m)~{}~{}~{}~{}~{}~{}~{}$ $\displaystyle+\int_{t(x)}^{t_{\rm h}}{\rm
d}t\,f_{\rm[M_{\rm h},t_{\rm
h}]}(t)\,\bar{F}_{\rm[M(t),t]}(1,m),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (8)

where $f_{\rm[M_{\rm h},t_{\rm h}]}(t)$ is the formation time probability
distribution function (PDF) of haloes with $M_{\rm h}$ at $t_{\rm h}$,
calculated within the CUSP formalism in Manrique et al. (1998) (see also Raig
et al. 2001 for a practical approximate expression in the extended Press-
Schechter (EPS) formalism; Press & Schechter 1974; Bower 1991; Bond et al.
1991; Lacey & Cole 1994). $\bar{F}_{\rm[M(t),t]}(x,m)]$ stands for the mean
fraction $F_{\rm[M(t),t]}(x,m)]$ inside $x$ of accreted subhaloes, but, given
the form of ${\cal N}^{\rm acc}(x,m)$ (eq. [3]), it coincides with the simple
radial average of $F_{\rm[M(t),t]}(x,m)$ inside $x$. This is why we denote it
with a bar. Note that the contribution on $F(x)$ from haloes formed after
$t(x)$ (the second term on the right) is averaged over their own formation
times. It thus takes into account the different weight of haloes formed at
$t(x)$ with previous different formations times (and so on so forth). On the
contrary, the contribution on $F(x)$ from haloes formed before $t(x)$ (the
first term) does not depend on their individual formation times ($F^{\rm PA}$
does not depend on any formation time) because, at the radius $x$, such haloes
are accreting material ex-novo.

To write equation (8) we have taken into account that, when a halo suffers a
major merger, its content is scrambled,444The scrambling must be complete,
otherwise major mergers would not cause haloes to fully lose the memory of
their past history as they do (Salvador-Solé & Manrique, 2021). so the
fraction $F_{\rm[M(t),t]}(x,m)$ at any radius $x$ equals its mean value within
the total radius at that moment, $\bar{F}_{\rm[M(t),t]}(1,m)$. In addition, we
have taken into account that, after the last major merger, haloes evolve by
pure accretion, so $M(t)/M_{\rm h}$ is the mass track of purely accreting
haloes with boundary condition $M(t_{\rm h})/M_{\rm h}=1$. Note that, even if
$F^{\rm PA}_{\rm[M_{\rm h},t_{\rm h}]}(x,m)$ (in scaled arguments) is
essentially universal, $F_{\rm[M_{\rm h},t_{\rm h}]}(x,m)$ will depend on
$M_{\rm h}$ and $t_{\rm h}$ through the explicit dependence on these
quantities of the halo formation time PDF (Fig. 4).

Figure 4: Formation time PDFs of haloes with the quoted masses (coloured
lines) at the present time $t_{\rm h}=t_{0}$. The coloured dashed vertical
lines mark the median formation time $t_{\rm med}$ of haloes of each mass.

(A colour version of this Figure is available in the online journal.)

Multiplying equation (8) by ${\cal N}^{\rm acc}(x,m)$ and integrating over $x$
out to 1, we are led, by partial integration and taking into account equation
(3), to

$\displaystyle\bar{F}_{\rm[M(t_{\rm h}),t_{\rm h}]}(1,m)=\overline{F^{\rm
PA}_{\rm[M_{\rm h},t_{\rm h}]}\,f_{\rm[M_{\rm h},t_{\rm h}]}^{\rm
c}}(1,m)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!+\\!\\!\int_{0}^{t_{\rm h}}\\!\\!{\rm d}t\,f_{\rm[M_{\rm
h},t_{\rm h}]}(t)\bar{F}_{\rm[M(t),t]}(1,m)\frac{M(t)}{M_{\rm
h}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (9)

where $f_{\rm[M_{\rm h},t_{\rm h}]}^{\rm c}(x)$ stands for the cumulative
formation time PDF of haloes with $M_{\rm h}$ at $t_{\rm h}$ up to $t(x)$.
Equation (9) is a Volterra equation of second kind for
$\bar{F}_{\rm[M(t),t]}(1,m)$ as a function of $t$. Note that, according to
equation (9), the mass average of $F$ out to 1, $\bar{F}_{\rm[M(t),t]}(1,m)$,
is different from $\bar{F}^{\rm PA}{\rm[M(t),t]}(1,m)$. The reason for this
will be seen below. Bringing the solution of this Volterra equation in the
integral on the right of equation (8), we arrive at the desired function
$F_{\rm[M(t_{\rm h}),t_{\rm h}]}(x,m)$ for any value of $x$ (and $m$). Note
also that, as $n_{\rm[M_{\rm h},t_{\rm h}]}^{\rm c}(x)$ does not depend on
$m$, equations (9) and (8) imply that $F_{\rm[M(t_{\rm h}),t_{\rm h}]}(x,m)$
would be independent of $m$ provided its counterpart $F^{\rm PA}$ were. (This
is approximately the case for the $\mu$ and $f_{\rm rel}$ functions; see
next.)

Figure 5: Mean truncated-to-original subhalo mass ratio profiles predicted
using the CUSP $M$–$c$ relation for ordinary haloes of several masses $M_{\rm
h}$ at the current time $t_{h}$ averaged over their formation times (solid
coloured lines). The results for different subhalo masses $M_{\rm s}^{\rm tr}$
(from $10^{6}$ M⊙ to $10^{9}$ M⊙) overlap when they are multiplied by $[M_{\rm
s}^{\rm tr}/10^{8}$ M⊙$)]^{-0.03}$. For comparison we plot the prediction for
purely accreting MW-mass haloes (solid black line). To illustrate the effects
of averaging over halo formation times we also plot in long-dashed and short-
dashed lines the profiles for two individual haloes of $10^{12}$ M⊙ formed at
a high- and low-redshift, respectively.

(A colour version of this Figure is available in the online journal.)

This procedure can be applied to the virtual fractions of accreted subhaloes
with $M_{\rm s}^{\rm tr}$ at $r$ converted into stripped subhaloes by direct
stripping and through the release of subsubhaloes we arrive at the $\mu(x,m)$
and $f_{\rm rel}(x,m)$ profiles in ordinary haloes of any mass averaged over
their formation times. In Figure 5 we depict the $\mu(x,m)$ profile of
ordinary haloes averaged over their formation times that is predicted for the
CUSP $M$–$c$ relation. To better realise the effect of that average, we also
plot the schematic (approximate) $\mu$ profiles of two individual haloes of
$10^{12}$ M⊙, one formed at $t_{\rm f}\sim 2$ Gyr (corresponding to
$r=0.2R_{h}$ in the final halo grown inside-out since that moment) and the
other one formed at $t_{\rm f}\sim 10$ Gyr (corresponding to $r=0.8R_{\rm
h}$). For the reason mentioned above when explaining the meaning of
$\bar{F}_{\rm[M(t),t]}(1,m)$ in the second integral on the right of equation
(8), the scrambling of the system at the last major merger yields a flat $\mu$
profile equal to its radial average within the formation radius. This is not
only the case just after the merger, but also long time after. Indeed, after
the major merger all subhaloes inside the scrambled region orbit and suffer
stripping during the same time interval $t_{\rm h}-t_{\rm f}$ regardless of
their past history, so $\mu$ in that region deepens, but it remain essentially
flat until the next merger or the final time. Strictly speaking, subhaloes at
different radii within that region suffer different tidal forces over their
orbits (the smaller the radius, the lower the concentration of the halo they
see). But, as mentioned in connection with Figure 1, this has a much less
marked effect on $\mu$ than the different time subhaloes suffer stripping. In
fact, if the ending $\mu$ profile within that region did show a substantial
dependence on radius, the accurate $\mu$ profile of haloes averaged over their
formation times commented below would show it, which is not the case (see Fig.
5). Likewise, the flat $\mu$ profile of those individual halos inside the
formation radius have been taken equal to the radial average inside that
radius of the $\mu^{\rm PA}$ at the ending time, while it is actually somewhat
lower due to the stripping suffered by subhaloes during the time elapsed since
the scrambling. Besides those simplifications, the important point to retain
from those schematic examples is that, since after the merger haloes evolve
inside-out by accretion, their (essentially flat) $\mu$ profile in the inner
region jumps at its edge to the $\mu^{\rm PA}$ profile of purely accreting
haloes with the same mass $M_{\rm h}$ at $t_{\rm h}$.

That behaviour of the $\mu$ profile in individual haloes of a given mass
translates into their formation-time average. This is the reason why the
formation-time-averaged profile is also flat at small radii and begins to
increase at some radius, dependent on the typical formation time of haloes of
that mass, towards the profile of purely accreting haloes. The difference is
that, while the profiles of individual haloes reach the purely accreting
solution right at the formation radius, the formation-time-averaged profile
only reaches it at $R_{\rm h}$ so that it stays systematically below that
solution at large radii. The reason for this difference is clear. Since the
$\mu$ profiles are outwards increasing, their mass average inside any
formation radius is always lower than the original value at that radius, so
the formation-time-averaged profile is also lower than the purely accreting
solution. And, as large radii contribute the most to the mass average of $\mu$
inside $R_{\rm h}$, that mass average in ordinary haloes averaged over their
formation times will always be somewhat smaller than that of purely accreting
haloes. In any event, the difference should be similar for haloes of any mass,
which implies that the mass average of $\mu$ in ordinary haloes will
essentially coincide with that in purely accreting ones arising, as mentioned
in Section 2, from the mass dependence of halo concentration. All these
conclusions referring to the $\mu$ profiles of ordinary haloes will translate
into their subhalo radial distributions and MFs, which will explain their
behaviour.

Figure 6: Scaled number density profiles of subhaloes of different masses
$M_{\rm s}$, which overlap with no added factor, in ordinary haloes of several
masses $M_{\rm h}$ (coloured lines), compared to the profile obtained in
purely accreting haloes of $10^{12}$ M⊙ (thick black line) and the scaled mass
density profile of such haloes (thin black line).

(A colour version of this Figure is available in the online journal.)

Once the functions $\mu$ and $f_{\rm rel}$ have been determined, we can
proceed to derive the subhalo abundance per infinitesimal mass and radius
around $r$ and $M_{\rm s}^{\rm tr}$, ${\cal N}^{\rm stp}(r,M_{\rm s}^{\rm
tr})$, in ordinary haloes by application of equation (4). Specifically, since
the abundance of accreted subhaloes does not depend on the halo formation time
and $f_{\rm rel}$ is negligible, the subhalo radial abundance in ordinary
haloes of a given mass averaged over their formation times is simply equal to
the formation-time-averaged $\mu$ profile times the abundance of accreted
subhaloes.

The first consequence of the form of the radial abundance of subhaloes in
ordinary haloes refers to the scaled subhalo number density profile (eq. [6]).
As can be seen in Figure 6, the profiles of subhaloes of different masses
overlap because their $\mu$ profiles differ, as mentioned, by the same
constant factor as in purely accreting haloes, which cancels with their
scaling. Also like in purely accreting haloes, the profiles for haloes of
different masses do not overlap because their $\mu$ profiles do not differ by
just a constant factor. But there is one interesting difference in comparison
with the case of purely accreting haloes: the scaled number density profiles
are substantially steeper now. In fact, at small enough radii they are
parallel to the halo density profile. (To avoid crowding, in Figure 6 we only
plot the mass density profile and the scaled subhalo number density profile
for haloes of $10^{12}$ M⊙, so that this comparison is only possible for such
haloes.) The reason for such a behaviour of the scaled number density profile
of ordinary haloes is that in the scrambled regions subhaloes of all masses
have been mixed up with the dDM, so their number density profiles become
proportional to the mass density profile of the halo there.

Figure 7: Cumulative subhalo MFs predicted for ordinary haloes of several
masses multiplied by $[M_{\rm h}/(10^{12}$ M⊙$)]^{-0.07}$. The universal MF of
purely accreting haloes (thick black line) is just slightly higher.

(A colour version of this Figure is available in the online journal.)

By integrating over $r$ the previous radial abundance of subhaloes, we are led
to the differential subhalo MF in ordinary haloes of different masses (eq.
[7]) and by integration over $M_{\rm s}^{\rm tr}$ the corresponding cumulative
MFs shown in Figure 7. For the reason mentioned when describing the mass
integral of $\mu$, the cumulative MF of ordinary haloes averaged over their
formation times shows the same dependence on halo mass as purely accreting
haloes. (Strictly speaking, they overlap now when they are multiplied by
$[M_{\rm h}/(10^{12}$ M⊙$)]^{-0.07}$. The slight difference in the power index
with respect to that of purely accreting haloes ($-0.08$) arises from the
distinct formation time PDFs of haloes with different masses.) In addition,
the MF of ordinary haloes averaged over their formation times is somewhat
smaller than that of purely accreting haloes, as expected from the discussion
on the values of $\mu$ near $R_{\rm h}$. On the contrary, the MF of individual
ordinary haloes fully overlap with the MF of purely accreting haloes of the
same mass, which explains the result in Figure 3. Only the MF of extreme late-
forming haloes is substantially lower because of the rapid decrease of their
$\mu$ profile when going away from $R_{\rm h}$. We will comeback to this
behaviour of the MF of the latest-forming haloes in next Section.

##  4 Ordinary Haloes of a Fixed mass and Different Formation Times

The fraction of accreted subhaloes satisfying any desired property in haloes
with a fixed mass $M_{\rm h}$ at $t_{\rm h}$ formed in an interval $\Delta$ of
time around any desired value $t_{\rm f}$ can be derived using equation (8)
with the formation time PDF restricted within that interval, that is

$f^{t_{\rm f},\Delta}_{\rm[M_{\rm h},t_{\rm h}]}(t)=A\,\Pi(t-t_{\rm
f},\Delta)\,f_{\rm[M_{\rm h},t_{\rm h}]}(t),$ (10)

where $\Pi(t-t_{\rm f},\Delta)$ is the top-hat function of width $\Delta$
around $t_{\rm f}$ and $A$ is the normalization factor. This way we are led to

$\displaystyle F_{\rm[M_{\rm h},t_{\rm h}]}^{t_{\rm
f},\Delta}(x,m)=\int_{0}^{t(x)}{\rm d}t\,f^{t_{\rm f},\Delta}_{\rm[M_{\rm
h},t_{\rm h}]}(t)\,F^{\rm PA}_{\rm[M_{\rm h},t_{\rm
h}]}(x,m)~{}~{}~{}~{}~{}~{}~{}$ $\displaystyle+\int_{t(x)}^{t_{\rm h}}{\rm
d}t\,f^{t_{\rm f},\Delta}_{\rm[M_{\rm h},t_{\rm
h}]}(t)\bar{F}_{\rm[M(t),t]}[1,m],~{}~{}~{}~{}~{}~{}~{}$ (11)

Thus, equation (11) with the function $\bar{F}_{\rm[M(t),t]}[1,m]$ solution of
the integral equation (9) leads to the desired fraction $F^{t_{\rm
f}\pm\Delta/2}_{\rm[M_{\rm h},t_{\rm h}]}(x,m)$.

This procedure can be applied to infer the $\mu$ and $f_{\rm rel}$ profiles in
ordinary haloes with formation times averaged inside any desired interval.
Following Bose et al. (2020), who studied MW-mass haloes assuming the WMAP7
cosmology like we do, we consider two extreme intervals: one with an upper
bound at $t/t_{\rm h}=0.28$ ($z=1.74$), corresponding to $r/R_{\rm h}=0.36$,
that embraces the 20% earliest-forming objects and another interval with a
lower bound at $t/t_{\rm h}=0.76$ ($z=0.29$), corresponding to $r/R_{\rm
h}=0.82$, that embraces the 20% latest-forming objects. (The previous figures
correspond to haloes with $M_{\rm h}=10^{12}$ M⊙; for haloes of different
masses we have chosen suited values to delimit the same kind of early- and
late-forming objects.) We remark that these two intervals coincide with those
used by Bose et al. (2020) despite the different halo formation time
definition they adopt: the time the halo reaches 50% of its final mass.
Indeed, haloes with $M_{\rm h}$ at $t_{\rm h}$ that follow the pure accretion
track reaches half the final mass at $z=3.5$. Thus, all haloes undergoing the
last major merger after that redshift automatically reach that pure accretion
track at the same moment and, hence, they are also seen to form there
according to the alternate formation time definition. While haloes undergoing
the last major merger before $z=3.5$, in the alternate formation time
definition will be seen to form some time after when they will reach $z=3.5$.
But, since this redshift is higher than the upper redshift of our interval of
the earliest-forming haloes, all these haloes will lie in that interval
according to both formation time definitions, even though their individual
formation times will differ in both cases.

The $\mu$ profiles for ordinary haloes formed in those two extreme intervals
are shown in Figure 8. Their shape in each case is readily understood from our
previous explanations of the behaviour of that profile in Section 3. It also
explains the behaviour of the corresponding scaled subhalo number density
profiles, $n(r,M_{\rm s})/\bar{n}(M_{\rm s})$ (eq. [6]), shown in Figure 9,
that follow from their radial abundances given by equation (4). As can be
seen, the mean scaled subhalo number density profiles of the earliest-forming
haloes of any mass are quite similar to each other: they are parallel to the
scaled mass density profile of their respective haloes until $r/R_{\rm h}\sim
0.3$, where such haloes typically formed, and then rapidly increase reaching
the $\mu$ profile of purely accreting haloes at a finite radius substantially
smaller than $R_{\rm h}$. On the contrary, the scaled profiles of the latest-
forming haloes keep their initial trend parallel to the density profile of the
respective haloes until a much larger radius, where they suddenly recover to
reach the $\mu$ profile of purely accreting haloes at $R_{\rm h}$. Thus, the
two scaled subhalo number density profiles are quite distinct, meaning that
this property would be a good tool for probing the halo formation time.
However, since these profiles are scaled to the total number of subhaloes of
each mass, they do not inform on the subhalo richness of haloes.

An alternative estimate of the radial distribution of subhaloes that is non-
scaled and easier to determine in simulations as well as observations because
using cumulative quantities is that put forward by Bose et al. (2020), namely
the profile of the total number of subhaloes with masses below some given
value that lies inside each radius $r$. In Figures 10 and 11 we plot two
versions of this integrated density profile: one dealing with “plain”
subhaloes as used so far, with masses in the range $10^{-6}M_{\rm h}<M_{\rm
s}^{\rm tr}<10^{-3}M_{\rm h}$ (the lower limit is supposed to account for the
typical mass resolution of simulations of MW-mass haloes), and another one
dealing with “luminous” subhaloes, i.e. subhaloes harbouring faint and ultra-
faint satellites (with stellar mass less than $3.2\times 10^{6}$ M⊙), ready to
be compared with the results of the simulations of MW-mass haloes carried by
Bose et al. (2020) (see their Fig. 5).

Figure 8: Mean truncated-to-original subhalo mass ratio profiles predicted
using the CUSP $M$–$c$ relation for the 20% earliest-forming haloes (long-
dashed lines) and the 20% latest-forming haloes (short-dashed lines) of
several masses. For comparison we plot the same profile for a purely accreting
halo of $10^{12}$ M⊙ (thick black line).

(A colour version of this Figure is available in the online journal.)

Figure 9: Same as Figure 8 (same lines) but for the scaled subhalo number
density profiles. Also plotted is the formation time-independent mass density
profile of haloes with $10^{12}$ M⊙ (thin black line) scaled so as to
facilitate the comparison.

(A colour version of this Figure is available in the online journal.)

As can be seen in Figure 10, the integrated radial number density profile of
plain subaloes in haloes of different masses essentially overlap, in both the
20% earliest- and 20% latest-forming haloes, when they are multiplied by
$[M_{\rm h}/(10^{12}$ M⊙)]-0.08; there is just a small dispersion at
intermediate radii due to the fact that the formation time intervals cannot be
scaled to $M_{\rm h}$. The fact that this mass dependence essentially
coincides with that of the subhalo MFs555It is again $0.08$ rather than $0.07$
due to the narrower formation time interval used here compared to that used in
the subhalo MF of ordinary haloes. is not unsurprising: both properties arise
from the integration over $r$ of the original radial abundance of stripped
subhaloes (eq. [4]). Regarding the radial distributions in haloes of any fixed
mass, we see that in the earliest-forming haloes it is slightly higher at
large radii ($r/R_{\rm h}\ga 0.6$) than in the latest-forming ones. Since
subhaloes at those radii contribute the most to the total number of objects in
the haloes, this means that the former haloes are richer than the latter ones
(see the MFs of the two kinds of haloes below). But this is not the only
difference: the profile of the earliest-forming haloes decreases towards the
halo centre much more steeply than that of the latest-forming haloes, which
essentially stay parallel to the integrated mass density profile of the halo.
As a consequence, both profiles cross at $r\sim 0.06R_{\rm h}$. Thus, the
integrated radial distribution of plain subhaloes would be a very good probe
of halo assembly history provided subhalo masses could be determined.

In Figure 11 we compare the theoretical radial distributions of luminous
subhaloes to those of satellites found by Bose et al. (2020) in the 20%
earliest- and latest-forming MW-mass haloes ($M_{200}\sim 1-1.3\times 10^{12}$
M⊙). In principle, to calculate such theoretical profiles we should find the
upper and lower mass limits of stripped subhaloes that correspond to the upper
and lower stellar mass limits of $3.2\times 10^{6}$ M⊙ and zero M⊙,
respectively, used by Bose et al. But this correspondence depends on the
radius in an unknown way, so that calculation is hard to achieve. Fortunately,
there is an alternative way to calculate the desired radial distributions. The
vast majority of satellites with stellar masses below $\sim 3.2\times 10^{6}$
M⊙ necessarily formed in haloes with masses between the minimum mass for star
formation, $M_{\rm s}\sim 1.4\times 10^{8}$ M⊙ (Bose et al., 2020), and the
maximum mass of haloes having formed $3.2\times 10^{6}$ M⊙ of stars, $\sim
5\times 10^{8}$ M⊙ . Certainly, a small fraction of satellites with a final
stellar mass slightly below $3.2\times 10^{6}$ M⊙ may have suffered very
strong tidal stripping affecting not only dark matter but also stars so that
the initial mass of their subhaloes may be slightly larger than $5\times
10^{8}$ M⊙. Consequently, the upper mass limit is actually somewhat fuzzy.
Nonetheless, since the radius of the stellar component is typically 1% of the
radius of their host haloes (Kravtsov, 2013), the fraction of satellites with
initial stellar mass larger than $\mu M_{\rm s}$ able to lose stars is very
small. In addition, any small variation in the upper mass limit has an
insignificant effect on the total number of haloes within the bracketed mass
range. (Only do variations in the lower mass limit substantially affect that
number due to the much higher abundance of objects at the low mass end.)
Therefore, the number of accreted subhaloes within that mass range is a very
robust estimate for the number of satellites with final masses below
$3.2\times 10^{6}$ M⊙. And, since low mass subhaloes, unaffected by dynamical
friction, stay at the same apocentric radius as when they were accreted, the
integrated radial distribution of satellites with masses in the included
sample should essentially coincide with that of accreted subhaloes in the
corresponding mass range independent of radius. We have confirmed the
robustness of this counting against small changes in the upper subhalo mass
limit.

Certainly, the previous reasoning seems to contradict the results of
simulations. As discussed in Paper I (see also Salvador-Solé et al. 2021c),
the number density profile of accreted subhaloes of any mass is independent of
the formation time of the host halo. Yet, Bose et al. (2020) found a
significant difference between the radial distributions of satellites in the
20% earliest- and 20% latest-forming haloes. The clue for that apparent
contradiction is the different formation time PDF of haloes of different
masses (Fig. 4). According to it, the 20% earliest- and 20% latest-forming
haloes should essentially coincide with the less and most massive objects,
respectively, in the sample of haloes with masses $M_{200}\sim 1-1.3\ \times
10^{12}$ M⊙ used by Bose et al. (2020). And, since the abundance of accreted
subhaloes within any given mass range is proportional to $M_{\rm h}^{-1}$ (see
Sec. 2), the 20% earliest-forming haloes should be about $1.3$ times richer
than the 20% latest-forming ones, where $1.3$ is the ratio of extreme halo
masses in the sample.

Figure 10: Integrated number density profiles of plain subhaloes with masses
$10^{-6}M_{\rm h}<M_{\rm s}^{\rm tr}<10^{-3}M_{\rm h}$ predicted by CUSP for
the 20% earliest-forming (long-dashed lines) and 20% latest-forming (short-
dashed lines) haloes of different masses $M_{\rm h}$ multiplied by $[M_{\rm
h}/(10^{12}$ M⊙)]-0.08. For comparison we plot the same profile for purely
accreting haloes of the MW-mass (thick black line) and the integrated mass
density profile of the (purely accreting or ordinary) halo (thin black line)
with suited zero-point so as to essentially overlap with the profiles of the
latest-forming haloes.

(A colour version of this Figure is available in the online journal.)

Figure 11: Integrated number density profiles of accreted subhaloes with
masses $1.4\times 10^{8}$ M⊙ $<M_{\rm s}<5\times 10^{8}$ M⊙ corresponding to
satellites with stellar masses $M_{\star}<10^{6}$ M⊙ in haloes with
$M_{200}=1\times 10^{12}$ M⊙ (long-dashed red line) and $M_{200}=1.3\times
10^{12}$ M⊙ (short-dashed red line) predicted by CUSP compared to the average
profiles of the 20% earliest-forming (long-dashed black line) and 20% latest-
forming (short-dashed black line) haloes found by Bose et al. (2020) for MW-
mass haloes ($M_{200}\sim 1-1.3\ \times 10^{12}$ M⊙).

(A colour version of this Figure is available in the online journal.)

In Figure 11 the integrated radial distributions of luminous subhaloes or,
more exactly, the corresponding accreted subhaloes predicted by CUSP for
haloes with the two extreme masses are compared to the integrated radial
distributions of satellites (including orphan objects) found by Bose et al.
(2020) for the 20 % earliest- and 20 % latest-forming haloes with $M_{200}\sim
1-1.3\times 10^{12}$ M⊙. These authors provide the results obtained from the
Copernicus Complexio Low Resolution (COLOR) simulation and the Copernicus
Complexio (COCO) suite of simulations, with substantially higher resolution
but much poorer statistics (COLOR has $\sim 20$ times more haloes than COCO).
In any case, the mass resolution in both simulations is sufficient to detect
star formation in all haloes with masses larger than $1.4\times 10^{8}$ M⊙, so
including orphan satellites in stripped subhaloes with ending masses below the
resolution mass should be sufficient to count all satellites. Yet, the higher
resolution of the COCO simulations yields twice more satellites of the
relevant masses than in the COLOR simulations (see Fig. 1 of Bose et al.
2020). Since all-level subhaloes are also twice more abundant than first-level
ones (Paper I), that result suggests that many satellites in the COCO
simulations lie in high-level subhaloes. This would explain why including one
orphan satellite per disrupted (first-level) subhalo in the COLOR simulation
is not enough to recover the satellite abundance found in the COCO
simulations. This conclusion is also supported by the fact that the satellite
abundance in the COCO simulations agrees with the (all-level) subhalo
abundance found in merger tree Monte Carlo simulations according to the EPS
formalism (Bose et al., 2020). But the idea that 50% of all satellites lies in
high-level subhaloes is little realistic if we think about MW satellites. We
have thus chosen to compare the radial distribution of (first-level) accreted
subhaloes predicted by CUSP to the radial distribution of satellites found in
the COLOR simulation.

As can be seen in Figure 11 the predicted profiles almost fully agree, indeed,
with those found by Bose et al. (2020). The only slight difference in the
profiles of the earliest-forming haloes is likely due to the effects of
dynamical friction, excluded from our treatment. The upper subhalo mass limit
of $\sim 5\times 10^{8}$ M⊙ is a factor $\sim 5$ higher than the minimum mass
of subhaloes suffering significant dynamical friction, so in early-forming
haloes subhaloes could indeed be slightly more concentrated towards the
centre. In any event, this agreement gives strong support to the conclusion
that the radial distributions of satellites in haloes with masses in a very
narrow range around any fixed value do not depend on their formation times. In
addition, it reinforces the idea that satellites in the COLOR simulation lie
in first-level subhaloes only, while those in the COCO simulations likely also
occupy higher-level subhaloes.

Figure 12: Same as Figure 10 but for the cumulative subhalo MFs (also
multiplied by $[M_{\rm h}/(10^{12}$ M⊙)]-0.08

. (A colour version of this Figure is available in the online journal.)

The inclusion of orphan satellites in simulations is crucial for the previous
comparison to be meaningful. Otherwise stripping would affect the radial
distribution of satellites through the varying number of disrupted subhaloes
produced (e.g. Graus et al. 2019) and we could not use accreted subhaloes to
count satellites. This does not mean, of course, that we can use the radial
profiles of non-orphan satellites to probe halo formation times. Even though
this is in principle possible in simulations, it is not in observations:
reaching higher magnitudes allows a better determination of the radial
distributions (and MFs) of satellites (e.g. Grand et al. 2021), but it does
not alter the fundamental result that the abundance of satellites coincides
with that of accreted subhaloes, independent of halo formation time. A better
determination of the radial distribution of ultra-faint satellites should be
useful, however, in connection with the “missing satellite problem” itself. By
comparing the radial distribution (or the MF) of ultra-faint satellites to the
theoretical one(s) of the corresponding low-mass accreted subhaloes (eq. 3,
with ${\cal N}^{\rm acc}(M_{\rm s})$ given by the well-normalised MF of
accreted subhaloes; see Paper I), it should be possible to confirm whether or
not there is a problem in the satellite abundance. Note that such a test,
unaffected by the complications arising from stripping and dynamical friction,
has the added advantage that it is independent of the assembly history of
haloes. To do that it would be advisable, however, to exclude the very central
region of haloes so as to avoid the effect of satellite destruction by central
galaxies (see Bose et al. 2020 and references therein).

Finally, integrating the radial abundances, ${\cal N}^{t_{\rm
f}\pm\Delta}_{\rm[M_{\rm h},t_{\rm h}]}(r,M_{\rm s})$ out to $R_{\rm h}$, we
are led to the differential MF of plain or luminous subhaloes in haloes of
different formation times in early-forming and late-forming haloes of any
fixed mass. In Figure 12 we depict the cumulative MFs of plain subhaloes for
haloes of different masses. As usual, all MFs depend on halo mass as $M_{\rm
h}^{0.08}$. Regarding their similarity with the respective MF of purely
accreting haloes, we see that, as expected, the MFs in the earliest-forming
haloes of any mass overlap with it, while those of the latest-forming haloes
are substantially lower. The difference is consistent with that found in the
radial distribution of subhaloes at large radii, while the shape of all these
MFs is always the same, essentially a power-law of index close to $-1$. In
this sense, subhalo MF is a poorer probe of halo formation time than the
radial distribution of subhaloes.

##  5 SUMMARY AND DISCUSSION

The present Paper is the last one of a series of three devoted to provide a
comprehensive insight on halo substructure. Making use of the CUSP formalism
we derived, in Paper I, the properties of accreted subhaloes from the
statistics of their seeds (nested peaks) in the primordial Gaussian random
density field. In Paper II, we developed a detailed stripping model and used
the properties of accreted subhaloes to derive those of stripped ones taking
into account the repetitive tidal truncation and shock heating they suffer as
they orbit within their hosts. As argued in Paper II, the properties of
stripped subhaloes depend on the particular assembly history of their host
halo. Nevertheless, their derivation was carried in the simplest, little
realistic case of purely accreting haloes. Here we have extended those results
to ordinary haloes, i.e. having undergone major mergers, with the aim to: 1)
see how the properties of substructure determined in Paper II are modified by
the action of major mergers, 2) determine their dependence on substructure on
halo mass and formation time (i.e. the time of the last major merger) and 3)
find the capabilities of substructure as a probe of halo assembly history.

In a formal plane, it has been shown that all the properties of substructure
in haloes of any mass and formation time are encoded in the mean truncated-to-
original subhalo mass ratio profile, $\mu(r,M_{\rm s}^{\rm tr})$, weakly
dependent, actually, on $M_{\rm s}^{\rm tr}$. This profile is flat inside the
radius reached by the object at its last major merger (inside which the system
was completely scrambled) and rapidly approaches, at larger radii, the $\mu$
profile of the purely accreting halo of the same mass, setting the
correspondence between $r$ and $t$ which can be used to determine the time of
the merger. This profile is at the base of all the remaining more practical
results.

Regarding our first goal, we have demonstrated that the agreement found in
Paper II between our predictions for purely accreting haloes and the results
of simulations dealing with ordinary ones was not casual. Even though major
mergers affect the properties of substructure in ordinary haloes with respect
to purely accreting ones, their subhalo MF remains essentially unaltered.
Major mergers only leave a substantial imprint in the radial distribution of
subhaloes.

These results have important repercussions on our second goal. On the one
hand, the dependence on halo mass of the subhalo MF is the same in ordinary
haloes as in purely accreting ones. Since in the latter the formation time
plays no role, that dependence can only arise from the mass dependence on mass
of halo concentration, as found in Paper II. This dependence on halo mass is
already seen, of course, in the (non-scaled) radial distribution of subhaloes
whose integral over the radius leads to the MF. But the radial distribution of
subhaloes in haloes of any given mass has the added interest that it also
harbours information on the formation time of the host halo. Indeed, the
scaled subhalo number density profiles in ordinary haloes, independent of
subhalo mass, is steeper than in purely accreting haloes, being proportional
to the density profile of the halo inside the scrambling radius reached at its
formation. This scaled version of the radial distribution of subhaloes does
not inform on the subhalo richness, but another non-scaled version of it, the
integrated subhalo number density profile, does. We have shown that this
latter profile is higher at large radii contributing the most to the total
subhalo abundance in early-forming haloes than in late-forming ones. Moreover,
the profiles in the two kinds of haloes are also very different as they cross
each other at some intermediate radii. Of course, the different richness of
haloes of a given mass is also reflected in their subhalo MF, which in
extremely late-forming haloes is found to be substantially lower than in all
the rest.

All these results lead to the following conclusion regarding the third goal:
the radial distribution of subhaloes is very useful to probe the halo
formation time, much better than the subhalo MF. Unfortunately, this
conclusion only holds for plain subhaloes, not for satellites. Indeed, the
properties of satellites do not depend on how their host halo stripped
subhaloes. They only depend on the masses of accreted subhaloes where they
formed, whose properties are independent of the halo formation time. Thus, the
properties of ultra-faint satellites do not inform on the formation time of
the host halo. The good news is that, by simply comparing their radial
distribution to that of very low mass subhaloes they should unambiguously tell
whether or not there really is a “missing satellite problem”, without
depending on the particular formation history of the host halo.

In its current form, our analytic treatment does not account for dynamical
friction. Of course, this is not a drawback when dealing with low enough mass
subhaloes (faint enough satellites). But, if we want to deal with more massive
subhaloes (more luminous satellites), we should account for that process. That
possibility would be very welcome because, even though the radial distribution
of satellites does not depend on subhalo stripping, it is sensitive to
dynamical friction. Thus, the properties of satellites could still inform on
the halo assembly history through the effect of dynamical friction. Work in
this line is currently in progress (Salvador-Solé et al., 2021d).

ACKNOWLEDGEMENTS

One of us, I.B., has benefited of a MEXT scholarship by the Japanese MECSST.
This work was funded by grants CEX2019-000918-M (Unidad de Excelencia ‘María
de Maeztu’) and PID2019-109361GB-100 (together with FEDER funds) by
MCIN/AEI/10.13039/501100011033 and by the grant 2017SGR643 funded by the
Catalan DEC.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Angulo et al. (2009) Angulo R. E., Lacey C. G., Baugh C. M., Frenk C. S., 2009, MNRAS, 399, 983
  * Applebaum et al. (2021) Applebaum E., Brooks A. M., Christensen C. R., Munshi F., Quinn T. R., Shen S., Tremmel M., 2021, ApJ, 906, 96
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15
  * Bennet et al. (2019) Bennet P., Sand D. J., Crnojević D., Spekkens K., Karunakaran A., Zaritsky D., Mutlu-Pakdil B., 2019, ApJ, 885, 153
  * Bennet et al. (2020) Bennet P., Sand D. J., Crnojević D., Spekkens K., Karunakaran A., Zaritsky D., Mutlu-Pakdil B., 2020, ApJL, 893, L9
  * Benson et al. (2013) Benson A. J., Farahi A., Cole S., et al., 2013, MNRAS, 428, 1774
  * Bose et al. (2016) Bose S., Hellwing W. A., Frenk C. S., Jenkins A., Lovell M. R., Helly J. C., Li B., 2016, MNRAS, 455, 318
  * Bose et al. (2020) Bose S., Deason A. J., Belokurov V., Frenk C. S., 2020, MNRAS, 495, 743
  * Bond et al. (1991) Bond J.R., Cole S., Efstathiou G., & Kaiser N., 1991, ApJ, 379, 440
  * Bower (1991) Bower R. G., 1991, MNRAS, 248, 332
  * Boylan-Kolchin et al. (2010) Boylan-Kolchin M., Springer V., White S. D. M., Jenkins A., 2010, MNRAS, 406, 896
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bullock & Boylan-Kolchin (2017) Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343
  * Carlsten et al. (2020) Carlsten S. G., Greene J. E., Peter A. H. G., Greco J. P., Beaton R. L., 2020, ApJ, 902, 124
  * Carlsten et al. (2021) Carlsten S. G., Greene J. E., Peter A. H. G., Beaton R. L., Greco J. P., 2021, ApJ, 908, 109
  * Cautun et al. (2014) Cautun M., Hellwing W. A., van de Weygaert R., Frenk C. S., Jones B. J. T., Sawala W., 2014, MNRAS, 445, 2049
  * Crnojević et al. (2019) Crnojević D., Sand D. J., Bennet P., Pasetto S., Spekkens K., Caldwell N., Guhathakurta P., et al., 2019, ApJ, 872, 80
  * Danieli et al. (2017) Danieli S., van Dokkum P., Merritt A., Abraham R., Zhang J., Karachentsev I. D., Makarova L. N., 2017, ApJ, 837, 136
  * Deason et al. (2013) Deason A. J., Belokurov V., Evans N. W., Johnston K. V., 2013, ApJ, 763, 113
  * Diemand et al. (2007) Diemand J., Kuhlen M., Madau P., 2007, ApJ, 657, 267
  * Elahi et al. (2009) Elahi P. J., Widrow L. M., Thacker R. J., 2009, Ph. Rev. D, 80, 123513
  * Engler et al. (2021) Engler C., Pillepich A., Pasquali A., Nelson D., Rodriguez-Gomez V., Chua K. T. E., Grebel E. K., et al., 2021, MNRAS, 507, 4211
  * Fakhouri, Ma, & Boylan-Kolchin (2010) Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267
  * Font et al. (2020) Font A. S., McCarthy I. G., Poole-Mckenzie R., Stafford S. G., Brown S. T., Schaye J., Crain R. A., et al., 2020, MNRAS, 498, 1765
  * Font, McCarthy, & Belokurov (2021) Font A. S., McCarthy I. G., Belokurov V., 2021, MNRAS, 505, 783
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown A. G. A., Vallenari A., Prusti T., de Bruijne J. H. J., Babusiaux C., Bailer-Jones C. A. L., et al., 2018, A&A, 616, A1
  * Gao et al. (2008) Gao L., Navarro J. F., Cole S., Frenk C. S., White S. D. M., Springel V., Jenkins A., et al., 2008, MNRAS, 387, 536
  * Gao et al. (2011) Gao L., Frenk C. S., Boylan-Kolchin M., Jenkins A., Springel V., White S. D. M., 2011, MNRAS, 410, 2309
  * Gao et al. (2012) Gao L., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 419, 1721
  * Geha et al. (2017) Geha M., Wechsler R. H., Mao Y.-Y., Tollerud E. J., Weiner B., Bernstein R., Hoyle B., et al., 2017, ApJ, 847, 4
  * Giocoli et al. (2008) Giocoli C., Tormen G., van den Bosch F. C., 2008, MNRAS, 386, 2135
  * Giocoli et al. (2010) Giocoli C., Tormen G., Sheth R. K., van den Bosch F. C., 2010, MNRAS, 404, 502
  * 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
  * Grand et al. (2021) Grand R. J. J., Marinacci F., Pakmor R., Simpson C. M., Kelly A. J., Gómez F. A., Jenkins A., et al., 2021, MNRAS.tmp
  * Graus et al. (2019) Graus A. S., Bullock J. S., Kelley T., Boylan-Kolchin M., Garrison-Kimmel S., Qi Y., 2019, MNRAS, 488, 4585
  * Griffen et al. (2016) Griffen B. F., Ji A. P., Dooley G. A., Gómez F. A., Vogelsberger M., O’Shea B. W., Frebel A., 2016, ApJ, 818, 10
  * Han et al. (2016) Han J., Cole S., Frenk C. S., Jing Y., 2016, MNRAS, 457, 1208
  * Han et al. (2018) Han J., Cole S., Frenk C. S., Benitez-Llambay A., Helly J., 2018, MNRAS, 474, 604
  * Hellwing et al. (2016) Hellwing W. A., Frenk C. S., Cautun M., Bose S., Helly J., Jenkins A., Sawala T., et al., 2016, MNRAS, 457, 3492
  * Henry (2000) Henry, J. P., 2000, ApJ, 534, 565
  * Ishiyama et al. (2020) Ishiyama T., Prada F., Klypin A. A., Sinha M., Metcalf R. B., Jullo E., Altieri B., et al., 2020, arXiv, arXiv:2007.14720
  * Jiang & van den Bosch (2016) Jiang F. & van den Bosch F. C. 2016, MNRAS, 458, 2848
  * Jiang et al. (2021) Jiang F., Dekel A., Freundlich J., van den Bosch F. C., Green S. B., Hopkins P. F., Benson A., et al., 2021, MNRAS, 502, 621
  * Juan et al. (2014a) Juan E., Salvador-Solé E., Domènec G., Manrique A., 2014, MNRAS, 439, 719
  * Juan et al. (2014b) Juan E., Salvador-Solé E., Domènech G., Manrique A., 2014, MNRAS, 439, 3156
  * Kampakoglou & Benson (2007) Kampakoglou M. & Benson A. J., 2007, MNRAS, 374, 775
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
  * Klypin et al. (2011) Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
  * Komatsu et al. (2011) Komatsu E., Smith K. M., Dunkley J., Bennet C. L., Gold B., Hinshaw G., Jarosik N., et al. others, 2011, ApJS, 192, 18
  * Kravtsov (2013) Kravtsov A. V., 2013, ApJL, 764, L31
  * Lacey & Cole (1994) Lacey, C., & Cole, S. 1994, MNRAS, 271, 671
  * Lancaster, Belokurov, & Evans (2019) Lancaster L., Belokurov V., Evans N. W., 2019, MNRAS, 484, 2556
  * Lee (2004) Lee J., 2004, ApJ, 604, L73
  * Lovell et al. (2014) Lovell M. R., Frenk C. S., Eke V. R., et al., 2014, MNRAS, 439, 300
  * Lovell et al. (2021) Lovell M. R., Cautun M., Frenk C. S., Hellwing W. A., Newton O., 2021, MNRAS, 507, 4826
  * Manrique & Salvador-Solé (1995) Manrique A. & Salvador-Solé E., 1995, ApJ, 453, 6
  * Manrique & Salvador-Solé (1996) Manrique A. & Salvador-Solé E., 1996, ApJ, 467, 504
  * Manrique et al. (1998) Manrique A., Raig A., Solanes J. M., González-Casado G., Stein, P., Salvador-Solé E., 1998, ApJ, 499, 548
  * Mao et al. (2020) Mao Y.-Y., Geha M., Wechsler R. H., Weiner B., Tollerud E. J., Nadler E. O., Kallivayalil N., 2020, arXiv, arXiv:2008.12783
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Oguri & Lee (2004) Oguri M. & Lee J. 2004, MNRAS, 355, 120
  * Onions et al. (2012) Onions J., Knebe A., Pearce F. R., Muldrew S. I., Lux H., Knollmann S. R., Ascasibar Y., et al., 2012, MNRAS, 423, 1200
  * Peñarrubia & Benson (2005) Peñarrubia J. & Benson A. J., 2005, MNRAS, 364, 977
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425
  * Pullen et al. (2014) Pullen A. R., Benson A. J., Moustakas L. A., 2014, ApJ, 792, 24
  * Raig et al. (2001) Raig, A., González-Casado, G., Salvador-Solé, E. 2001, MNRAS, 327, 939
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893
  * Richings et al. (2020) Richings J., Frenk C., Jenkins A., Robertson A., Fattahi A., Grand R. J. J., Navarro J., et al., 2020, MNRAS, 492, 5780
  * Ruchti et al. (2015) Ruchti G. R., Read J. I., Feltzing S., Serenelli A. M., McMillan P., Lind K., Bensby T., et al., 2015, MNRAS, 450, 2874
  * Salvador-Solé, Solanes, & Manrique (1998) Salvador-Solé E., Solanes J. M., Manrique A., 1998, ApJ, 499, 542
  * Salvador-Solé & Manrique (2021) Salvador-Solé E., Manrique A., 2021, ApJ, 914,141
  * Salvador-Solé et al. (2012b) Salvador-Solé E., Serra S., Manrique A., González-Casado G., 2012b, MNRAS, 424, 3129 (SSMG)
  * Salvador-Solé et al. (2012a) Salvador-Solé E., Viñas J., Manrique A., Serra S., 2012a, MNRAS, 423, 2190 (SVMS)
  * Salvador-Solé et al. (2021a) Salvador-Solé E., Manrique A., Botella I., 2021a, in press in MNRAS(Paper I)
  * Salvador-Solé et al. (2021b) Salvador-Solé E., Manrique A., Botella I., 2021b, in press in MNRAS(Paper II)
  * Salvador-Solé et al. (2021c) Salvador-Solé E., Canales D., Manrique A., Juan E., Botella I., 2021c, in preparation
  * Salvador-Solé et al. (2021d) Salvador-Solé E., Manrique A., Canales D., Botella I., 2021d, in preparation
  * Samuel et al. (2020) Samuel J., Wetzel A., Tollerud E., Garrison-Kimmel S., Loebman S., El-Badry K., Hopkins P. F., et al., 2020, MNRAS, 491, 1471
  * Sheth (2003) Sheth R. K., 2003, MNRAS, 345, 1200
  * Smercina et al. (2018) Smercina A., Bell E. F., Price P. A., D’Souza R., Slater C. T., Bailin J., Monachesi A., et al., 2018, ApJ, 863, 152
  * Springel et al. (2008) Springel V., Wang J., Vogelsberger M., et al., 2008a, MNRAS, 391, 1685
  * Sugiyama (1995) Sugiyama, N. 1995, ApJS, 100, 281
  * Taylor & Babul (2001) Taylor J. E. & Babul A., 2001, ApJ, 559, 716
  * Taylor & Babul (2004) Taylor J. E. & Babul A., 2004, MNRAS, 348, 811
  * Tollerud, Boylan-Kolchin, & Bullock (2014) Tollerud E. J., Boylan-Kolchin M., Bullock J. S., 2014, MNRAS, 440, 3511
  * van den Bosch et al. (2005) van den Bosch, F. C., Tormen G., Giocoli C., 2005, MNRAS, 359, 1029
  * van den Bosch & Jiang (2016) van den Bosch F. C. & Jiang F., 2016, MNRAS, 458, 2870
  * van den Bosch et al. (2018) van den Bosch F. C., Ogiya G., Hahn O., Burkert A., 2018, MNRAS, 474, 3043
  * Wang et al. (2011) Wang J., Navarro J. F., Frenk C. S., et al., 2011, MNRAS, 413, 1373
  * Wyse (2001) Wyse R. F. G., 2001, ASPC, 230, 71
  * Zentner & Bullock (2003) Zentner A. R., Bullock, J. S., 2003, ApJ, 598, 49
  * Zentner et al. (2005) Zentner A. R., Berlin A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H., 2005, ApJ, 624, 505