# An Improved and Physically-Motivated Scheme for Matching Galaxies with Dark
Matter Halos

Stephanie Tonnesen  Flatiron Institute, CCA, 162 5th Avenue, New York, NY
10010 USA<EMAIL_ADDRESS>(ST) Jeremiah P. Ostriker  Flatiron Institute,
CCA, 162 5th Avenue, New York, NY 10010 USA  Princeton University Observatory,
Ivy Lane, Princeton, NJ 08544 USA  Columbia Astrophysics Laboratory, 538 W
120th St, New York, NY 10027 USA

(January 2021)

###### Abstract

The simplest scheme for predicting real galaxy properties after performing a
dark matter simulation is to rank order the real systems by stellar mass and
the simulated systems by halo mass and then simply assume monotonicity - that
the more massive halos host the more massive galaxies. This has had some
success, but we study here if a better motivated and more accurate matching
scheme is easily constructed by looking carefully at how well one could
predict the simulated IllustrisTNG galaxy sample from its dark matter
computations. We find that using the dark matter rotation curve peak velocity,
$v_{max}$, for normal galaxies reduces the error of the prediction by 30% (18%
for central galaxies and 60% for satellite systems) - following expectations
from the physics of monolithic collapse. For massive systems with halo mass
$>$ 1012.5 M⊙ hierarchical merger driven formation is the better model and
dark matter halo mass remains the best single metric. Using a new single
variable that combines these effects, $\phi$ $=$ $v_{max}$/$v_{max,12.7}$ \+
Mpeak/(1012.7 M⊙) allows further improvement and reduces the error, as
compared to ranking by dark matter mass at $z=0$ by another 6% from $v_{max}$
ranking. Two parameter fits – including environmental effects produce only
minimal further impact.

##  1 Introduction

There has been a great deal of progress in recent years in our understanding
and simulations of galaxy formation. The initial conditions seem to be well
specified by LCDM cosmological models and their variants, and hydrodynamic
codes now have the capacity and resolution to include many of the physical
processes relevant to galaxy formation and evolution. These include dark
matter gravitational collapse, gas cooling, star formation and the mechanical
and radiative feedback from stars and central, massive black holes (BHs).
There is no widely accepted mechanism for the formation of these massive BHs,
but, assuming their formation as seed BHs, their evolution and effects on the
surrounding galaxies are now reasonably well modeled. Recent summaries of the
status of existing work on galaxy formation and evolution are presented in
Somerville & Davé (2015) and Naab & Ostriker (2017). For more detailed
treatments one can consult the EAGLE (Schaye et al., 2015; Crain et al.,
2015), FIRE (Hopkins et al., 2018a, b), IllustrisTNG (Vogelsberger et al.,
2014; Genel et al., 2014; Weinberger et al., 2017; Pillepich et al., 2018),
MUFASA (Davé et al., 2016), NIHAO (Wang et al., 2015; Blank et al., 2019) and
other simulations.

