# Semi-analytical frameworks for subhalos from the smallest to the largest
scale

Nagisa Hiroshima, 1,2 Shin’ichiro Ando,3,4 and Tomoaki Ishiyama5
1Department of Physics, University of Toyama, 3190 Gofuku, Toyama 930-8555,
Japan
2RIKEN iTHEMS, Wako, Saitama 351-0198, Japan
3GRAPPA Institute, Institute of Physics, University of Amsterdam, 1098 XH
Amsterdam, The Netherlands
4Kavli Institute for the Physics and Mathematics of the Universe (WPI),
University of Tokyo, Kashiwa 277-8583, Japan
5Institute of Management and Information Technologies, Chiba University, Chiba
263-8522, Japan

###### Abstract

Substructures of dark matter halo, called subhalos, provide important clues to
understand the nature of dark matter. We construct a useful model to describe
the properties of subhalo mass functions based on the well-known analytical
prescriptions, the extended Press-Schechter theory. The unevolved subhalo mass
functions at arbitrary mass scales become describable without introducing free
parameters. The different host halo evolution histories are directly recast to
their subhalo mass functions. As applications, we quantify the effects from
(i) the Poisson fluctuation, (ii) the host mass scatter, and the (iii)
different tidal evolution models on observables in the current Universe with
this scheme. The Poisson fluctuation dominates in the number count of the mass
ratio to the host of $\sim{\cal O}(10^{-2})$, where the intrinsic scatter is
smaller by a factor of a few. The host-mass scatter around its mean does not
affect the subhalo mass function. Different models of the tidal evolution
predict a factor of $\sim 2$ difference in numbers of subhalos with $\la{\cal
O}(10^{-5})$, while the dependence of the Poisson fluctuation on the tidal
evolution models is subtle. The scheme provides a new tool for investigating
the smallest-scale structures of our Universe which are to be observed in near
future experiments.

###### keywords:

dark matter — cosmology: theory — galaxies: haloes

††pagerange: Semi-analytical frameworks for subhalos from the smallest to the
largest scale–Semi-analytical frameworks for subhalos from the smallest to the
largest scale

##  1 Introduction

The nature of dark matter (DM) is still elusive, while the evidence for its
existence is strong (e.g. rotation curves of galaxies (van Albada et al.,
1985), bullet clusters (Clowe et al., 2004), etc.). One of the strongest, and
the original motivations for DM comes from the structure formation history of
our Universe (Peebles, 1982; Aghanim et al., 2020). New particles beyond the
Standard Model as well as non-particle solutions, such as primordial black
holes, which source the gravitational potential for baryonic structures to
evolve are required to explain the current structures. Various possibilities
are considered for solving this long-standing problem (see Bertone et al.
(2005); Bertone & Tait (2018); Carr et al. (2021) for reviews).

The hierarchical structure of halos that form through the gravitational
collapse of matter overdensities is a key to access the nature of DM. The
structure formation proceeds through mergers, accretion, and tidal stripping.
Its statistical property reflects the nature of particle DM. For example, the
minimum halo mass of the halo ranges from ${\cal O}(10^{-12})$–${\cal
O}(10^{-3})M_{\odot}$ depending on parameters of the weakly interacting
massive particles (Diamanti et al., 2015; Bringmann, 2009; Hofmann et al.,
2001; Green et al., 2004; Loeb & Zaldarriaga, 2005; Profumo et al., 2006).
Introduction of the self-interaction in DM sector or velocity dispersion could
erase the small-scale structures (Spergel & Steinhardt, 2000; Buckley et al.,
2014; Nadler et al., 2020; Bond et al., 1980; Hogan & Dalcanton, 2000). The
importance of the clustering of the small-scale halos in its host halo appears
in another context of the indirect searches of DM. The annihilation signals of
DM should be enhanced by a factor of several due to the clustering effects,
which are known as the boost factor (SÃ¡nchez-Conde & Prada, 2014; Ishiyama,
2014; Moline et al., 2017; Hiroshima et al., 2018; Ishiyama & Ando, 2020). The
determination of the slope of the subhalo mass function as well as the minimum
halo mass is crucial in its estimation (see Ando et al. (2019) for a recent
review).

Numerical (Gao et al., 2004; Diemand et al., 2007; Giocoli et al., 2008a;
Ishiyama et al., 2021; Moline et al., 2021), semi-analytical (Penarrubia &
Benson, 2005), and analytical (Giocoli et al., 2008b; van den Bosch et al.,
2005; Yang et al., 2011) studies have investigated the properties of subhalos
such as their mass functions. Numerical simulations are advantageous in
treating all effects determining the properties, while they are limited by
numerical resolutions and computational costs. Semi-analytical models cure
such difficulties by approximating several parts of the physical processes.
The merger tree based on the extended Press-Schechter (EPS) theory (Bond et
al., 1991; Lacey & Cole, 1993), which is a theoretically-motivated extension
of the original Press-Schechter model (Press & Schechter, 1974), takes an
essential role in the semi-analytical approaches.

Previous studies have developed schemes to reconstruct the so-called unevolved
mass function, which is the mass function corresponding to that before the
tidal evolution, in semi-analytical ways. Based on the EPS prescription, Zhao
et al. (2009) has derived a fitting formula to follow the main branch of the
merger tree calibrated against the numerical simulation. The main branch of a
halo is obtained by repeatedly connecting the most massive progenitor of the
most massive progenitor. The work by Zhao et al. (2009) is extended in Yang et
al. (2011) to include the scatter around the main branch of the tree, showing
that the unevolved subhalo mass function is well-accessible in this scheme. On
the other hand, Correa et al. (2015a, b) have derived analytical formulae
describing the evolution of the main branch halos in a way that does not
require the fitting to the simulation results. It is achieved by examining the
growth rate of the initial density perturbations.

In order to access the astrophysical observables such as the number counts of
the satellite galaxies, those models must be implemented with the prescription
for tidal evolution after the subhalo accretion, as pointed out in Giocoli et
al. (2008b). van den Bosch et al. (2005) has proposed to adopt the orbit-
averaged mass-loss rate to treat the tidal evolution in a semi-analytical way.
Jiang & van den Bosch (2016) calculate the mass-loss rate by segmenting the
time steps for treating the host evolution. The free parameters of the mass-
loss model are tuned to reproduce the evolved subhalo mass function in
numerical simulations in their succeeding works (Jiang & van den Bosch, 2014,
2016). The normalization of the subhalo mass function before the tidal
evolution had to be fitted using numerical results as those in previous works
(Giocoli et al., 2008a; Li & Mo, 2009).

