††thanks<EMAIL_ADDRESS>
# A Differentiable Model of the Assembly of Individual and Populations of Dark
Matter Halos

Andrew P. Hearin1,⋆ Jonás Chaves-Montero1,2 Matthew R. Becker1 Alex Alarcon1
1HEP Division, Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL
60439, USA  2Donostia International Physics Centre, Paseo Manuel de Lardizabal
4, 20018 Donostia-San Sebastian, Spain

###### Abstract

We present a new empirical model for the mass assembly of dark matter halos.
We approximate the growth of individual halos as a simple power-law function
of time, where the power-law index smoothly decreases as the halo transitions
from the fast-accretion regime at early times, to the slow-accretion regime at
late times. Using large samples of halo merger trees taken from high-
resolution cosmological simulations, we demonstrate that our 3-parameter
model, Diffmah, can approximate halo growth with a typical accuracy of 0.1 dex
for $t\gtrsim 1\,{\rm Gyr}$ for all halos of present-day mass $M_{\rm
halo}\gtrsim 10^{11}{\rm M}_{\odot},$ including subhalos and host halos in
gravity-only simulations, as well as in the IllustrisTNG hydrodynamical
simulation. We additionally present a new model for the assembly of halo
populations, DiffmahPop, which not only reproduces average mass growth across
time, but also faithfully captures the diversity with which halos assemble
their mass. Our python implementation is based on the autodiff library JAX,
and so our model self-consistently captures the mean and variance of halo mass
accretion rate across cosmic time. We show that the connection between halo
assembly and the large-scale density field, known as halo assembly bias, is
accurately captured by Diffmah, and that residual errors in our approximations
to halo assembly history exhibit a negligible residual correlation with the
density field. Our publicly available source code can be used to generate
Monte Carlo realizations of cosmologically representative halo histories; our
differentiable implementation facilitates the incorporation of our model into
existing analytical halo model frameworks.

###### keywords:

Cosmology: large-scale structure of Universe; methods: N-body simulations

##  1 Introduction

In the standard cosmological model, the matter content of the Universe is
dominated by Cold Dark Matter (CDM), and gravitationally self-bound objects
referred to as dark matter halos are the fundamental building blocks of
structure formation (see Mo et al., 2010, for a comprehensive review). Dark
matter halos are the natural sites of galaxy formation (White & Rees, 1978;
Blumenthal et al., 1984), and so a detailed understanding of the buildup and
evolution of halos is a key ingredient of any theory of structure growth.

The basic physical picture of halo evolution is now relatively well
understood. At early times, halo assembly is characterized by a “fast-
accretion phase” in which the growth of halo mass is rapid and major mergers
are common; mass accretion rates diminish considerably at later times in a
“slow-accretion phase”, during which time major mergers are comparatively rare
(Bullock et al., 2001; Wechsler et al., 2002; Tasitsiomi et al., 2004). The
mass density profile of dark matter halos is well-described by a double power
law known as the NFW profile (Navarro et al., 1997), an approximation that
remains reasonably accurate across most of cosmic time (although see Ludlow et
al., 2013, for shortcomings of this approximation). Although the concentration
parameter that defines the NFW profile exhibits a well-known dependence upon
total mass (e.g., Diemer & Kravtsov, 2015; Child et al., 2018), all halos in
the fast-accretion phase tend to have similar values of concentration
$c\approx 3-4$ (Zhao et al., 2003a), and as halos transition to the slow-
accretion phase, their concentrations steadily increase as they predominantly
pile up mass onto their outskirts (see Wang et al., 2020, for a detailed
physical picture of the relationship between halo assembly and NFW
concentration).

Many additional characteristics of the internal structure of dark matter halos
are closely related to their assembly histories, including substructure
abundance (Gao et al., 2004; Wechsler et al., 2006; Giocoli et al., 2010), the
“splashback” feature in the outer profile (Diemer & Kravtsov, 2014), and
ellipsoidal shape (Chen et al., 2020a; Lau et al., 2021). The growth history
of a dark matter halo is also tightly connected to the larger-scale
environment in which it evolves, a phenomenon known as halo assembly bias
(e.g., Gao et al., 2005; Mao et al., 2018; Mansfield & Kravtsov, 2020).

The mass assembly history of dark matter halos is widely used in theoretical
models of galaxy formation as a fundamental quantity regulating the star
formation rate of the halo’s resident galaxy, including early formulations of
traditional semi-analytic modeling (Kauffmann et al., 1993; Baugh et al.,
1996; Avila-Reese et al., 1998, SAMs, hereafter), practically all contemporary
SAMs (e.g., Bower et al., 2006; Somerville et al., 2008; Henriques et al.,
2015; Lagos et al., 2018), simple empirical models (Watson et al., 2015;
Behroozi & Silk, 2015), and recent empirical approaches (Becker, 2015; Moster
et al., 2018; Behroozi et al., 2019). Strong motivation for treating halo
assembly as the backbone of galaxy formation models also comes directly from
hydrodynamical simulations, which show that the accretion rate of baryonic
mass closely tracks that of dark matter (Wetzel & Nagai, 2015).

The assembly of a halo is fully characterized by its merger tree, which
describes the mass and accretion time of the “progenitors” that merged into
the main halo over time, as well as the progenitors of those progenitors, etc.
The most common theoretical technique for modeling halo merger trees is the
(extended) Press-Schechter formalism (Press & Schechter, 1974; Bond et al.,
1991; Bower, 1991; Lacey & Cole, 1993, EPS, hereafter). The EPS framework
yields considerable insight into many aspects of the structure and buildup of
dark matter halos (e.g., Somerville & Kolatt, 1999; Nusser & Sheth, 1999; Cole
et al., 2008), and has been shown to yield reasonably accurate predictions for
the average growth of halos over time (van den Bosch, 2002), the mass function
of the progenitor halos (Parkinson et al., 2008), and the statistical
distribution of halo formation times (Li et al., 2007; Giocoli et al., 2012).
Merger trees predicted by EPS furthermore enable suitably implemented SAMs to
generate predictions for galaxy evolution without the need for high-resolution
N-body simulations (e.g., Somerville et al., 2008; Benson, 2012). We refer the
reader to Zentner (2007) for a pedagogical overview of EPS, and to Jiang & van
den Bosch (2014) for a comparison of modern implementations of EPS predictions
of merger trees.

Empirically formulated fitting functions calibrated against cosmological
simulations offer a simpler route to describing the evolutionary history of
halo mass. The one-parameter exponential model developed in Wechsler et al.
(2002) is one of the first of such models, and a number of improvements have
followed this early empirical approach that capture richer phenomena with
greater accuracy. For example, the fitting formula in Tasitsiomi et al. (2004)
is a two-parameter generalization of the exponential growth model with
improved performance for a greater diversity of cluster-mass halos. In
Fakhouri et al. (2010), the authors present a simple and accurate fitting
formula for $\langle dM_{\rm halo}(t)/dt|M_{0}\rangle,$ the average mass
accretion rate of halos as a function of cosmic time and present-day mass
$M_{0},$ and they furthermore find that they can use their model to accurately
predict $\langle M_{\rm halo}(t)|M_{0}\rangle$ by numerically integrating
their parameterized mass accretion rate histories. In McBride et al. (2009),
the authors introduce a different two-parameter fitting function for
individual halo mass growth, and also include additional modeling ingredients
for the dependence of how the fitting parameters depend on present-day halo
mass.

The program to develop empirical approximations to dark matter halo assembly
is faced with a number of challenges. Each of the efforts summarized above
were conducted at fixed cosmology, and also for gravity-only N-body
simulations, and the robustness of these findings to effects of cosmology and
hydrodynamics remains relatively coarsely characterized (although see van den
Bosch et al., 2014, discussed further in §5 below). Efforts to improve
flexibility beyond the one-parameter exponential model also face a significant
difficulty of ensuring physically consistent behavior across the full range of
mass and redshift; for example, as pointed out in Appendix H of Behroozi et
al. (2013a), predictions of the model presented in McBride et al. (2009)
exhibit unphysical behavior at high redshift due to the crossing of average
assembly histories of halos of different present-day mass. Of special concern
to the present paper, empirical fitting functions historically only reproduce
average halo growth across time, but do not enjoy capability to predict the
diversity of halo assembly, a shortcoming made especially glaring through
comparison to EPS-based techniques.

In this paper, we seek to improve upon this situation by developing an
empirical approach to modeling the growth of both individual halos, as well as
the assembly of cosmologically representative populations. In §3, we will
introduce a new parameterization for individual halo growth, Diffmah,
validating that our formulation is sufficiently flexible to capture the
assembly of halos in both the gravity-only and full-physics hydro simulations
described in §2. In §4, we outline DiffmahPop, our statistical model for the
assembly of halo populations, relegating exposition on mathematical and
computational details to the appendices. We discuss our results and potential
applications in §5, and conclude in §6 with a brief summary of our primary
results.

##  2 Simulations and Merger Trees

We have built our model for the mass assembly history (MAH hereafter) of dark
matter halos to approximate results taken from the merger trees of halos
identified in both gravity-only and hydrodynamical simulations. For our
gravity-only simulations, we use a combination of the Bolshoi Planck
simulation (BPL, Klypin et al., 2011) and MultiDark Planck 2 (MDPL2, Klypin et
al., 2016). The BPL simulation was carried out using the ART code (Kravtsov et
al., 1997) by evolving $2048^{3}$ dark-matter particles of mass $m_{\rm
p}=1.55\times 10^{8}{\rm M}_{\odot}$ on a simulation box of $250\,{\rm Mpc}$
on a side; the MDPL2 simulation was carried out using the L-GADGET-2 code
(Springel, 2005) with $3840^{3}$ dark-matter particles of mass $m_{\rm
p}=1.51\times 10^{9}M_{\odot}$ on a simulation box of $1000\,{\rm Mpc}$ on a
side; both simulations were run under cosmological parameters closely matching
Planck Collaboration et al. (2014). For both BPL and MDPL2, we used publicly
available111https://www.peterbehroozi.com/data.html merger trees that were
identified with Rockstar and ConsistentTrees (Behroozi et al., 2013b, c;
Rodríguez-Puebla et al., 2016).

We additionally studied the assembly of halos in the IllustrisTNG
hydrodynamical simulations. IllustrisTNG is a suite of hydro simulations that
models the joint evolution of dark matter, gas, stars, and supermassive black
holes, and includes treatment of radiative gas cooling, star formation,
galactic winds, and AGN feedback (Weinberger et al., 2017; Pillepich et al.,
2018). We use publicly available data from the largest hydrodynamical
simulation of the suite, TNG300-1 (Nelson et al., 2018). The TNG300-1
simulation was carried out using the moving-mesh code Arepo (Springel, 2010)
$2500^{3}$ with gas tracers together with the same number of dark matter
particles in a simulation box of $205$ Mpc on a side under a cosmology very
similar to Planck Collaboration et al. (2014). For TNG300-1 the corresponding
mass resolution is $4.1\times 10^{7}{\rm M}_{\odot}$ and $7.7\times 10^{6}{\rm
M}_{\odot}$ for dark matter and gas, respectively. Halos and subhalos in
IllustrisTNG were identified with the SUBFIND algorithm (Springel et al.,
2001), and the merger trees we use were constructed with SUBLINK (Rodriguez-
Gomez et al., 2015).

For some of the results in the paper, we calculate cross-correlation functions
between simulated halos and the density field, using a random downsampling of
particles from the appropriate snapshot. For BPL we use a random downsampling
of dark matter particles provided by halotools (Hearin et al., 2017). For
IllustrisTNG, we use a random downsampling of TNG300-1 particles at $z=0$
kindly provided by Benedikt Diemer.

All results in the main body of the paper pertain to the assembly histories of
present-day host halos (i.e., upid=-1 for Rockstar, and the “main halo” for
SUBFIND). Comparable results for the case of subhalos can be found in Appendix
D. Throughout the paper, including the present section, values of mass and
distance are quoted assuming $h=1.$ For example, when writing $M_{\rm
peak}=10^{12}{\rm M}_{\odot},$ we suppress the ${\rm M}_{\odot}/h$ notation
and write the units as ${\rm M}_{\odot}.$

Figure 1: Physical interpretation of the parameters of the Diffmah model for
individual halo assembly. Each panel shows the assembly history of a dark
matter halo with present-day mass $M_{0}=10^{12}{\rm M}_{\odot};$ the top row
of panels shows the history of cumulative peak mass, $M_{\rm peak}(t);$ the
bottom row of panels shows the history of mass accretion rate, ${\rm d}M_{\rm
peak}/{\rm d}t.$ Our model approximates halo assembly as a power-law function
of time with rolling index, $M_{\rm peak}(t)\propto t^{\alpha(t)},$ as in Eq.
1. The left columns show the influence of $\alpha_{\rm early},$ the power-law
index at early times; the middle columns show the influence of $\alpha_{\rm
late},$ the power-law index at late times; the right columns show $\tau_{\rm
c},$ which controls the transition time between the fast- and slow-accretion
regimes. Smaller values of each parameter correspond to halos with earlier
formation times.

