# The mass accretion history of dark matter haloes down to Earth mass

Yizhou Liu*  |1,2, Liang Gao1,2,3,4, Sownak Bose*  |4, Carlos S. Frenk4,
Adrian Jenkins*  |4, Volker Springel*  |5, Jie Wang1,2,3, Simon D. M. White*
|5, Haonan Zheng*  |1,2,4
1National Astronomical Observatories, Chinese Academy of Sciences, Beijing,
100101, China
2School of Astronomy and Space Science, University of Chinese Academy of
Sciences, Beijing 100049, China
3Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal
University, Beijing 102206, China
4Institute for Computational Cosmology, Department of Physics, Durham
University, Science Laboratories, South Road, Durham DH1 3LE
5Max-Planck Institute for Astrophysics, Karl-Schwarzschild Str. 1, D-85748,
Garching, Germany
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We take advantage of the unprecedented dynamical range provided by the
"Cosmic-Zoom" project to study the mass accretion history (MAH) of present-day
dark matter haloes over the entire mass range present in the $\Lambda$CDM
paradigm when the dark matter is made of weakly interacting massive particles
of mass $100\ \mathrm{GeV}$. In particular, we complement previous studies by
exploring the MAHs of haloes with mass from $10^{8}\ h^{-1}\mathrm{M_{\odot}}$
down to Earth mass, $10^{-6}\ h^{-1}\mathrm{M_{\odot}}$. The formation
redshift of low-mass haloes anti-correlates weakly with mass, peaking at $z=3$
for haloes of mass $10^{-4}\ h^{-1}\mathrm{M_{\odot}}$. Even lower masses are
affected by the free-streaming cutoff in the primordial spectrum of density
fluctuations and form at lower redshift. We compare MAHs in our simulations
with predictions from two analytical models based on the extended Press-
Schechter theory (EPS), and three empirical models derived by fitting and
extrapolating either results from cosmological $N$-body simulations or Monte
Carlo realizations of halo growth. All models fit our simulations reasonably
well over the mass range for which they were calibrated. While the empirical
models match better for more massive haloes, $M>10^{10}\
h^{-1}\mathrm{M_{\odot}}$, the analytical models do better when extrapolated
down to Earth mass. At the higher masses, we explore the correlation between
local environment density and MAH, finding that biases are relatively weak,
with typical MAHs for haloes in extremely low-density and in typical regions
differing by less than $20$ percent at high redshift. We conclude that EPS
theory predicts the hierarchical build-up of dark matter haloes quite well
over the entire mass range.

###### keywords:

methods: analytical – methods: numerical – galaxies: haloes

††pubyear: 2022††pagerange: The mass accretion history of dark matter haloes
down to Earth mass–A

##  1 Introduction

In the standard cold dark matter theory, dark matter haloes are the building
blocks of the large-scale structure of the Universe. They provide the
framework within which gas cools and condenses to form galaxies (White & Rees,
1978). Understanding the growth of dark matter haloes is therefore of
fundamental importance to all aspects of physical cosmology (Mo et al., 2010).

Extensive efforts have been made to build analytical models for the mass
accretion histories (MAHs) of dark matter haloes. Most previous models are
empirical, derived from merging trees of dark matter haloes, which can be
extracted either from cosmological simulations (Kauffmann & White, 1993;
Wechsler et al., 2002; Neistein & Dekel, 2008; Zhao et al., 2009; McBride et
al., 2009; Genel et al., 2010; Fakhouri et al., 2010; Giocoli et al., 2012;
Ludlow et al., 2013) or from random Monte Carlo realizations of a statistical
model (Kauffmann & White, 1993; Somerville & Kolatt, 1999; van den Bosch,
2002). These models are very successful at predicting MAHs for dark matter
haloes in certain mass and redshift ranges. For example, one of the latest
models based on this approach (van den Bosch et al., 2014) matches well the
Bolshoi simulation (Klypin et al., 2011) in the halo mass range
$[10^{11},10^{14}]\ h^{-1}\mathrm{M_{\odot}}$ up to redshift $z=3$.

An alternative analytic approach to building a MAH model directly utilizes the
extended Press-Schechter theory (EPS: Bond et al., 1991; Bower, 1991; Lacey &
Cole, 1993); Gao et al. (2005b, hereafter G05) was perhaps the first such
model. Neistein et al. (2006) also developed an analytic model based on EPS
and used a number of approximations to find a relatively simple solution. As
shown by Correa et al. (2015a), with some further assumptions and
approximations, this model agrees well with dark-matter-only simulations from
the OWLS project (Schaye et al., 2010) for $10^{11}\ \mathrm{M_{\odot}}$
haloes at redshift $z=1,2,3,4$ (Correa et al., 2015b).

In this paper, we take advantage of the unprecedented dynamical range of the
Cosmic-Zoom simulation suite (Wang et al., 2020), to explore the MAHs of dark
matter haloes over the entire halo mass range populated at $z=0$ if dark
matter is made of 100GeV WIMPs (see the review of Roszkowski et al., 2018). We
compare the results with a variety of models that have previously been tested
only over a restricted part of this halo mass range.

This paper is organized as follows: in Section 2, we introduce the simulations
used for this study. In Section 3, we briefly review the EPS formalism based
on the spherical collapse model, together with five previously published MAH
models. A detailed comparison with the numerical results is presented and
discussed in Section 4, while in Section 4.2 we we use the Millennium-II
simulation (Boylan-Kolchin et al., 2009) to explore the impact of local
environment on the MAH. We summarize our main findings in Section 5.

##  2 Simulations