In this work, we propose an EPS-based analytical scheme to access the
statistics of subhalos making full use of its advantage. The construction of
the unevolved subhalo mass function in this scheme requires no parameters to
be fitted. The evolved subhalo mass function is obtained by implementing the
orbit-averaged tidal evolution model to the unevolved mass functions, which is
a prescription proposed in van den Bosch et al. (2005); Jiang & van den Bosch
(2014, 2016). Unlike the earlier work that focused on constructing full merger
trees including second and third order branches, we only follow the main
branch of the merger tree. This approach enables us to construct the subhalo
mass functions without being limited by the numerical resolution. We can fully
adopt the analytical formula of the subhalo accretion rates combined with the
host evolution, both of which are described by the EPS theory. The downside of
this approach is incapability of tracing higher order subhalos: i.e., subhalos
in a parent subhalo. However, our interest in this work is the properties of
subhalos at the first order. Those at higher orders are still very difficult
to resolve with numerical simulations or to identify observationally. As we
explain below, we evaluate the subhalo statistics using 500 of host
realizations, which reasonably treats the halo-to-halo scatter. In our
calculation, the computational cost is highly reduced and an unlimited range
of the subhalo-to-halo mass ratio can be covered simultaneously. We quantify
the numerical or observational effects of host-mass scatter and the Poisson
fluctuation, which we need to manage in extracting the indication about the
intrinsic properties of DM as well as the dependence of the observables on the
tidal-evolution model based on this work.

This work improves the host-subhalo connection from our previous works
(Hiroshima et al., 2018; Ando et al., 2019) which provides an access to the
subhalo observables adopting semi-analytical method. In the following, we
directly connect the different host evolution histories to those of subhalo
mass functions. This reduces the needs to rely on fitting formulae calibrated
against numerical simulations and the accessibility to the statistical
quantities is enhanced. Also, the scheme paves the way for halo mass functions
in much general cosmological setups.

The structure of this article is as follows. Sec. 2 is devoted for the
simulation we have performed to compare with the analytical calculations.
After summarizing the host evolution in the EPS framework in Sec. 3, we show
the prediction about the subhalo mass function in Sec. 4. The scatter of the
subhalo mass function due to the extrinsic effects is investigated in Sec. 5.
The dependence on the tidal-evolution model is also discussed in this section.
We finally conclude in Sec. 6. Throughout the paper, we use “$\ln$" and
“$\log$" as notations for natural and 10-based logarithmic quantity. We assume
cosmological parameters of Ade et al. (2016) (Table 4, “TT+lowP+lensing”)
unless mentioned specifically.

##  2 Numerical simulation

In this section, we explain our numerical simulations used for deriving the
host evolution and subhalo mass function to compare with the EPS-based
calculations. Numerical results provided in this section appear in Secs. 3, 4,
and 5.

We use a high-resolution cosmological $N$-body simulation, Shin-Uchuu
(Ishiyama et al., 2021). The Shin-Uchuu simulation consists of $6400^{3}$ dark
matter particles with a mass of $8.97\times 10^{5}\,h^{-1}M_{\odot}$ in a
comoving cubic box with a side length of 140$\,h^{-1}\rm Mpc$. The
gravitational softening length is 0.4$\,h^{-1}\rm kpc$. The cosmological
parameters are $\Omega_{0}=0.3089$, $\Omega_{\rm b}=0.0486$,
$\lambda_{0}=0.6911$, $h=0.6774$, $n_{\rm s}=0.9667$, and $\sigma_{8}=0.8159$,
which are consistent with observational results of cosmic microwave background
by the Planck satellite (Aghanim et al., 2020).

rockstar halo/subhalo finder (Behroozi et al., 2013a) was used to find
gravitational bound halos and subhalos. Then consistent trees merger tree code
(Behroozi et al., 2013b) was used to construct their merger trees. Those data
are available on Skies & Universes site. Further details of the Shin-Uchuu
simulation and its basic halo and subhalo properties are presented in the
literature (Ishiyama et al., 2021; Moline et al., 2021).

##  3 Host halo evolution

Our analytical calculation of the host evolution is based on the prescription
in Yang et al. (2011). We summarize the scheme in this section. Focusing on
the Milky-Way-like host at $z=0$, we follow the evolution path of the halo
backward in time up to $z=10$. Concretely, we follow the steps below.

  1. 1.

We set the host mass at $z=0$ and prepare $N_{\rm host}$ of host halos. We use
$N_{\rm host}=500$ halo realizations for each single mass host in our
calculations. We set our fiducial host mass as $M_{\rm host}(z=0)=1.30\times
10^{12}M_{\odot}$, adopting the convention of the $M_{200}$ that is the mass
enclosed in a radius $r_{200}$ within which the average density is 200 times
the critical density. This convention matches those in the literature we refer
to in our calculations.

  2. 2.

We prepare the steps of the redshift to sample the host mass. We take a
constant interval in the $\ln(1+z)$ space taking $N_{z}=800$ points from $z=0$
to $z=10$. In the following, the step size is denoted with $\Delta z$.

  3. 3.

The distribution of the accreting mass is determined using the EPS formula by
specifying the redshift and the host mass. Following Yang et al. (2011), the
number distribution of the subhalos that accreted in the redshift interval
between $z+\Delta z$ and $z$ is given as

$\displaystyle F(m)=$ $\displaystyle f(s(m),\delta(z+\Delta
z)|S(M),\delta(z))$ (1) $\displaystyle\times
G\left(\frac{\sqrt{S(M)}}{\sqrt{s(m)}},\frac{\delta(z+\Delta
z)}{\sqrt{s(m)}}\right)\left|\frac{ds}{dm}(m)\right|\left(\frac{M}{m}\right)$

where

$\displaystyle f(s(m),\delta(z+\Delta z)|S(M),\delta(z))$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\frac{\delta(z+\Delta
z)-\delta(z)}{\left[s(m)-S(M)\right]^{3/2}}\exp\left[-\frac{\left(\delta(z+\Delta
z)-\delta(z)\right)^{2}}{2(s(m)-S(M))}\right].$ (2)

Eq. 1 gives the number of the subhalo in the mass range of $[m,\ m+dm]$ that
have accreted in the redshift interval of $[z+\Delta z,\ z]$, i.e. the time
interval corresponding to the host evolution from $M(z+\Delta z)$ to $M(z)$.
The factor $M/m$ appears in the last part of the equation to obtain the number
fraction of mass $m$ from its mass fraction. Each quantity in Eq. 1 is:

     * •

$\delta(z)=1.686/D(z)$: the linear growth factor.

     * •

$S(M)$ and $s(m)$: the variance of the linear density field, $\sigma^{2}(M)$
and $\sigma^{2}(m)$ at $z=0$. Each of the mass scale $M$ and $m$ corresponds
to that of the host and accreting halo at that redshift, respectively.

     * •

$G(x,y)=G_{0}x^{\gamma_{1}}y^{\gamma_{2}}$: a correction function for the EPS
formula tuned against the numerical simulations.

We adopt the parametrized form of the $D(z)$ and $\sigma(m)$ provided in
Ludlow et al. (2016). The parameter values of the correction function $G$ are
set to $(G_{0},\ \gamma_{1},\ \gamma_{2})=(0.57,\ 0.38,\ -0.01)$, which we
have taken from Parkinson et al. (2008) and Cole et al. (2008). Those works
have tuned the parameters with the Millennium $N$-body simulation (Springel et
al., 2005) which covers the halo mass in the range of $[\sim
10^{10}h^{-1}M_{\odot},\ \sim 10^{16}h^{-1}M_{\odot}]$ and redshift up to
$z\sim 10$, and suitable for generating the merger tree (Jiang & van den
Bosch, 2014) .