##  3 Modeling Individual Halo Assembly with Diffmah

Figure 2: Example fit to individual halo growth. The blue curve in each panel
shows the assembly history of the same dark matter halo taken directly from
the merger trees of the IllustrisTNG simulation. The top panel shows the
history of the halo’s mass accretion rate, and the bottom panel shows the
history of its cumulative peak mass. In each panel, the orange curve shows the
approximate halo history based on our differentiable model.

In this section, we present results for Diffmah, our model of the assembly of
individual dark matter halos. In §3.1 we describe the functional form of the
model and its differentiable implementation, and in §3.2 we show the accuracy
with which the model can approximate the formation history of simulated halos.

###  3.1 Diffmah Model formulation

The Diffmah model for the growth of individual halos is designed to capture
the evolution of cumulative peak halo mass, $M_{\rm peak}(t),$ defined as the
largest mass the main progenitor halo has ever attained up until the time $t.$
Thus as a matter of definition, $M_{\rm peak}(t)$ is a non-decreasing function
for all $t,$ so that $M_{\rm peak}(t)$ will remain constant if a halo
experiences a period of mass loss.

We model $M_{\rm peak}(t)$ with a power-law function of cosmic time with a
rolling index,

$\displaystyle M_{\rm peak}(t)=M_{0}(t/t_{0})^{\alpha(t)},$ (1)

where $t_{0}$ is the present-day age of the universe, $M_{0}\equiv M_{\rm
peak}(t_{0}),$ and $\alpha(t)$ is given by a sigmoid function defined as
follows:

$\displaystyle\alpha(t;\tau_{\rm c},k,\alpha_{\rm early},\alpha_{\rm
late})\equiv\alpha_{\rm early}+\frac{\alpha_{\rm late}-\alpha_{\rm
early}}{1+\exp(-k(t-\tau_{\rm c}))}.$ (2)

The parameters $\alpha_{\rm early}$ and $\alpha_{\rm late}$ define the
asymptotic value of the power-law index at early and late times, respectively,
and $\tau_{\rm c}$ is the transition time between the early- and late-time
indices. The parameter $k$ controls the speed of the transition between the
two regimes; as described in Appendix A, for all results in the paper we hold
$k$ fixed to a constant value of $3.5.$

In Figure 1, we give a visual illustration of the physical interpretation of
the three free parameters of our model: $\alpha_{\rm early},$ $\alpha_{\rm
late},$ and $\tau_{\rm c}.$ Each of the six panels shows the assembly history
of halos with present-day mass $M_{0}=10^{12}{\rm M}_{\odot},$ with variations
of the MAH parameters color-coded according to the legend. For each of the
three MAH parameters, $\alpha_{\rm early},$ $\alpha_{\rm late},$ and
$\tau_{\rm c},$ smaller values of the parameter corresponds to a halo that was
accreting mass more rapidly at early times and more slowly today; thus a halo
with a smaller value of a MAH parameter has an earlier formation time relative
to a halo with a larger value of that parameter. For curves labeled “fiducial
model”, we have $(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm
c})=(2.5,0.3,1.25{\rm[Gyr]});$ for parameters set to their smaller values, we
have $(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm
c})=(1.25,0.05,0.6{\rm[Gyr]});$ for larger values, we have $(\alpha_{\rm
early},\alpha_{\rm late},\tau_{\rm c})=(5.0,0.6,2.5{\rm[Gyr]}).$ The top row
of panels shows the history of cumulative mass, $M_{\rm peak}(t),$ which is
parameterized directly according to Eq. 1. The bottom row of panels shows
${\rm d}M_{\rm peak}/{\rm d}t,$ the history of mass accretion rate. While
${\rm d}M_{\rm peak}/{\rm d}t$ is analytically calculable by symbolically
differentiating Eq. 1, we have computed each curve in the bottom panels of
Figure 1 via automatic differentiation (autodiff, see Baydin et al., 2015, for
a recent review), as our source code is implemented based on the JAX autodiff
library (Bradbury et al., 2018).

###  3.2 Diffmah Model performance

Figure 3: Residuals of Diffmah fits to individual halos. Using on a large
collection of fits to the assembly histories of simulated halos, each panel
shows the logarithmic difference between the simulated and approximate values
of cumulative peak mass, $M_{\rm peak}(t),$ as a function of time; results for
samples of halos of different present-day mass are shown in different panels
as indicated by the in-panel annotation. The black curve in each panel shows
the average residual difference, and the shaded band shows the $1\sigma$
scatter; residuals for gravity-only simulations are shown with the solid black
curve and purple shaded region; residuals for IllustrisTNG are shown with the
dashed black curve and green shaded region. The figure shows that the rolling
power-law model defined by Eq. 1 gives an unbiased approximation to halo mass
assembly with a typical error of 0.1 dex for $t\gtrsim 1$ Gyr.

In the previous section, we illustrated the basic features of the Diffmah
model for dark matter halo assembly. Here we show the accuracy with which
simulated halo histories can be approximated by Diffmah. When fitting the mass
growth of an individual dark matter halo with the fitting function in Equation
1, we allow the parameters $\alpha_{\rm early},$ $\alpha_{\rm late},$ and
$\tau_{\rm c}$ to vary freely, and use gradient descent to minimize the
logarithmic difference between the simulated and predicted values of $M_{\rm
peak}(t),$ with $k$ held fixed to a constant value of $3.5$ for all halos, and
$M_{0}$ held fixed to $M_{\rm peak}(t_{0}).$ We refer the reader to Appendix A
for a detailed description of our optimization technique.

####  3.2.1 Capturing individual halo growth

Figure 2 shows a typical example of the ability of Diffmah to approximate
individual dark matter halo growth. Both panels show the assembly history of
the same dark matter halo taken from the IllustrisTNG simulation; the top
panel shows the history of the halo’s mass accretion rate, and the bottom
panel shows the history of its cumulative peak mass. The blue curve in each
panel shows the assembly history of the halo as taken directly from the
simulated merger tree, and the orange curve shows the differentiable
approximation. As described in Appendix A, our differentiable approximation
comes from fitting $M_{\rm peak}(t)$ of the simulated halo, and our
approximation for ${\rm d}M_{\rm peak}/{\rm d}t$ is a derived quantity.

Using the optimization techniques detailed in Appendix A, we have identified a
set of best-fitting parameters for every halo in the samples described in §2.
Figure 3 displays the quality of these fits for both gravity-only N-body
simulations and IllustrisTNG. Each panel shows the logarithmic difference
between the simulated and best-fitting MAH plotted as a function of cosmic
time; results for samples of halos of different present-day mass are shown in
different panels as indicated by the in-panel annotation. The black curves in
each panel show the average residual difference, and the shaded band shows the
$1\sigma$ scatter; residuals for gravity-only simulations are shown with the
solid black curve and purple shaded region; residuals for IllustrisTNG are
shown with the dashed black curve and green shaded region.

Figure 3 shows that the rolling power-law model defined by Eq. 1 gives an
unbiased approximation to halo mass assembly with a typical error of 0.1 dex
for $t\gtrsim 1$ Gyr; residual errors are slightly larger at higher halo mass
due to the increased prominence of major mergers. In Appendix D, we show that
analogous results hold for the case of present-day host halos that have
experienced a flyby event at some point in their past history (Figure D13),
for present-day subhalos (Figure D14), and also for subhalos that have either
merged or become disrupted prior to $z=0$ (Figure D15). In Figure 3, we have
restricted attention to $t>1$ Gyr for the sake of ensuring good mass
resolution across the mass range of all our samples; this cut in cosmic time
is driven by our focus on the redshifts $z\lesssim 5$ that are most relevant
to cosmological surveys, and by the lack of a publicly-available, higher-
resolution simulation with homogeneously processed merger trees (see §A for
further discussion).

####  3.2.2 Capturing halo assembly bias

Even though Figure 3 shows that simulated MAHs are well approximated by the
Diffmah model, of course the agreement is not perfect, and a question that
naturally arises is whether the residual errors are correlated with the cosmic
density field. As mentioned in §1, at fixed present-day mass, the large-scale
density of halos exhibits a significant dependence upon assembly history, a
phenomenon referred to as halo assembly bias. This effect is commonly
quantified in terms of the two-point cross-correlation between halos and dark
matter particles, $\xi_{\rm\delta m}(r).$ For example, for a sample of dark
matter halos with the same value of $M_{0}\equiv M_{\rm peak}(t_{0}),$ the
quantity $\xi_{\rm\delta m}(r)$ varies with the formation time of the halo,
$t_{\rm form},$ defined as the first time the halo attains some specified
fraction of $M_{0}.$ The non-zero residual errors shown in Fig. 3 indicate
that each halo’s approximate and actual values of $t_{\rm form}$ will differ,
and it is plausible that the difference correlates with the density field in a
way that alters the $t_{\rm form}$-dependence of $\xi_{\rm\delta m}(r).$ For
example, halos with a merger-rich history tend to reside in denser regions
relative to halos with a more quiescent assembly, and the prevalence of
mergers may be correlated with the residual errors of our best-fitting
approximation.

To investigate this question, we select a sample of host halos in the BPL
simulation with $10^{11.9}{\rm M}_{\odot}\leq M_{0}\leq 10^{12.1}{\rm
M}_{\odot},$ divide the sample in half according to the median value of
$t_{\rm form}$ for the sample, and use halotools (Hearin et al., 2017) to
compute $\xi_{\rm\delta m}(r)$ for each subsample. We have repeated this
exercise using the value of $t_{\rm form}$ computed directly from each halo’s
simulated merger tree, and compared the clustering to the case when $t_{\rm
form}$ is defined instead based on the differentiable approximation to each
individual halo’s MAH. We display our results in Figure 4; the top panel shows
$t_{\rm form}=t_{4\%},$ and the bottom panel shows $t_{\rm form}=t_{50\%}.$ We
find a similar level of agreement for halos of different mass in both gravity-
only simulations, and Figure E17 in Appendix E shows the analogous result for
the case of IllustrisTNG. We have also verified that we obtain a similar level
of agreement between simulation- and model-based clustering when dividing our
halo samples based on quartiles of formation time rather than medians.

In summary, the results of §3 and Appendices D-E demonstrate that the
3-parameter fitting function defined by Eq. 1 is sufficiently flexible to
approximate $M_{\rm peak}(t)$ with an accuracy of better than 0.1 dex across
nearly all of cosmic time; that this level of accuracy is robust to variations
in halo growth due to hydrodynamics and baryonic feedback, and applies equally
well to subhalos and splashback centrals; and that the small residual errors
of the approximation exhibit an essentially negligible correlation with the
large-scale density field.

Figure 4: Assembly bias in the Diffmah model for individual halo growth. Using
a sample of host halos in the BPL simulation with $M_{\rm peak}(t_{0})\approx
10^{12}{\rm M}_{\odot},$ we divide the sample in half according to $t_{\rm
form},$ the median value of halo formation time for the sample, and for each
subsample we compute $\xi_{\rm\delta m}(r),$ the cross-correlation between
halos and dark matter particles. Red curves show $\xi_{\rm\delta m}(r)$ for
early-forming halos, blue curves for late-forming halos; solid curves show
results for the case where $t_{\rm form}$ is computed directly from each
halo’s simulated merger tree; dashed curves show results for the case where
$t_{\rm form}$ is defined by the differentiable approximation to each halo’s
assembly history. The top and bottom panels show results for different
definitions of $t_{\rm form}$ as indicated in the legend. The figure
demonstrates that the correlation between halo formation time and the density
field is retained when simulated merger trees are approximated with our
differentiable model.

##  4 Modeling the Assembly of Halo Populations with DiffmahPop

In the previous section, we presented Diffmah, a model for the assembly
history of individual dark matter halos defined by Eq. 1. Our model for
$M_{\rm peak}(t)$ is characterized by three parameters: $\tau_{\rm
c},\alpha_{\rm early},$ and $\alpha_{\rm late}.$ In this section, we present
DiffmahPop, a model for $P(M_{\rm peak}(t)|M_{0}),$ the statistical
distribution of assembly histories for halos of present-day peak mass,
$M_{0}.$ Having shown that individual simulated MAHs can be accurately
approximated by our parametric fitting function, our approach to modeling
$P(M_{\rm peak}(t)|M_{0})$ is to construct a suitably accurate statistical
model for $P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}|M_{0}),$ the
distribution of our MAH model parameters as a function of $M_{0}.$