But there are virtues to considering simpler treatments that do not need to
rely on sub-grid modeling and which are easily adaptable to analyzing large
data sets. Of these, the simplest, perhaps, is the abundance matching scheme
which is based on the fact that, in all variants of the LCDM modeling,
galaxies live in more massive dark matter (DM) halos, quasi-spherical lumps of
dark matter, which grow via gravitational instabilities from very low
amplitude ( 10-5), gaussian perturbations imparted at very early times in a
roughly power law distribution by unknown processes (thought to be related to
inflation) ((Vale & Ostriker, 2004, 2006, 2008, e.g.).

There is good agreement on how to compute the formation of these dark matter
halos with several computational codes now able to make moderately high-
resolution cosmic scale volumes containing accurate distributions of DM halos
having well defined properties. Early analyses by Navarro et al. (1997) showed
that these could be represented to a reasonable approximation by three
numbers, a mass, a virial radius and a core radius, with the ratio of the
latter two numbers represented as the concentration. To zeroth order observed
galaxies can be represented by their stellar masses, or alternatively by their
stellar luminosities.

Thus, the simplest possible scheme for populating a volume in the universe
with galaxies would be to populate it first with DM halos, and then make a
rank ordered list of these with the most massive first. Next, one could take
the same volume from the real universe and rank order the observed galaxies by
mass (or luminosity), and then simply assume that the more massive halos
hosted more massive galaxies, putting in each halo the corresponding galaxy.

This, almost ridiculously simple, scheme was pursued by Vale & Ostriker in
three papers (Vale & Ostriker, 2004, 2006, 2008; Kravtsov et al., 2004;
Tasitsiomi et al., 2004, see also). Using this scheme (or a variation
thereof), one can take different variants of the LCDM model, compute the halo
distribution, populate the halos with galaxies using the simple abundance
matching scheme and then compare to observations. By construction, the
luminosity functions must come out to be correct, but correlation functions
(Conroy et al., 2006; Marín et al., 2008; Guo et al., 2010; Trujillo-Gomez et
al., 2011), pair counts (Berrier et al., 2006), magnitude gap statistics
(Hearin et al., 2013; Ostriker et al., 2019), and galaxy-galaxy lensing
(Hearin & Watson, 2013) can be usefully compared to observations.

Two immediate questions arise. First, is there a better zeroth order scheme
than ranking by mass and matching. Several authors have considered this
question. For example, while the DM in subhalos is quickly stripped once they
are accreted onto halos, stripping of the more centralized galaxy starts later
(Nagai & Kravtsov, 2005), and in fact their optical sizes are observed to grow
with cosmic time (van Dokkum et al., 2010, cf), presumably due to the
accretion of smaller satellite systems. Thus, Conroy et al. (2006) improved
subhalo abundance matching by matching galaxies to halos at the time at which
they are accreted onto a central halo. Even earlier, Kravtsov et al. (2004)
proposed using the maximum circular velocity of (sub)halos, $v_{max}$, which
is more stable than halo mass to stripping. Reddick et al. (2013) tested
abundance matching models using several halo properties, and found that only
vpeak (defined as the peak value of $v_{max}$ over the history of the halo) or
a combination of $v_{max}$ for central galaxies and vpeak for satellite
galaxies is able to reproduce observations of galaxy clustering. Indeed,
Zentner et al. (2014) argues that abundance matching using $v_{max}$ matches
several observed galaxy statistics (Conroy et al., 2006; Hearin et al., 2013;
Hearin & Watson, 2013; Reddick et al., 2013) because halo mass alone does not
determine the halo velocity profile.

Xu & Zheng (2018) confirmed that for the central galaxies in the original
Illustris simulations, M∗ is more tightly correlated with vpeak than with halo
mass. They also find that at fixed vpeak, the correlation between M∗ and other
halo properties is removed. In He (2020), the author uses subhalo abundance
matching and finds that vpeak correlates best with the stellar mass at the
epoch of vpeak in both central and satellite galaxies in EAGLE, Illustris, and
IllustrisTNG. Stellar mass stripping of satellite galaxies results in
increased scatter in the $z=0$ M∗ to vpeak relation. Chaves-Montero et al.
(2016) find that $v_{relax}$, defined as the maximum of the circular velocity
of a dark matter structure while it fulfils a relaxation criterion, as
evaluated along its entire history, correlates most strongly with M∗.

Second, does there exist any first order refinements of the mass-matching
scheme that could be implemented, which would be easy to apply and would
significantly increase its accuracy. In fact, Lehmann et al. (2017) point out
that while ranking by $v_{max}$ is similar to ranking by halo mass, at fixed
halo mass more concentrated halos have higher $v_{max}$ (Klypin et al., 2011).
They therefore use a parameterization from Mao et al. (2015) that includes
both halo mass and concentration:

$V_{\alpha}\equiv V_{vir}(\frac{V_{max}}{V_{vir}})^{\alpha}$ (1)

where $v_{max}$ is the maximal circular velocity of the halo and

$V_{vir}\equiv(\frac{GM_{vir}}{R_{vir}})^{1/2}$ (2)

They find that an $\alpha$ $\sim$ 0.57, with a scatter of 0.17 dex is the best
fit to the SDSS clustering measurements, indicating that both the halo mass
and the maximum circular velocity impact the stellar mass of galaxies.

Identifying the variables causing the scatter in an abundance matching model
also leads to a better understanding of the important physical processes
influencing galaxy formation. For example, Matthee et al. (2017) matched
galaxies in the hydrodynamic EAGLE simulation to those in the dark matter only
simulations, and found that much of the scatter is due to halo concentration.
Specifically, at a constant halo mass, higher concentration is correlated to
higher stellar mass, likely because higher concentrations imply earlier
formation times. However, they do not find a halo property than can explain
the remaining scatter. Martizzi et al. (2020) considered the influence of
formation time and environment on the scatter in the stellar mass to halo mass
relation using IllustrisTNG. Sorting by the current subhalo mass, they find
that the scatter in the relation for central galaxies is correlated more
strongly with formation time, while the scatter for satellite galaxies is
correlated more strongly with environment. On the other hand, Dragomir et al.
(2018) compare the galaxy luminosity functions and the galaxy stellar mass
function at different environmental densities for SDSS observations and a
subhalo abundance matching model applied to the Bolshoi-Planck simulation and
find that the model predictions agree well with observations.

In this paper we will try to both physically motivate and improve the simplest
matching scheme. We will use IllustrisTNG to determine whether complicating
the most simple form of subhalo abundance matching, ie matching stellar mass
to a single halo property, reduces scatter in the assignment of galaxies to
halos. We first verify the halo property that produces the least scatter in
the relation, $v_{max}$, considering the central and satellite populations
separately and selecting different mass ranges even in this initial step. We
then provide a physical motivation for the parameters used in the optimal
matching scheme. Then, similarly to Martizzi et al. (2020), we calculate how
the scatter is reduced when we fold in a second halo feature. However, unlike
these recent works, we test a wide range of possible parameters and possible
combinations thereof.

In Section 2 we describe our sample of galaxies and the variables we measure
for each galaxy. In Section 3.1 we outline a theoretical basis for selecting
the dark matter halo property on which to rank, and in section 3.2 we try out
the simplest one parameter schemes and show that the physically motivated
focus on peak velocity dispersion is best for normal galaxies but that total
halo mass remains best for first brightest, massive, central systems as would
be expected from the presented physical arguments. Section 4 broadens the
treatment to include multiple variables - including environment - and then in
Sections 5 and 6 we present an overall discussion of the results and possible
tests of the scheme. Finally, in Section 7 we summarize our findings and
conclusions.

##  2 Methods

###  2.1 IllustrisTNG

The IllustrisTNG100 (public data release: Nelson et al., 2019)111www.tng-
project.org is part of a suite of cosmological simulations run using the AREPO
moving mesh code (Springel, 2010). TNG100 has a volume of 110.7 Mpc3 and a
mass resolution of $7.5\times 10^{6}\mathrm{M}_{\odot}$ and $1.4\times
10^{6}\mathrm{M}_{\odot}$ for dark matter and baryons, respectively. The TNG
suite implements upgraded subgrid models compared to the Illustris simulation
(Vogelsberger et al., 2014; Genel et al., 2014); specifically, a modified
black hole accretion and feedback model (Weinberger et al., 2017), updated
galactic winds (Pillepich et al., 2018). TNG also includes
magnetohydrodynamics (Pakmor et al., 2011).

###  2.2 Galaxy Selection

We use galaxy populations from the IllustrisTNG 100 simulation described
above. We consider galaxies at the $z=0$ output with dark matter masses of
$\geq$1011 M${}_{\odot}/h$ in the dark-matter only run (DMO) that are matched
in the full hydrodynamical simulation with galaxies whose stellar mass is
greater than 109 M${}_{\odot}/h$.

In detail, we first selected all galaxies in the DMO simulation with a dark
matter mass greater than 1011 M⊙/$h$. We then used the publicly available
matching data to find the corresponding galaxy in the full hydro simulation.
We use all galaxies identified with masses above 5 $\times$ 106 M⊙/$h$, which
clearly includes galaxies that are underresolved in the simulation. However,
at our minimum dark matter halo mass, the lowest stellar mass of any galaxy in
our sample is 7 $\times$ 107 M⊙/$h$. In order to only include well-resolved
galaxies, that are more likely to be in an observational sample, our final
analysis only includes galaxies with stellar masses above 109 M⊙/$h$.

###  2.3 Galaxy Environmental Measures

We use nearby galaxies to measure the local environment. We include all
galaxies with dark matter masses above 109 M${}_{\odot}/h$ within 1 Mpc, 2
Mpc, 5 Mpc, 8 Mpc or 15 Mpc of each galaxy. In order to have a more physical
measure of the local mass density, we summed the total mass of all these
galaxies.

We also were able to separate galaxies into satellites or centrals using the
GroupFirstSub identifier in the IllustrisTNG DMO simulation. This allowed us
to perform our fits for the entire sample and for satellites and centrals
separately, and, as we shall see, the two categories are significantly
different in their properties. This gives us three samples: “all", “centrals"
and “satellites". The fourth sample is labeled as “mix", which is the combined
sample in which satellites and centrals are fit separately.

###  2.4 Concentration

We use three measures of the concentration. First, from Bose et al. (2019) we
use:

$c_{v}\equiv\frac{V_{max}}{H_{o}R_{max}}$ (3)

where $v_{max}$ is the maximum velocity of the simulated rotation curve and
Rmax is the radius at which Vc is maximal. Bose et al. (2019) show that this
is equivalent to the concentration calculated using all the particles in a
halo and assuming an NFW profile (Moliné et al., 2017, see also).

Second, we use the ratio $c_{h}$$\equiv$$v_{max}$/Vhalfmass, where Vhalfmass
is the circular velocity at the half mass radius of the dark matter halo,
calculated as $\sqrt{GM_{halfmass}/R_{halfmass}}$.

Finally, we use the ratio $c_{R}$$\equiv$Rmax/Rhalfmass.

###  2.5 Percent Error

We define the error as:

$Error\equiv\frac{\sum\limits_{N}\mid log(M_{true}/M_{prediction})\mid}{N}$
(4)

so that for small errors our definition is equivalent to 0.43 times the
average fractional error.

##  3 Rank Ordering

In this section we discuss using rank ordering to match halos and galaxies.
Specifically, we first present a straightforward theoretical scaffold for
selecting the halo property best suited for rank ordering. Then, using
IllustrisTNG we confirm our derivation.

###  3.1 A Simple Theoretical Basis for Selecting the Ordering Halo Property

We can begin with an assumption that stellar mass is related to the baryonic
mass scaled to the dark matter mass, corrected by the fraction of matter that
cools and forms stars in the center of the halo:

$M_{*}\propto
M_{peak}\frac{\Omega_{b}}{\Omega_{d}}\frac{t_{form}}{t_{cool,form}}$ (5)

Here $t_{form}$ is the formation time of the halo and $t_{cool,form}$ is the
time required for the baryons to cool and condense into a galaxy. This is
simply putting in the form of an equation the classical idea of “monolithic
collapse" first proposed by Eggen et al. (1962).

We can relate the mean density of the galaxy to its mass and radius:

$\rho_{max}\equiv\frac{M_{max}}{\frac{4}{3}\pi r_{max}^{3}}$ (6)

In which $\rho_{max}$, $M_{max}$, and $r_{max}$ are the density, mass, and
radius at which the circular velocity reaches $v_{max}$, where $v_{max}^{2}$
$=$ $\frac{GM_{max}}{r_{max}}$.

We can relate tform to the halo density assuming standard gravitational
collapse (Gunn & Gott 1972):

$G<\rho>\equiv t_{form}^{-2}$ (7)

We can also relate tcool,form to density using energy conservation and the
standard cooling equations:

$\frac{\frac{3}{2}kT_{max}}{m}\equiv\frac{GM_{max}}{r_{max}}$ (8)

The above equation defines Tmax. We subsequently can define tcool,form as:

$\Lambda(T_{max})\rho_{max}^{2}\equiv\frac{\frac{3}{2}\rho_{max}kT_{max}}{t_{cool,form}}$
(9)

Thus, tcool,form $\propto$ $\rho_{max}^{-1}f^{-1}$ where
$f$$\equiv$$\Lambda(T_{max})/T_{max}$, and tform $\propto$
$\rho_{max}^{-\frac{1}{2}}$. If we use these relations in Equation 5, we find
that

$M_{*}\propto
M_{peak}\rho_{max}^{\frac{1}{2}}f\propto(\frac{M_{max}}{r_{max}})^{\frac{3}{2}}f\propto
v_{max}^{3}f$ (10)

For low mass halos for which $f$ is nearly proportional to Tmax this indicates
a steep dependence on M∗ on $v_{max}$, which flattens for higher mass galaxies
with higher Tmax.

We highlight that this derivation is based on the assumption of spherical
collapse of the halo and pure radiative cooling. We do not consider any
complicating processes that we know affect galaxies in the universe, such as
mergers or feedback from star formation or AGN. In fact, we might expect this
relation between M∗ and $v_{max}^{3}f$ to break down more often for higher
mass galaxies, as they have been found to have later growth times where these
assumptions clearly breakdown (Behroozi et al., 2013). For first brightest
systems, sitting in massive halos from which they can accrete satellites, one
would expect Mpeak to be more relevant, and, as we have noted, both
observations (van Dokkum et al., 2010, e.g.) and LCDM theory argue that
hierarchical accretion is the dominant process for first brightest galaxies.

Indeed, we can go a step farther and ask the basis for and the value of the
transition mass above which “normal" growth of the stellar component from a
cooling collapse becomes difficult. This was addressed in a paper by Rees &
Ostriker (1977)(eqn 20) in an elementary discussion of the maximum mass of
cosmic gas that can cool and collapse in a dynamical time. They did not
include the important effects of dark matter in their treatment and obtained a
mass of $[(\frac{\hbar c}{Gm_{p}^{2}})^{2}(\frac{e^{2}}{\hbar
c})^{5}(\frac{m_{p}}{m_{e}})^{\frac{1}{2}}]m_{p}$$\sim$1012 M⊙ in baryons. Had
the effects of dark matter been included the value of the baryonic, transition
mass would have been reduced somewhat, but the corresponding dark matter mass
would have approximated 1012.5 M⊙. In fact, the mass function of galaxies in
the standard Press-Schechter parameterization declines exponentially above a
certain critical mass, the stellar mass being roughly 1011 M⊙ and the
corresponding halo mass being roughly 1012.5 M⊙. Consequently, we have both an
observational and a physical basis for expecting that galaxies above some
critical mass will grow primarily by accreting satellites and cannot be formed
easily by a monolithic collapse. Thus, while matching based on $v_{max}$ will
be best for normal systems, we can expect that, for first brightest galaxies
in massive clusters, Mpeak should be the best metric.

###  3.2 Rank Ordering in IllustrisTNG

Here, we test these theoretical predictions using the IllustrisTNG simulation.
As described in Section 2.2, we use a sample of galaxies with dark matter mass
greater than 1011 M⊙/$h$ in the DMO simulation and stellar masses above 109
M⊙/$h$ in the full hydrodynamical simulation.

We rank-ordered our selected galaxies by total mass and the stellar mass
separately for each of our samples: "all" (11927), "satellites" (2337), and
"centrals" (9590). We have used three simple proxies for dark matter halo mass
in our ranking schemes: the current dark matter mass, MDM, the peak dark
matter mass, Mpeak, and the current $v_{max}$. These are shown in order from
the top to bottom panels in Figure 1. We show the total mass proxy and stellar
mass for each of our galaxies as orange “o". The blue lines show the predicted
stellar mass using the rank-ordering method for each of our samples.

We highlight that M∗ is the current stellar mass at $z$ $=$ 0, including any
mass loss from star particles due to stellar evolution (Vogelsberger et al.,
2013; Leitner & Kravtsov, 2011; Wiersma et al., 2009). Within $\sim$3 Gyr of
formation, 45% of a star particle’s mass can be lost due to stellar evolution
(Leitner & Kravtsov, 2011), inserting a factor of $\sim$2 into Equation 10.
However, we cannot simply include a constant term into the relation between
current stellar mass and halo mass (or $v_{max}$) because of the time-
dependence of stellar evolution, and the fact that mass recycled from earlier
stellar populations may be required to form later generations of stars. The
time-dependence of stellar mass loss would move galaxies with earlier
formation times to lower M∗ as they have had more time to recycle stellar mass
back to gas mass. However, we note that earlier formation times are correlated
with higher stellar mass at a constant halo mass (Matthee et al., 2017, e.g.),
so stellar evolution likely only flattens this relation.

We see that the scatter decreases as we move from using MDM to Mpeak, although
the M∗ $\propto$ M${}_{halo}^{(3/4)}$ for higher masses holds for both
variables. We also highlight that the rank-order line for satellite galaxies
is much closer to that for centrals when we use Mpeak than when we use MDM.
The scatter continues to decrease when we use $v_{max}$ as our dark matter
halo mass proxy, particularly at lower masses (lower $v_{max}$). To guide the
eye we have overplotted simple power-law relations between M∗ and $v_{max}$.

Figure 1: The stellar mass of galaxies versus possible variables to use for
ranking. Blue lines show the predicted stellar mass using the rank-ordering
method for each of our samples. Top panel: Ranking using the current dark
matter halo mass has the most scatter. Middle panel: Using Mpeak for ranking
reduces scatter, and the rank ordering predictions for satellite and central
galaxies is much closer. Bottom panel: Ranking using $v_{max}$ reduces scatter
dramatically, particularly for lower mass (lower $v_{max}$) halos, as
quantified in Table 1.

The results shown in Figure 1 are quantified using the percent error as
described in Section 2.5 (eqn 4), with the results shown in Table 1. We see
that while using the current MDM in the matching scheme is reasonably accurate
for central galaxies, it is much less accurate for satellite systems.
Therefore, we also consider the peak mass of the halo, Mpeak. Using Mpeak
should correct for mass loss from satellite galaxies due to tidal stripping.
Because dark matter is distributed to a larger radius than the stars in a
galaxy, it will be more strongly stripped.

Therefore, while we expect Mpeak to be very similar to MDM for central
galaxies, it can vary by a considerable amount for satellites. Indeed, we see
in Table 1 that the improvement for centrals is very small when using Mpeak
rather than MDM, but it is dramatic for satellite galaxies. We also consider
$v_{max}$, as tidal stripping is found to have little effect on this property,
likely because the maximum rotational velocity is reached at relatively low
radii. Using $v_{max}$ for our variable we find that the error for central
galaxies has improved by more than 15%, although the correlation between
$v_{max}$ and stellar mass for satellites is somewhat weaker than the
correlation between Mpeak and stellar mass. However, because most of our
galaxies are centrals, $v_{max}$ remains the best single variable for rank-
ordering our galaxy sample.

We stress that, because of the shape of the mass function, any relation
between stellar mass and halo mass will be dominated by the lowest-mass
galaxies. Therefore, we also consider separately only galaxies whose mass in
the DMO simulation is greater than 1012 M${}_{\odot}/h$ in order to remove the
bulk of low mass galaxies while still retaining a sample with $\sim$200
satellite galaxies.

Unsurprisingly, we find that, as in the full sample, ordering using MDM is
reasonably accurate for central galaxies, but much less so for satellite
systems. Again, we find a large improvement in the ranking scheme for
satellite galaxies using Mpeak.

However, unlike in the full sample, $v_{max}$ is the worst ranking variable
for central galaxies with halo masses above 1012 M${}_{\odot}/h$. This agrees
well with our theoretical argument that at large masses Mpeak will be the best
ranking variable due to merging.

With this empirical support for the trends predicted in our model, we also
develop a straightforward variable, using the physical intuition from above,
that $v_{max}$ will be the best ranking variable for low mass galaxies and
Mpeak will be the best ranking variable for high mass galaxies (Section 3.1).
For this variable we normalize both $v_{max}$ and Mpeak to their values at a
“pivot mass" of Mpeak $=$ 1012.7 M⊙. We call these variables $v_{norm}$
$\equiv$ $v_{max}/v_{max,12.7}$ and $m_{norm}\equiv M_{peak}/10^{12.7}$. We
then rank order our galaxies using the parameter based on these normalized
values:

$\phi\equiv v_{norm}+m_{norm}$ (11)

Using this parameter, low mass galaxies depend more strongly on $v_{max}$,
while high mass galaxies depend on Mpeak. Both the exact value of the pivot
mass and the powers of $v_{norm}$ and $m_{norm}$ were selected to minimize
error while fleshing out our theoretical scaffold.

As shown in Table 1, using this parameter $\phi$ gives some improvement on the
fit to the central galaxies in our sample, and dramatically reduces the error
for the satellite galaxies. Using this variable for the mix of all galaxies
reduces the error by a substantial 33% when compared with rank ordering by
MDM.

Number of galaxies (MDM $>$ 1011 M⊙) | 11927 | 9590 | 2337 | 11927 | 11927
---|---|---|---|---|---
Galaxy Sample | All | Centrals | Satellites | Mix | % Improvement
Rank Ordering using MDM | 0.198 | 0.130 | 0.279 | 0.159 | –
Rank Ordering using Mpeak | 0.136 | 0.127 | 0.133 | 0.128 | 19
Rank Ordering using $v_{max}$ | 0.116 | 0.106 | 0.137 | 0.112 | 30
Rank Ordering using $\phi$ $\equiv$ $v_{norm}$ \+ $m_{norm}$ | 0.111 | 0.101 | 0.119 | 0.105 | 34

Table 1: The percent error (eqn 4) of the ranking method using a single
variable chosen to be either MDM (at $z=0$), Mpeak, $v_{max}$ (at $z=0$), and
$\phi$ $\equiv$ $v_{norm}$ \+ $m_{norm}$ (eqn 11) for the dark matter mass for
the different samples (a fit to all galaxies, only centrals, only satellites,
and mix of all galaxies fitting the centrals and satellites separately). We
use a galaxy sample with dark matter mass in the DMO simulation greater than
1011M${}_{\odot}/h$ that is matched to any galaxy in the hydrodynamical run
with stellar mass greater than 109M${}_{\odot}/h$. The final column shows the
percent improvement of ranking by the selected variable compared to MDM (at
$z=0$) on the “mix" sample. We see that ranking using the single variable
$\phi$ reduces the error 34% compared to matching by MDM (at $z=0$).

Number of galaxies (MDM $>$ 1011 M⊙) | 11927 | 9590 | 2337 | 11927 | 11927
---|---|---|---|---|---
Galaxy Sample | All | Centrals | Satellites | Mix | % Improvement
$\phi$ \+ $v_{disp}$ | 0.112 | 0.102 | 0.117 | 0.105 | 0
$\phi$ \+ $v_{max}$ | 0.111 | 0.101 | 0.117 | 0.104 | 1
$\phi$ \+ MDM | 0.105 | 0.101 | 0.110 | 0.103 | 2
$\phi$ \+ Mpeak | 0.111 | 0.101 | 0.117 | 0.104 | 1
$\phi$ \+ rmax | 0.111 | 0.101 | 0.118 | 0.105 | 0
$\phi$ \+ rDM | 0.105 | 0.100 | 0.114 | 0.103 | 2
$\phi$ \+ $c_{v}$ | 0.111 | 0.101 | 0.118 | 0.105 | 0
$\phi$ \+ $c_{h}$ | 0.109 | 0.101 | 0.117 | 0.104 | 1
$\phi$ \+ $c_{r}$ | 0.109 | 0.101 | 0.116 | 0.104 | 1
$\phi$ \+ $t_{peak}$ | 0.105 | 0.101 | 0.110 | 0.103 | 2
$\phi$ \+ $t_{50}$ | 0.106 | 0.099 | 0.116 | 0.102 | 3
$\phi$ \+ $t_{85}$ | 0.104 | 0.099 | 0.111 | 0.101 | 4
$\phi$ \+ MDM,r<1Mpc | 0.104 | 0.100 | 0.115 | 0.103 | 2
$\phi$ \+ MDM,r<2Mpc | 0.103 | 0.099 | 0.113 | 0.102 | 3
$\phi$ \+ MDM,r<5Mpc | 0.105 | 0.099 | 0.115 | 0.102 | 3
$\phi$ \+ MDM,r<8Mpc | 0.107 | 0.100 | 0.116 | 0.103 | 2
$\phi$ \+ MDM,r<15Mpc | 0.109 | 0.100 | 0.117 | 0.104 | 1
$\phi$ \+ MDM,r<2Mpc \+ $t_{85}$ | 0.101 | 0.097 | 0.108 | 0.099 | 6

Table 2: The percent error using two variables: the $\phi$ ranking method plus
the listed corrections. We use the galaxy sample with dark matter mass greater
than 1011M${}_{\odot}/h$ that is matched to any galaxy in the hydro run with
stellar mass greater than 109M${}_{\odot}/h$. Note that MDM,r<XMpc is the mass
of all halos within that radius from the DMO simulation with MDM $>$ 109M⊙,
including the mass of the halo from which the measurement originates. All halo
properties are measured using the DMO simulation. Here the final column shows
the percent improvement on the “mix" sample of using the correction variable
in addition to rank-ordering by $\phi$ in comparison to only rank-ordering by
$\phi$.

##  4 Using Secondary Variables to Improve Rank Ordering

We now attempt to minimize the scatter in the $\phi$ \- M∗ relation using
other features of dark matter halos. These features are listed in Table 2. We
have roughly grouped the halo properties into those related to the halo mass
(MDM, Mpeak, $v_{disp}$$\equiv$dark matter velocity dispersion, and
$v_{max}$), size (rmax$\equiv$ $v_{max}$ radius and rDM$\equiv$ dark matter
half mass radius), shape (concentration using the three methods described in
Section 2.4), formation time (the lookback time to Mpeak, to when the halo
reaches 50% of its z=0 mass, and to when the halo reaches 85% of its z=0
mass), and surrounding environment (using the mass of dark matter halos within
various radii).

In this section we first describe the method we use to include secondary
features, and then discuss the results.

###  4.1 Method of Correction

Figure 2: The top and bottom panels show the first and second steps used to
include a secondary halo feature to reduce the scatter in rank ordering halos
(Section 4.1). Here we use the MDM,r<2Mpc environmental measure, written as
M2Mpc. Top: First we plot this variable as a function of $\phi$, our rank-
ordering variable. Here we show the total sample (“all") as well as the
centrals and satellites ("sats") separately. The points are color-coded as
centrals and satellites. Bottom: Using the scatter from a rolling median, we
can find that the ratio of the true stellar mass of the galaxy to the rank-
ordered assigned mass has a dependence on MDM,r<2Mpc. The points are not
color-coded for satellites and centrals, as for the “all" fit we use all of
the galaxies in the sample. We can then correct our stellar mass using this
dependence. Notice that the rolling median for the total sample is similar to
that for the separated satellite and central samples.

We first plot our feature as a function of our best single variable $\phi$,
and find the running median of the feature using a window size of 50 galaxies.
We have tested using other window sizes (25 and 100 galaxies as well as a
constant $\Delta log(\phi)$ bin of $\pm$0.3 around each galaxy) and find
similar results in our percent improvement. The top panel of Figure 2 shows
this plot using the environmental density MDM,r<2Mpc variable. Clearly there
is a trend of increasing environmental density as a function of $\phi$, and it
differs for satellites and centrals.

Because we use the rolling median, we need to remove the first and last 25
values, so we are left with an "all" sample of 11877, a "centrals" sample of
9540, and a "satellites" sample of 2287 galaxies. Removing these galaxies has
little impact on the percent errors using the rank ordering method for each
sample (a change of less than 1%).

We then plot Mtrue/Mrank as a function of $\Delta$log(feature), which is the
difference between the log(feature) for each dark matter halo and the
log(featurerollingmedian) found at each $\phi$. In the bottom panel of Figure
2 we show how Mtrue/Mrank is related to the scatter in MDM,r<2Mpc. This
relation is fit using first, second and third order polynomial fits. Finally,
we correct our prediction for the stellar mass using our chosen fit as below
(the second order polynomial fit shown tends to give the best results):

$log(M_{*,pred})=log(M_{*,rank})+\\\ \alpha\Delta log(feature)^{2}+\beta\Delta
log(feature)+\gamma$ (12)

Finally, we calculate the percent error of the new prediction. This value for
each feature and galaxy population is shown in Table 2.

###  4.2 Results

All of our quantitative results are shown in Table 2. The most glaring result
is that most corrections do not result in a large improvement of the percent
error from ranking using $\phi$.

Using random resampling of 70% of our data sets (“all",“central", and
“satellites") 60 times, we find a distribution of errors with means matching
the values listed for $\phi$ of the complete sample in Table 1, and standard
deviations of 0.001, 0.001, 0.002, and 0.0009 for “all", “centrals",
“satellites" and “mix" samples, respectively. With this in mind we can look
more closely at the improvement when adding a second feature to our matching
scheme.

In some more detail, it is not surprising that all of the halo features
describing halo mass do not improve the fit to central galaxies at all. These
are well-fit by our $\phi$ variable. However, interestingly, the error is
reduced for the satellite sample when we include a MDM correction. This may be
because we have largely ignored satellite galaxy evolution by choosing
$v_{max}$ and Mpeak as the components of $\phi$. Including MDM may start to
include the later evolution of these galaxies.

We find universally small improvement when considering our variables
describing halo size (rmax and rDM) and shape (concentration). Although
previous work has shown that the scatter in abundance matching is related to
concentration (Matthee et al., 2017, e.g.), it is not surprising that
concentration does not improve the scatter when ranking using our $\phi$
variable. This is because $\phi$ includes $v_{max}$, which is already related
to concentration (Klypin et al., 2011).

Interestingly, there is some improvement in the error when folding formation
time into the stellar mass estimate. For example, $t_{85}$ is the halo feature
that results in the smallest percent errors across all of our samples: “all"
galaxies, “centrals", “satellites", and the “mix" sample. As discussed above,
although we find that earlier formation times are correlated with higher
stellar mass at a constant $\phi$ (in agreement with previous work), stellar
evolution causing mass loss over several Gyr may flatten this relationship.

Finally, using environment to correct for the stellar mass also has a small
impact on the overall error. Despite this, we note that including the mass
from galaxies within 2-5 Mpc seems to produce a slightly better correction
than smaller or larger environment windows.

####  4.2.1 Correcting Using A Combination of Environment and Formation Time

Finally, we use our fits for each of our strongest individual corrections,
$t_{85}$ and MDM,r<2Mpc, to create a combined correction on the rank-ordering
technique.

$\begin{array}[]{l}M_{*,pred}=log(M_{*,rank})+\\\ (\alpha_{M2}\Delta
log(M_{DM,r<2Mpc})^{2}+\beta_{M2}\Delta log(M_{DM,r<2Mpc})+\\\
(\alpha_{t85}\Delta log(t_{85})^{2}+\beta_{t85}\Delta
log(t_{85})+\gamma\end{array}$ (13)

We use the curvefit module in scipy to perform a least-squares fit to the
above equation, and find that we can reduce the error using both MDM,r<2Mpc
and $t_{85}$ as shown in the final line of Table 2.

| All | Centrals | Satellites | Mix
---|---|---|---|---
$\phi$ | 0.111,0.111 | 0.102,0.101 | 0.120,0.117 | 0.105,0.104
$\phi$ \+ MDM,r<2Mpc | 0.104,0.104 | 0.099,0.099 | 0.114,0.112 | 0.102,0.101
$\phi$ \+ $t_{85}$ | 0.104,0.104 | 0.099,0.100 | 0.111,0.109 | 0.101,0.102
$\phi$ \+ MDM,r<2Mpc \+ $t_{85}$ | 0.101,0.100 | 0.097,0.097 | 0.109,0.104 | 0.099,0.097

Table 3: The median percent error on ten iterations using a training and test
sample. For each sample and fitting method we list the training then test
sample values.

###  4.3 Verifying our Results

Here we use two methods to verify our results on the improvement using
multiple halo features to determine stellar mass.

####  4.3.1 Random Forest Regression

Now that we have gained insight into the level of improvement that can be
gained by using more than one feature of dark matter halos in the abundance
matching technique, we turn to machine learning to provide an independent
check of our modeling and ranking scheme.

Using Random Forest Regression (RFR) allows us to rank the features according
to their effect on the model output, and has the additional benefit of
expanding the space of available models beyond polynomial fitting. For this
work we use scikit-learn (Pedregosa et al., 2011).

First, we are able to reproduce the percent error on the entire sample using
only our defined $\phi$ feature (0.111), and in a two-feature setting where we
add a central/satellite galaxy label (0.105). We check the rest of our ranking
parameters from Table 1 and verify that $\phi$ produces the best ranking
variable to match DM halos to galaxies. Also, we confirm that our $\phi$
variable produces lower error values than the combination of $v_{max}$ and
Mpeak.

We also use our four selected halo features that we found produced the best
match between the DMO and hydrodynamical simulations, $\phi$, MDM,r<2Mpc,
$t_{85}$, and the central/satellite label. Using an optimized RF regressor,
the expected test set error is 0.098 with a standard deviation of 0.0013,
quite similar to the percent error we find ranking the satellites and centrals
separately using $\phi$ and applying our analytic correction using MDM,r<2Mpc
and $t_{85}$. This is reassuring because it shows that our results are only
very mildly dependent on the modeling assumptions.

We can also include all the features and use a parameter optimization
technique to find the minimum possible error of a Random Forest Regression. We
find a minimum error of 0.092 using eight randomly selected features, creating
100 trees (nestimators) with a maximum depth of 14 branches (maxdepth).
However, there are more than 30 combinations within one standard deviation
(0.0015), including one using only 4 features. We can conclude that there may
be many similarly relevant predictors in our feature list. This supports our
analytic reasoning that several of our halo features are reasonable proxies
for halo mass, and we have already noted that our other halo features can be
separated into only a few types of variables (halo size, concentration,
formation time, and environment).

Indeed, if we optimize the RFR including one feature of each type we can reach
an error of 0.095 ($\phi$, MDM,r<2Mpc, $t_{50}$, rDM, and $c_{v}$). This is
within two standard deviations of the four halo parameters we use in our
analytic model, and so does not indicate a dramatic improvement.

Comparing our results to the errors found using the RFR machine learning
technique gives us assurance that our analytic method for including extra halo
features is reasonable, and that our conclusions are not strongly model-
dependent. While continuing to add features can reduce the error on the
matching scheme, we do not find other clear DM halo features that dramatically
improve upon our analytic method.

####  4.3.2 Cross-Validation

In order to obtain another view on whether increasing the number of halo
features improves our estimate of stellar mass we can use cross-validation.
This can be used to determine how meaningful our derived improvements are when
used to predict the stellar mass of galaxies. Cross-validation is specifically
designed to trade off over- and under-fitting to give the highest prediction
accuracy. For this, we randomly select 80% of our sample as our test set, on
which we perform the fitting processes as described. We use the remaining 20%
as our test set to determine if the percent error on the stellar mass
prediction improves when including more features. Specifically, we select 80%
of our total sample for the “all" fits, and then 80% of the central and
satellite samples, in order to determine the “central", “satellite", and “mix"
fits.

We performed this cross-validation routine ten times using ten different
random subsets of the data, and universally find improvement in both the
training and test sets when using MDM,r<2Mpc, $t_{85}$, or their combination.
In Table 3, we list the median percent error values for the ten sets of
training and test samples. We can conclude that we have not yet overfit using
these halo features, and our improvement in predicting stellar masses from
halo masses is real, albeit small.

##  5 Discussion

What have we learned from this exercise? The zeroth order conclusion is that a
matching scheme based on the maximum velocity in a dark matter halo is a good
single predictor of the final stellar mass for normal galaxies, whether they
are central galaxies or satellites. The typical error in the prediction (in
the IllustrisTNG100 simulations) is 11.6 percent in log(M∗) compared to 19.8
using MDM, and the dependence of stellar mass on $v_{max}$ is unsurprisingly
log (M∗) $\sim$ (3.8 $\pm$ 0.02) log($v_{max}$) (using bootstrap resampling
with 70% of the data set), close to the Faber-Jackson relation (Faber &
Jackson, 1976). This result is just what one would have expected from the
simplest physical argument that estimates the amount of gas that can be turned
into stars in the standard Gunn & Gott (1972) collapse of a dark matter halo.

But, for high mass systems comparable to the first brightest galaxies in
clusters living in halos more massive than 1012.7 M⊙, the accretion of
satellite systems will significantly increase the stellar mass and the most
relevant halo parameter is simply the peak dark matter mass, Mpeak. Using a
single variable, $\phi$ (Equation (11)), which incorporates both features
reduces the error to 10.5% when satellites and centrals are ranked separately.

These prescriptions should be easy to implement and can replace the simplest,
halo mass based initial matching schemes when estimating the expected galaxy
stellar masses given a dark matter simulation.

If one wants to go farther and improve the best zeroth order scheme by first
order corrections then we have found that a roughly 6% improvement is
possible. Interestingly, environmental considerations that we considered did
not lead to significant improvement even in satellite galaxies, and the best
single variable for improvement was $t_{85}$, the time at which a halo reached
85% of its peak mass.

However, an almost mindless combination of the two variables ($v_{max}$,
Mpeak) worked best. We further found that a simple linear combination based on
these two variables enables predictions to a typical accuracy of 10.5 percent
error in log(M∗).

##  6 Tests

All of these results are based on simulated data and it is important to test
them in the real world. We have been able to think of two tests that might be
applied to help determine whether the proposed matching scheme provides a
significant improvement over the simplest matching scheme. First one
constructs a standard LCDM, dark matter only simulation and, using a standard
halo finding algorithm, makes a catalog of dark matter halos labeling each of
them with the final dark matter mass, MDM, the peak dark matter mass Mpeak
over the history of the halo and the current halo maximum circular velocity
$v_{max}$. Then, to test the classic matching scheme (as has been done before
– Conroy et al., 2006), one takes a representative volume and rank orders the
halos by MDM, takes catalog values for a comparable volume (from, say, the
Sloan Digital Sky Survey) and rank orders the observed galaxies by (for
example) $g$ or $r$ magnitudes and then identifies each DM halo with the
matched by ranking, real galaxy. This gives one an artificial catalog of
galaxies each tagged with a position, a velocity and a $g$ or $r$ optical
magnitude.

Then one would “observe” this synthesized catalog and construct two spatial
distribution functions, a galaxy-galaxy spatial correlation function (Conroy
et al., 2006; Hearin et al., 2013; Hearin & Watson, 2013; Reddick et al.,
2013) and a void distribution function (Walsh & Tinker, 2019, e.g.). These
could then be compared to the known galaxy-galaxy spatial correlation
functions and the known void distribution functions with both one parameter
functions specified as a function of magnitude. The magnitude distribution
itself is of course correct by construction. Then comparing – say – the
autocorrelation length as a function of galaxy magnitude between the real and
synthesized data sets allows one to determine the fractional error as a
function of galaxy magnitude.

Then one would go back to the original DM halo catalog and, using
(Mpeak,$v_{max}$), construct for each halo the value of $\phi$ =
($v_{max}$/$v_{1}$) + (Mpeak/(M1)) (Equation 11), where ($v_{1}$, M1) are the
values of ($v_{max}$, Mpeak) for the average halo of mass 1012.7 M⊙. Now, with
each halo tagged with its value of $\phi$, one can rank order the synthetic
sample by $\phi$ and attach visual magnitudes to each galaxy by the same
method as was done using MDM. Now one has a new catalog to observe with
respect to spatial distribution metrics and can again find the fractional
error in – for example – the spatial autocorrelation length as a function of
visual brightness and compute the error by comparing to real observed data.

This procedure would give us a quantitative estimate as to how well the
matching scheme was working compared to both reality and the previous simpler
matching scheme which has had considerable success. And, unlike the exercises
in this paper, the tests would not be dependent on the accuracy of our current
galaxy formation algorithms, which, while well tested, suffer from the
“confirmation bias” inevitable when uncertain modelling parameters are
adjusted to fit observations. We look forward to pursuing these independent
tests in future work.

##  7 Conclusions

In this paper we have examined schemes to populate a synthesized dark matter
only set of cosmological simulations with galaxies to see if we could devise a
simple and accurate scheme. We took as our starting point a matching scheme
(Vale & Ostriker, 2004, 2006, 2008) which, while almost naively simple has had
some success. In that scheme, one rank orders DM halos by final mass and rank
orders real galaxies in a similar cosmic volume by luminosity and attaches to
the kth ranked halo the kth ranked galaxy. Table 1 represents of one parameter
efforts which we compared to the computed luminosities in the IllustrisTNG
simulated galaxy catalog. Table 2 summarizes our results with two parameter
fits where we used combinations of velocity dispersion, mass and environmental
density. We did not find that adding an environmental variable produced a
significant improvement over simpler schemes nor did we find that any of the
two parameter fits that we investigated were statistically significantly
superior to the one parameter fits. What we did discover was that a new single
variable, $\phi$ (cf equation 12), which combines information from both mass
and velocity variables, provides a quite significant improvement over the
basic ranking scheme using final dark matter mass, the error being reduced by
about 33% percent. In our examination of the physical basis for the success of
this new variable we examined simple arguments starting with the over half
century old paper by Rees & Ostriker (1977).

There is a critical mass for galaxies – the mass above which it cannot cool by
normal radiative processes in roughly a free fall time. That mass corresponds
roughly to 1012.7 M⊙ which we designate as M1. Below this mass there is a
simple analytic argument that asks if a gaseous object can cool in its own
free fall time and is equivalent to M∗ $\sim$ $v_{max}^{3}f$. For masses above
M1 growth only occurs by accretion of satellites and that is proportional to
M. So, we designed a metric, $\phi$, which is dominated by velocity for low
mass objects and dominated by mass for high mass objects more massive than M1.
This single variable, based on the physical motivation given above, seems to
provide a matching scheme superior to others which we have tested. We did try
other combinations of (Mpeak,$v_{max}$) and found none superior to the simple
variable, $\phi$, that we had tested. So our bottom line is that the variable,
$\phi$ (Equation 11), is the best single variable to use in predicting the
stellar mass of galaxies, given their halo properties.

We would like to thank the referee for helpful comments that improved the
paper. ST would like to thank Claire Kopenhafer and Tjitske Starkenburg for
their help and scripts in reading in and analyzing TNG outputs, Viviana
Aquaviva for her Machine Learning class and comments on the draft, and Dan
Foreman-Mackey for discussions and comments on cross-validation. ST gratefully
acknowledges support from the Center for Computational Astrophysics at the
Flatiron Institute, which is supported by the Simons Foundation. The data used
in this work were hosted on facilities supported by the Scientific Computing
Core at the Flatiron Institute, a division of the Simons Foundation.

## References

  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 762, L31, doi: 10.1088/2041-8205/762/2/L31
  * Berrier et al. (2006) Berrier, J. C., Bullock, J. S., Barton, E. J., et al. 2006, ApJ, 652, 56, doi: 10.1086/507573
  * Blank et al. (2019) Blank, M., Macciò, A. V., Dutton, A. A., & Obreja, A. 2019, MNRAS, 487, 5476, doi: 10.1093/mnras/stz1688
  * Bose et al. (2019) Bose, S., Eisenstein, D. J., Hernquist, L., et al. 2019, MNRAS, 490, 5693, doi: 10.1093/mnras/stz2546
  * Chaves-Montero et al. (2016) Chaves-Montero, J., Angulo, R. E., Schaye, J., et al. 2016, MNRAS, 460, 3100, doi: 10.1093/mnras/stw1225
  * Conroy et al. (2006) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201, doi: 10.1086/503602
  * Crain et al. (2015) Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937, doi: 10.1093/mnras/stv725
  * Davé et al. (2016) Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265, doi: 10.1093/mnras/stw1862
  * Dragomir et al. (2018) Dragomir, R., Rodríguez-Puebla, A., Primack, J. R., & Lee, C. T. 2018, MNRAS, 476, 741, doi: 10.1093/mnras/sty283
  * Eggen et al. (1962) Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748, doi: 10.1086/147433
  * Faber & Jackson (1976) Faber, S. M., & Jackson, R. E. 1976, ApJ, 204, 668, doi: 10.1086/154215
  * Genel et al. (2014) Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175, doi: 10.1093/mnras/stu1654
  * Gunn & Gott (1972) Gunn, J. E., & Gott, J. Richard, I. 1972, ApJ, 176, 1, doi: 10.1086/151605
  * Guo et al. (2010) Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111, doi: 10.1111/j.1365-2966.2010.16341.x
  * He (2020) He, J.-h. 2020, MNRAS, 493, 4453, doi: 10.1093/mnras/staa620
  * Hearin & Watson (2013) Hearin, A. P., & Watson, D. F. 2013, MNRAS, 435, 1313, doi: 10.1093/mnras/stt1374
  * Hearin et al. (2013) Hearin, A. P., Zentner, A. R., Berlind, A. A., & Newman, J. A. 2013, MNRAS, 433, 659, doi: 10.1093/mnras/stt755
  * Hopkins et al. (2018a) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018a, MNRAS, 477, 1578, doi: 10.1093/mnras/sty674
  * Hopkins et al. (2018b) —. 2018b, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * 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
  * Lehmann et al. (2017) Lehmann, B. V., Mao, Y.-Y., Becker, M. R., Skillman, S. W., & Wechsler, R. H. 2017, ApJ, 834, 37, doi: 10.3847/1538-4357/834/1/37
  * Leitner & Kravtsov (2011) Leitner, S. N., & Kravtsov, A. V. 2011, ApJ, 734, 48, doi: 10.1088/0004-637X/734/1/48
  * Mao et al. (2015) Mao, Y.-Y., Williamson, M., & Wechsler, R. H. 2015, ApJ, 810, 21, doi: 10.1088/0004-637X/810/1/21
  * Marín et al. (2008) Marín, F. A., Wechsler, R. H., Frieman, J. A., & Nichol, R. C. 2008, ApJ, 672, 849, doi: 10.1086/523628
  * Martizzi et al. (2020) Martizzi, D., Vogelsberger, M., Torrey, P., et al. 2020, MNRAS, 491, 5747, doi: 10.1093/mnras/stz3418
  * Matthee et al. (2017) Matthee, J., Schaye, J., Crain, R. A., et al. 2017, MNRAS, 465, 2381, doi: 10.1093/mnras/stw2884
  * Moliné et al. (2017) Moliné, Á., Sánchez-Conde, M. A., Palomares-Ruiz, S., & Prada, F. 2017, MNRAS, 466, 4974, doi: 10.1093/mnras/stx026
  * Naab & Ostriker (2017) Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59, doi: 10.1146/annurev-astro-081913-040019
  * Nagai & Kravtsov (2005) Nagai, D., & Kravtsov, A. V. 2005, ApJ, 618, 557, doi: 10.1086/426016
  * 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. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Ostriker et al. (2019) Ostriker, J. P., Choi, E., Chow, A., & Guha, K. 2019, ApJ, 885, 97, doi: 10.3847/1538-4357/ab3288
  * Pakmor et al. (2011) Pakmor, R., Bauer, A., & Springel, V. 2011, MNRAS, 418, 1392, doi: 10.1111/j.1365-2966.2011.19591.x
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Reddick et al. (2013) Reddick, R. M., Wechsler, R. H., Tinker, J. L., & Behroozi, P. S. 2013, ApJ, 771, 30, doi: 10.1088/0004-637X/771/1/30
  * Rees & Ostriker (1977) Rees, M. J., & Ostriker, J. P. 1977, MNRAS, 179, 541, doi: 10.1093/mnras/179.4.541
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Tasitsiomi et al. (2004) Tasitsiomi, A., Kravtsov, A. V., Wechsler, R. H., & Primack, J. R. 2004, ApJ, 614, 533, doi: 10.1086/423784
  * Trujillo-Gomez et al. (2011) Trujillo-Gomez, S., Klypin, A., Primack, J., & Romanowsky, A. J. 2011, ApJ, 742, 16, doi: 10.1088/0004-637X/742/1/16
  * Vale & Ostriker (2004) Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189, doi: 10.1111/j.1365-2966.2004.08059.x
  * Vale & Ostriker (2006) —. 2006, MNRAS, 371, 1173, doi: 10.1111/j.1365-2966.2006.10605.x
  * Vale & Ostriker (2008) —. 2008, MNRAS, 383, 355, doi: 10.1111/j.1365-2966.2007.12544.x
  * van Dokkum et al. (2010) van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018, doi: 10.1088/0004-637X/709/2/1018
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031, doi: 10.1093/mnras/stt1789
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * Walsh & Tinker (2019) Walsh, K., & Tinker, J. 2019, MNRAS, 488, 470, doi: 10.1093/mnras/stz1351
  * Wang et al. (2015) Wang, L., Dutton, A. A., Stinson, G. S., et al. 2015, MNRAS, 454, 83, doi: 10.1093/mnras/stv1937
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009, MNRAS, 399, 574, doi: 10.1111/j.1365-2966.2009.15331.x
  * Xu & Zheng (2018) Xu, X., & Zheng, Z. 2018, MNRAS, 479, 1579, doi: 10.1093/mnras/sty1547
  * Zentner et al. (2014) Zentner, A. R., Hearin, A. P., & van den Bosch, F. C. 2014, MNRAS, 443, 3044, doi: 10.1093/mnras/stu1383