We regard the most massive subhalo in the redshift interval $[z,z+\Delta z]$
as the progenitor of the host of $M(z)$ at $z$. When the accreting halo mass
$m$ of Eq. 1 satisfies $m>M(z)/2$, the hierarchy of halo assures that the halo
to be the unique progenitor. We determine the evolution history by following
the progenitor adopting the inverse transform method.

  4. 4.

Prepare subhalos in the mass range of $0.5M(z)\leq m<M(z)$ taking $10^{5}$
points in linear space.

  5. 5.

Construct the cumulative distribution function by integrating Eq. 1 in the
above mass range, then normalize it.

  6. 6.

Generate a random number following the uniform distribution in the interval of
$[0,1]$, which are assigned to a point in the normalized cumulative
distribution function in the above calculation step. We identify the
progenitor at $z+\Delta z$ as that at the corresponding mass in the normalized
cumulative mass function.

  7. 7.

Repeat the same procedure succeedingly to redshift $z=10$.

In Fig. 1, we show the resulting evolution of our fiducial hosts. The solid
line is the median of the $N_{\rm host}=500$ realizations and thin lines are
the histories of 100 host realizations. We also show the 16% (2.5%) and 84%
(97.5%) quantiles with dashed (dotted) lines. Square points with error bars
are the results obtained in numerical simulations which we have explained in
Sec. 2. As it is seen in the figure, analytical and numerical calculations
agree within the range of the error bar in both the median and scatter while
the median in the analytical model is slightly smaller. The host evolution
history also agrees with the fitting function in Correa et al. (2015a) in the
regime where the function is tuned against simulations. In this regime, the
functional form is similar to that predicted in Wechsler et al. (2002) and
similar results are obtained also in previous works ( Bosch (2002); Zhao et
al. (2009); Giocoli et al. (2012).) All the $N_{\rm host}=500$ realization of
the host evolution history is used to derive the subhalo mass function in the
next section.

Figure 1: Evolution history determined in the EPS formalism for hosts of
$M=1.30\times 10^{12}M_{\odot}$ at $z=0$. The blue solid line is the median.
Dashed (dotted) lines correspond to 16% and 84% (2.5% and 97.5%) quantiles.
Each realization of the host is shown in thin (cyan) lines. We only plots 100
of the total 500 host realizations in order for the visibility. Numerical
results are shown in square points, of which thick (thin) error bars represent
the quantiles same as dashed (dotted) lines for EPS-based analytical
calculations. The analytical calculation slightly underestimates the median of
the history while it agrees with the simulation within the range of the error
bar.

##  4 Subhalo mass function

In this section, we explain our self-consistent framework to derive subhalo
mass functions reflecting the different host evolution histories. We start
from constructing unevolved mass functions based on the host evolution history
(see Sec. 3) in Sec. 4.1. Assuming that the increment of the host mass is
solely determined by accretions of smaller halos, the normalization of the
unevolved subhalo mass function at each redshift step is determined without
introducing any free parameters. The unevolved mass function at $z=0$ is then
obtained by summing up contributions from accreting halos in all redshift
steps. The evolved subhalo mass function, which directly relates to the
observables in the current Universe, is calculated by combining the unevolved
mass function with the tidal-evolution model in Hiroshima et al. (2018). Since
we hold the mass of the host in each redshift step and realization, the host-
evolution dependent mass-loss becomes calculable, as it is explained in Sec.
4.2.

###  4.1 Unevolved mass function

  1. 1.

Set a single value for the host halo $M_{\rm host}$ at $z=0$.

  2. 2.

Generate the subhalo sample of which mass ranges of $[m_{\rm min},\ m_{\rm
max}]=[10^{-18},\ 1]\times M_{\rm host}(z=0)$ for each step of the redshift
and host realization. The minimum mass of the sample subhalo can be
arbitrarily chosen and we set the value so that to match the prediction in
neutralino DM scenarios (Green et al., 2004). We take $N_{\rm sh}=300$
sampling points for the accreting subhalos. The interval of the accreting mass
sampling is constant in the $\ln(m_{\rm 200})$ space. The virial-quantities
are used for calculating the evolved mass function as will explain in Sec.
4.2.

  3. 3.

Determine the total accreted mass during the interval of $[z+\Delta z,z]$ as
$M_{\rm acc}(z)=M(z)-M(z+\Delta z)$ by directly using the host evolution data.

  4. 4.

Count the number of subhalo in a similar manner with that explained in Sec. 3.
Taking the $\delta(z)$ and $s(m)=\sigma^{2}(m,z=0)$ as the variables for the
redshift and mass, the distribution function is

$\displaystyle\left.\frac{dN}{dm}\right|_{z}=$ $\displaystyle
f(s(m),\delta(z+\Delta z)|S(M),\delta(z))$ (3) $\displaystyle\times
G\left(\frac{\sqrt{S(M)}}{\sqrt{s(m)}},\frac{\delta(z+\Delta
z)}{\sqrt{s(m)}}\right)\frac{1}{m}\left|\frac{ds}{dm}\right|\frac{dM}{(1+z)},$

where $dM=M(z)-M(z+\Delta z)=M_{\rm acc}(z)$. Eq. 3 gives the distribution of
the accreting subhalo at the redshift step $z$. As it is seen from the
functional form of Eq. 3, the distribution function diverges at $m\sim M(z)$.
We discard the regime where the distribution function increases with $m$,
regarding that such a regime is not physical.111Our treatment is close to that
proposed as the Model III in Yang et al. (2011) in practice, while we do not
need to introduce extensional parameters.

  5. 5.

Normalize the distribution $dN/dm$ by applying the condition

$M_{\rm acc}(z)=\sum m(z)dm\left.\frac{dN}{dm}\right|_{z}.$ (4)

The sum is taken for subhalos remaining after the cut in the previous step.

  6. 6.

Integrate the mass function over the redshift interval of $z=[10,0]$.

The above procedure of the calculation enables us to obtain the unevolved mass
functions reflecting the different history of host evolutions. There are no
needs for tuning parameters in the above calculation procedure.

Figure 2 shows the unevolved mass function at $z=0$ derived in the above
procedure. In the left panel, we show the differential mass function which
covers a wider range for the subhalo of
$m=[10^{-4}M_{\odot},10^{12}M_{\odot}]$, while the cumulative mass function in
a narrower range of $m/M_{\rm host}(z=0)=[10^{-4},1]$ in the right panel. We
take 50 and 100 mass bins for each mass function, respectively. The wavy
features in the panels are caused by the binning effect. The solid (dashed)
line corresponds to the average ($1\sigma$ scatter) of 500 host realizations
in each panel. The scatter of the mass function which originates from the
different host evolution histories is small and its amount is comparable to
those obtained in numerical simulations. The comparison with the numerical
subhalo mass functions is shown in Sec. 5.2.