The simulations used in this study comprise a suite of high resolution multi-
zoom cosmological simulations from the “Cosmic-Zoom” project and a high
resolution, large volume cosmological simulation – the Millennium-II (MS-II)
(Boylan-Kolchin et al., 2009), which we now briefly describe.

###  2.1 The Cosmic-Zoom simulation

Table 1: Basic parameters of the “Cosmic-Zoom” simulations. From left to right the columns give: the name of the zoom level; the radii of the high-resolution regions; the number of high-resolution particles; the softening length of the high-resolution regions; the mass of the high-resolution particles; the rms linear overdensity extrapolated to $z=0$ in the top-hat filter of $M_{\rm tot}=n_{p}m_{p}$; the ratio between the mean mass density in the high-resolution regions at $z=0$ and the cosmic mean. level | $R_{high}[h^{-1}\mathrm{Mpc}]$ | $n_{p}$ | $\epsilon[h^{-1}\mathrm{kpc}]$ | $m_{p}[h^{-1}\mathrm{M_{\odot}}]$ | $\sigma(M_{tot},z=0)$ | $\bar{\rho}/\rho_{mean}$
---|---|---|---|---|---|---
L0 | 500 | $1.0\times 10^{10}$ | 5.0 | $1.0\times 10^{9}$ |  | 1.0
L1 | 35 | $1.0\times 10^{10}$ | $3.0\times 10^{-1}$ | $5.0\times 10^{5}$ | 0.34 | 0.39
L2 | 6.0 | $5.4\times 10^{9}$ | $3.8\times 10^{-2}$ | $1.0\times 10^{3}$ | 1.66 | 0.082
L3 | 0.68 | $1.8\times 10^{9}$ | $5.6\times 10^{-3}$ | $1.9$ | 4.22 | 0.036
L4 | 0.18 | $2.0\times 10^{9}$ | $6.8\times 10^{-4}$ | $3.7\times 10^{-3}$ | 6.96 | 0.026
L5 | 0.024 | $1.5\times 10^{9}$ | $1.5\times 10^{-4}$ | $3.9\times 10^{-5}$ | 9.36 | 0.024
L6 | 0.0045 | $1.7\times 10^{9}$ | $2.6\times 10^{-5}$ | $1.8\times 10^{-7}$ | 12.12 | 0.014
L7c | 0.00075 | $2.5\times 10^{9}$ | $3.6\times 10^{-6}$ | $5.8\times 10^{-10}$ | 15.06 | 0.016
L8c | 0.00016 | $1.5\times 10^{9}$ | $9.5\times 10^{-7}$ | $1.1\times 10^{-11}$ | 17.60 | 0.028

Using a multi-zoom technique, the “Cosmic-Zoom” project simulated the
formation and evolution of dark matter haloes over the full mass range
assuming the dark matter is a weakly interacting massive particle (WIMP) of
mass about $100\ \mathrm{GeV}$ (Wang et al., 2020). The Cosmic-Zoom simulation
suite contains 1 base simulation L0, and 8 individual “zoom-in” resimulations
L1 to L8c carried out with the GADGET-4 simulation code (Springel et al.,
2021). The strategy adopted to make these computations feasible consists of
successively re-simulating low density zoomed-in regions of the parent
simulation until the lowest mass halo is resolved with about 1 million
particles. The base simulation L0 is a Planck Cosmological version of the
Millennium simulation (Springel et al., 2005). These simulations adopt the
Planck cosmological parameters: $\Omega_{\mathrm{m}}=0.307$,
$\Omega_{\mathrm{b}}=0.04825$, $H_{0}=67.77\ \mathrm{km\ s^{-1}Mpc^{-1}}$,
$\sigma_{8}=0.8288$ and $n_{\mathrm{s}}=0.9611$ (Planck Collaboration et al.,
2014). We list the basic parameters of the Cosmic-Zoom simulations in Table 1
and refer the readers to Wang et al. (2020) for details.

###  2.2 The Millennium-II simulation

The Cosmic-Zoom simulations are extreme in that they follow structure
formation and evolution in low density regions. Gao & White (2007)
demonstrated that statistically halo properties not only depend on the halo
mass but also on its large scale environment, the “assembly bias” effect.
Thus, results from the Cosmic-Zoom simulations may be biased relative to
results for regions of average density. We use the Millennium-II simulation
(MS-II) to examine the impact of the environment on the MAHs. The cosmological
parameters of the MS-II are consistent with the WMAP-1 model:
$\Omega_{\mathrm{m}}=0.25$, $\Omega_{\mathrm{b}}=0.045$, $H_{0}=73.0\
\mathrm{km\ s^{-1}Mpc^{-1}}$, $\sigma_{8}=0.90$ and $n_{\mathrm{s}}=1.0$. The
cube size of the Millennium-II is $L_{\mathrm{box}}=100\ h^{-1}\mathrm{Mpc}$
and the total particle number is $N_{\mathrm{p}}=1024^{3}$. The simulation has
a mass resolution of $6.89\times 10^{6}\ h^{-1}\mathrm{M_{\odot}}$ and a
spatial resolution of $1\ h^{-1}\mathrm{kpc}$, respectively. While the MS-II
adopts different cosmological parameters from the Cosmic-Zoom simulations,
this will not qualitatively change our main results.

The dark matter haloes in these simulations are identified by means of the
friends-of-friends (FoF) algorithm, linking together particles within $0.2$
times the mean separation (Davis et al., 1985). The mass of a dark matter halo
is defined as its virial mass, $M_{200}$, namely the mass within the radius
$r_{200}$ of the sphere enclosing a mean density of 200 times the critical
density (see halo mass conversion in Diemer, 2018). Dark matter halo merging
trees are constructed with the algorithm described in Boylan-Kolchin et al.
(2009).