In order to motivate the form of our model for $P(\alpha_{\rm
early},\alpha_{\rm late},\tau_{\rm c}|M_{0}),$ in Figure 5 we show two
different cross sections of this multi-dimensional distribution, focusing on
BPL halos with $M_{0}=10^{12}{\rm M}_{\odot}.$ The top panel shows a histogram
of $P(\tau_{\rm c}|M_{0}),$ which shows a roughly log-normal shape with a
modest bimodality in which two distinct sub-populations emerge: earlier-
forming halos with $\tau_{\rm c}\lesssim 2$ Gyr, and later-forming halos with
$\tau_{\rm c}\gtrsim 2$ Gyr. In the bottom panel of Figure 5, we show a
scatter plot of the two-dimensional distribution of $P(\beta_{\rm
e},\beta_{\ell}|M_{0}),$ where we define $\beta_{\ell}=U(\alpha_{\rm late}),$
the variable $\beta_{\rm e}=U(\alpha_{\rm early}-\alpha_{\rm late}),$ and the
function $U(s)\equiv\ln(\exp(s)-1).$ As described in Appendix A, the variables
$\beta_{\rm e}$ and $\beta_{\ell}$ are the actual quantities that we
programmatically vary when seeking best-fitting approximations to the growth
of individual halos, as this enforces physical constraints in the
approximation to each halo’s assembly. In the bottom panel of Figure 5, points
with $\tau_{\rm c}>2$ Gyr are shown in purple, and points with $\tau_{\rm
c}<2$ Gyr are shown in orange. The color-coding in the bottom panel brings the
bimodal structure of $P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm
c}|M_{0})$ into sharper relief; while the two-dimensional marginal
distribution $P(\alpha_{\rm early},\alpha_{\rm late}|M_{0})$ is rather
complex, the same two-dimensional distribution becomes significantly simpler
when transforming to the variables $\beta_{\rm e}$ and $\beta_{\ell},$ and
conditioning the distribution upon whether $\tau_{\rm c}<2$ Gyr. In the
remainder of this section, we will outline how we have leveraged this apparent
bimodality to build a simple two-population model of halo assembly; we refer
the reader to Appendix B for a more expansive discussion of the full three-
dimensional distribution, $P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm
c}|M_{0}).$

Figure 5: Bimodal distribution of Diffmah parameters. The figure shows the
distribution of best-fitting parameters of our model for individual halo
growth for a sample of host halos in the BPL simulation with $M_{\rm
peak}(t_{0})\approx 10^{12}{\rm M}_{\odot}.$ The top panel shows the
distribution of $\tau_{\rm c},$ which presents a modest bimodality of early-
forming halos ($\tau_{\rm c}<2$ Gyr) and late-forming halos ($\tau_{\rm c}>2$
Gyr). The bottom panel shows a scatter plot of $U(\alpha_{\rm early}),$ and
$U(\alpha_{\rm early}-\alpha_{\rm late}),$ where $U(s)\equiv\ln(\exp(s)-1),$
with points color-coded according to their value of $\tau_{\rm c}.$ The two
panels together highlight the bimodal structure of the three-dimensional
probability distribution, $P(\tau_{\rm c},\alpha_{\rm early},\alpha_{\rm
late}|M_{0}).$ We leverage this bimodality when building our model for the
growth of halo populations.

We model $P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}|M_{0})$ as a
two-population normal distribution of three variables: $\log_{10}\tau_{\rm
c},$ $\beta_{\ell}=U(\alpha_{\rm late}),$ and $\beta_{\rm e}=U(\alpha_{\rm
early}-\alpha_{\rm late}).$ For BPL, MDPL2, and IllustrisTNG, we find that
relative abundance of the two populations is roughly $50\%$, but exhibits a
shallow mass-dependence. Each sub-population is described by an independently-
defined mean, $\mu_{\rm i},$ and covariance matrix, $\Sigma_{\rm i},$ each of
which also have a shallow, monotonic mass-dependence, so that $\mu_{\rm
i}=\mu_{\rm i}(M_{0})$ and $\Sigma_{\rm i}=\Sigma_{\rm i}(M_{0}).$ In
calibrating our model parameters for this mass dependence, our principal goal
is to faithfully describe $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm
peak}/{\rm d}t|M_{0}).$ Thus when optimizing our model parameters, for our
target data we use measurements of the first two moments of $P(M_{\rm
peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm d}t|M_{0}),$ rather than
attempting to directly recover the distribution $P(\alpha_{\rm
early},\alpha_{\rm late},\tau_{\rm c}|M_{0}).$ Briefly, we tabulate $P(M_{\rm
peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm d}t|M_{0})$ using several
narrow bins of halo mass, trimming $3\sigma$ outliers before computing mean
and variance; we optimize the parameters of our model for $\mu_{\rm i}(M_{0})$
and $\Sigma_{\rm i}(M_{0})$ by minimizing the logarithmic difference between
the mean and variance of simulated and predicted values for $P(M_{\rm
peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm d}t|M_{0}).$ Our
differentiable formulation in JAX enables us to optimize our parameters with
the gradient-based BFGS algorithm implemented in scipy (Broyden, 1970;
Fletcher, 1970; Goldfarb, 1970; Shanno, 1970). We refer the reader to Appendix
C for a detailed description of our application of these techniques.

Figure 6: Average assembly history in the DiffmahPop model for halo
populations. The top panel shows the average cumulative peak mass as a
function of cosmic time, and the bottom panel shows average accretion rates.
The histories of halos with different present-day mass are color-coded as
indicated in the legend. Solid curves show the average histories of simulated
halos; dashed curves show the corresponding quantity approximated by our model
for the assembly of halo populations. Figure 7: Diversity of assembly
histories in the DiffmahPop model for halo populations. In the left column of
panels, we plot the history of cumulative mass of simulated halos, $M_{\rm
peak}(t),$ showing its average value with the solid black curve, and the
$1\sigma$ standard deviation with the gray band. The right hand column of
panels shows the mean and standard deviation of halo mass accretion rates. The
history of halos of different present-day mass are shown in different rows of
panels as indicated by the in-panel annotation. In all panels, the dashed
green curves show our model for the assembly of halo populations.

Figure 8: Example comparisons of halo assembly distributions. In the top
panel, we focus on halos with present-day mass $M_{0}=10^{14}{\rm M}_{\odot},$
and show the distribution of masses at higher redshift, $M_{\rm peak}(z),$
with results at different redshift color-coded as indicated in the legend. In
the bottom panel, for halos with $M_{0}=10^{12}{\rm M}_{\odot},$ we show the
distribution of formation times, $t_{\rm form},$ using several different
definitions as indicated in the legend. Dashed black curves in each panel show
the corresponding predictions based on our model for the assembly of halo
populations. Plotted quantities shown in this figure are distinct from those
shown in Fig. 7 in that here we show comparisons of the full probability
distributions of quantities computed directly from simulated merger trees
without trimming $3\sigma$ outliers.

In Figure 6, we compare the average assembly of halos in gravity-only
simulations to our best-fitting model; the top panel shows the average
cumulative peak mass as a function of cosmic time, and the bottom panel shows
average accretion rates. Assembly histories for halos of different present-day
mass are color-coded as shown in the legend, with dashed curves showing the
model, and solid curves showing the corresponding quantity in the simulations;
we use BPL for halos with $M_{\rm peak}\leq 10^{13.5}{\rm M}_{\odot},$ and
MDPL2 for more massive halos.

Several of the basic physical principles of CDM structure formation are on
plain display in Figure 6. In the top panels, we can see that for halos of all
mass, the slope of $M_{\rm peak}(t)$ is steeper at early times relative to
late times; the flattening of the slope of $M_{\rm peak}(t)$ is a
manifestation of the evolution of halos from the fast-accretion regime to the
slow-accretion regime. Structure growth in CDM is said to be “hierarchical”,
such that the formation of smaller halos precedes that of more massive halos.
The hierarchical nature of halo assembly is visible in both panels, as we can
see that the slope of the bluer curves reaches an inflection point at earlier
times relative to the redder curves. In the bottom panel we can see the
influence of dark energy on halo mass accretion rates: at low redshift when
the cosmic acceleration becomes significant, the value of ${\rm d}M_{\rm
peak}/{\rm d}t$ is declining for halos of all mass, with lower-mass halos
becoming subject to this trend at earlier epochs.

In Figure 7, we show that the DiffmahPop model for the assembly of halo
populations captures both the average growth of halos, as well as the
diversity of halo assembly. In the left column of panels of Fig. 7, we plot
$\langle M_{\rm peak}(t)|M_{0}\rangle,$ and in the right column of panels we
plot $\langle{\rm d}M_{\rm peak}/{\rm d}t|M_{0}\rangle;$ in all panels, solid
black curves show the quantity taken from the simulation, gray shaded bands
show the $1\sigma$ standard deviation of the simulated quantity, and the
dashed green curves show the results of our best-fitting model. Each row of
panels corresponds to the history of halos with different present-day mass, as
indicated by the in-panel annotation.

As mentioned above, the plotted data in Figure 7 are the first and second
moments of $P(M_{\rm peak}(t)|M_{0}),$ which we used directly as target data
to calibrate our model parameters. In Figure 8, we show comparisons between
simulations and our model predictions for the full distribution, $P(M_{\rm
peak}(t)|M_{0}).$ In the top panel of Fig. 8, we focus on a sample of host
halos in MDPL2 with present-day mass of $M_{0}=10^{14}{\rm M}_{\odot},$ and
show the distribution of halo masses at several different redshifts as
indicated in the legend. In the bottom panel of Fig. 8, we focus on host halos
in BPL with present-day mass of $M_{0}=10^{12}{\rm M}_{\odot},$ and show
distributions of halo formation time, $t_{\rm form},$ defined as the first
time the halo attained some specified fraction of its present-day mass;
different shaded histograms show different choices for the definition of
$t_{\rm form}$ as indicated in the legend. In both panels, each shaded
histogram is accompanied by a dashed black curve showing the predictions of
our differentiable model of halo populations.

Fig. 8 does not display the same level of agreement between simulations and
model as illustrated in Figure 7, because in Fig. 8 we show the full
probability distributions taken directly from the simulated merger trees
without trimming $3\sigma$ outliers. We refer the reader to §C.2 for a
detailed description of our outlier exclusion algorithm, which has the effect
of excluding 2-3% of halos at each mass that have recently experienced an
unusually large major merger. Thus the shaded histograms in Fig. 8 include
contributions from short-term fluctuations associated with recent mergers,
whereas we optimized our model for halo populations using a sample that
excludes such contributions. This difference can be seen in two different
manifestations in Fig. 8. First, the PDFs shown in the shaded histograms have
slightly broader support relative to our model, reflecting the increased
variance associated with the inclusion of outliers. Second, the shaded
histograms tend to be shifted towards slightly earlier redshifts relative to
the model; this second effect is due to our election to model the evolution of
peak halo mass, $M_{\rm peak}(t),$ whose associated formation times are
asymmetrically shifted to earlier epochs by major mergers. Both of these
effects are visible in Fig. 8. The full distributions for Milky Way mass halos
shown in the bottom panel display a somewhat better level of agreement
relative to the full distributions for cluster-mass halos shown in the top
panel; this is sensible since the assembly of more massive halos tends to be
more rich in mergers.

##  5 Discussion & Future Work

We have presented a new model for the assembly history of individual dark
matter halos, Diffmah. In our model, halo mass is a power-law function of
cosmic time with a rolling index, $M_{\rm peak}(t)\propto t^{\alpha(t)},$
where $\alpha(t)$ is a sigmoid function that smoothly transitions halo growth
from the early-time, fast-accretion regime, to the late-time, slow-accretion
regime (see Eqs. 1-2). We note that we have chosen to model peak halo mass,
$M_{\rm peak}(t),$ the largest mass attained by the halo up until the time
$t,$ and so our model is not intended to capture the instantaneous halo mass,
$M_{\rm halo}(t),$ which is subject to additional physical effects whose
treatment is beyond the scope of the present paper.

The results in §3 in the main body of the paper demonstrate that for host
halos identified with Rockstar over a wide range of present-day mass,
$M_{0}\equiv M_{\rm peak}(z=0)\gtrsim 10^{11}{\rm M}_{\odot},$ the Diffmah
model faithfully describes the mass assembly of individual halos over nearly
all of cosmic time, $t\gtrsim 1$ Gyr, with a typical accuracy of $\sim 0.1$
dex or better. In Appendix D, we show that our model for $M_{\rm peak}(t)$
enjoys the same level of accuracy regardless of whether the host halo
experiences a splashback event in its past history, or whether the halo
eventually becomes a subhalo of a more massive host, including cases where the
subhalo is disrupted or merges with its host. In Appendix E (together with
Fig. 3), we show that our model achieves the same level of accuracy in its
description of $M_{\rm peak}(t)$ for both host halos and subhalos in
IllustrisTNG, a full-physics hydrodynamical simulation whose halos and merger
trees were identified with entirely different algorithms from those used in
the analysis of the gravity-only simulations we consider.