---
Figure 2: Differential (left) and cumulative (right) unevolved mass function
at $z=0$ for host of $M(z=0)=1.30\times 10^{12}M_{\odot}$. The solid (dashed)
line shows the average ($1\sigma$ scatter) for 500 realizations of the host.
Non-smooth features are induced by the binning effect and can be removed when
we change the mass bins for plotting.

###  4.2 Evolved mass function

Unevolved mass functions obtained in Sec. 4.1 is then combined with the
analytical models of the tidal evolution to determine the evolved mass
functions. The model of Hiroshima et al. (2018) is based on the orbit-averaged
mass-loss prescriptions (Jiang & van den Bosch, 2014, 2016) and capable of
treating the dependence of the mass-loss rate on the host mass and the
redshift. We starts the calculation from the data of the host evolution and
the accreting mass distribution generated in the previous sections. The data
is in the form of the $(N_{\rm sh},N_{z})$ matrix, tracking the history from
high to low redshifts in the following manner:

  1. 1.

Convert the $M_{200}$ mass, which is adopted in previous sections, to the
virial one222Functions of $\delta(z)$ and $\sigma(M)$ in Ludlow et al. (2016)
are given in terms of $M_{200}$ while the tidal prescription in this section
is based on the virial quantities. The transformation from $M_{200}$ to
$M_{\rm vir}$ keeps the consistency of the formalism.. The concentration-mass
relation of Correa et al. (2015a, b) are applied to determine the virial mass
of host and subhalos. Assuming the NFW halo profile (Navarro et al., 1996), we
use the following quantities

$\displaystyle r_{\rm vir}$ $\displaystyle=$ $\displaystyle\left(3m_{\rm
vir}\times\left(4\pi\rho_{c}(z)\varDelta(z)\right)^{-1}\right)^{1/3},$ (5)
$\displaystyle c_{\rm vir}$ $\displaystyle=$ $\displaystyle
c_{200}\frac{r_{\rm vir}}{r_{200}},$ (6) $\displaystyle r_{\rm s}$
$\displaystyle=$ $\displaystyle r_{\rm vir}/c_{\rm vir},$ (7)
$\displaystyle\rho_{\rm s}$ $\displaystyle=$ $\displaystyle\frac{m_{\rm
vir}}{4\pi r_{\rm s}^{3}f(c_{\rm vir})},$ (8) $\displaystyle r_{\rm max}$
$\displaystyle=$ $\displaystyle 2.163r_{\rm s},$ (9) $\displaystyle V_{\rm
max}$ $\displaystyle=$ $\displaystyle\sqrt{\frac{4\pi G\rho_{\rm
s}}{4.625}}r_{\rm s}.$ (10)

The $\varDelta(z)$ which describes the overdensity corresponding to the virial
mass is taken from Bryan & Norman (1998). By assuming the NFW model, the halo
profile is described with two parameters: the characteristic density
$\rho_{\rm s}$ and the scale radius $r_{\rm s}$. The set of parameters
$(\rho_{\rm s},r_{\rm s})$ has one-to-one correspondence with the set $(V_{\rm
max},r_{\rm max})$. The function $f(c)=\ln(1+c)-c/(1+c)$ relates the NFW
profile parameter and the halo mass in a direct way.

  2. 2.

Track the evolution of accreted halos starting from $z=10$. The part of the
mass matrix of subhalo corresponding to those have ever accreted before the
redshift step of the interest enters the calculation. We calculate the tidal
mass loss during the interval $[z+\Delta z,z]$ as

$\Delta m(z)=\frac{\Delta z}{(1+z)^{2}H(z)}\times\left[A(M_{\rm
host},z)\frac{m}{\tau_{\rm dyn}}\left(\frac{m}{M_{\rm
host}(z)}\right)^{\zeta(M_{\rm host},z)}\right].$ (11)

The host mass is evaluated at each redshift of the calculation by directly
referring to the host evolution data. We take the time-dependent mass-loss
model of Hiroshima et al. (2018) as our fiducial one. The factor
$1/(1+z)/H(z)$ appears to use the redshift as the time variable for the mass
loss. We discuss the influence of the uncertainty in the tidal evolution model
in Sec. 5.3.

  3. 3.

Replace the subhalo mass matrix with that of after the tidal mass-loss, i.e.,
the matrix of $m(z)$ is replaced with $m(z)-\Delta m$. If $\Delta m(z)\geq
m(z)$, we set the mass of that halo to zero.

  4. 4.

Go to the next redshift step of $z-\Delta z$. The treatment is the same as
that of the one-step before, but the mass matrix describing subhalos those
have accreted by redshift $z$ is already replaced with those after the tidal
mass-loss down to that redshift step. Also, a new part of the mass matrix
corresponding to those accrete between $[z,\ z-\Delta z]$ enters the
calculation. The process is repeated to $z=0$. This treatment enables us to
calculate the mass-loss of subhalos including the information about their
accretion redshift.

After finishing the calculation, we obtain the mass matrix of all the subhalos
which have accreted between $z=[10,0]$ and evolved. The $V_{\rm max}$ and
$r_{\rm max}$ for NFW halo before and after the tidal stripping are related as
(Penarrubia et al., 2010)

$\displaystyle\frac{V_{\rm max}(z=0)}{V_{\rm max}(z_{\rm acc})}$
$\displaystyle=$ $\displaystyle\frac{2^{0.4}x^{0.3}}{(1+x)^{0.4}},$ (12)
$\displaystyle\frac{r_{\rm max}(z=0)}{r_{\rm max}(z_{\rm acc})}$
$\displaystyle=$ $\displaystyle\frac{2^{-0.3}x^{0.4}}{(1+x)^{-0.3}},$ (13)

where $x=m(z=0)/m(z_{\rm acc})$. The relationships for $(V_{\rm max},r_{\rm
max})$ before and after the tidal stripping are recasted to those for
$(\rho_{\rm s},r_{\rm s})$. Using $m$, $r_{\rm s}$, and $\rho_{\rm s}$ at
$z=0$, we evaluate the truncation radius $r_{t}=c_{t}r_{\rm s}$ as

$m=4\pi\rho_{\rm s}r_{t}^{3}f(c_{t}).$ (14)

If $c_{t}<0.77$, we omit those halos regarding them to be fully disrupted
(Hayashi et al., 2003). The mass matrix of the evolved subhalo is integrated
over the redshift with the same weight as that for unevolved mass function
after this cut.

We show the evolved subhalo mass functions in the virial mass unit in Fig. 3.
The host mass in the virial unit is $M_{\rm vir}(z=0)=1.5\times
10^{12}M_{\odot}$, which is the same as that for unevolved mass function of
$M_{\rm host}(z=0)=1.30\times 10^{12}M_{\odot}$ in the $M_{200}$ unit. The
number distribution of the accreting subhalo is the same as that for unevolved
mass function. We just re-count the subhalos in each mass bin at $z=0$ using
the same assignment for the mass bins with that of Fig. 2. The evolved mass
function is composed of all the survived subhalos which have been accreted and
evolved between $z=10$ and $=0$. While we only consider the subhalos which
have accreted in $z\leq 10$, the extension of the calculation to a higher
redshift does not change the results. This is due to the fact that subhalos
that have accreted in such a high redshift are fully disrupted and do not
contribute much to the mass function.