##  3 MAH models

In this section we briefly introduce the EPS theory and review existing models
selected for comparison with our simulations, among which two are analytical,
based on EPS theory (Gao et al., 2005b; Correa et al., 2015a), and three are
empirical, based on N-body or Monte Carlo simulations (van den Bosch, 2002;
Zhao et al., 2009; Ludlow et al., 2013; van den Bosch et al., 2014).

###  3.1 Analytical models

  * i)

The model of Gao et al. (2005)

Perhaps the first analytical MAH model was proposed by Gao et al. (2005b).
This is based upon the conditional mass function from EPS theory,

$\displaystyle n(M_{1},z_{1}|M_{0},z_{0}){\rm d}M_{1}=$ (1)
$\displaystyle\frac{M_{0}}{M_{1}}f(S(M_{1}),\delta_{c}(z_{1})|S(M_{0}),\delta_{c}(z_{0}))\Bigg{|}\frac{\mathrm{d}S(M_{1})}{\mathrm{d}M_{1}}\Bigg{|}\mathrm{d}M_{1},$

where $S(M)$ is the mass variance of the present density field smoothed with a
top-hat filter of radius $R=[M/(\frac{4}{3}\pi\bar{\rho})]^{1/3}$, and
$\delta_{c}(z)$ is the critical density for spherical collapse, linearly
extrapolated to the present time. Eqn. (1) describes the average number of
progenitors at redshift $z_{1}$, in the mass interval
$(M_{1},M_{1}+\mathrm{d}M_{1})$, which have merged into a descendent halo of
mass, $M_{0}$, at time $z_{0}$. The most massive progenitor (hereafter MMP) of
a halo at redshift $z$ is derived by requiring the average number of
progenitors above a given mass $M^{\prime}$ to be exactly unity:

$\int^{M_{0}}_{M^{\prime}}n(M,z|M_{0},0)\mathrm{d}M=1,$ (2)

We take the mean progenitor mass as our estimate of the MMP; this can be
derived by solving the following equation,

$\displaystyle M(z)$
$\displaystyle=\frac{\int^{M_{0}}_{M^{\prime}}Mn(M,z|M_{0},0)\mathrm{d}M}{\int^{M_{0}}_{M^{\prime}}n(M,z|M_{0},0)\mathrm{d}M}$
(3)
$\displaystyle=M_{0}\left[1-\mathrm{erf}\left(\frac{\delta_{c}(z)-\delta_{c}(0)}{\sqrt{2(S(M^{\prime})-S(M_{0}))}}\right)\right].$

While this model is simple, it quite accurately predicts the MAHs of simulated
rich galaxy clusters out to redshift 40, when such clusters were only
$1/10^{9}$ of the mass at the present day (Gao et al., 2005b). The model was
later applied to predict MAHs of a suite of zoom-in high-resolution
simulations of individual Milky Way-size dark matter haloes from the Aquarius
simulations (Gao et al., 2010).

  * ii)

The model of Correa et al. (2015)

Using a similar approach to that of G05, Neistein et al. (2006) developed an
analytical MAH model based on the EPS theory. In contrast with G05, this model
describes the growth of a dark matter halo by following the main branch of its
merger tree. In this model, the mass growth of the main progenitor is
described by the differential equation,

$\frac{\mathrm{d}M(z)}{\mathrm{d}\delta_{c}(z)}=-\sqrt{\frac{2}{\pi}}\frac{M(z)}{\sqrt{S(M(z)/q)-S(M(z))}},$
(4)

where $M/q$ is equivalent to the lowest limit, $M^{\prime}$, of the integral.
Neistein et al. (2006) further found that $q$ is within the range $[2.1,2.3]$,
with a weak dependence on $\Omega_{\mathrm{m}}$ and halo mass. Using the fit
to $S(M)$ given by van den Bosch (2002) and assuming $q=2.2$, Neistein et al.
(2006) obtained a fitting formula for the MAHs of dark matter haloes in the
mass range, $[10^{6},10^{15}]\ h^{-1}\mathrm{M_{\odot}}$, and found good
agreement between this model and the Monte Carlo result of Somerville & Kolatt
(1999).

McBride et al. (2009) found that the MAHs of dark matter haloes can be well
fit by

$M(z)=M_{0}(1+z)^{\alpha}\mathrm{e}^{\beta z},$ (5)

where $\alpha$, $\beta$ are free parameters that are correlated with each
other (Wong & Taylor, 2012). Correa et al. (2015a) substituted this solution
in Eqn. (4) and found that, with some approximations, $\alpha$ and $\beta$ can
be analytically determined. We refer to this model as C15 and point the reader
to the original paper for further details of this model.

###  3.2 $N$-body or Monte Carlo simulations based fitting models

The MAHs of dark matter haloes can also be derived from merger trees extracted
either from N-body simulations or from Monte Carlo realizations. We consider
three popular models.

  * i)

van den Bosch 2002 & Ludlow et al. 2014

van den Bosch (2002) used the N-branch Monte Carlo algorithm (Somerville &
Kolatt, 1999) to construct a number of random MAH realizations of dark matter
halos in a wide range of masses. The average MAH of dark matter haloes was
fitted with the formula,

$\log\frac{M(z)}{M_{0}}=-0.301\left[\frac{\log(1+z)}{\log(1+z_{f})}\right]^{\nu},$
(6)

where $z_{f}$ and $\nu$ are free, strongly correlated parameters. The value of
$z_{f}$ is obtained directly by solving the following equation,

$\delta_{c}(z_{f})=\delta_{c}(0)+0.477\sqrt{2[\sigma^{2}(fM_{0})-\sigma^{2}(M_{0})]},$
(7)