We have additionally presented a new model for the assembly of halo
populations, DiffmahPop. In typical empirical approaches to modeling halo
assembly, average halo growth is directly parameterized with a fitting
function, and variance in halo assembly is modeled simply as random scatter.
Our alternative approach is formulated by using the Diffmah model for
individual halo growth as its fundamental basis, and by building a model for
the statistical distribution of the Diffmah parameters that regulate
individual halo assembly. As a result, DiffmahPop not only reproduces $\langle
M_{\rm peak}(t)|M_{0}\rangle,$ but also a cosmologically realistic diversity
of individual halo trajectories across time, $P(M_{\rm peak}(t)|M_{0}).$ The
DiffmahPop model for the growth of halo populations is differentiable courtesy
of its implementation in the autodiff library JAX, and accurately captures the
average mass accretion rate across cosmic time, $\langle{\rm d}M_{\rm
peak}/{\rm d}t|M_{0}\rangle,$ as well as the diversity of accretion rates
$P({\rm d}M_{\rm peak}/{\rm d}t|M_{0}).$ As a convenience, our publicly
available python code includes functionality to generate a Monte Carlo
realization of the diversity of halos with the same $M_{0}.$

The accuracy of our approximation may not be surprising as previous results
have shown that halo mass is predominantly assembled via a combination of
smooth accretion and mergers with a very small mass ratio (e.g., Genel et al.,
2010). However, our model could naturally be used as the basis of empirical
techniques to predict merger rates and/or the subhalo mass function at the
time of accretion (also referred to as the unevolved subhalo mass function, or
USHMF hereafter). For example, in Yang et al. (2011), the authors developed an
analytical model that accurately predicts the USHMF; as the basis of their
model, they used the fitting function of Zhao et al. (2003b) for the median
assembly history of dark matter halos, supplementing this with a simple log-
normal distribution to describe random scatter in the main-branch mass as a
function of time. The techniques in the present work could be used to improve
upon such an approach, as our model provides an accurate description of
$P(M_{\rm peak}(t)|M_{0}),$ as shown in Figure 7. Such an extension could also
prove fruitful when used in concert with subhalo orbital evolution codes
(e.g., Zentner et al., 2005; Jiang et al., 2021). Recent progress in modeling
subhalo orbital evolution has demonstrated the importance of accounting for
the evolution of the host halo potential (Ogiya et al., 2021), and so
deploying our model in this context could be leveraged to capture physically
realistic accretion-rate correlations in the halo-to-halo variance of
substructure abundance and orbits across cosmic time (see, e.g., Jiang & van
den Bosch, 2017). Along similar lines, it has recently been shown that scatter
in the thermal Sunyaev-Zel’dovich effect is largely driven by variance in the
assembly histories of cluster-mass halos (Green et al., 2020), and so our
model may also help improve the predictive power of models for the mass-
observable relation of galaxy clusters.

While previous work has demonstrated that the cosmology-dependence of halo
assembly is relatively simple (e.g., Zhao et al., 2009), all results in the
present work are limited to a fixed set of cosmological parameters that
closely match those in Planck Collaboration et al. (2014). In van den Bosch et
al. (2014), it was shown that halo assembly histories have a universal form
controlled entirely by the cosmology-dependence of cosmic time, and so it may
be relatively straightforward to extend our model for $P(M_{\rm
peak}(t)|M_{0})$ to capture the influence of cosmological parameters on halo
growth. Such an extension could be directly incorporated into frameworks that
model cosmology-dependence via physically-motivated rescalings (e.g., Angulo &
White, 2010; Renneby et al., 2018; Aricò et al., 2020). In principle, our
model could also be used in applications with emulators built upon
cosmological suites of simulations (e.g., Heitmann et al., 2016; DeRose et
al., 2019; Nishimichi et al., 2019), although this would require suites with
high mass-resolution that include merger trees, which remains uncommon
(although see Contreras et al., 2020, for recent progress).

Our restriction to cosmological parameters similar to Planck Collaboration et
al. (2014) is just one example of how in the present work, we have not
attempted to provide a comprehensive study of $P(M_{\rm peak}(t)|M_{0}).$ As
another example, even though we have provided separate calibrations of our
model for halos in gravity-only simulations and in IllustrisTNG, we have not
attempted a detailed characterization of the differences between the assembly
histories of halos in the two simulations. For such a purpose, a simulation
suite such as CAMELS (Villaescusa-Navarro et al., 2020) that spans a range of
baryonic feedback parameters would enable a more rigorous study than an
uncontrolled comparison between IllustrisTNG and N-body halos. Furthermore,
when calibrating our models for $P(M_{\rm peak}(t)|M_{0}),$ we have focused on
the mass range $10^{11.5}{\rm M}_{\odot}\leq M_{\rm peak}\leq 10^{14.5}{\rm
M}_{\odot}$ as this range of halo masses contains thousands of halos that are
resolved by over $2000$ particles in the N-body simulations we use; however, a
series of nested-resolution boxes would be considerably more effective towards
the goal of extending the mass range and conducting rigorous resolution tests.
As discussed in §A, we have similarly focused on the range of cosmic time
$t>1$Gyr to avoid potential resolution issues in our simulations, but our work
could be extended to redshifts $z>5$ with a simulation suite that would permit
such resolution tests. Improvements such as these would require a dedicated
effort along the lines of what has been done to calibrate models of the halo
mass function (e.g., Jenkins et al., 2001; Tinker et al., 2008; McClintock et
al., 2019; Bocquet et al., 2020); we thus consider these improvements beyond
our present scope, as here we merely demonstrate that a modeling approach such
as ours has capability to deliver a precision tool.

When fitting the growth history of large samples of simulated halos, we find
that the distribution of best-fitting parameters exhibits a clear bimodality,
with earlier-forming and later-forming halos occupying distinct regions of the
three-dimensional parameter space of $\alpha_{\rm early},\alpha_{\rm late},$
and $\tau_{\rm c}.$ Bimodal halo growth has been reported previously in Shi et
al. (2018) in terms of the distribution of halo formation times, although this
previous finding applied only to infall halos that eventually became subhalos
of some more massive host. The bimodality we find pertains to host halos at
$z=0,$ applies to halos of all mass, persists even when excluding halos that
have previously experienced a flyby event, and is present with comparable
strength in both gravity-only simulations and IllustrisTNG. We have leveraged
this bimodal structure in building our analytical model for $P(M_{\rm
peak}(t)|M_{0}),$ as capturing this feature considerably simplifies the task
of achieving the level of agreement with simulations shown in Fig. 7,
particularly for the case of the second moment.

The two-point clustering of dark matter halos of the same mass has well-known
dependence upon halo formation time, $t_{\rm form},$ a phenomenon known as
halo assembly bias. In Figure 4, we showed that the $t_{\rm form}$-dependence
of halo clustering remains the same whether one uses the actual formation time
taken directly from simulated merger trees, or the value of $t_{\rm form}$
taken from our best-fitting differentiable model. This result implies that our
model can approximate halo growth in such a way that residual errors in the
fit are uncorrelated with the density field on large scales. This finding is
especially interesting in light of the relationship between assembly bias and
halo concentration (Wechsler et al., 2006), the latter of which is closely
connected to assembly history (Wechsler et al., 2002; Chen et al., 2020b), and
is known to be strongly influenced by mergers (e.g., Wang et al., 2020). We
will investigate this issue systematically in follow-up work to the present
paper.

As discussed in §1, halo assembly history plays a fundamental role in
contemporary semi-analytic models of galaxy formation (SAMs), since the star
formation rate of a main sequence galaxy is commonly assumed to be
proportional to the mass accretion rate of its parent halo. We have formulated
our model for the assembly of halo populations with numerous applications to
SAMs squarely in mind. For example, our model has potential to considerably
accelerate efforts to explore the portions of SAM parameter space that
regulate $\epsilon_{\rm eff},$ the efficiency with which galaxies transform
accreted mass into stars, since running SAMs on main-branch histories
alleviates the need to solve galaxy formation ODEs along every branch of the
merger tree. On the one hand, this computational speedup comes at the cost of
neglecting the role of mergers on SAM-predicted galaxy properties. On the
other hand, the capability of our model to generate a physically realistic
diversity of smooth, individual halo assembly histories enables a rather
comprehensive investigation of the specific role of mergers in shaping the
properties of the galaxy population. We intend to pursue these directions in
future work based on the Galacticus semi-analytic model (Benson, 2012);
Galacticus already includes an independent implementation of our model in its
publicly available source
code222https://github.com/galacticusorg/galacticus/wiki, which may be useful
in its own right for researchers more accustomed to programming in Fortran
rather than Python.

As the principal aim of our future work, we are using our model for halo
assembly as the foundation of a new approach to forward modeling cosmological
structure formation, Surrogate Modeling the Baryonic Universe (SBU). The basis
of SBU is a minimally flexible parameterized family of solutions to the
differential equations of galaxy formation; we make predictions for galaxy
populations by building a statistical model for a cosmologically realistic
diversity of individual galaxy histories (Alarcon et al., 2021). The approach
taken in the present work to modeling halo populations essentially provides a
template for this framework.

This forward modeling effort is already well underway. In Chaves-Montero &
Hearin (2020), we studied the influence of star formation history (SFH) on the
broadband photometry of galaxies and found that, despite the complexity of the
processes that impact optical color, physically-motivated variations in SFH
correspond to changes in a unique direction in color-color space. In a follow-
up paper (Chaves-Montero & Hearin, 2021), we studied how the presence of
unresolved SFH fluctuations impact model predictions for broad-band
photometry, finding that the statistical distributions of the broad-band
colors of a galaxy population are quite insensitive to short-term variability.
These results taken together indicate that the influence of halo assembly upon
galaxy photometry is likely to be simple, and that the absence of an explicit
treatment of halo mergers in our model may not pose a problem for the ability
of SBU models to make accurate predictions for the color distributions
observed by imaging surveys such as the Dark Energy
Survey333https://www.darkenergysurvey.org/ (DES, Dark Energy Survey
Collaboration, 2016) and the Legacy Survey of Space and
Time444http://www.lsst.org/ (LSST, Ivezić et al., 2019; LSST Science
Collaboration et al., 2009). Furthermore, the results shown in Figure 4
indicate that SBU predictions for color-dependent clustering are also likely
to be unbiased by the smooth nature of our model for halo assembly.

We expect that the differentiable formulation of our halo assembly model based
on the JAX autodiff library will play a critical role in this forward-modeling
program. Modern analyses of large-scale structure rely upon a considerable
number of “nuisance parameters” in order to ensure that cosmological inference
is robust to systematic uncertainty. Conventional Bayesian techniques such as
MCMC scale poorly with the dimension of the model parameter space, and since
the statistical precision of cosmological surveys will dramatically improve in
the 2020s, the parameter space of models used to analyze large-scale structure
measurements will only continue to increase. The capability to make
differentiable predictions addresses this growing problem, since gradient-
based optimization techniques such as Adam are extremely efficient even in
very high dimensions (Kingma & Ba, 2014), as are gradient-based inference
methods such as Hamiltonian Monte Carlo (Hoffman & Gelman, 2014). In the
present work, we have shown that not only is our model for individual halo
growth differentiable, but using the weighted sampling techniques described in
Appendix C, so are one-point function predictions based on our model for the
assembly of halo populations. In a closely related paper to the present work
(Hearin et al., 2021), we will introduce a new theoretical framework for
making differentiable predictions for two-point functions such as galaxy
clustering and lensing, using abundance matching (Kravtsov et al., 2004) as a
toy model for the galaxy-halo connection. The halo assembly model presented
here will provide a core ingredient to our program to make differentiable,
multi-wavelength predictions for large-scale structure that are physically
self-consistent across cosmic time.

##  6 Summary

We conclude by summarizing our primary results:

  1. 1.

We have introduced Diffmah, a new fitting function for $M_{\rm peak}(t),$ the
evolution of cumulative peak mass of individual dark matter halos; our model
approximates halo assembly as a power-law function of cosmic time with rolling
index, $M_{\rm peak}(t)\propto t^{\alpha(t)}$ (see Eq. 1). Using both gravity-
only simulations and IllustrisTNG, we have demonstrated that the Diffmah model
can approximate halo growth with an accuracy of 0.1 dex over the range
$t\gtrsim 1\,{\rm Gyr}$ for halos of present-day mass $M_{0}\gtrsim
10^{11}{\rm M}_{\odot}.$

  2. 2.

We find that the connection between halo assembly and the large-scale density
field, known as halo assembly bias, is entirely captured by Diffmah, and that
residual errors of our differentiable approximation to halo growth exhibit a
negligible correlation with the density field on large scales.

  3. 3.

We have introduced a new model for the growth of halo populations, DiffmahPop,
and shown that our model faithfully reproduces the evolution of average halo
mass, $\langle M_{\rm peak}(t)|M_{0}\rangle,$ and average mass accretion rate,
$\langle{\rm d}M_{\rm peak}/{\rm d}t|M_{0}\rangle,$ for cosmic times $t\gtrsim
1\,{\rm Gyr}.$ The DiffmahPop model additionally captures the diversity of
halo mass assembly, $P(M_{\rm peak}(t)|M_{0}),$ as well as the diversity of
accretion rates, $P({\rm d}M_{\rm peak}/{\rm d}t|M_{0}).$

  4. 4.