---
Figure 3: The same as Fig. 2 but for evolved subhalos. The scatter becomes
apparent reflecting the different host evolution histories. In this figure,
the mass is measured in the virial mass for both the host and subhalos.

After the tidal evolution, the total number of the subhalo decreases. Also,
different from the unveolved cases, the scatter of the subhalo mass function
becomes apparent, reflecting the different host evolution histories. The
height of the subhalo mass function is similar to those in previous works
(Jiang & van den Bosch, 2017; Hiroshima et al., 2018).

##  5 Quantifying the source of the scatter

In the previous section, we have explained our basic formalism to derive the
unevolved- and evolved- mass functions. Those mass functions depend only on
the physical properties of DM halos, which are described in the EPS theory.
However, those in either of the observations or simulations include extrinsic
contributions such as the Poisson fluctuation or the host-mass scatters. We
first discuss the Poisson fluctuation and the host-mass scatter effects
separately to quantify each contribution. Then we show the mass function
including both effects, which should be the one compared to those in
observations or simulations. This could be a new scheme to extract intrinsic
property from observables. The uncertainties from the treatments of the tidal
disruption are also discussed in this section.

###  5.1 The effect of the Poisson fluctuation

The Poisson fluctuation in the number count of the subhalo at the massive end
is inevitable for both in the simulations and observations. By removing the
contributions from the Poisson fluctuation, we can retrieve the nature of DM
halo from observables. We provide a mapping between the subhalo mass function
with and without the Poisson fluctuation by implementing this effect to our
mass function derived in Sec. 4 in the following way:

  1. 1.

Construct the cumulative subhalo function for each realization (see Sec. 4).
The cumulative mass function is normalized by dividing it by the total number
of the subhalo, $N_{\rm tot}$.

  2. 2.

Generate a random number following the Poisson distribution of which
expectation value is $N_{\rm tot}$. We denote this generated number as $N_{\rm
tot,P}$.

  3. 3.

Generate $N_{\rm tot,P}$ random numbers assuming a uniform distribution in the
interval of $[0,1]$, then convert the random numbers to the halo mass using
inverse transform method. This step corresponds to generate a new subhalo
distribution of which total number is $N_{\rm tot,P}$ and the distribution
function is the same as that determines the $N_{\rm tot}$.

  4. 4.

Construct the mass function using the re-sampled $N_{\rm tot,P}$ subhalos.

  5. 5.

Repeat the above procedure for $N_{\rm MC}=500$ times for each $N_{\rm tot}$
(i.e. the host realization).

We show the cumulative mass function including the Poisson fluctuation in Fig.
4. The host mass at $z=0$ is the same as previous figures. The left and right
panel corresponds to the unevolved and evolved mass function, respectively. We
limit the regime of the subhalo mass in this figure to clearly see the Poisson
fluctuation effects. As it is seen in Fig. 4, the $1\sigma$ scatter of the
number counts becomes larger by introducing the effect. The Poisson
fluctuation dominates the scatter at $N_{\rm tot}\lesssim 10$. In terms of the
mass ratio, this corresponds to $m/M\gtrsim{\cal O}(10^{-2})$ for evolved
subhalos, which is in a similar range pointed out in the previous work based
on numerical simulations (Boylan-Kolchin et al., 2010). In Fig. 4, we also
plot the results from the numerical simulations of Sec. 2 with square points.
The error bar denotes the $1\sigma$ scatter. As it is seen in the figure, the
analytical subhalo mass functions show good agreements with those in numerical
simulations when the Poisson fluctuation effects are included. However, the
analytical subhalo mass function shows a sharp cutoff at $m/M\sim 0.5M_{\rm
host}$. This feature originates from our treatment allowing only the subhalos
of $m/M(z)$ at each accretion redshift. The prediction from the analytical
calculation is slightly larger compared to the numerical ones. Our assumption
that to determine the normalization of the accreting halo distribution would
be responsible for this difference. In this picture, smooth accretion of the
mass could not be included and hence the mass fraction to be assigned in
subhalos should get larger.

---
Figure 4: The cumulative mass function of unevolved (left) and evolved (right)
subhalos at $z=0$ including the Poisson fluctuation effects. The blue solid
lines show the average, while the dashed lines show the $1\sigma$ range. We
also plot the mass function without the Poisson effects with black lines, of
which $1\sigma$ range is visible in the dotted curve for the evolved case,
while it degenerates with the average for the unevolved one. The average
without the Poisson effect underlies the solid line. Square points with error
bars are results from the $N$-body simulations of Sec. 2.

###  5.2 The effect of the host-mass scatter

The scatter of the host-mass estimate also contributes to the observed subhalo
mass function. In our calculation, the host mass can be fixed to an arbitrary
single value. However, simulated halos distribute in a finite mass range so
careful treatments are required. The mass estimate from observations contains
errors in a certain amount (e.g. Karukes et al. (2020)) which should propagate
to the subhalo mass functions. We take the contributions from the uncertainty
in the host-mass estimate.

In addition to our fiducial host of $M_{200}(z=0)=1.3\times 10^{12}M_{\odot}$,
we calculate the subhalo mass functions for host of $M_{\rm host}=4.43\times
10^{11}M_{\odot}$ and $2.16\times 10^{12}M_{\odot}$. The results are shown in
Fig. 5. The left and right panels corresponds to the unevolved and evolved
cases, respectively. A factor of $\sim 5$ difference in the host mass estimate
affects the cumulative mass function negligibly when it is plotted in terms of
the mass ratio, $m/M$. The behavior is reasonable since we do not put
assumptions that could significantly change the physics of the host and/or
subhalo in the variation scales of the host mass. The results quantitatively
agree with the previous works. In the figure, we plot the fitting function
provided in Jiang & van den Bosch (2016) with red-dotted lines for comparison.
The halo formation time which appears in the fitting function is evaluated as
the redshift where the average host mass is the half of the current value. The
relatively sharp cutoff compared to the fitting formula originates from our
formulation that does not account for mergers of halos with comparable mass
(Moreno & Sheth (2008) and references therin).

---
Figure 5: Unevolved (left) and evolved (right) subhalo mass function for
$M_{200}(z=0)=4.43\times 10^{11}M_{\odot}$ (black, solid), $1.30\times
10^{12}M_{\odot}$ (blue, dashed), and $2.16\times 10^{12}M_{\odot}$ (green,
dot-dashed). Each line shows the average of the 500 host realizations of that
host value. The mass function with different host mass overlaps each other.
The red-dotted line in each panel corresponds to the fitting formula of Jiang
& van den Bosch (2016) assuming $M_{200}=1.3\times 10^{12}M_{\odot}$. For the
formula of the evolved mass function, we estimate the halo formation time as
that the average host mass is the half of the current value.