where $f$ is a free parameter. The strong correlation between $z_{f}$ and
$\nu$ can be approximated as

$\displaystyle\nu$
$\displaystyle=1.211+1.858\log(1+z_{f})+0.308\Omega_{\Lambda}^{2}$ (8)
$\displaystyle-0.032\log[M_{0}/(10^{11}h^{-1}\mathrm{M_{\odot}})].$

van den Bosch (2002) found that $f=0.254$ results in a good fit to the average
MAHs of the Monte Carlo realizations. However, the model does not fit
cosmological simulations well. Ludlow et al. (2014) improved the model by
calibrating it against the Millennium simulation, and found that $f=0.068$
produces a good fit to the average MAHs in the Millennium simulation. We refer
to this improved model as VL14 model.

  * ii)

van den Bosch 2014

Based on the Bolshoi simulation (Klypin et al., 2011) and the binary Monte
Carlo merger tree method of Parkinson et al. (2008), van den Bosch et al.
(2014) developed another model describing the MAH of the main branch of dark
matter halo merger trees.

Some studies (van den Bosch et al., 2005; Giocoli et al., 2008; Li & Mo, 2009;
Yang et al., 2011) indicated that the unevolved subhalo mass functions of dark
matter halos of different masses are almost universal (but see Xie & Gao,
2015), meaning that, when expressed in terms of the normalized mass,
$M_{p}/M_{0}$ (where $M_{p}$ is progenitor mass and $M_{0}$ is final halo
mass), the accretion mass function of dark matter halos is independent of halo
mass, $M_{0}$, but the assembly time of progenitors is different for different
haloes. Motivated by this, van den Bosch et al. (2014) speculated that the MAH
has a universal form when written as a function, $\psi(\widetilde{\omega})$,
of a ’universal time’ $\widetilde{\omega}$, where $\psi=M/M_{0}$. They found
that this universal time can be expressed as

$\widetilde{\omega}(z|M,M_{0},z_{0})=\frac{\delta_{c}(z)-\delta_{c}(z_{0})}{\sqrt{\sigma^{2}(M)-\sigma^{2}(M_{0})}}G^{0.4}(M,M_{0},z_{0}),$
(9)

where $G(M,M_{0},z_{0})$ is a perturbing function proposed by Parkinson et al.
(2008) to improve the accuracy of the progenitor mass functions derived from
Monte Carlo realizations according to the Millennium simulation.

van den Bosch et al. (2014) fitted the inverse function,
$\widetilde{\omega}(\psi)$, of this universal form with the equation

$F(\psi)=a_{1}[1-a_{2}\log\psi]^{a_{3}}(1-\psi^{a_{4}})^{a_{5}}$ (10)

to the Bolshoi simulation and Parkinson merger trees, with the latter used to
complement the simulation data below the resolution limit. The average or
median MAH for a halo of any mass in any cosmology can be computed by solving
this equation numerically,

$F(\psi)=\widetilde{\omega}(\psi,z).$ (11)

We call this model vdB14. For the specific implementation details, we refer
the reader to the original paper.

  * iii)

Zhao et al. 2009

Based on a set of $N$-body simulations in different cosmologies (scale-free,
standard CDM, open CDM and $\Lambda$CDM - see Jenkins et al. 1998 for the
definitions), Zhao et al. (2009) developed an empirical model to describe the
median MAH of the main branch of dark matter halos,

$\mathrm{d}\lg\sigma(M)/\mathrm{d}\lg\delta_{c}(z)=\frac{\omega(z,M)-p(z,z_{\rm
obs},M_{\rm obs})}{5.85},$ (12)

where $\omega(z,M)\equiv\delta_{c}(z)/s(M)$ and $s(M)\equiv\sigma(M)\times
10^{\mathrm{d}\lg\sigma/\mathrm{d}\lg m|_{M}}$. The function $p(z,z_{\rm
obs},M_{\rm obs})$ is given by

$\displaystyle p(z,z_{\rm obs},M_{\rm obs})=$ (13) $\displaystyle p(z_{\rm
obs},z_{\rm obs},M_{\rm
obs})\times\mathrm{max}\left[0,1-\frac{\lg\delta_{c}(z)-\lg\delta_{c}(z_{\rm
obs})}{0.272/\omega(z_{\rm obs},M_{\rm obs})}\right],$

where

$\displaystyle p(z_{\rm obs},z_{\rm obs},M_{\rm obs})=$ (14)
$\displaystyle\frac{1}{1+[\omega(z_{\rm obs},M_{\rm
obs})/4]^{6}}\frac{\omega(z_{\rm obs},M_{\rm obs})}{2},$

and $M_{\rm obs}$ and $z_{\rm obs}$ are the mass and redshift of the final
halo.

Given a cosmological model and its power spectrum, the median MAH of a halo of
mass, $M_{\rm obs}$, at $z_{\rm obs}$ can be constructed step by step using
Eqn. (12). We call this model Z09 and refer the readers to the original paper
for details of the model.

##  4 Results

###  4.1 MAHs of dark matter haloes over the entire mass range: simulation
results vs theoretical predictions