Our python code, diffmah, is publicly available and can be installed with pip
or conda. The repository for our source code is on GitHub,
https://github.com/ArgonneCPAC/diffmah, and includes Jupyter notebooks
providing demonstrations of typical use cases. A parallelized script in the
diffmah repository can be used to fit the assembly histories of individual
simulated halos with the Diffmah parameters. The diffmah code provides a
differentiable description of $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm
peak}/{\rm d}t|M_{0})$ for both gravity-only simulations and IllustrisTNG, and
also includes a convenience function that can be used to generate Monte Carlo
realizations of cosmologically realistic populations of halo assembly
histories. Precomputed fits for hundreds of thousands of halos in the BPL,
MDPL2, and IllustrisTNG simulations are available at
https://portal.nersc.gov/project/hacc/aphearin/diffmah_data/.

## Acknowledgements

This manuscript has been revised since the version originally published by the
Open Journal of Astrophysics (OJA). In this revised version, a change to the
nomenclature reflects an explicit distinction between the Diffmah model of
individual halo growth, and the DiffmahPop model of the growth of halo
populations.

We thank Raul Angulo, Andrew Benson, Joe DeRose, Alexie Leauthaud, and Daisuke
Nagai for useful discussions. Thanks also to Michael Buehlmann and Aurora
Cossairt for feedback on an early release of our source code, and Benedikt
Diemer for kindly providing downsampled particles for IllustrisTNG. APH thanks
Lauren Mastro for keeping the WWOZ Sunday Gospel Show going throughout the
pandemic.

We thank the developers of NumPy (Van Der Walt et al., 2011), SciPy (Jones et
al., 2001-2016), Jupyter (Ragan-Kelley et al., 2014), IPython (Pérez &
Granger, 2007), scikit-learn (Pedregosa et al., 2011), JAX (Bradbury et al.,
2018), and Matplotlib (Hunter, 2007) for their extremely useful free software.
While writing this paper we made extensive use of the Astrophysics Data
Service (ADS) and arXiv preprint repository. The authors gratefully
acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and
the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu)
for funding the MultiDark simulation project by providing computing time on
the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ,
www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi
project of the University of California High-Performance AstroComputing Center
(UC-HiPACC) and were run at the NASA Ames Research Center.

Work done at Argonne National Laboratory was supported under the DOE contract
DE-AC02-06CH11357. We gratefully acknowledge use of the Bebop cluster in the
Laboratory Computing Resource Center at Argonne National Laboratory.

## References

  * Alarcon et al. (2021) Alarcon, A., Hearin, A., Becker, M., & Chaves-Montero, J. 2021, in prep
  * Angulo & White (2010) Angulo, R. E., & White, S. D. M. 2010, MNRAS, 405, 143, doi: 10.1111/j.1365-2966.2010.16459.x
  * Aricò et al. (2020) Aricò, G., Angulo, R. E., Hernández-Monteagudo, C., et al. 2020, MNRAS, 495, 4800, doi: 10.1093/mnras/staa1478
  * Avila-Reese et al. (1998) Avila-Reese, V., Firmani, C., & Hernández, X. 1998, ApJ, 505, 37, doi: 10.1086/306136
  * Baugh et al. (1996) Baugh, C. M., Cole, S., & Frenk, C. S. 1996, MNRAS, 283, 1361, doi: 10.1093/mnras/283.4.1361
  * Baydin et al. (2015) Baydin, A. G., Pearlmutter, B. A., Radul, A. A., & Siskind, J. M. 2015, JMLR, 18, 1
  * Becker (2015) Becker, M. R. 2015, arXiv:1507.03605, arXiv:1507.03605.  https://arxiv.org/abs/1507.03605
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Behroozi & Silk (2015) Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32, doi: 10.1088/0004-637X/799/1/32
  * Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013a, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013b, ApJ, 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Behroozi et al. (2013c) Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013c, ApJ, 763, 18, doi: 10.1088/0004-637X/763/1/18
  * Benson (2012) Benson, A. J. 2012, New A, 17, 175, doi: 10.1016/j.newast.2011.07.004
  * Blumenthal et al. (1984) Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517, doi: 10.1038/311517a0
  * Bocquet et al. (2020) Bocquet, S., Heitmann, K., Habib, S., et al. 2020, ApJ, 901, 5, doi: 10.3847/1538-4357/abac5c
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440, doi: 10.1086/170520
  * Bower (1991) Bower, R. G. 1991, MNRAS, 248, 332, doi: 10.1093/mnras/248.2.332
  * Bower et al. (2006) Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645, doi: 10.1111/j.1365-2966.2006.10519.x
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018
  * Broyden (1970) Broyden, C. G. 1970, IMA Journal of Applied Mathematics, 6, 76, doi: 10.1093/imamat/6.1.76
  * Bullock et al. (2001) Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559, doi: 10.1046/j.1365-8711.2001.04068.x
  * Chaves-Montero & Hearin (2020) Chaves-Montero, J., & Hearin, A. 2020, MNRAS, 495, 2088, doi: 10.1093/mnras/staa1230
  * Chaves-Montero & Hearin (2021) —. 2021, MNRAS, doi: 10.1093/mnras/stab1831
  * Chen et al. (2020a) Chen, Y., Mo, H. J., Li, C., et al. 2020a, ApJ, 899, 81, doi: 10.3847/1538-4357/aba597
  * Chen et al. (2020b) —. 2020b, ApJ, 899, 81, doi: 10.3847/1538-4357/aba597
  * Child et al. (2018) Child, H. L., Habib, S., Heitmann, K., et al. 2018, ApJ, 859, 55, doi: 10.3847/1538-4357/aabf95
  * Cole et al. (2008) Cole, S., Helly, J., Frenk, C. S., & Parkinson, H. 2008, MNRAS, 383, 546, doi: 10.1111/j.1365-2966.2007.12516.x
  * Contreras et al. (2020) Contreras, S., Angulo, R. E., Zennaro, M., Aricò, G., & Pellejero-Ibañez, M. 2020, MNRAS, 499, 4905, doi: 10.1093/mnras/staa3117
  * Dark Energy Survey Collaboration (2016) Dark Energy Survey Collaboration. 2016, MNRAS, 460, 1270, doi: 10.1093/mnras/stw641
  * DeRose et al. (2019) DeRose, J., Wechsler, R. H., Tinker, J. L., et al. 2019, ApJ, 875, 69, doi: 10.3847/1538-4357/ab1085
  * Diemer & Kravtsov (2014) Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1, doi: 10.1088/0004-637X/789/1/1
  * Diemer & Kravtsov (2015) —. 2015, ApJ, 799, 108, doi: 10.1088/0004-637X/799/1/108
  * Fakhouri et al. (2010) Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267, doi: 10.1111/j.1365-2966.2010.16859.x
  * Fletcher (1970) Fletcher, R. 1970, The Computer Journal, 13, 317, doi: 10.1093/comjnl/13.3.317
  * Gao et al. (2005) Gao, L., Springel, V., & White, S. D. M. 2005, MNRAS, 363, L66, doi: 10.1111/j.1745-3933.2005.00084.x
  * Gao et al. (2004) Gao, L., White, S. D. M., Jenkins, A., Stoehr, F., & Springel, V. 2004, MNRAS, 355, 819, doi: 10.1111/j.1365-2966.2004.08360.x
  * Genel et al. (2010) Genel, S., Bouché, N., Naab, T., Sternberg, A., & Genzel, R. 2010, ApJ, 719, 229, doi: 10.1088/0004-637X/719/1/229
  * Giocoli et al. (2012) Giocoli, C., Tormen, G., & Sheth, R. K. 2012, MNRAS, 422, 185, doi: 10.1111/j.1365-2966.2012.20594.x
  * Giocoli et al. (2010) Giocoli, C., Tormen, G., Sheth, R. K., & van den Bosch, F. C. 2010, MNRAS, 404, 502, doi: 10.1111/j.1365-2966.2010.16311.x
  * Goldfarb (1970) Goldfarb, D. 1970, Math. Comp., 23, doi: 10.1090/S0025-5718-1970-0258249-6
  * Green et al. (2020) Green, S. B., Aung, H., Nagai, D., & van den Bosch, F. C. 2020, MNRAS, 496, 2743, doi: 10.1093/mnras/staa1712
  * Hearin et al. (2021) Hearin, A., Ramachandra, N., Becker, M., & DeRose, J. 2021, in prep
  * Hearin et al. (2017) Hearin, A. P., Campbell, D., Tollerud, E., et al. 2017, AJ, 154, 190, doi: 10.3847/1538-3881/aa859f
  * Heitmann et al. (2016) Heitmann, K., Bingham, D., Lawrence, E., et al. 2016, ApJ, 820, 108, doi: 10.3847/0004-637X/820/2/108
  * Henriques et al. (2015) Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663, doi: 10.1093/mnras/stv705
  * Higham (2009) Higham, N. 2009, Wiley Interdisciplinary Reviews: Computational Statistics, 1, 251, doi: 10.1002/wics.18
  * Hoffman & Gelman (2014) Hoffman, M., & Gelman, A. 2014, Journal of Machine Learning Research, 15, 1593
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jenkins et al. (2001) Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372, doi: 10.1046/j.1365-8711.2001.04029.x
  * Jiang et al. (2021) Jiang, F., Dekel, A., Freundlich, J., et al. 2021, MNRAS, 502, 621, doi: 10.1093/mnras/staa4034
  * Jiang & van den Bosch (2014) Jiang, F., & van den Bosch, F. C. 2014, MNRAS, 440, 193, doi: 10.1093/mnras/stu280
  * Jiang & van den Bosch (2017) —. 2017, MNRAS, 472, 657, doi: 10.1093/mnras/stx1979
  * Jones et al. (2001-2016) Jones, E., Oliphant, T., Peterson, P., et al. 2001-2016, http://www.scipy.org
  * Kauffmann et al. (1993) Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201, doi: 10.1093/mnras/264.1.201
  * Kingma & Ba (2014) Kingma, D. P., & Ba, J. 2014, arXiv:1412.6980, arXiv:1412.6980.  https://arxiv.org/abs/1412.6980
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340, doi: 10.1093/mnras/stw248
  * Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102, doi: 10.1088/0004-637X/740/2/102
  * Kravtsov et al. (2004) Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35, doi: 10.1086/420959
  * Kravtsov et al. (1997) Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73, doi: 10.1086/313015
  * Lacey & Cole (1993) Lacey, C., & Cole, S. 1993, MNRAS, 262, 627, doi: 10.1093/mnras/262.3.627
  * Lagos et al. (2018) Lagos, C. d. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573, doi: 10.1093/mnras/sty2440
  * Lau et al. (2021) Lau, E. T., Hearin, A. P., Nagai, D., & Cappelluti, N. 2021, MNRAS, 500, 1029, doi: 10.1093/mnras/staa3313
  * Li et al. (2007) Li, Y., Mo, H. J., van den Bosch, F. C., & Lin, W. P. 2007, MNRAS, 379, 689, doi: 10.1111/j.1365-2966.2007.11942.x
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv:0912.0201, arXiv:0912.0201
  * Ludlow et al. (2013) Ludlow, A. D., Navarro, J. F., Boylan-Kolchin, M., et al. 2013, MNRAS, 432, 1103, doi: 10.1093/mnras/stt526
  * Mansfield & Kravtsov (2020) Mansfield, P., & Kravtsov, A. V. 2020, MNRAS, 493, 4763, doi: 10.1093/mnras/staa430
  * Mao et al. (2018) Mao, Y.-Y., Zentner, A. R., & Wechsler, R. H. 2018, MNRAS, 474, 5143, doi: 10.1093/mnras/stx3111
  * McBride et al. (2009) McBride, J., Fakhouri, O., & Ma, C.-P. 2009, MNRAS, 398, 1858, doi: 10.1111/j.1365-2966.2009.15329.x
  * McClintock et al. (2019) McClintock, T., Rozo, E., Becker, M. R., et al. 2019, ApJ, 872, 53, doi: 10.3847/1538-4357/aaf568
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (United Kingdom: Cambridge University Press)
  * Moster et al. (2018) Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822, doi: 10.1093/mnras/sty655
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624, doi: 10.1093/mnras/stx3040
  * Nishimichi et al. (2019) Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29, doi: 10.3847/1538-4357/ab3719
  * Nusser & Sheth (1999) Nusser, A., & Sheth, R. K. 1999, MNRAS, 303, 685, doi: 10.1046/j.1365-8711.1999.02197.x
  * Ogiya et al. (2021) Ogiya, G., Taylor, J., & Hudson, M. 2021, MNRAS, 503, 1233, doi: 10.1093/mnras/stab361
  * Ogiya et al. (2019) Ogiya, G., van den Bosch, F. C., Hahn, O., et al. 2019, MNRAS, 485, 189, doi: 10.1093/mnras/stz375
  * Parkinson et al. (2008) Parkinson, H., Cole, S., & Helly, J. 2008, MNRAS, 383, 557, doi: 10.1111/j.1365-2966.2007.12517.x
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16, doi: 10.1051/0004-6361/201321591
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425, doi: 10.1086/152650
  * Ragan-Kelley et al. (2014) Ragan-Kelley, M., Perez, F., Granger, B., et al. 2014, in American Geophysical Union Fall Meeting Abstracts, Vol. D7
  * Renneby et al. (2018) Renneby, M., Hilbert, S., & Angulo, R. E. 2018, MNRAS, 479, 1100, doi: 10.1093/mnras/sty1332
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49, doi: 10.1093/mnras/stv264
  * Rodríguez-Puebla et al. (2016) Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893, doi: 10.1093/mnras/stw1705
  * Romberg (1955) Romberg, W. 1955, Norske Vid. Selsk. Forh., Trondheim, 28, 30
  * Shanno (1970) Shanno, D. 1970, Math. Comp., 647, doi: 10.1090/S0025-5718-1970-0274029-X
  * Shi et al. (2018) Shi, J., Wang, H., Mo, H. J., et al. 2018, ApJ, 857, 127, doi: 10.3847/1538-4357/aab775
  * Somerville et al. (2008) Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481, doi: 10.1111/j.1365-2966.2008.13805.x
  * Somerville & Kolatt (1999) Somerville, R. S., & Kolatt, T. S. 1999, MNRAS, 305, 1, doi: 10.1046/j.1365-8711.1999.02154.x
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105, doi: 10.1111/j.1365-2966.2005.09655.x
  * Springel (2010) —. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726, doi: 10.1046/j.1365-8711.2001.04912.x
  * Tasitsiomi et al. (2004) Tasitsiomi, A., Kravtsov, A. V., Gottlöber, S., & Klypin, A. A. 2004, ApJ, 607, 125, doi: 10.1086/383219
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709, doi: 10.1086/591439
  * van den Bosch (2002) van den Bosch, F. C. 2002, MNRAS, 331, 98, doi: 10.1046/j.1365-8711.2002.05171.x
  * van den Bosch et al. (2014) van den Bosch, F. C., Jiang, F., Hearin, A., et al. 2014, MNRAS, 445, 1713, doi: 10.1093/mnras/stu1872
  * van den Bosch & Ogiya (2018) van den Bosch, F. C., & Ogiya, G. 2018, MNRAS, 475, 4066, doi: 10.1093/mnras/sty084
  * van den Bosch et al. (2018) van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043, doi: 10.1093/mnras/stx2956
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, ArXiv:1102.1523
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2020, arXiv:2010.00619, arXiv:2010.00619.  https://arxiv.org/abs/2010.00619
  * Wang et al. (2020) Wang, K., Mao, Y.-Y., Zentner, A. R., et al. 2020, MNRAS, 498, 4450, doi: 10.1093/mnras/staa2733
  * Watson et al. (2015) Watson, D. F., Hearin, A. P., Berlind, A. A., et al. 2015, MNRAS, 446, 651, doi: 10.1093/mnras/stu2065
  * Wechsler et al. (2002) Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52, doi: 10.1086/338765
  * Wechsler et al. (2006) Wechsler, R. H., Zentner, A. R., Bullock, J. S., Kravtsov, A. V., & Allgood, B. 2006, ApJ, 652, 71, doi: 10.1086/507120
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * Wetzel & Nagai (2015) Wetzel, A. R., & Nagai, D. 2015, ApJ, 808, 40, doi: 10.1088/0004-637X/808/1/40
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341
  * Yang et al. (2011) Yang, X., Mo, H. J., Zhang, Y., & van den Bosch, F. C. 2011, ApJ, 741, 13, doi: 10.1088/0004-637X/741/1/13
  * Zentner (2007) Zentner, A. R. 2007, International Journal of Modern Physics D, 16, 763, doi: 10.1142/S0218271807010511
  * Zentner et al. (2005) Zentner, A. R., Berlind, A. A., Bullock, J. S., Kravtsov, A. V., & Wechsler, R. H. 2005, ApJ, 624, 505, doi: 10.1086/428898
  * Zhao et al. (2009) Zhao, D. H., Jing, Y. P., Mo, H. J., & Börner, G. 2009, ApJ, 707, 354, doi: 10.1088/0004-637X/707/1/354
  * Zhao et al. (2003a) Zhao, D. H., Mo, H. J., Jing, Y. P., & Börner, G. 2003a, MNRAS, 339, 12, doi: 10.1046/j.1365-8711.2003.06135.x
  * Zhao et al. (2003b) —. 2003b, MNRAS, 339, 12, doi: 10.1046/j.1365-8711.2003.06135.x