We also examine how the uncertainty in the fiducial value of the host mass
affect the subhalo mass function, by comparing the subhalo mass functions for
different host masses in terms of the subhalo mass. As it is shown in Fig. 6,
subhalo mass function scales as the mass of the host. The host-mass dependence
in obtained in our scheme is similar to that in previous works. The fitting
function of Jiang & van den Bosch (2016) plotted in this figure explicitly
shows the matching between the results. The maximum mass allowed for the
subhalo increases for more massive hosts, hence the peak of the subhalo mass
shifts with the host mass.

---

Figure 6: Differential mass function of unevolved (left) and evolved (right)
halos. The black solid line corresponds to the case of the $M_{\rm
host}=1.30\times 10^{12}M_{\odot}$, which is the center value we have adopted.
The blue-dashed is for $M_{\rm host}=4.43\times 10^{11}M_{\odot}$ and green
dot-dashed is for $M_{\rm host}=2.16\times 10^{12}M_{\odot}$. The peak
position of the subhalo mass is determined by the host mass. The host
difference between the subhalo mass function roughly scales as that of the
host mass. The red-dotted lines correspond to the formulae of Jiang & van den
Bosch (2016) calculated in the same way with those of Fig. 5

.

###  5.3 Uncertainties from the tidal model

In the above calculations, we have adopted the tidal evolution model for
subhalo, which is characterized with two parameters: the amplitude $A$ and the
dependence to the mass-ratio to the host $\zeta$ (see Eq. 11). The tidal
evolution are still under debate (van den Bosch & Ogiya, 2018; Green & van den
Bosch, 2019), while it is essential for relating the subhalo mass functions
constructed in the EPS framework to those for observables. We quantify how the
tidal model affects the predictions of the subhalo mass function at $z=0$ in
this section.

For this purpose, we compare the evolved subhalo mass function assuming a
time-independent parameter model of Jiang & van den Bosch (2016, 2017) with
those predicted in the time-dependent parameter model (Hiroshima et al., 2018)
which we have shown in previous sections. Jiang & van den Bosch (2016, 2017)
have calibrated the model parameter $A$ and $\zeta$ with simulation (Klypin et
al., 2011) so that the the average and the scatter of cumulative mass function
in the regime $m/M=[10^{-4},1]$ agrees, obtaining $(A,\zeta)=(0.86,0.07)$. The
result is shown in Fig. 7. The difference between the time-dependent and time-
independent model predictions is subtle for relatively massive subhalos and
both models agree well with the simulation at $m/M>{\cal O}(10^{-4})$.

Figure 7: Evolved mass function with a different model assuming time-
independent parameters. In the range of $m/M\gtrsim{\cal O}(10^{-4})$, the
difference from the result assuming time-dependent parameters is subtle. Both
models show good agreement with numerical results which are shown as square
with error bars.

The difference between two models becomes apparent when we compare in a wider
mass range of $m/M\geq 10^{-10}$, as we show in Fig. 8. In this figure, we
take the ratio of the mean and variance of the time-independent tidal-
parameter model, which is denoted with “const" in the figure, to those of the
time-dependent model. A single mass host of $M_{\rm host}=1.30\times
10^{12}M_{\odot}$ is assumed for both cases. The time-dependent parameter
model predicts a smaller number of subhalos at $m/M\lesssim{\cal O}(10^{-6})$.
The amount of the scatter is slightly larger for the time-independent
parameter model at $m/M\gtrsim{\cal O}(10^{-3})$ while it decreases at
$m/M\lesssim{\cal O}(10^{-4})$. On the other hand, the scatter of the subhalo
number count remains at small mass-ratio regimes of $m/M\lesssim{\cal
O}(10^{-4})$ for time-dependent parameter model.

This reflects the fact that the trend of the model dependence is controlled by
the mass-loss strength. The smaller amplitude of the tidal model parameter in
this regime leads to a smaller scatter. The amplitude of the mass-loss rate is
larger at $m/M\gtrsim 10^{-2}$ for the time-independent parameter model hence
the mean is smaller while the scatter, which is determined by the amount of
the mass-loss, becomes larger. The quantitative agreement between models at
$m/M\gtrsim 10^{-4}$ could be an outcome for the time-independent parameter
model as they have explained the tuning against the simulation in the Sec.3 of
Jiang & van den Bosch (2016).

The effects from the mass-ratio dependence, which is characterized with the
parametar $\zeta$, becomes apparent in the regime of $m/M\lesssim{\cal
O}(10^{-4})$ where it becomes difficult to compare with simulations. The mass-
ratio dependence is weaker for the time-dependent parameter model, then the
tidal power of the less massive hosts becomes larger compared to the case
assuming the time-independent parameter. This means that the relative strength
of the tidal force is larger for the former at high redshifts, i.e. in the
younger hosts. Subhalos those have accreted in such a stage of the host
dominates in the low-mass regime of the mass function. As a result, the
scatter remains at $m/M\lesssim 10^{-4}$ and the prediction of the number
count could be smaller by a few factors at $m/M\sim{\cal O}(10^{-10})$ for the
time-dependent parameter model.

Figure 8: Comparison between the cumulative mass function with two different
tidal evolution models. We take the time-dependent parameter model as the
reference and denoting the time-independent model with the suffix “const" in
this figure. The increasing number of the subhalo number count towards the low
mass regime reflects a stronger mass-ratio dependence of the time-independent
parameter model. The stronger tidal power of the low-mass host in the time-
dependent parameter model generates a larger scatter in the low subhalo mass
end.

In Fig. 9, we quantify the dependence of the Poisson fluctuation on tidal
models. In this figure, we take the ratio between cases with and without
Poisson fluctuations. Again we consider a single mass host of $M_{\rm
host}=1.30\times 10^{12}M_{\odot}$, ignoring the host-mass scatter. The
Poisson fluctuation does not affect the mean for both tidal models, while the
variance at the massive end is apparently changed. The variance is enhanced by
several factors at $m/M\gtrsim{\cal O}(10^{-1})$, where the regime corresponds
to $N_{\rm tot}\lesssim$ a few. The amount of the Poisson fluctuation is not
sensitive to the tidal-evolution models. This indicates that the estimates of
the intrinsic subhalo statistics from observations or simulations with this
relationship is robust against the uncertainty in the tidal evolution models.

Figure 9: The ratio of the mean and the variance between the predictions with-
and without- the Poisson-fluctuation effects. The “base" corresponds to the
case without Poisson fluctuation. In the upper(lower) panel we show the mean
(1$\sigma$ scatter). The mean is not affected in either of the tidal evolution
models. The Poisson fluctuation affects the 1$\sigma$ scatter at
$m/M\gtrsim{\cal O}(10^{-1})$. For both models, the effect is most apparent at
the massive end of the subhalo mass as predicted.

##  6 Conclusion

In this article, we propose a new scheme to construct subhalo mass functions
in an analytical way, which bases on the well-known EPS theory. Subhalo mass
functions at arbitrary mass scales and redshifts are accessible with their
statistical properties in this model. No additional parameters to be
calibrated with numerical simulations are required for unvevolved mass
functions while the model reconstructs the numerical results well in the
regime where we can compare with. The model could be a new probe for subhalo
mass functions at the smallest-scales, which are difficult to investigate in
numerical simulations.