Figure 1: MAHs of dark matter haloes from the “Cosmic-Zoom” simulations and $5$ empirical or analytical models over the mass range $[10^{-6},10^{14}]\ h^{-1}\mathrm{M_{\odot}}$ at $z=0$. In the left and right panels, the black dots indicate the MAHs extracted from the Cosmic-Zoom simulations, while lines of different colours and styles are predictions from G05, Z09, VL14, vdB14, C15, as indicated by the labels. Figure 2: Same as the figure 1, but showing the residuals. Figure 3: Left panel: MAHs of dark matter haloes, scaled to the final halo mass, for different masses. Dots connected by dashed lines, coloured according to halo mass, show results from our simulations; solid lines, also coloured according to halo mass, show predictions from the G05 model. Right panel: halo mass-formation time relation predicted by the G05 model. Here, the formation time is defined as the time when the main progenitor acquired $50\%$ (red), $10\%$ (blue) and $1\%$ (yellow) of its present mass, respectively. The black dots show results for the definition with $50\%$ from the simulation. Table 2: Basic information of the halo samples from the “Cosmic-Zoom” simulations. halo mass ($h^{-1}\mathrm{M_{\odot}}$) | number of samples | level
---|---|---
$10^{14}$ | 2588 | L0
$10^{13}$ | 18981 | L0
$10^{12}$ | 203 | L1
$10^{11}$ | 1288 | L1
$10^{8}$ | 883 | L2
$10^{7}$ | 8980 | L2
$10^{5}$ | 379 | L3
$10^{4}$ | 4137 | L3
$10^{2}$ | 837 | L4
$10^{1}$ | 8780 | L4
$10^{0}$ | 645 | L5
$10^{-1}$ | 6319 | L5
$10^{-3}$ | 2317 | L6
$10^{-4}$ | 24769 | L6
$10^{-5}$ | 706 | L7c
$10^{-6}$ | 2784 | L7c

In Figure 1 we present the MAHs predicted by the five models described above
and derived from the Cosmic-Zoom simulations for dark matter haloes in 16
different mass bins. The masses of these haloes cover almost the entire mass
range, $[10^{-6},10^{14}]\ h^{-1}\mathrm{M_{\odot}}$, for a SUSY dark matter
model. The number of haloes in different mass bins varies from a few hundred
to a few tens of thousands and thus provides good statistics to investigate
the MAHs of haloes of different masses. The exact number of haloes in each
mass bin and the corresponding simulation are listed in Table 2. Since most of
the MAH models introduced in the last Section predict MAHs along the main
branch (MPM) of the merging tree, except G05 which predicts the MAHs of the
MMPs along the entire branch, we present our simulation results for MAHs along
the main branch for easy comparison. Note that statistically, MMP only differs
from the MPM at very high redshift (see Appendix A). In practice, we trace
back each present-day dark matter halo along its main branch step by step
until it first emerges as a halo (i.e. with 20 dark matter particles); the
median values for the whole halo sample in each mass bin are then obtained.

In this figure, we show MAHs from the simulations and models as dots and
lines, respectively. The predictions from the different models are
distinguished with different line types, as indicated in the label. Figure 2
displays the residuals between the models and the simulations, using the same
line styles as in Figure 1. Only the G05 and VL14 models are able to predict
the MAHs over the full mass range considered in this paper. The predictions of
the other models are only available above some mass limits, and as a result
some predictions are missing in some of the panels.

Overall, all models predict our simulation results fairly well given the fact
that the dynamical ranges both in time and mass are huge; the extent of the
agreement amongst the models, however, varies with halo mass and redshift. For
the relatively high mass haloes, $[10^{11},10^{14}]\
h^{-1}\mathrm{M_{\odot}}$, the best-performing model is vdB14, especially at
relatively low redshifts $z<3$. The agreement between the model and the
simulation results is almost perfect for cluster size dark matter haloes,
while the agreement becomes slightly worse for decreasing halo mass and
increasing redshift. This good agreement is not surprising, since the vdB14
model was calibrated with N-body merger trees of haloes in this mass range.
The predictions from the other four models are also acceptable, with
deviations within $25$ percent at most redshifts except the highest. For the
MAHs of low mass haloes, the two analytical models of G05 and C15 perform
extremely well over most of the redshift range we explore here; deviations
from the simulation results are within a few percent. The Z09 model also
performs well at redshifts $z<2$ but tends to underestimate the MAHs beyond
that. The VL14 model is in reasonable agreement with our data with deviations
within $25$ percent in most redshift ranges. For the least massive haloes, the
G05 model still performs well except for Earth mass haloes. Note that G05
actually predicts MMP which gives slightly larger values of the MAHs at high
redshift than the MPM MAHs (see Appendix A). The C15 model also performs well
for the range of halo masses for which results from this model are available.
Apart from the Earth mass haloes, the performance of the VL14 model is also
acceptable.

For easier identification of the relative differences in the MAHs for dark
matter haloes of different masses, we plot together, in the left panel of
Figure 3, the MAHs of the same haloes shown above. Here the solid dots
connected by dashed lines display the Cosmic-Zoom results and the lines,
coloured according to halo mass, show predictions from the G05 model, which we
have shown gives the best results for the MAHs over the entire mass range.
Clearly, more massive haloes increase their masses more rapidly except for the
least massive haloes which are affected by the free-streaming cutoff in the
power spectrum. In the right panel of the figure, we show the formation times
of dark matter haloes over the entire mass range. Here, the formation times
are defined in several ways, namely when the main progenitor of the halo had
accumulated $50\%$, $10\%$ and $1\%$ of its present-day mass. Overall, the
halo mass-formation time relation is monotonic, with more massive haloes
having later formation times, except for the least massive haloes. The
logarithmic slope of the relation varies with halo mass: it is steeper at the
high-mass end and flattens towards the low-mass end. The most massive dark
matter haloes, of mass $10^{15}\ h^{-1}\mathrm{M_{\odot}}$, acquire half of
their present-day mass at $z\sim 0.5$, while those of mass $10^{-4}\
h^{-1}\mathrm{M_{\odot}}$ assemble the same mass fraction at $z\sim 3$.