##  Appendix A Fitting the Assembly History of Individual Halos with Diffmah

In this appendix, we give a detailed description of Diffmah, our parametric
model for the mass assembly of individual halos. As outlined in §3, the
Diffmah model for individual halo growth is based on Eq. 1, which describes a
power-law scaling between peak halo mass and cosmic time, $M_{\rm
peak}=M_{0}(t/t_{0})^{\alpha(t)},$ where the power-law index $\alpha(t)$
smoothly rolls between asymptotic values $\alpha_{\rm early}$ and $\alpha_{\rm
late}$ according to the sigmoid function in Eq. 2. Thus there are a total of
six numbers that fully characterize our model for individual halo growth: the
normalization, $M_{0},$ the present-day cosmic time, $t_{0},$ and the four
sigmoid parameters, $\alpha_{\rm early},\alpha_{\rm late},k,$ and $\tau_{\rm
c}.$ For all halos in each simulation, we hold $t_{0}$ fixed to the present-
day age of the universe in the simulation, and $M_{0}$ fixed to the value
$M_{\rm peak}(t_{0})$ in the simulation. All results in the paper also hold
$k$ fixed to 3.5; to determine this particular value, we ran the optimization
algorithm described below while allowing $k$ to be a free parameter, and then
observed no appreciable changes in the quality of the fits when holding $k$
fixed to any value in its typical best-fitting range, $2\lesssim k\lesssim 5.$
Thus the Diffmah model for individual halo growth has a total of 3 free
parameters: $\alpha_{\rm early},\alpha_{\rm late},$ and $\tau_{\rm c}.$

When evaluating the sigmoid behavior of the power-law index, we implement this
scaling relation in logarithmic time, $x\equiv\log_{10}t,$ rather than linear
time. Thus in practice, we calculate the power-law index $\alpha(t)$ via:

$\displaystyle\alpha(t)=\alpha(10^{x})=\alpha_{\rm early}+\frac{\alpha_{\rm
late}-\alpha_{\rm early}}{1+\exp(-k(x-x_{0}))}.$ (3)

In order to identify an optimal set of parameters for each halo, we searched
our three-dimensional model parameter space, ${\theta},$ for the combination
of $\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}$ that minimizes the
quantity $\mathcal{L}_{\rm MSE},$ defined as

$\displaystyle\mathcal{L}_{\rm
MSE}\equiv\frac{1}{N}\sum_{i}\left(w(t_{i};{\theta})-v(t_{i})\right)^{2},$ (4)

where $w$ and $v$ are the base-10 logarithm of the predicted and simulated
values of $M_{\rm peak}$ evaluated at a set of $N$ control points, $t_{i}.$
For the control points of each halo, we use the collection of snapshots in the
halo’s merger tree after a time $t_{\rm min},$ defined as

$\displaystyle t_{\rm min}\equiv{\rm max}(t_{\rm cut},\ t_{\rm
thresh},t_{\rm\delta m}),$ (5)

where we set $t_{\rm cut}=1$ Gyr, and the quantity $t_{\rm thresh}$ is the
most-recent simulated snapshot where the halo mass falls below a simulation-
dependent threshold mass, $M_{\rm thresh};$ for the MDPL2 simulation, we use
$M_{\rm thresh}=10^{11.25}{\rm M}_{\odot},$ and for the Bolshoi-Planck and TNG
simulations we use $M_{\rm thresh}=10^{10}{\rm M}_{\odot}.$

Figure A9: Gradients of the Diffmah parameters. The figure shows logarithmic
derivatives of $M_{\rm peak}(t)$ with respect to the three parameters of our
model for individual halo growth. This figure is similar to Fig. 1, only here
we show the behavior of the autodiff-computed gradients that we use to
determine the best-fitting parameters describing the assembly of each
simulated halo.

In Eq. 5, the quantity $t_{\rm\delta m}$ refers to the most-recent snapshot
where the halo mass falls below some fraction $f_{\rm\delta M}$ of its
present-day mass; for all simulations, we set $f_{\rm\delta M}=10^{2.5}.$ In
each of the three simulations we studied, for a surprisingly large number of
halos, we find that at early times, $t\lesssim 3$Gyr, the behavior of $M_{\rm
peak}(t)$ suddenly drops like a stone by an order of magnitude or more from
one snapshot to the next, even for halos that should nominally be resolved by
hundreds of particles. We found through experimentation that the inclusion of
the term $t_{\rm\delta m}$ in Eq. 5 is reasonably effective at limiting the
influence of this phenomenon on the best-fitting parameters describing the
halo MAH. We suspect that this behavior is a reflection of the difficulty of
halo-finding when the high-redshift density field resembles a sea of shallow
peaks, however, the root cause of this phenomenon is unclear, and other
origins of this phenomenon are certainly plausible. For example, recent work
on structure formation has uncovered surprisingly stringent resolution demands
of N-body simulations for science applications that require robust tracking of
the evolution of halo substructure (van den Bosch et al., 2018; van den Bosch
& Ogiya, 2018; Ogiya et al., 2019). It could be the case that analogous
results apply to the simulation requirements associated with reliable tracking
of halo mass accretion history in the first Gyr of cosmic time; alternatively,
halo assembly at high redshift may simply require an even more flexible
fitting function than the smoothly rolling power-law model adopted here. As
discussed further in §5, the most reliable way to address this and related
issues is through a dedicated resolution-requirement study that uses
homogeneously-processed, nested simulation boxes spanning a wide range of
resolution. Since such a simulation suite is not currently publicly available,
and since our own science applications are focused on the redshift range
$z\lesssim 5$ that is most relevant to present-day and near-future
cosmological surveys, we have opted to relegate further study of this issue to
future work.

To minimize $\mathcal{L}_{\rm MSE}$ for each halo, we use the JAX
implementation of the Adam algorithm (Kingma & Ba, 2014), which is a gradient
descent technique with an adaptive learning rate. Our use of Adam requires
calculating the gradients $\partial\mathcal{L}_{\rm MSE}/\partial\theta,$
which is facilitated by our JAX implementation, allowing us to efficiently
compute these gradients to machine precision without reliance upon finite
differencing methods. We show example gradients of our model predictions for
$M_{\rm peak}(t)$ in Figure A9. To tune the performance of the optimization
beyond the default settings recommended in Kingma & Ba (2014), we use 2
successive burn-in cycles with a relatively large step-size parameter of 0.25
for $\sim 50$ updates, followed by a final cycle for an additional $\sim 200$
updates using a step-size parameter of 0.01.

###  A.1 Imposing physical constraints on individual halo growth

When minimizing $\mathcal{L}_{\rm MSE}$, the parameters we actually vary in
the gradient descent are $\beta_{\rm e},\,\beta_{\ell},$ and $x_{0},$ where
the relationships between these parameters and the quantities that directly
appear in Eq. 1, $\alpha_{\rm early},\,\alpha_{\rm late},$ and $\tau_{\rm c},$
are defined as follows:

$\displaystyle x_{0}$ $\displaystyle=$ $\displaystyle\log_{10}\tau_{\rm c}$
$\displaystyle\alpha_{\rm late}$ $\displaystyle=$
$\displaystyle\mathcal{S}(\beta_{\ell})$ (6) $\displaystyle\alpha_{\rm early}$
$\displaystyle=$ $\displaystyle\alpha_{\rm late}+\mathcal{S}(\beta_{\rm e})$
$\displaystyle\mathcal{S}(z)$ $\displaystyle\equiv$
$\displaystyle\ln(1+\exp(z))$

The function $\mathcal{S}(z)$ is the softplus function, which is strictly
positive across the real line, exhibits approximately linear behavior for
$s\gtrsim 1,$ and asymptotically approaches zero as $s\rightarrow-\infty.$ For
any positively-valued $z,$ the inverse of the softplus function is well-
defined and computed as $\mathcal{S}^{-1}(z)\equiv\ln(\exp(z)-1),$ which we
have written as $U(z)$ in the main body of the paper and plotted on the axes
of Figure 5, as well as on the axes of the figures in Appendix C. Figure A10
gives a visual illustration of the behavior of the softplus function we use to
enforce these constraints for the case of $\alpha_{\rm late}.$