In our scheme, we track the evolution of the main branch of the host merger
tree and different evolution histories of the hosts are directly reflected to
their subahlo mass functions. As a demonstration, we investigate the scatter
in the subhalo number counts considering the Milky-Way-like galaxies of
$M_{\rm host}=1.30\times 10^{12}M_{\odot}$ at $z=0$. Our findings are:

  1. 1.

The Poisson fluctuation dominates at $m/M_{\rm host}\gtrsim{\cal O}(10^{-2})$.
The cumulative number of subhalo exceeds one in the range of $m/M\sim[3\times
10^{-2},\ 5\times 10^{-2}]$ if the Poisson fluctuation effect is removed,
while it extends to a much wider range of $m/M\sim[10^{-2},\ 10^{-1}]$. The
amount of the enhancement in the scatter due to this effect is about a factor
of a few at its peak.

  2. 2.

The host mass variation in a finite range of several factors does not affect
the cumulative subhalo mass function when it is described in terms of the mass
ratio $m/M$. However, the determination of the center value of the host mass
itself does affect the prediction of the halo mass function. The subhalo mass
function scales as the host mass difference.

  3. 3.

Comparing two models for the tidal evolution, the difference between the tidal
models becomes apparent in the regime $m/M\lesssim{\cal O}(10^{-4})$ where the
tuning of the tidal parameter with numerical simulations becomes difficult.
The expected number count of the subhalo of $m/M\gtrsim{\cal O}(10^{-10})$ in
the Milky-Way-like galaxies differs by a factor of $\sim 2$ between the models
we have tested in this paper. The scatter in the small mass-ratio regime
becomes larger when the time-dependent tidal parameter model is adopted. This
feature originates from a relatively larger tidal power of the model for the
hosts in their low-mass stages.

The above-mentioned point 3 could cause an impact on the small-scale halo
search in this era of the precise astrophysics. The difference between two
models becomes clearer at $m/M\lesssim 10^{-6}$, where the regime
gravitational lensing observations are sensitive to (Gilman et al., 2019;
Hsueh et al., 2020). Also, the smallest structure of DM halo is to be probed
in future pulsar timing array experiment (Ishiyama et al., 2010; Clark et al.,
2016; Kashiyama & Oguri, 2018; Lee et al., 2021b; Delos & Linden, 2021) hence
our result could affect those predictions. Another impact should appear in
constraining the small-scale power spectrum based on the DM annihilation
signals in ultra-compact minihalos (Bringmann et al., 2012; Yang et al., 2013;
Nakama et al., 2018; Delos et al., 2018; Lee et al., 2021a). The dependence on
the models of the tidal evolution shifts the such constraints by factors. The
precise models for tidal evolution are needed to obtain much robust picture of
the small-scale structure of our Universe.

## acknowledgement

This work is supported by JSPS KAKENHI Grants No. JP19K23446 (N. H.), MEXT
KAKENHI Grants No. JP20H05852 (N. H.), JP20H05850 (S. A.), JP20H05861 (S. A.),
JP20H05245 (T. I.), and JP21H01122 (T. I.). T. I. has been supported by MEXT
as “Program for Promoting Researches on the Supercomputer Fugaku”
(JPMXP1020200109), JICFuS, and IAAR Research Support Program, Chiba
University, Japan.

We thank Instituto de Astrofísica de Andalucía (IAA-CSIC), Centro de
Supercomputación de Galicia (CESGA) and the Spanish academic and research
network (RedIRIS) in Spain for hosting Uchuu Data Release one in the Skies &
Universes site for cosmological simulations. The Uchuu simulations were
carried out on Aterui II supercomputer at Center for Computational
Astrophysics, CfCA, of National Astronomical Observatory of Japan, and the K
computer at the RIKEN Advanced Institute for Computational Science. The
numerical analysis were partially carried out on XC40 at the Yukawa Institute
Computer Facility in Kyoto University. The Uchuu DR1 effort has made use of
the skun@IAA_RedIRIS and skun6@IAA computer facilities managed by the IAA-CSIC
in Spain (MICINN EU-Feder grant EQC2018-004366-P).

## Data Availability

The calculation code underlying this article are available in
https://github.com/shinichiroando/sashimi-c at
https://doi.org/10.3390/galaxies7030068 and
https://doi.org/10.1103/PhysRevD.97.123002. The data are available in the
article and in http://skiesanduniverses.org/ at
https://doi.org/10.1093/mnras/stab1755.