###  4.2 Impact of local density on MAHs of dark matter haloes

Figure 5: The correlation between halo formation time and large-scale environment density in the MS-II. The left panel shows the mean (grey), the median (black), the $20\%$ (blue) and the $80\%$ (red) points of the environment density distribution as a function of formation time, while the right panel shows the corresponding points of the formation time distribution as a function of environment density. Table 3: Properties of the halo samples from the MS-II simulation. halo mass $h^{-1}\mathrm{M_{\odot}}$ |  tophat smoothing scale $h^{-1}\mathrm{Mpc}$ |  $a\%$ lowest/highest density regions |  overdensity range for lowest samples |  overdensity range for highest samples |  number of haloes
---|---|---|---|---|---
$10^{10}$ |  2 |  $2\%$ |  $-0.916\sim-0.798$ |  $41.152\sim 302.291$ |  1251
$10\%$ |  $-0.916\sim-0.645$ |  $11.702\sim 302.291$ |  6256
$30\%$ |  $-0.916\sim-0.250$ |  $2.307\sim 302.291$ |  18769
5 |  $2\%$ |  $-0.923\sim-0.781$ |  $12.498\sim 27.788$ |  1251
$10\%$ |  $-0.923\sim-0.613$ |  $4.501\sim 27.788$ |  6256
$30\%$ |  $-0.923\sim-0.277$ |  $1.236\sim 27.788$ |  18769
10 |  $2\%$ |  $-0.883\sim-0.692$ |  $3.907\sim 5.689$ |  1251
$10\%$ |  $-0.883\sim-0.523$ |  $2.002\sim 5.688$ |  6256
$30\%$ |  $-0.883\sim-0.218$ |  $0.648\sim 5.688$ |  18769
$10^{11}$ |  2 |  $10\%$ |  $-0.830\sim-0.468$ |  $9.491\sim 175.844$ |  1429
$30\%$ |  $-0.830\sim-0.035$ |  $2.231\sim 175.844$ |  4287
5 |  $10\%$ |  $-0.878\sim-0.551$ |  $4.153\sim 27.758$ |  1429
$30\%$ |  $-0.878\sim-0.229$ |  $1.228\sim 27.758$ |  4287
10 |  $10\%$ |  $-0.858\sim-0.499$ |  $1.962\sim 5.684$ |  1429
$30\%$ |  $-0.858\sim-0.198$ |  $0.645\sim 5.684$ |  4287

By construction, the Cosmic-Zoom simulations follow the evolution of cosmic
structure in extremely low density regions, with most of them having densities
of less than 3 percent of the cosmic mean. We list the overdensities of the
L1-L8c high-resolution regions in Table 1. While the Markov process of the
overdensity trajectories on which EPS theory is based implies that the mass
growth of a halo is unaffected by its large-scale environment, Gao & White
(2007) showed that dark matter haloes “know” about their large-scale
environment, a result confirmed by many other studies (Sheth & Tormen, 2004;
Gao et al., 2005a; Harker et al., 2006; Wechsler et al., 2006; Li et al.,
2008). It is therefore necessary to explore the impact of the local density on
the MAHs.

We make use of the MS-II to examine the impact of local density on the MAHs of
dark matter haloes. The MS-II only resolves dark matter haloes more massive
than $1.38\times 10^{8}\ h^{-1}\mathrm{M_{\odot}}$ with $20$ particles. In
order to have reliable results, we only consider dark matter haloes more
massive than $10^{10}\,h^{-1}\mathrm{M_{\odot}}$, corresponding to more than
about $1400$ dark matter particles. In order to have good statistics, we
select 2 halo samples in two different narrow mass bins, $10^{10}$ and
$10^{11}\,h^{-1}\mathrm{M_{\odot}}$, respectively. For each dark matter halo
in the sample, we calculate its local top-hat environment density on 3
different scales, $2$, $5$ and $10\,h^{-1}\mathrm{Mpc}$, respectively. In
Table 3, we list the basic properties of our halo samples.

In Figure LABEL:fig:Impact_of_Environment, we present MAHs of dark matter
haloes of different top-hat overdensities defined at the three scales
mentioned above for our two halo samples, and show the residuals relative to
the median MAHs of each sample as a whole. Since the lowest halo mass sample
has very good statistics, for this sample we also show the $2$ percent highest
and lowest top-hat overdensities at $2\,h^{-1}\mathrm{Mpc}$. Panel (a) shows
results for our least massive halo sample. The local density does affect the
MAHs. For top-hat densities calculated within $2\,h^{-1}\mathrm{Mpc}$, the
MAHs of overdense halo samples correlate strongly with their top-hat
densities. At redshift $z=3$, the $2\%$, $10\%$ and $30\%$ highest
overdensities deviate from the median of the whole sample by about $60\%$,
$40\%$ and $20\%$, respectively. Interestingly, while the MAHs of the
underdense halo samples are qualitatively the inverse of the overdense samples
as expected, the largest deviations from the median are only about $20\%$,
much smaller than their overdense counterparts. In addition, the MAHs of
underdense haloes seem to correlate with their top-hat densities quite weakly.
The results for other density definitions are shown in the middle and right
panels, which are very similar to the left panel except of the impact becomes
weaker with increasing smoothing scale for the overdense samples. Panel (b)
shows that the results for the $10^{11}\,h^{-1}\mathrm{M_{\odot}}$ sample are
very similar to the results for haloes with mass an order of magnitude lower,
suggesting that the impact of the local density on the MAHs does not correlate
strongly with halo mass.