Through our use of these variable transformations, we ensure that the best-
fitting parameters determined by our application of gradient descent will
always respect $0<\alpha_{\rm late}<\alpha_{\rm early}.$ Mathematically, this
guarantees that in the best-fitting approximation to each individual halo, for
all $t>0$ the quantity $M_{\rm peak}(t)$ will increase monotonically, and
${\rm d}M_{\rm peak}/{\rm d}t$ will strictly slow down as $M_{\rm peak}(t)$
eventually attains the value $M_{0}$ at the time $t_{0}.$ These mathematical
constraints embody the physical expectation that dark matter halo growth
begins with a fast-accretion phase at early times, and transitions to a slow-
accretion phase at late times. One can imagine employing a more flexible
fitting function such as one based directly on a machine-learning algorithm
that would not hard-wire this expectation into the approximation; such an
approach has potential to capture short-term fluctuations that our model does
not. We note, however, that even when these constraints are imposed on the
fits as described above, dark matter halo growth can be approximated with a
typical accuracy of 0.1dex or better for all $t\gtrsim 1$ Gyr.

Figure A10: Illustration of physical constraints on $\alpha_{\rm late}.$ The
figure shows the behavior of the softplus function that controls the
relationship between $\alpha_{\rm late}$ and $\beta_{\ell},$ the actual
variable that we vary in the gradient descent algorithm used to fit the
assembly history of each individual halo. We use a similar technique to ensure
$0<\alpha_{\rm late}<\alpha_{\rm early},$ and a simple logarithmic
transformation to enforce $\tau_{\rm c}>0.$ These transformations guarantee
that the growth history of each model halo will always increase monotonically
and transition from a fast-accretion phase to a slow-accretion phase, even for
outlier halos with unusual merger trees. See Eq. A.1 and associated discussion
for details.

In carrying out this investigation, we studied a wide range of alternative
implementations in which alternative variables were used in the optimization.
For example, rather than using $\beta_{\rm e}$ and $\beta_{\ell},$ we have
also explored the use of logarithmic variables, and an alternative formulation
in which the varied parameters were restricted to finite bounds via a sigmoid
function. The specific choice for which parameters are varied in the
optimization amounts, in effect, to a choice for a prior distribution on the
parameters. We found through considerable experimentation that this choice can
have a significant impact on the distribution of best-fitting parameters, even
though the quality of the individual fits themselves was typically unaffected
by this choice. Since one of our primary goals was to build a simple model for
the distribution of best-fitting parameters, then the choice for the effective
priors warranted special attention. Ultimately, our decision on the
parameterization was driven by two considerations: first, that the
differentiable approximation to each simulated halo respects the physical
constraints defined above, and second, that the distribution of best-fitting
parameters admits a reasonably simple empirical description.

Amongst all the variations we explored for the effective priors, as well as
amongst several variations in the actual functional form used to fit the
growth of individual halos, we found a clear bimodality in the halo
population. The manifestation of this bimodality of course changed from
variation to variation, but in all cases, we consistently found two separate
sub-populations within the distribution of best-fitting parameters, roughly
corresponding to early-forming and late-forming regions of parameter space.
While these two sub-populations are easily identified in terms of their
clustering within the distribution of best-fitting parameters, the
distribution of the actual values of the formation time of the halos, $t_{\rm
form},$ shows heavy overlap between the two groupings. There is also no strong
dependence of the NFW concentration, $c,$ on sub-population membership, even
though $c$ is strongly correlated with $t_{\rm form}.$ Moreover, when
splitting a sample of halos of the same $M_{0}$ according to sub-population
membership, the two-point clustering of two samples hardly differs, even for
halo masses where the actual $t_{\rm form}$-dependence of $\xi_{\rm
m\delta}(r)$ is strong. Thus even though this bimodality considerably
simplified our task of modeling the distribution of best-fitting parameters
(as described in detail in Appendix B), the sub-population membership of a
halo may not have special physical significance.

##  Appendix B Modeling the Assembly of Halo Populations with DiffmahPop

In this appendix, we describe our model for the diversity of assembly
histories of halo populations. The goal of our halo population model is to
capture $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm
d}t|M_{0}).$ Our approach to this problem is to use our model for individual
halo growth as a basis, and to build a model of the statistical connection
between present-day halo mass and the distribution of parameters describing
individual halo MAHs, $P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm
c}|M_{0}).$

As outlined in §4, we model the distribution $P(\alpha_{\rm early},\alpha_{\rm
late},\tau_{\rm c}|M_{0})$ as a two-population normal distribution of 3
parameters. We model this distribution in terms of the same parameter
transformations used in Appendix A to define the variables $\beta_{\ell}\equiv
U(\alpha_{\rm late}),$ $\beta_{\rm e}\equiv U(\alpha_{\rm early}-\alpha_{\rm
late}),$ and $x_{0}\equiv\log_{10}\tau_{\rm c},$ where
$U(s)\equiv\ln(\exp(s)-1)$ is the inverse of the softplus function (see Eq.
A.1). Thus each sub-population is described by an independently-defined mean,
$\mu_{\rm i}\equiv(\beta_{\rm e},\beta_{\ell},x_{0}),$ and covariance matrix,
$\Sigma_{\rm i},$ defined as

$\displaystyle\Sigma_{\rm i}$ $\displaystyle\equiv$
$\displaystyle\begin{bmatrix}\sigma(\beta_{\rm e},\beta_{\rm
e})&\sigma(\beta_{\rm e},\beta_{\ell})&\sigma(\beta_{\rm e},x_{0})\\\
\sigma(\beta_{\ell},\beta_{\rm
e})&\sigma(\beta_{\ell},\beta_{\ell})&\sigma(\beta_{\ell},x_{0})\\\
\sigma(x_{0},\beta_{\rm
e})&\sigma(x_{0},\beta_{\ell})&\sigma(x_{0},x_{0})\end{bmatrix}$ (7)

By defining our halo population model in terms of the transformed variables
$\beta_{\rm e},\,\beta_{\ell},$ and $x_{0},$ we ensure that all values of
$\alpha_{\rm early},\,\alpha_{\rm late},$ and $\tau_{\rm c}$ in the support of
$P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}|M_{0})$ will respect the
physical constraints described in §A.1, even in the tails of the distribution.

The relative weight of the two populations, $F_{\rm late},$ has shallow,
monotonic mass-dependence, so that $F_{\rm late}=F_{\rm late}(M_{0}).$ Each of
the three components of $\mu_{\rm i},$ and each of the six components of
$\Sigma_{\rm i},$ also have a shallow, monotonic mass-dependence, so that
$\mu_{\rm i}=\mu_{\rm i}(M_{0}),$ and $\Sigma_{\rm i}=\Sigma_{\rm i}(M_{0}).$
For the remainder of this appendix, we will describe our models for $F_{\rm
late}(M_{0}),$ $\mu_{\rm i}(M_{0})$ and $\Sigma_{\rm i}(M_{0})$ in turn.

We model the mass-dependence of $F_{\rm late}(M_{0})$ as a power-law function
of present-day mass with rolling index. We capture this behavior using the
same sigmoid functional form defined in Eq. 2, so that $F_{\rm late}$ has
sigmoid-type dependence upon $\log_{10}M_{0}.$ We calibrate the specific
parameters of the sigmoid using the techniques described in detail in Appendix
C.

We model the mass-dependence of each component of $\mu_{\rm i}$ in the same
way we modeled $F_{\rm late},$ i.e., as a power-law function of $M_{0}$ with
rolling index, implemented based on a sigmoid function of $\log_{10}M_{0}.$
Again we relegate discussion of our determination of the specific parameters
of these sigmoid functions to Appendix C.

In order to model the components of each $\Sigma_{\rm i},$ we make use of the
Cholesky matrix decomposition. Briefly, for any real-valued, symmetric,
positive-definite matrix, $\Sigma,$ there exists a unique, real-valued, lower-
triangular matrix, $L,$ that defines the Cholesky decomposition, written as
$\Sigma=L\cdot L^{\rm T}.$ The diagonal entries of $L$ are strictly positive
with a product equal to the determinant of $\Sigma.$ We use JAX to calculate
$L$ so that our computations will be differentiable with autodiff, but many
modern linear algebra libraries have efficient implementations of the Cholesky
decomposition (for a modern review, see Higham, 2009).

In modeling each of the two $3\times 3$ matrices, $\Sigma_{\rm i},$ the
quantities we parameterize are the 6 entries of the associated Cholesky
matrix:

$\displaystyle L_{\rm i}$ $\displaystyle\equiv$
$\displaystyle\begin{bmatrix}a_{\rm i}&0&0\\\ d_{\rm i}&b_{\rm i}&0\\\ e_{\rm
i}&f_{\rm i}&c_{\rm i}\end{bmatrix}$ (8)

In fitting our model for $P(M_{\rm peak}(t)|M_{0})$ as described in Appendix
C, we formulate our parameterization in terms of $\log_{10}a_{\rm i},$
$\log_{10}b_{\rm i},$ and $\log_{10}c_{\rm i}$ to ensure that $\Sigma=L\cdot
L^{\rm T}$ will always be positive definite, but the parameters $d_{\rm
i},e_{\rm i},$ and $f_{\rm i}$ can take on any value on the real line and
still produce a valid covariance matrix, and so we use linear variables for
our parameterization of the off-diagonal entries of $L.$

##  Appendix C Fitting the Parameters of DiffmahPop

Our goal in §4 is to identify a point in the parameter space of our model that
can be used to generate a realistic distribution of individual, smooth
trajectories of halo assembly, and at the same time results in an accurate
reproduction of $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm
d}t|M_{0}).$ To achieve these goals, our model for $P(M_{\rm peak}(t)|M_{0})$
has numerous free parameters that required tuning in order to achieve the
level of agreement shown in Figures 7 & 8. In C.1, we review the general
techniques we used for making differentiable predictions for the mean and
variance of $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm
d}t|M_{0}),$ and in C.2 we describe how we optimized the parameters of our
best-fitting models.

###  C.1 Differentiable predictions for halo population assembly

The first moment of our model prediction for $P(M_{\rm peak}(t)|M_{0})$ is
given by

$\displaystyle\langle M_{\rm peak}(t)|M_{0}\rangle$ $\displaystyle\equiv$
$\displaystyle\int\int\int{\rm d}\tau_{\rm c}{\rm d}\alpha_{\rm early}{\rm
d}\alpha_{\rm late}M_{\rm peak}(t)$ (9) $\displaystyle\times$ $\displaystyle
P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}|M_{0}),$

where the PDF in the integrand is given by555Note that in Eq. 9 we have
written $M_{\rm peak}(t|\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c})$
simply as $M_{\rm peak}(t),$ and in Eq. 10 we have suppressed the conditioning
upon $M_{0},$ e.g., $F_{\rm late}(M_{0})$ is written simply as $F_{\rm late}.$

$\displaystyle P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c})$
$\displaystyle=$ $\displaystyle(1-F_{\rm late})\cdot P_{\rm early}(\alpha_{\rm
early},\alpha_{\rm late},\tau_{\rm c})$ (10) $\displaystyle+$ $\displaystyle
F_{\rm late}\cdot P_{\rm late}(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm
c}),$

where $P_{\rm early}$ and $P_{\rm late}$ are the two, separately-defined
three-dimensional Gaussian distributions described in Appendix B. The
calculation of the second moment of $P(M_{\rm peak}(t)|M_{0})$ is similar to
Eq. 9, only rather than $M_{\rm peak}(t)$ appearing in the integrand, the PDF-
weighted quantity is instead $\left[M_{\rm peak}(t)-\langle M_{\rm
peak}(t)\rangle\right]^{2}.$

There are a number of approaches one could take to computing these nested
integrals. First, we note that our model for $P(\alpha_{\rm early},\alpha_{\rm
late},\tau_{\rm c}|M_{0})$ is completely analytic and smooth, without any
complex oscillatory behavior or sharp turns in the parameter space, and so
commonly-used integration routines (e.g., Romberg, 1955) would have no trouble
giving a high-accuracy result. Alternatively, each of the component
ingredients of our model for the PDF in the integrand of Eq. 9 have well-
established numerical algorithms for generating stochastic realizations of the
distributions, and so it would be straightforward to use a Monte Carlo
approach to compute these integrals via moments of randomly generated samples.
Either calculation could also be conveniently conducted in python using
implementations in the scipy library. However, we have instead opted for a
third method based on weighted sampling that simplifies gradient calculations
with autodiff.

Our autodiff-friendly method for calculating the integrals in Eq. 9 can be
easily understood in terms of a simple toy example of a generic PDF
convolution:

$\displaystyle y=\int_{-\infty}^{\infty}{\rm d}xP(x)f(x)\approx\sum_{x_{i}\in
D}\delta x_{i}P(x_{i})f(x_{i}),$ (11)