## References

  * Ade et al. (2016) Ade P. A. R., et al., 2016, Astron. Astrophys., 594, A13
  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astron. Astrophys., 641, A6
  * Ando et al. (2019) Ando S., Ishiyama T., Hiroshima N., 2019, Galaxies, 7, 68
  * Behroozi et al. (2013a) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, Astrophys. J., 762, 109
  * Behroozi et al. (2013b) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013b, Astrophys. J., 763, 18
  * Bertone & Tait (2018) Bertone G., Tait M. P. T., 2018, Nature, 562, 51
  * Bertone et al. (2005) Bertone G., Hooper D., Silk J., 2005, Phys. Rept., 405, 279
  * Bond et al. (1980) Bond J. R., Efstathiou G., Silk J., 1980, Phys. Rev. Lett., 45, 1980
  * Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, Astrophys. J., 379, 440
  * Bosch (2002) Bosch F. C. v. d., 2002, Mon. Not. Roy. Astron. Soc., 331, 98
  * Boylan-Kolchin et al. (2010) Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., 2010, Mon. Not. Roy. Astron. Soc., 406, 896
  * Bringmann (2009) Bringmann T., 2009, New J. Phys., 11, 105027
  * Bringmann et al. (2012) Bringmann T., Scott P., Akrami Y., 2012, Phys. Rev. D, 85, 125027
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, Astrophys. J., 495, 80
  * Buckley et al. (2014) Buckley M. R., Zavala J., Cyr-Racine F.-Y., Sigurdson K., Vogelsberger M., 2014, Phys. Rev. D, 90, 043524
  * Carr et al. (2021) Carr B., Kohri K., Sendouda Y., Yokoyama J., 2021, Rept. Prog. Phys., 84, 116902
  * Clark et al. (2016) Clark H. A., Lewis G. F., Scott P., 2016, Mon. Not. Roy. Astron. Soc., 456, 1394
  * Clowe et al. (2004) Clowe D., Gonzalez A., Markevitch M., 2004, Astrophys. J., 604, 596
  * Cole et al. (2008) Cole S., Helly J., Frenk C. S., Parkinson H., 2008, Mon. Not. Roy. Astron. Soc., 383, 546
  * Correa et al. (2015a) Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015a, Mon. Not. Roy. Astron. Soc., 450, 1514
  * Correa et al. (2015b) Correa C. A., Wyithe J. B., Schaye J., Duffy A. R., 2015b, Mon. Not. Roy. Astron. Soc., 452, 1217
  * Delos & Linden (2021) Delos M. S., Linden T., 2021
  * Delos et al. (2018) Delos M. S., Erickcek A. L., Bailey A. P., Alvarez M. A., 2018, Phys. Rev. D, 98, 063527
  * Diamanti et al. (2015) Diamanti R., Catalan M. E. C., Ando S., 2015, Phys. Rev. D, 92, 065029
  * Diemand et al. (2007) Diemand J., Kuhlen M., Madau P., 2007, Astrophys. J., 667, 859
  * Gao et al. (2004) Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, Mon. Not. Roy. Astron. Soc., 355, 819
  * Gilman et al. (2019) Gilman D., Birrer S., Treu T., Nierenberg A., Benson A., 2019, Mon. Not. Roy. Astron. Soc., 487, 5721
  * Giocoli et al. (2008a) Giocoli C., Tormen G., Bosch F. C. v. d., 2008a, Mon. Not. Roy. Astron. Soc., 386, 2135
  * Giocoli et al. (2008b) Giocoli C., Pieri L., Tormen G., 2008b, Mon. Not. Roy. Astron. Soc., 387, 689
  * Giocoli et al. (2012) Giocoli C., Tormen G., Sheth R. K., 2012, Mon. Not. Roy. Astron. Soc., 422, 185
  * Green & van den Bosch (2019) Green S. B., van den Bosch F. C., 2019, Mon. Not. Roy. Astron. Soc., 490, 2091
  * Green et al. (2004) Green A. M., Hofmann S., Schwarz D. J., 2004, Mon. Not. Roy. Astron. Soc., 353, L23
  * Hayashi et al. (2003) Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T. R., 2003, Astrophys. J., 584, 541
  * Hiroshima et al. (2018) Hiroshima N., Ando S., Ishiyama T., 2018, Phys. Rev. D, 97, 123002
  * Hofmann et al. (2001) Hofmann S., Schwarz D. J., Stoecker H., 2001, Phys. Rev. D, 64, 083507
  * Hogan & Dalcanton (2000) Hogan C. J., Dalcanton J. J., 2000, Phys. Rev. D, 62, 063511
  * Hsueh et al. (2020) Hsueh J.-W., Enzi W., Vegetti S., Auger M., Fassnacht C. D., Despali G., Koopmans L. V. E., McKean J. P., 2020, Mon. Not. Roy. Astron. Soc., 492, 3047
  * Ishiyama (2014) Ishiyama T., 2014, Astrophys. J., 788, 27
  * Ishiyama & Ando (2020) Ishiyama T., Ando S., 2020, Mon. Not. Roy. Astron. Soc., 492, 3662
  * Ishiyama et al. (2010) Ishiyama T., Makino J., Ebisuzaki T., 2010, Astrophys. J. Lett., 723, L195
  * Ishiyama et al. (2021) Ishiyama T., et al., 2021, Mon. Not. Roy. Astron. Soc., 506, 4210
  * Jiang & van den Bosch (2014) Jiang F., van den Bosch F. C., 2014, Mon. Not. Roy. Astron. Soc., 440, 193
  * Jiang & van den Bosch (2016) Jiang F., van den Bosch F. C., 2016, Mon. Not. Roy. Astron. Soc., 458, 2848
  * Jiang & van den Bosch (2017) Jiang F., van den Bosch F. C., 2017, Mon. Not. Roy. Astron. Soc., 472, 657
  * Karukes et al. (2020) Karukes E. V., Benito M., Iocco F., Trotta R., Geringer-Sameth A., 2020, JCAP, 05, 033
  * Kashiyama & Oguri (2018) Kashiyama K., Oguri M., 2018, doi:10.48550/arXiv.1801.07847
  * Klypin et al. (2011) Klypin A., Trujillo-Gomez S., Primack J., 2011, Astrophys. J., 740, 102
  * Lacey & Cole (1993) Lacey C. G., Cole S., 1993, Mon. Not. Roy. Astron. Soc., 262, 627
  * Lee et al. (2021a) Lee V. S. H., Mitridate A., Trickle T., Zurek K. M., 2021a, JHEP, 06, 028
  * Lee et al. (2021b) Lee V. S. H., Taylor S. R., Trickle T., Zurek K. M., 2021b, JCAP, 08, 025
  * Li & Mo (2009) Li Y., Mo H., 2009, doi:10.48550/arXiv.0908.0301
  * Loeb & Zaldarriaga (2005) Loeb A., Zaldarriaga M., 2005, Phys. Rev. D, 71, 103520
  * Ludlow et al. (2016) Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., Navarro J. F., Cole S., Frenk C. S., 2016, Mon. Not. Roy. Astron. Soc., 460, 1214
  * Moline et al. (2017) Moline A., Sanchez-Conde M. A., Palomares-Ruiz S., Prada F., 2017, Mon. Not. Roy. Astron. Soc., 466, 4974
  * Moline et al. (2021) Moline A., et al., 2021, doi:10.48550/arXiv.2110.02097
  * Moreno & Sheth (2008) Moreno J., Sheth R. K., 2008, Mon. Not. Roy. Astron. Soc., 391, 1729
  * Nadler et al. (2020) Nadler E. O., Banerjee A., Adhikari S., Mao Y.-Y., Wechsler R. H., 2020, Astrophys. J., 896, 112
  * Nakama et al. (2018) Nakama T., Suyama T., Kohri K., Hiroshima N., 2018, Phys. Rev. D, 97, 023539
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, Astrophys. J., 462, 563
  * Parkinson et al. (2008) Parkinson H., Cole S., Helly J., 2008, Mon. Not. Roy. Astron. Soc., 383, 557
  * Peebles (1982) Peebles P. J. E., 1982, Astrophys. J. Lett., 263, L1
  * Penarrubia & Benson (2005) Penarrubia J., Benson A. J., 2005, Mon. Not. Roy. Astron. Soc., 364, 977
  * Penarrubia et al. (2010) Penarrubia J., Benson A. J., Walker M. G., Gilmore G., McConnachie A., Mayer L., 2010, Mon. Not. Roy. Astron. Soc., 406, 1290
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, Astrophys. J., 187, 425
  * Profumo et al. (2006) Profumo S., Sigurdson K., Kamionkowski M., 2006, Phys. Rev. Lett., 97, 031301
  * SÃ¡nchez-Conde & Prada (2014) SÃ¡nchez-Conde M. A., Prada F., 2014, Mon. Not. Roy. Astron. Soc., 442, 2271
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, Astrophys. J., 568, 52
  * Yang et al. (2011) Yang X., Mo H., Zhang Y., Bosch F. C. d., 2011, Astrophys. J., 741, 13
  * Yang et al. (2013) Yang Y., Yang G., Huang X., Chen X., Lu T., Zong H., 2013, Phys. Rev. D, 87, 083519
  * Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Boerner G., 2009, Astrophys. J., 707, 354
  * van Albada et al. (1985) van Albada T. S., Bahcall J. N., Begeman K., Sancisi R., 1985, Astrophys. J., 295, 305
  * van den Bosch & Ogiya (2018) van den Bosch F. C., Ogiya G., 2018, Mon. Not. Roy. Astron. Soc., 475, 4066
  * van den Bosch et al. (2005) van den Bosch F. C., Tormen G., Giocoli C., 2005, Mon. Not. Roy. Astron. Soc., 359, 1029