In order to confirm that our results are consistent with previous papers (Gao
et al., 2005a; Boylan-Kolchin et al., 2009), we investigate the correlation
between halo formation time and large-scale environment in the MS-II. We
consider halo samples selected in the mass bin $\sim
10^{10}\mathrm{M_{\odot}}$ for good statistics. The formation time of a halo
is defined as the epoch when its main progenitor first reaches $50\%$ of its
present-day mass, as noted above in Section 4.1. Here, however, the large-
scale environment density is estimated within a thick spherical shell with
inner and outer radii of 5 and 10 $h^{-1}\mathrm{Mpc}$, to provide a
definition as close as possible to that used in earlier measurements of
assembly bias. The correlation between formation redshift and this environment
density is presented in Fig. 5, where the halo samples are divided into ten
equally sized sub-samples in either formation time (left panel) or environment
density (right panel) in order to display the mean, median, and $20\%$ and
$80\%$ points of the distribution of each quantity as a function of the other.

Fig. 5 shows that there is substantial scatter in the relationship between
formation time and environment density, but that there is nevertheless a clear
correlation between the two. In the left panel, the ratio of the mean
overdensity around the 20% earliest forming haloes to that around the 20%
latest forming haloes is consistent with the ratio $\sim 3$ of the assembly
bias of these two samples as shown in Figure 12 of Boylan-Kolchin et al.
(2009). In both panels, the scatter in the dependent variable is roughly
lognormal, is independent of the value of the independent variable and is
considerably larger than the trend in the median of the distribution.

If the results above can be extrapolated to lower halo masses, we tentatively
conclude that the impact of local density on the results presented in the
previous Section is, at most, mild. The Cosmic-Zoom simulations only
underestimate the MAHs of halos in regions of typical density at the level of
$20$ percent at most, and hence the conclusions drawn from the last Section
are unaffected.

##  5 Conclusion and Discussion

In this paper, we take advantage of the dynamical range of the Cosmic-Zoom
simulations, to study halo MAHs across the full mass range populated if the
dark matter is made of 100GeV WIMPs. In particular, we present MAHs for $z=0$
dark matter haloes in the previously unexplored mass range $[10^{-6},10^{9}]\
h^{-1}\mathrm{M_{\odot}}$. As expected from studies of more massive haloes,
the formation redshifts of low-mass haloes anti-correlate with mass except at
the extremely low masses affected by the free streaming cutoff in the initial
power spectrum. In this regime, the least massive haloes assemble latest. The
formation redshift–mass relation for dark matter haloes thus has a peak at
$z\sim 3$ for halo mass $\sim 10^{-4}\ h^{-1}\mathrm{M_{\odot}}$.

We compare the MAHs of our simulated dark matter haloes with five
representative analytical models, two of which were derived directly from EPS
theory, while the rest were based on merging trees extracted either from Monte
Carlo realizations or from cosmological N-body simulations. Overall, all these
models are in reasonable agreement with our simulations, but the level of
agreement varies with halo mass and redshift. For relatively massive haloes,
$M>10^{10}\ h^{-1}\mathrm{M_{\odot}}$ the empirical models generally agree
better with our simulation data than the analytic models based on EPS theory.
This is not perhaps surprising since the empirical models were calibrated
using N-body simulations which only resolve relatively massive haloes. At
lower mass, the analytic models are superior and match our numerical data
accurately except for masses $\sim 10^{-6}\ h^{-1}\mathrm{M_{\odot}}$, the
least massive sample we considered.

We further examine the impact of local density on the MAHs of dark matter
haloes using a large-volume N-body cosmological simulation – the Millennium-
II. As expected from the assembly bias effect (Gao et al., 2005a; Gao & White,
2007), we find that local environment density does affect the MAHs of dark
matter haloes: for haloes of given mass halo, those residing in high-density
regions assemble systematically earlier than their counterparts residing in
low-density regions. However, the impact is only strong for high density
regions and when the local density is defined within small radii. For low-
density regions, differences of $<20\%$ are found between the mean MAH of
haloes in regions in the lowest $10\%$ of the overdensity distribution and
those of haloes in regions of typical overdensity. The impact of environment
overdensity on the MAHs is similar for our $10^{10}$ and
$10^{11}\,h^{-1}\mathrm{M_{\odot}}$ halo samples, suggesting a weak
correlation with halo mass consistent with that already known for the assembly
bias effect (e.g. Boylan-Kolchin et al., 2009). Furthermore, we do not see
stronger density effects for the MAHs of haloes in the lowest $2\%$, rather
than the lowest $10\%$, of the overdensity distribution. MAHs depend at most
weakly on environment density in the extreme low-density tail. Taken together,
we conclude that the impact of local density on MAHs is mild, except possibly
in substantially overdense regions. Hence, a combination of existing
analytical and empirical MAH models can provide reliable predictions over the
entire halo mass range that is populated at $z=0$ if the dark matter is made
up of WIMPs. Depending on the mass range of interest, one can choose the best
model according to the results presented in this paper.

## Acknowledgements

We acknowledge the support from the National Natural Science Foundation of
China (Grant No. 11988101) and K.C.Wong Education Foundation. CSF acknowledges
support from the European Research Council (ERC) through Advanced Investigator
grant DMIDAS (GA 786910) and from the STFC Consolidated Grant ST/T000244/1. HZ
acknowledges support from the China Scholarships Council (No. 202104910325).
SB is supported by the UK Research and Innovation (UKRI) Future Leaders
Fellowship [grant number MR/V023381/1]. This work used the DiRAC@Durham
facility managed by the Institute for Computational Cosmology on behalf of the
STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS
capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1,
ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant
ST/R000832/1. DiRAC is part of the UK National e-Infrastructure.

## Data Availability