where we will think of the integrand as some smooth target function, $f(x),$
weighted by a probability distribution, $P(x),$ and where in the second
equation we have replaced the continuous integration with a finite sampling of
points spanning the domain of support, $D.$ In other words, to calculate the
integral in Eq. 11, rather than using an iterative algorithm such as Romberg
(1955) that is challenging to implement in an autodiff library, instead one
can use simple Gaussian integration with a sufficiently dense sampling grid.
This toy problem is of the same form as the problem at hand in the evaluation
of Eq. 9 for cases where $P(x)=P(x|\theta)$ and $f(x)=f(x|\theta),$ only with
this formulation, we can readily see how simple it is to compute gradients of
$y$ with respect to model parameters, $\theta:$

$\displaystyle\partial y/\partial\theta=\sum_{x_{i}\in D}\delta
x_{i}\frac{\partial}{\partial\theta}P(x_{i}|\theta)f(x_{i}|\theta),$ (12)

By implementing the model components described in Appendix B in an autodiff
library such as JAX, evaluating expressions such as 12 is entirely
straightforward and does not require laboriously working out the analytical
result for each contribution to the summations. One need only ensure that the
integration domain $D$ is sufficient to span the range of non-negligible
support for $P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}|M_{0}).$ We
found this flexibility especially helpful during development phases of our
model, when the need to conduct a range of modeling experiments required many
variations of the underlying functional forms. In the following section, we
describe how we carried out these calculations for the particular case of the
PDF convolutions that arise from optimizing our model predictions for the
first and second moments of $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm
peak}/{\rm d}t|M_{0}).$

###  C.2 Optimizing DiffmahPop predictions for halo population assembly

Figure C11: Example comparison of the distribution of Diffmah parameters for
individual halo assembly. For halos of present-day peak mass
$M_{0}=10^{12}{\rm M}_{\odot},$ the figure shows cross-sections of
$P(\alpha_{\rm early},\alpha_{\rm late},\tau_{\rm c}|M_{0}).$ For each
individual BPL halo in the mass bin, we approximate the assembly history of
the halo with the rolling power-law model described in §3, and plot the
distribution of parameters $U(\alpha_{\rm early}),$ $U(\alpha_{\rm
early}-\alpha_{\rm late}),$ and $\log_{10}\tau_{\rm c},$ where
$U(s)\equiv\ln(\exp(s)-1).$ The orange ellipsoids show the earlier-forming
sub-population $(\tau_{\rm c}\lesssim 2$ Gyr) of the two-component Gaussian
distributions identified by scikit-learn, and the purple ellipses show the
later-forming sub-population $(\tau_{\rm c}\gtrsim 2$ Gyr). The colored
ellipses are defined by $\bar{\mu}_{\rm i}$ and $\bar{\Sigma}_{\rm i}$ that we
use to define the initial guess for the parameters regulating our model
prediction for $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm
d}t|M_{0}).$ See text in C.2 for details. Figure C12: Model for the mass-
dependence of the distribution of the Diffmah parameters for individual halo
assembly. The axes in each row of panels are the same as those described in
Figure C11. Different colored ellipses correspond to normal distributions of
halos with different present-day mass as indicated in the legend. The top row
(bottom row) of panels shows distributions for the early-forming (later-
forming) population of halos, where the mean and variance of each population
was determined according to the optimization procedure described in C.2.

In this section, we describe how we optimized the parameters of our model
predictions for $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm
d}t|M_{0}),$ where the target distributions come from the merger trees of
halos in gravity-only N-body simulations. To generate our target data for a
particular bin of present-day halo mass, we begin by selecting host halos in
bin of $M_{\rm peak}(z=0)$ with width $0.1$dex. One issue that arises in
making accurate predictions for the variance of halo histories is that in our
model, the trajectory of each halo has exactly the same final mass, $M_{0},$
by construction, whereas the bin of simulated halos will have a source of
unwanted variance due to the finite width of the bin. To account for this
effect, we rescale the history of each simulated halo according to its actual
value of $M_{\rm peak}(z=0)$ in the simulation before computing the target
mean and variance.

As mentioned in §4, when defining the target data for a particular bin of
present-day halo mass, $M_{0},$ we limit the impact of wild outliers in the
merger trees by trimming $3\sigma$ outliers from the halo sample. We conduct
this outlier exclusion procedure as follows. For each snapshot of the
simulation, $t,$ we compute the trimmed mean, $\mu_{\rm p}(t|M_{0}),$ and
trimmed variance, $\sigma_{\rm p}(t|M_{0}),$ using
scipy.stats.mstats.trimmed_mean and scipy.stats.mstats.trimmed_variance,
respectively; within the scipy library, these quantities are computed by rank-
ordering the objects in the sample, computing the percentile of each object,
$p,$ excluding the objects lying outside the range $(p,100-p),$ and computing
the mean and variance from the resulting population. At each snapshot, and for
each bin of halo mass, we use $\mu_{\rm p}$ and $\sigma_{\rm p}$ to compute
the $z$-score,

$z_{\rm p}(t|M_{0})\equiv\frac{\log_{10}M_{\rm peak}(t|M_{0})-\mu_{\rm
p}(t|M_{0})}{\sigma_{\rm p}(t|M_{0})}.$

We use $p=0.01$ to define each snapshot’s trimmed mean and variance, and at
each snapshot over the redshift range $0<z<3,$ we flag any halo with $|z_{\rm
p}(t|M_{0})|>3.$ We exclude any flagged halo from the sample whose mean and
variance define the target data we use to optimize our model for halo
populations; this outlier exclusion procedure typically removes 2-3% of the
halos in each mass bin; by construction, this excludes halos that experience
an unusually large major merger after $z<3;$ the function implementing this
algorithm can be found in the measure_mahs.py module of the publicly available
diffmah source code. To compute the target data defined in this fashion, we
use the BPL simulation for bins of present-day halo mass centered at
$M_{0}\leq 10^{13.5}{\rm M}_{\odot},$ and MDPL2 for larger halo masses,
choosing bin centers spanning the range
$\log_{10}M_{0}\in\left[11.75,14.5\right],$ with 0.25dex separation.

In building our model for $P(M_{\rm peak}(t)|M_{0}),$ our goal is to
faithfully recover the diversity of smooth assembly histories for a population
of halos with the same present-day mass, $M_{0}.$ As our model for individual
halo assembly does not explicitly account for contributions from mergers, for
our target data we use each individual halo’s best-fitting $M_{\rm peak}(t)$
and ${\rm d}M_{\rm peak}/{\rm d}t$ to define the target mean and variance of
$P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm d}t|M_{0}).$ If we
were to instead use the directly simulated histories to define the target one-
point functions, then our model for $P(\alpha_{\rm early},\alpha_{\rm
late},\tau_{\rm c}|M_{0})$ would instead converge to a result with
unrealistically large variance in the smooth trajectories traced by real
halos. Using the directly simulated values of the MAHs as target data would
require building an additional model component describing the variance from
mergers about the smooth trajectories, which is beyond our intended scope.

The predictions of our model are controlled by the parametrized behavior of
$F_{\rm late}(M_{0}),$ $\mu_{\rm i}(M_{0})$ and $\Sigma_{\rm i}(M_{0}),$ and
optimizing these functions requires beginning with an initial guess for the
parameters. In order to determine this initial guess, we directly inspect the
distribution of best-fitting parameters describing the individual trajectories
of the halos in each target bin, $M_{0}.$ We use the Gaussian mixture model
implemented in scikit-learn to decompose the binned halos into two weighted
sub-populations defined by $\beta_{\rm e},\beta_{\ell},$ and $x_{0};$ we thus
use scikit-learn to supply an estimate for the fraction of later-forming halos
in the bin, $\bar{F}_{\rm late}(M_{0}),$ and the mean and covariance of each
sub-population, $\bar{\mu}_{\rm i}(M_{0})$ and $\bar{\Sigma}_{\rm i}(M_{0}),$
respectively. We proceed in this fashion and record each result for each bin
of $M_{0}.$ Figure C11 shows an example comparison between initial guess for
the two-component model and the true distribution of parameters for halos with
$M_{0}=10^{12}{\rm M}_{\odot}.$

Using the collection of $\bar{F}_{\rm late}(M_{0}),$ $\bar{\mu}_{\rm
i}(M_{0})$ and $\bar{\Sigma}_{\rm i}(M_{0})$ for each of our target bins of
present-day halo mass, we find that each of these quantities exhibits a
shallow, monotonic mass-dependence that is well-approximated by the
4-parameter sigmoid function defined in Eq. 2. To define our initial guess for
the best-fitting parameters of our halo population model, the parameters of
each sigmoid can simply be hand-tuned to give a reasonable approximation to
each quantity, and in all cases we found acceptable approximations by holding
the transition-mass parameter fixed to $x_{0}=13.5,$ and the transition-speed
parameter fixed to $k=0.5.$

Using the procedure described above to define the initial guess, we use the
L-BFGS-B algorithm implemented in scipy to optimize our model parameters. The
loss function we minimize in this case is analogous to $\mathcal{L}_{\rm MSE}$
used in Appendix A (see Eq. 4), only here we minimize the sum of the
logarithmic differences between the collection of means and variances of
$P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm d}t|M_{0})$
tabulated for each of our target mass bins. In searching for the best-fitting
point, we do not vary all four parameters of each sigmoid function of each
component, but only the two parameters controlling the two asymptotic bounds
on each sigmoid, holding the $k$ and $x_{0}$ parameters fixed to their hand-
tuned values. Moreover, we do not allow the parameters to vary arbitrarily,
but rather we only permit the L-BFGS-B algorithm to search within a relatively
narrow range of $\sim 30\%$ of the initial guess. We refer the reader to our
publicly available source code for further quotidian details on the
implementation of this optimization procedure. Figure C12 displays the results
for the mass-dependence of the two-component normal distribution whose
predictions for $P(M_{\rm peak}(t)|M_{0})$ and $P({\rm d}M_{\rm peak}/{\rm
d}t|M_{0})$ are shown in the main body of the paper. We refer the reader to
the rockstar_pdf_model.py of our source code for the specific best-fitting
values that result from this optimization exercise.

##  Appendix D Assembly Histories of Subhalos, Orphans, and Splashback-
Centrals

In the publicly available Rockstar catalogs described in §2, the column
a_first_infall indicates the scale factor of the first snapshot where the
(sub)halo was contained within the virial radius of some more massive halo.
For a present-day host halo (as defined by upid=-1), the value of
a_first_infall will be smaller than unity if and only if the halo experienced
a flyby event at some point in its past history. For the results in the main
body of the paper, we defined our host halo samples by requiring upid=-1,
without regard for the a_first_infall column. Figure D13 in this appendix
shows the residual errors of our model for individual growth for samples of
flyby host halos defined by upid=-1 and a_first_infall<1, using narrow bins of
present-day peak halo mass as indicated by the in-panel annotation. Figure D14
in this appendix shows the analogous results for subhalos in the publicly
available Rockstar catalogs in which upid$\neq-1.$ Finally, Figure D15 shows
the analogous results for samples of subhalos that were either merged or
disrupted at some point prior to $z=0.$ These subhalos do not appear in the
publicly available Rockstar catalogs, but their histories can be extracted
from the publicly available merger trees using the extract_orphan_info.c
program in the UniverseMachine source
code666https://bitbucket.org/pbehroozi/universemachine/src/main/.

Figure D13: Same as Fig. 3, but for host halos in BPL that experienced a flyby
event in their past history. Figure D14: Same as Fig. 3, but for subhalos in
BPL that survive to $z=0.$ Figure D15: Same as Fig. 3, but for subhalos in BPL
that either merged or were disrupted prior to $z=0.$

##  Appendix E Assembly Histories of Halos in IllustrisTNG

In this appendix, we illustrate the results of calculations described in the
main body of the paper, only here shown for halos in IllustrisTNG. When
analyzing the IllustrisTNG merger trees, we compute $M_{\rm peak}(t)$ using
the total halo mass defined by the sum of dark matter, stars and gas. We
remind the reader that halos and merger trees in IllustrisTNG were identified
by SUBFIND and SUBLINK, respectively, which are entirely different algorithms
than the Rockstar and Consistent Trees codes used on the gravity-only
simulations studied in this paper.

Figure E16 is analogous to Fig. 3 in the main body of the paper, but for the
case of present-day subhalos in IllustrisTNG. Figure E17 demonstrates that the
residual errors in approximating halo growth are uncorrelated with the large-
scale density field in IllustrisTNG. For IllustrisTNG, we have also repeated
the halo population model calibration exercise described in detail in Appendix
C. Making no changes to the form of any model components, and varying the same
set of parameters, we simply optimized the model according to a different set
of target data. Figures E18-E19 show that our model for halo population growth
achieves a comparable level of accuracy for IllustrisTNG as for gravity-only
simulations.

Figure E16: Same as Figure 3, but for subhalos in the IllustrisTNG simulation.
Figure E17: Same as Figure 4, but for the IllustrisTNG simulation. Figure E18:
Same as Figure 6, but for the IllustrisTNG simulation. Figure E19: Same as
Figure 7, but for the IllustrisTNG simulation.