The data used in this study can be made available upon reasonable request to
the corresponding author.

## References

  * Bond et al. (1991) Bond J., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  * Bower (1991) Bower R. G., 1991, MNRAS, 248, 332
  * Boylan-Kolchin et al. (2009) Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150
  * Correa et al. (2015a) Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015a, MNRAS, 450, 1514
  * Correa et al. (2015b) Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015b, MNRAS, 452, 1217
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * Diemer (2018) Diemer B., 2018, ApJS, 239, 35
  * Fakhouri et al. (2010) Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267
  * Gao & White (2007) Gao L., White S. D. M., 2007, MNRAS, 377, L5
  * Gao et al. (2005a) Gao L., Springel V., White S. D. M., 2005a, MNRAS, 363, L66
  * Gao et al. (2005b) Gao L., White S. D. M., Jenkins A., Frenk C., Springel V., 2005b, MNRAS, 363, 379
  * Gao et al. (2010) Gao L., Theuns T., Frenk C., Jenkins A., Helly J., Navarro J., Springel V., White S. D. M., 2010, MNRAS, 403, 1283
  * Genel et al. (2010) Genel S., Bouché N., Naab T., Sternberg A., Genzel R., 2010, ApJ, 719, 229
  * Giocoli et al. (2008) Giocoli C., Tormen G., van den Bosch F. C., 2008, MNRAS, 386, 2135
  * Giocoli et al. (2012) Giocoli C., Tormen G., Sheth R. K., 2012, MNRAS, 422, 185
  * Harker et al. (2006) Harker G., Cole S., Helly J., Frenk C., Jenkins A., 2006, MNRAS, 367, 1039
  * Jenkins et al. (1998) Jenkins A., et al., 1998, ApJ, 499, 20
  * Kauffmann & White (1993) Kauffmann G., White S. D. M., 1993, MNRAS, 261, 921
  * Klypin et al. (2011) Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
  * Lacey & Cole (1993) Lacey C., Cole S., 1993, MNRAS, 262, 627
  * Li & Mo (2009) Li Y., Mo H., 2009, arXiv preprint arXiv:0908.0301
  * Li et al. (2008) Li Y., Mo H., Gao L., 2008, MNRAS, 389, 1419
  * Ludlow et al. (2013) Ludlow A. D., et al., 2013, MNRAS, 432, 1103
  * Ludlow et al. (2014) Ludlow A. D., Navarro J. F., Angulo R. E., Boylan-Kolchin M., Springel V., Frenk C., White S. D. M., 2014, MNRAS, 441, 378
  * McBride et al. (2009) McBride J., Fakhouri O., Ma C.-P., 2009, MNRAS, 398, 1858
  * Mo et al. (2010) Mo H., van den Bosch F. C., White S. D. M., 2010, Galaxy formation and evolution. Cambridge University Press
  * Neistein & Dekel (2008) Neistein E., Dekel A., 2008, MNRAS, 383, 615
  * Neistein et al. (2006) Neistein E., van den Bosch F. C., Dekel A., 2006, MNRAS, 372, 933
  * Parkinson et al. (2008) Parkinson H., Cole S., Helly J., 2008, MNRAS, 383, 557
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
  * Roszkowski et al. (2018) Roszkowski L., Sessolo E. M., Trojanowski S., 2018, Reports on Progress in Physics, 81, 066201
  * Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  * Sheth & Tormen (2004) Sheth R. K., Tormen G., 2004, MNRAS, 350, 1385
  * Somerville & Kolatt (1999) Somerville R. S., Kolatt T. S., 1999, MNRAS, 305, 1
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * van den Bosch (2002) van den Bosch F. C., 2002, MNRAS, 331, 98
  * van den Bosch et al. (2005) van den Bosch F. C., Tormen G., Giocoli C., 2005, MNRAS, 359, 1029
  * van den Bosch et al. (2014) van den Bosch F. C., Jiang F., Hearin A., Campbell D., Watson D., Padmanabhan N., 2014, MNRAS, 445, 1713
  * Wang et al. (2020) Wang J., Bose S., Frenk C. S., Gao L., Jenkins A., Springel V., White S. D. M., 2020, Nature, 585, 39
  * Wechsler et al. (2002) Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002, ApJ, 568, 52
  * Wechsler et al. (2006) Wechsler R. H., Zentner A. R., Bullock J. S., Kravtsov A. V., Allgood B., 2006, ApJ, 652, 71
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Wong & Taylor (2012) Wong A. W., Taylor J. E., 2012, ApJ, 757, 102
  * Xie & Gao (2015) Xie L., Gao L., 2015, MNRAS, 454, 1697
  * Yang et al. (2011) Yang X., Mo H., Zhang Y., van den Bosch F. C., 2011, ApJ, 741, 13
  * Zhao et al. (2009) Zhao D., Jing Y., Mo H., Börner G., 2009, ApJ, 707, 354

##  Appendix A Mass history definition

As discussed in the main text, the MAH of a dark matter halo can be defined in
two ways. The most popular is to use mass of the main progenitor by walking
the main branch of the merging tree step by step. The main progenitor may not
be the most massive progenitor at an earlier time. The other definition,
adopted by G05, is to use the mass of the most massive progenitor (MMP) at
each time. In Figure LABEL:fig:MAH_and_MMP we compare the median of the two
definitions of MAHs for 7 halo samples of different mass from the Cosmic-Zoom
simulation. We overplot predictions from G05. The right panel shows the
residuals relative to the MAHs for the main branch. At low redshift, the
results for the two definitions are very similar, but at higher redshift, they
begin to diverge more and more. The G05 predictions agree better with the
simulation results when the MMP definition is used, as it should.