# High-redshift halo-galaxy connection via constrained simulations

Adi Nusser<EMAIL_ADDRESS>Department of Physics and the Asher
Space Research Institute
Israel Institute of Technology Technion, Haifa 32000, Israel

###### Abstract

The evolution of halos with masses around $M_{\textrm{h}}\approx
10^{11}\;\textrm{M}_{\odot}$ and $M_{\textrm{h}}\approx
10^{12}\;\textrm{M}_{\odot}$ at redshifts $z>9$ is examined using constrained
N-body simulations. The specific mass accretion rates,
$\dot{M}_{\textrm{h}}/M_{\textrm{h}}$, exhibit minimal mass dependence and
agree with existing literature. Approximately one-third of simulations reveal
an increase in $\dot{M}_{\textrm{h}}$ around $z\approx 13$, possibly implying
a dual-age stellar population. Comparing simulated halos with observed
galaxies having spectroscopic redshifts, we find that for galaxies at
$z\gtrsim 9$, the ratio between observed star formation rate (SFR) and
$\dot{M}_{\textrm{h}}$ is approximately $2\%$. This ratio remains consistent
for the stellar-to-halo mass ratio (SHMR) but only for $z\gtrsim 10$. At
$z\simeq 9$, the SHMR is notably lower by a factor of a few. At $z\gtrsim 10$,
there is an agreement between specific star formation rates (sSFRs) and
$\dot{M}_{\textrm{h}}/M_{\textrm{h}}$ of halos. However, at $z\simeq 9$,
observed sSFRs exceed simulated values by a factor of two. To explain the
relatively high star formation efficiencies in high-$z$ halos with
$M_{\textrm{h}}\approx 10^{11}M_{\odot}$, a simplified model is proposed,
assuming the applicability of the local Kennicutt-Schmidt law. The enhanced
efficiency relative to low-$z$ is mainly driven by the reduced effectiveness
of stellar feedback due to deeper gravitational potential for halos of a fixed
mass.

cosmology: large-scale structure – galaxies: formation – galaxies: high
redshift – galaxies: ISM – galaxies: luminosity function

##  1 Introduction

The standard $\Lambda$CDM cosmological model, incorporating a cosmological
constant, $\Lambda$, and cold dark matter (DM), has been remarkably successful
in interpreting and predicting fundamental properties of the large-scale
structure of the Universe. Despite potential tensions (e.g. Asgari et al.,
2021; Riess et al., 2023), this success extends to temperature anisotropies of
the cosmic microwave background (CMB), clustering of the distribution of
galaxies, and deviations of galaxy motions from a purely Hubble flow (e.g.
Bennett et al., 2013; Tegmark et al., 2004; Cole et al., 2005; Eisenstein et
al., 2007; Davis et al., 2011; Carrick et al., 2015; Planck Collaboration et
al., 2018; Lilow et al., 2021). On galactic scales, predicting the properties
of the galaxy population and its evolution with redshift has been less
straightforward. This complexity arises from the intricate nature of baryonic
physics involved in star formation processes, including gas dynamics, heating
and cooling mechanisms, and notably, the energetic feedback from supernovae
(SN) and active galactic nuclei (AGNs) (e.g. Larson, 1974; Dekel & Silk, 1986;
White & Frenk, 1991a; Silk & Rees, 1998; Dekel et al., 2009b; Okalidis et al.,
2021; Krumholz et al., 2018; Nusser & Silk, 2022).

Prior to the era of the James Webb Space Telescope (JWST) (Gardner et al.,
2023), significant efforts have been invested in fine-tuning models of galaxy
formation to adequately describe observations at low and moderately high
redshifts ($z~{}\hbox to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}10$)
(e.g. Finkelstein et al., 2022).

Observations obtained with the JWST have significantly deepened our view of
the universe, revealing galaxies as far back as a couple of hundred million
years near the Big Bang. However, the JWST has also detected an unexpected
excess of luminous galaxies at higher redshifts. While the initial findings
from the JWST posed serious challenges for the standard $\Lambda$CDM model of
structure formation, the severity of these discrepancies were largely
alleviated with more precise calibration and the availability of spectroscopic
redshifts (cf. Yung et al., 2024, for an overview).

It should be emphasized that the star formation rates (SFRs) in high-redshift
JWST galaxies are not particularly unusual in themselves (e.g. Robertson et
al., 2023; Harikane et al., 2024). These galaxies exhibit SFRs that can be
adequately sustained by cosmological gas accretion onto halos (Mason et al.,
2023). Matching the abundance of halos to the observed distribution of UV
magnitudes (used as proxies for the SFRs) of galaxies at $z\mathrel{\hbox
to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10$
implies that these galaxies should be hosted in halos of mass
$M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$ (Boylan-Kolchin, 2023; Mason
et al., 2023). For such halos, the star formation efficiency $f_{*}$, i.e. the
fraction of accreting gas needed account for the SFRs, is $\gtrsim 0.15$ (see
3.2.1 below). At low redshifts ($z<4$), the stellar-to-halo mass ratio (SHMR)
inferred from abundance matching is in the range $0.001-0.01$ for
$M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$ halos (e.g. Moster et al.,
2013; Rodríguez-Puebla et al., 2017; Girelli et al., 2020; Fu et al., 2022).
Assuming a global gas fraction $f_{b}=0.157$ in galaxies (Planck Collaboration
et al., 2018), this implies an average star formation efficiency,
$f_{*}\approx 0.06-0.006$, which is at least a factor of two lower than the
inferred value at $z>10$.

Additionally, uncertainties in estimated stellar masses remain large even for
best observed objects. Different assumptions in the star formation history
(SFH) in SED fitting can lead to non-negligible differences of the stellar
masses Whitler et al. (2023); Pacifici et al. (2023); Wang et al. (2023). This
is the case for example for the galaxy GS-z12 with a spectroscopic redshift
$z=12.48$ (D’Eugenio et al., 2023). Using the PROSPECTOR SED fitting (Johnson
et al., 2021), H24 found a stellar mass $M_{\ast}=4.3^{+1.8}_{-2}\times
10^{8}\,{\rm M}_{\odot}$ for this galaxy, while analysis using the BEAGLE core
(Chevallard & Charlot, 2016) performed by D’Eugenio et al. (2023) lead to
$M_{\ast}=4.36^{+1.8}_{-1.27}\times 10^{7}$ for the same object.

The evolution of the SHMR inferred through abundance matching techniques is
sensitive to various factors, including the shape of the galaxy stellar mass
function (e.g. Yang et al., 2012) As highlighted by Fu et al. (2022),
uncertainties in the observed stellar mass functions at redshifts $z<4$ can
lead to differing interpretations regarding the evolution of the SHMR from
redshift $z=0$ to $z=4$. Depending on specific assumptions made about the
stellar mass function, the SHMR could exhibit either a decrease or an increase
over this redshift range. For instance, for a halo with $M_{\textrm{h}}\approx
10^{11}\,{\rm M}_{\odot}$, the SHMR could vary by two orders of magnitude.
Seen in this light, it entirely plausible that the inferred star formation
efficiency varies between low and high redshifts. The key question is whether
galaxy formation can account for the boost in efficiency at high redshifts.

Numerous models aim to understand the formation of luminous high-redshift
galaxies. Harikane et al. (2023) suggested that the SFR inferred from UV
luminosities might be overestimated because of a prevalence of a top-heavy
initial stellar mass function (IMF) at high redshift. Such a top-heavy IMF
could arise naturally in the low metallicity environment characteristic of
Population III stars. Furthermore Yung et al. (2024) emphasized that a top-
heavy IMF could account for a factor of 4 boost in the UV luminosities, which
is needed to bring their semi-analytic models to agree with observations.

Another model predicts a feedback-free star formation at high redshifts. Dekel
et al. (2023) present a set of conditions on the gas density and surface
density that allow for a feedback-free star formation in gas clouds of
$\approx 10^{7}\,{\rm M}_{\odot}$ at high redshifts. At $z\approx 10$, the
conditions are satisfied in halos of mass $10^{11}\,{\rm M}_{\odot}$. The high
SFR efficiency depletes most the gas in the galaxy, effectively halting star
formation until enough new gas mass is accumulated by accretion and the next
episode of star formation is initiated.

However, Silk et al. (2024) invoke AGN positive feedback to explain the
elevated SFR. They propose that a short-lived AGN activity triggers vigorous
star formation activity in high redshift galaxies via the generation momentum-
conserving outflows. The predict that a late transition to energy-conserving
flows at $z\approx 6$ is associated with gas depletion and hence quenched star
formation at low redshifts.

According to Ferrara et al. (2023) the high abundance of $z\gtrsim 10$
galaxies could be explained if dust attenuation is decreased at high relative
to lower redshifts, in such a way that it compensates for the reduced
abundance of the host halos.

Moreover, there is a growing evidence for a highly variable star formation
history at high redshifts (Cole et al., 2023; Dressler et al., 2024). The
variability could be attributed to a combination of factors, including
mergers, interactions, and environmental conditions. The variability could
bias the inferred UV luminosity function as well as the stellar masses (Sun et
al., 2023a, b). The reason is that galaxies residing in low mass halos are
likely to be preferentially detected during phases of increased star
formation, which are associated with heightened production of UV photons.
Consequently, these galaxies might seem less numerous than they actually are.
This also could overestimate the inferred stellar mass. However, the required
variability for explaining the observations seems to exceed the predictions of
simulations (Pallottini & Ferrara, 2023), but see (Sun et al., 2023b).

In vanilla semi-analytic galaxy formation models (SAMs) (e.g. Kauffmann et
al., 1993; Benson, 2010; Cora, 2006; Henriques et al., 2015; Somerville &
Davé, 2015), the SHMR versus halo mass is almost independent of redshift
(Mitchell et al., 2018; Yung et al., 2024), suggesting that efficiency with
which gas converts into stars relative to the mass of their associated dark
matter halos remains relatively constant throughout the cosmic history.
According to these models, the expected star formation efficiency stars in a
$10^{11}\,{\rm M}_{\odot}$ halo is approximately $10^{-3}$ (e.g. Mitchell et
al., 2018; Yung et al., 2024), which is an order of magnitude lower that the
efficiency required at $z\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10$.

Another potential challenge is that the observed galaxy sizes tend to be
larger than expected in models for structure formation (Baggen et al., 2023;
Robertson et al., 2023). We argue in §5 that if gas resides in clumpy disks
then angular momentum transfer from the disk to the DM halo can be achieved by
dynamical friction (Chandrasekhar, 1949; El-Zant et al., 2001; Dekel et al.,
2009a).

Numerical studies of individual objects typically rely on the methodology of
zoom-in simulations (e.g Katz et al., 1993; Libeskind et al., 2020; Grand et
al., 2021; Sun et al., 2023b; Pallottini & Ferrara, 2023). In this type of
simulations high resolution is employed only in s small region allowing a
detailed study of its small scale dynamics while simultaneously capturing the
interaction within the larger cosmic environment . In this paper we invoke an
alternative approach of constrained simulations to model the high redshift
$z>9$ accretion history of halos above $M_{\textrm{h}}\approx 10^{11}\,{\rm
M}_{\odot}$. We utilize the Hoffman & Ribak (1991) method of constrained
random realizations to generate initial conditions that are guaranteed to
contain a halo in specified mass range, when evolved forward to a specified
redshift, $z$.

The technique of constrained simulations is very useful in large scale
structure studies, especially for assessing uncertainties in parameter
estimations realistically and in mitigating cosmic variance (Hellwing et al.,
2017). This approach has been used to derive simulation initial conditions
from the observed peculiar velocities (Hoffman et al., 2015) and from the 2MRS
redshift survey (Lilow et al., 2021). Initial conditions based on the galaxy
distribution in the Sloan Digital Sky Survey (SDSS) survey have also been
generated to run a constrained simulation to study the local Universe (Wang et
al., 2016). The same simulation has also been utilised to study the $z\simeq
0$ descendants of galaxies at $z\approx 8-9$ galaxies (Chen et al., 2023).

The structure of the paper is outlined as follows. In §2, we assess the
abundance of halos and compare is with the luminosity function of galaxies at
high redshifts. We underscore the advantage of constrained simulations based
on the expected abundance of the host halos. The simulations are detailed in
§3, which includes the method for generating suitable initial conditions.
Additionally, this section contrasts the accretion history of halos in the
simulations with observational data. In §4, a straightforward model for star
formation is introduced. The model produces a significantly higher ratio of
stellar to halo mass at high redshifts compared to low redshift. Further, in
§5 we present arguments that dynamical friction of galactic disks at high
redshift could lead to a reduction of the effective radii of galaxies by
decelerating the rotation of the disk. A summary and discussion are provided
in §6.

We adopt the standard flat $\Lambda$CDM cosmological model (Planck
Collaboration et al., 2018) with a total mass density parameter
$\Omega_{m}=0.311$, baryonic density $\Omega_{b}=0.049$, Hubble constant
$H_{0}=67.7\,{\rm km}\,{\rm s}^{-1}\,\,\rm Mpc^{-1}$ and normalization
$\sigma_{8}=0.81$.

##  2 Halo abundance at high redshift

Figure 1 displays the abundance of DM halos per logarithmic bins per $\,\rm
Mpc^{3}$ at different redshifts. We employ a halo definition in terms of a
spherical overdensity, where he halo virial radius $r_{\textrm{h}}(t)$, at any
given time $t$ is determined such that the mean density within this radius
equals $200$ times the critical density of the universe,
$\rho_{c}(t)=3H(t)^{2}/8\pi G$. Therefore,

$\displaystyle r_{\textrm{h}}(t)$ $\displaystyle=0.1H(t)^{-1}V_{c}\;,$
$\displaystyle M_{\textrm{h}}(t)$
$\displaystyle=0.1G^{-1}H(t)^{-1}V_{c}^{3}\;,$ (1)

where the circular velocity
$V_{\textrm{c}}=\sqrt{GM_{\textrm{h}}/r_{\textrm{h}}}$. At high redshifts
where $H\sim 1/t$, these relations yield

$\displaystyle r_{\textrm{h}}(t)$ $\displaystyle=13.1\,{\rm
kpc}\left(\frac{11}{1+z}\right)^{3/2}\frac{V_{c}}{181\,{\rm km}\,{\rm
s}^{-1}}\;,$ $\displaystyle M_{\textrm{h}}(t)$ $\displaystyle=10^{11}\,{\rm
M}_{\odot}\left(\frac{11}{1+z}\right)^{3/2}\left(\frac{V_{c}}{181\,{\rm
km}\,{\rm s}^{-1}}\right)^{3}\;,$ (2)

The plots are generated using the mass function outlined by Tinker et al.
(2008) for the $\Lambda$CDM cosmological model, as incorporated within the
COLOSSUS cosmology Python package (Diemer, 2018). The curves corresponding to
halos within the mass range $M_{\textrm{h}}\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10^{11}\,{\rm
M}_{\odot}$ at redshifts $z\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}9$ exhibit comparable abundance to
large groups and clusters at $z=0$. At the upper end of the mass range, the
dependence on mass steepens significantly at such high redshifts. For
instance, at $z=9$, the abundance of halos with $M_{\textrm{h}}=10^{12}\,{\rm
M}_{\odot}$ is four orders of magnitudes lower that than halos with
$M_{\textrm{h}}=10^{11}\,{\rm M}_{\odot}$. The dashed lines represent the
following simple approximation to the mass function

$\frac{{\rm d}n(z,M_{\textrm{h}})}{{\rm d}\log
M_{\textrm{h}}}=\frac{10^{-5}{\,\rm
Mpc}^{-3}}{\left[\left(\frac{M_{\textrm{h}}}{M_{-5}(z)}\right)^{2}+0.25\left(\frac{M_{\textrm{h}}}{M_{{-5}}(z)}\right)^{0.8}(z)\right]^{2}}\;.$
(3)

The approximation is valid in the redshift range $10<z<20$ with $M_{-5}(z)$ is
defined as,

$\log\left(M_{-5}/10^{11}\,{\rm M}_{\odot}\right)=0.185(z-10)^{1.03}\;,$ (4)

and is equal to the halo mass where ${\rm d}n/{\rm d}\log M=10^{-5}\,\rm
Mpc^{-3}$. For $M\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}M_{-5}(z)$, we have the steep
dependence ${\rm d}n/{\rm d}\log M\sim M^{-4}$.

Various fitting formulae for the mass function are available in the literature
(e.g. Press et al., 1974; Sheth & Tormen, 1999; Jenkins et al., 2001; Angulo
et al., 2012; Seppi et al., 2021). These formulas are primarily calibrated
using low redshift outputs from simulations, although exceptions exist (see
Behroozi et al. (2013)). Therefore, it is important to investigate whether
discrepancies among these mass functions might be particularly notable when
applied to high redshifts.

In fig. 2, compare the Tinker et al. (2008) with the widely used Despali et
al. (2016) mass function by examining the ratio of the latter to the former as
a function of halo mass. The depicted curves represent two specific redshift
values, as indicated in the figure. The disparity between the two mass
functions grows with both mass and redshift, albeit remaining within a factor
of a few. This discrepancy remains modest compared to the uncertainty in the
abundance of observed galaxies at $z\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10$. Indeed, we have
conducted numerous analyses using the Despali et al. (2016) mass function and
found negligible changes in the relevant results.

Figure 1: Halo abundance versus mass at different redshifts as denoted in the
figure. The solid lines are derived from the Tinker et al. (2008) expression
for the mass function in the Planck $\Lambda$CDM cosmology (see text). Dashed
lines are obtain from the fit eq. 3. At the high mass end, the curves are
$\sim M_{\textrm{h}}^{-4}$ and thus the cumulative number density
$n(>M_{\textrm{h}})\approx(1/4){\rm d}n/{\rm d}\log M_{\textrm{h}}$
Figure 2:  Ratio of Despali et al. (2016) to Tinker et al. (2008) mass
functions

Galaxies observed at high redshifts are linked with halos having masses
$M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$ (e.g., Boylan-Kolchin, 2023;
Chen et al., 2023; Yung et al., 2024). This connection is established through
abundance matching, where the galaxy UV luminosity function is mapped to the
halo mass function. Specifically, this involves equating the cumulative number
density of galaxies with UV magnitudes above a given value $M_{\textrm{UV}}$
to the number density of halos with masses greater than $M_{\textrm{h}}$.

Harikane et al. (2024) (hereafter H24) constrain the UV luminosity function of
high-redshift galaxies using 25 galaxies with spectroscopic redshifts spanning
$z\approx 8.61-13.20$. Their constraints align with various luminosity
distributions derived from photometric redshifts (e.g., Pérez-González et al.,
2023; Donnan et al., 2023; Harikane et al., 2023; Bouwens et al., 2023). H24
show that the observed UV luminosity distribution can be effectively modeled
by a double power-law function, denoted as $\Phi_{\textrm{UV}}$. We have
verified (not shown here) that this fit also agrees with the observed
luminosity functions in Yung et al. (2024).

Using this double power-law fit, we conduct a straightforward abundance
matching to associate galaxies with a UV magnitude, $M_{\textrm{UV}}$, to
halos of mass $M_{\textrm{h}}$

The results of the abundance matching are summarized in fig. 3 for three
redshift values. The width of the shaded area corresponds to variation in the
normalization of the double power law fit, spanning a factor in between $0.2$
and $5$. While this should provisionally reflect the uncertainty in the
measured UV luminosity function of H24, it is important to note that the
uncertainty in the observed luminosity in a single bin of one UV magnitude
width could be be as large as two orders of magnitude.

At $M_{\textrm{UV}}=-21$, the halo mass is in the range
$M_{\textrm{h}}=6\times 10^{10}-1.5\times 10^{11}\,{\rm M}_{\odot}$ at $z=10$,
with a cumulative halo abundance (indicated by log the number density in the
shaded areas) similar to massive groups and clusters at $z=0$. Although
galaxies with a fixed $M_{\textrm{UV}}$ correspond to lower halo masses as we
move from low high redshifts, the decrease in halo mass is insufficient to
maintain the same abundance. In fact, halos corresponding to a fixed
$M_{\textrm{UV}}$ become rarer.

At $z\approx 10$, on average only a single halo with $M_{\textrm{h}}=2\times
10^{11}\,{\rm M}_{\odot}$ is expected in a box a $100\,\rm Mpc$. While it is
possible to employ zoom-in techniques for simulating massive halos at
$z\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}10$, this necessitates significant computational resources.
Instead, we adopt a computationally friendlier approach, utilizing constrained
random realizations to generate initial conditions that are guaranteed to
contain a massive halo when evolved forward to a specified redshift, $z$.

Figure 3: Halo mass versus UV magnitude obtained by abundance matching of the
Tinker et al. (2008) halo mass function and the power law fit of H24 to UV
luminosity function, $\Phi_{\textrm{UV}}$. Each shaded area corresponds to a
single redshift and is bounded by $M_{\textrm{h}}$ curves obtained using
$5\Phi_{\textrm{UV}}$ (yielding the lower boundary of the shaded area) and by
$0.2\Phi_{\textrm{UV}}$ (leading to the upper boundary). The numbers in each
shaded area represent the log of the cumulative halo abundance (in
$\textrm{Mpc}^{-3}$) expected at the corresponding $M_{\textrm{h}}$.

##  3 Simulations

In 3.1, we outline the method for generating constrained initial conditions.
Using these initial conditions, simulations of dark matter particles were
conducted in periodic boxes using the SWIFT cosmological code (Schaller et
al., 2023). All simulations started at redshift $z_{i}=80$ and concluded at
$z=9$.

We performed nine simulations constrained to have halo masses of approximately
$M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$ in boxes of $L=17.1\,\rm
Mpc$, along with one unconstrained simulation in a box of the same size.
Additionally, two simulations were conducted with initial conditions
constrained to include a halo mass of $M_{\textrm{h}}\approx 10^{12}\,{\rm
M}_{\odot}$ at $z=9$, positioned at the center of cubic boxes of $L=36.9\,\rm
Mpc$.

Each simulation consisted of $512^{3}$ equal-mass particles, resulting in
particle masses of $10^{7}\,{\rm M}_{\odot}$ and $10^{6}\,{\rm M}_{\odot}$ in
the large and small boxes, respectively. The maximum physical softening used
in the simulations was $100\ {\rm pc}$. Throughout each simulation run, the
output of positions and velocities of all particles was retained at 18
different redshifts spanning from $z=20$ to $z=9$. Halos were identified from
the simulation outputs utilizing the VELOCIraptor halo finder (Elahi et al.,
2019). This halo finder provides halo masses according to several definitions.
Here we use VELOCIraptor masses that match the definition in section 2.

###  3.1 Constrained Initial Conditions

We formulate the condition (constraint) for the presence of a halo of a given
mass as follows. Let $\delta(\mathbf{x},z)$ be the linearly evolved density at
redshift $z$ and $\delta_{R}(\mathbf{x},z)$ be the convolution of $\delta$
with a top-hat window of comoving radius $R$. We associate a halo of mass
$M_{\textrm{h}}$ with a comoving Lagrangian radius
$R_{L}=(3M_{\textrm{h}}/4\pi\rho_{\textrm{m}})^{1/3}$, where
$\rho_{\textrm{m}}$ is the background density in comoving coordinates.

The formation of a halo of mass $M_{\textrm{h}}$ at redshift $z$, located at
position at $\mathbf{x}_{0}$ is determined by the condition
$\delta_{R=R_{L}}(\mathbf{x}_{0},z)=\delta_{c}\simeq 1.68$ (e.g. Press &
Schechter, 1974; Peebles, 1980). Here, $\delta_{c}$ serves as a critical
threshold indicative of the collapse or virialization of dark matter halos. At
$z=10$, halos with masses $M_{\textrm{h}}=10^{11}\,{\rm M}_{\odot}$ and
$M_{\textrm{h}}=10^{12}\,{\rm M}_{\odot}$ correspond to Lagrangian comoving
radii $R_{L}=0.84,\textrm{Mpc}$ and $1.82,\textrm{Mpc}$, respectively. In
terms of the ratio $\delta_{c}/\sigma_{R_{L}}$, these values equate to $4.8$
and $6.5$. This ratio provides a measure of the likelihood of halo formation
under the specified conditions.

This formulation is only approximate as the superposition of generic
fluctuations on all scales and non-linear evolution will lead to deviations
from the desired halo mass and position. Nonetheless, the prescription is
reasonable for massive (rare) halos (Robertson et al., 2009; Ludlow et al.,
2014), as is the case in the current study.

We adopt the methodology of Hoffman & Ribak (1991) to generate gaussian random
fields that satisfy the aforementioned condition. This elegant method
expresses the constrained random field, denoted by $\delta_{1}$, in terms of
an unconstrained random gaussian field, $\delta^{\textrm{unc}}$, as follows

$\delta_{1}(\mathbf{x},z)=\delta^{\textrm{unc}}(\mathbf{x},z)-\langle\delta(\mathbf{x},z)|{\delta_{R_{L}}^{\textrm{unc}}}\rangle+\langle\delta(\mathbf{x},z)|{\delta_{c}}\rangle\;$
(5)

where ${\delta_{R_{L}}^{\textrm{unc}}}$ is the value of the filtered
unconstrained field at $\mathbf{x}_{0}$.

The ensemble average of all $\delta$ fields satisfying the constraint
$\delta_{R_{L}}(\mathbf{x}_{0},z)=C$ is given by

$\langle\delta(\mathbf{x})|C\rangle=\zeta(|\mathbf{x}-\mathbf{x}_{0}|)\frac{C}{\sigma_{R_{L}}^{2}}\;,$
(6)

where $\sigma^{2}_{R_{L}}$ is the variance of the smoothed density,
$\delta_{R_{L}}$.

Figure 4: Projected 2D density field computed from one of the constrained
simulations with $L=17.1\,\rm Mpc$ at six distinct redshifts. Circles denote
halos classified by mass ranges: black ($M_{\textrm{h}}/\,{\rm
M}_{\odot}>10^{11}$), white ($5\times 10^{10}-10^{11}$), blue
($10^{10}-5\times 10^{10}$), red ($5\times 10^{9}-10^{10}$), and yellow
($10^{9}-5\times 10^{9}$), with size logarithmically scaled by mass. The color
bar indicates density in units relative to the mean density within the
simulation box.
Figure 5: Halo abundance computed from the simulations output at $z=9$. The
blue and orange histograms represent the constrained simulations with halo
masses of $10^{11},M_{\odot}$ and $10^{12},M_{\odot}$, respectively. The green
curve depicts the Tinker et al. (2008) mass function.

###  3.2 Results

We begin by visually inspecting the projected (2D) density field in one of the
constrained simulations. Figure 4 illustrates the 2D density (in units of the
mean 2D density value) for one of the $M_{\textrm{h}}=10^{11}\,{\rm
M}_{\odot}$ constrained simulations. The overlaid circles represent identified
halos in different mass ranges, as described in the figure caption. Only halos
above $10^{9}\,{\rm M}_{\odot}$ are marked.

The displayed region of the box is focused around the center. As expected, the
most massive halo (MMH) forms near the center. Additional massive halos are
associated with the growth of the MMH, however, their masses are significantly
lower than the MMH. Halos in the mass range $5\times
10^{9}<M_{\textrm{h}}/\,{\rm M}_{\odot}<10^{10}$ (red circles) are present at
$z\simeq 12$. These halos are not considered rare, as their expected abundance
is $10^{-2}-10^{-3}\textrm{dex}^{-1}\,\rm Mpc^{-3}$ (see fig. 1). Therefore,
we expect to identify a few such halos even in our simulation box of
$L=17.1\,\rm Mpc$.

Inspection of the panels at $z=11$ and $z=10$ reveals a major merger event
with the two halos marked by two blue circles ($10^{10}<M_{\textrm{h}}/\,{\rm
M}_{\odot}<5\times 10^{10}$) at $z=11$ merging to form the larger halo
indicated by the white circle ($5\times 10^{10}<M_{\textrm{h}}/\,{\rm
M}_{\odot}<10^{11}$) at $z=10$. There is a clear tendency of increasing halo
mass as we move nearer to the MMH, particularly at the lowest redshifts.

Figure 5 plots the abundance of halos from the simulations as a function of
$M_{\textrm{h}}$. At the high mass end, the simulated abundance is
significantly above the predictions from the Tinker et al. (2008) fitting
formula (green curve). In both $M_{\textrm{h}}=10^{11}\,{\rm M}_{\odot}$ and
$M_{\textrm{h}}=10^{12}\,{\rm M}_{\odot}$ constrained simulations, the MMH is
associated with other relatively rare massive halos with abundance well above
the green curve.

####  3.2.1 Halo Accretion History

The halo accretion history is directly linked to the availability of gas for
star formation. Newly accreted gas replenishes the reservoir, which is
subsequently converted into stars and may escape the galaxy through processes
like supernovae (SN) and active galactic nuclei (AGN) feedback.

Figure 6 displays curves of the accretion histories of the MMH in all
simulations. The top, middle and bottom panels correspond, respectively, to
the halo mass versus time $M_{\textrm{h}}(t)$, the accretion rate
${\dot{M}}_{\textrm{h}}$, and the specific accretion rate
${\dot{M}}_{\textrm{h}}/M_{\textrm{h}}$. Squares and circles represent
quantities inferred from observed data, as detailed below.

The curves of $M_{\textrm{h}}(t)$ (top panel) reveal that half of the
simulated halos, including the unconstrained halo (dotted), have acquired 80%
of their final masses in the last $150\,{\rm Myr}$. Only the dashed curves,
corresponding to constrained simulations with $M_{\textrm{h}}\approx
10^{12}\,{\rm M}_{\odot}$ and one of the nine solid curves, have acquired a
mass $\gtrsim 10^{10}\,{\rm M}_{\odot}$ by $z=15$. Significant variations
exist between individual halos, evident from the curves of
${\dot{M}}_{\textrm{h}}$ (middle panel), where some halos exhibit fluctuating
$M_{\textrm{h}}$ while others (e.g., dashed and a few solid curves) display
smoother behavior, albeit with fluctuations at a $~{}\hbox
to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}100\,{\rm Myr}$ scale.

The cyan circles in the top panel represent halo masses inferred via abundance
matching from the observed UV magnitudes, $M_{\textrm{UV}}$, of galaxies
listed in table 1 of H24. We refrain from attaching error bars to these
estimates due to uncertainties in abundance matching, which are relatively
small compared to the scatter of the points. The cyan points fall within the
range of solid curves corresponding to simulations constrained to contain a
$M_{\textrm{h}}\approx 10^{11}$ halo. While this is expected given that the
choice of this mass motivated by the observed abundance of galaxies, the
result is not entirely trivial. At a given redshift, the spread in halo masses
between different simulations is large, ranging from $M_{\textrm{h}}\approx
10^{10}\,{\rm M}_{\odot}$ to $10^{11}\,{\rm M}_{\odot}$ at $z=10$, as
indicated by the solid curves. This mass range is associated with more than a
two-order-of-magnitude difference in the abundance of halos, as shown in fig.
1.

Using SED fitting, H24 provide SFRs and stellar masses for a subset of the
galaxies in table 1 of their paper. The magenta squares in the top panel are
halo masses estimated from these stellar masses (table 3 in H24) assuming

$M_{\textrm{h}}=f_{t}^{-1}M_{\ast}\;,$ (7)

where $f_{t}=f_{*}f_{b}$ is assigned the value $f_{t}=0.02$ and the error bars
in the figure are estimated from the uncertainties on $M_{\ast}$. The open
square symbols refer to the galaxy GS-z12 for which different values of
redshift, SFR and $M_{\ast}$ are reported in D’Eugenio et al. (2023) and H24.
Both sets of data are shown with D’Eugenio et al. (2023) being the point with
the lower redshift ($z=12.43$).

At $z\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}10$, halo masses inferred from observed stellar masses
(magenta squares) by eq. 7 are consistent with the constrained simulations
with $M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$ and the results from
abundance matching (cyan points). However at $z\simeq 9$, the estimate in eq.
7 falls below the simulation and abundance matching. Decreasing the parameter
$f_{t}$ by a factor of $\approx 5-10$ would bring the $z\simeq 9$ galaxies to
agree with the simulations and abundance matching. This is in line with the
finding of Yung et al. (2024) that the luminosity function predicted by their
semi-analytic models is consistent the observations well at $z~{}\hbox
to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}9$, while it underestimates the
observed abundances at $z\approx 11$.

The UV magnitudes, $M_{\textrm{UV}}$, of galaxies in table 1 of H24 can also
be used to deduce the SFRs using ${\dot{M}}_{\ast}=1.15\times 10^{-28}L_{UV}$.
The time derivative of eq. 7 relates the SFR, $\dot{M}_{\ast}$, to the halo
accretion rate as $\dot{M}_{\textrm{h}}=\dot{M}_{\ast}/f_{t}$. The results for
$f_{t}=0.02$ are plotted as cyan points in the middle panel.

In addition to $M_{\textrm{UV}}$, the SED-fitted SFRs in table 3 of H24 are
used here to obtain ${\dot{M}}_{\textrm{h}}$, also for $f_{t}=0.02$. The
results are represented as the magenta squares in the middle panel. The UV
magnitude and SED-based value of ${\dot{M}}_{\textrm{h}}$ agree with each
other, but are not identical, as the SED fitting involves a more detailed
modeling of the SFR than the UV magnitudes alone. Curves of accretion rates
from the $M_{\textrm{h}}=10^{11}\,{\rm M}_{\odot}$ constrained simulations,
are consistent with the observations via $M_{\textrm{UV}}$ (cyan) and SED-SFRs
(magenta). Compared to the top panel, halo masses inferred from the observed
stellar masses (magenta square) are below those predicted through abundance
matching (cyan). It is important to note that while both the cyan and magenta
points in the top and middle panels rely on $M_{\textrm{UV}}$, their
underlying methodologies differ - abundance matching for the top panel versus
an empirical SFR-$M_{\textrm{UV}}$ relationship for the middle panel.

The results at $z\simeq 9$ are intriguing. While the choice of $f_{t}=0.02$
for inferring $M_{\textrm{h}}$ from $M_{\ast}$ (top panel) and
${\dot{M}}_{\textrm{h}}$ from ${\dot{M}}_{\ast}$ (middle panel) matches well
the $M_{\textrm{h}}=10^{11}\,{\rm M}_{\odot}$ simulations at $z\gtrsim 10$, a
lower $f_{t}$ is needed for consistency at $z\simeq 9$. However, the
deviations in $M_{\textrm{h}}$ derived from simulations at $z\simeq 9$, when
compared to observational inferences, are about $2-4$ times larger than those
seen in ${\dot{M}}_{\textrm{h}}$ at the same redshift. This discrepancy is
clearly illustrated in the specific accretion rates shown in the bottom panel.

The bottom panel shows the specific halo accretion rates,
${\dot{M}}_{\textrm{h}}/M_{\textrm{h}}$, from simulations. These rates hover
around a simple fit denoted by a black line, represented by the equation
${{\rm d}{\rm ln}M_{\textrm{h}}}/{{\rm d}t}=3.15t_{\textrm{Gyr}}^{-4/3}\,{\rm
Gyr}^{-1}$. This fit approximates the mean accretion rate for halos of mass
$M_{\textrm{h}}=10^{11}\,{\rm M}_{\odot}$ as proposed by Fakhouri et al.
(2010).

The magenta squares in the bottom panel indicate the sSFR of galaxies as
catalogued in table 3 of H24. There is a reasonable match between the observed
sSFRs (magenta squares) and the specific halo accretion rates from the
simulations at $z>10$. This consistency with the simulations corroborates the
findings of Bouwens et al. (2023), who noted that the sSFRs tend to follow the
scaling $(1+z)^{2.5}$ proposed by Fakhouri et al. (2010) for the specific halo
accretion rate at high redshift. However, it is important to note that the
normalization of the specific halo accretion rate in Bouwens et al. (2023) is
higher by a factor of a few compared to the fit by Fakhouri et al. (2010).

According to equation eq. 7, the specific accretion rate is equal to the sSFR,
i.e. ${\dot{M}}_{\textrm{h}}/M_{\textrm{h}}=\dot{M}_{\ast}/M_{\ast}$,
independently of $f_{t}$. Note that $M_{\ast}\propto M_{\textrm{h}}^{\alpha}$
($\alpha=const$), then
${\dot{M}}_{\textrm{h}}/M_{\textrm{h}}=\alpha^{-1}\dot{M}_{\ast}/M_{\ast}$
(Behroozi & Silk, 2015). Therefore, also for $\alpha\neq 1$ curves of the sSFR
and the specific halo accretion rate should trace each other, with a constant
ratio between them.

Figure 6: Accretion histories of the main halos in the simulations. The dotted
line represents the unconstrained simulation, while solid and dashed lines
correspond to constrained simulations with $M_{h}=10^{11}\,{\rm M}_{\odot}$
and $10^{12}\,{\rm M}_{\odot}$, respectively. Top: mass as a function of time.
Magenta squares are based on the observed $M_{\ast}$ using eq. 7. Cyan circles
represent $M_{\textrm{h}}$ obtained by abundance matching from observed
$M_{\textrm{UV}}$. Middle: The accretion rate $\dot{M}_{h}$. Squares represent
$\dot{M}_{h}=f_{t}^{-1}\dot{M}_{\ast}$ derived from observed $\dot{M}_{\ast}$.
The cyan circles show results obtained from the $\dot{M}_{\ast}$ estimated
from UV magnitudes directly. Bottom: Squares show the sSFR,
$\dot{M}_{\ast}/M_{\ast}$, of observed galaxies. Open squares refer to the
galaxy GS-z12 (see text). The black curve is a fit $\sim t^{-4/3}$.

##  4 A simplified Model for Star Formation

We present a straightforward model that offers insights into certain aspects
of galaxy formation at high redshifts. Our primary argument is that star
formation processes at both high and low redshifts may not significantly
differ. Given a halo with of mass $M_{\textrm{h}}(t_{z})$ at a time $t_{z}$
corresponding to a redshift $z$, we develop a model to calculate the stellar
mass, $M_{\ast}$, and the the SFR $\dot{M}_{\ast}$, within this halo. The main
assumption of our model is that star formation at high redshift occurs in a
rotationally supported disk, governed by the Schmidt-Kennicutt (SK) law, as
deduced from observations in the local universe (Schmidt, 1959; Kennicutt,
1998). At a specific time $t$ and for a disk radius $R$, let $\Sigma_{g}(t,R)$
and $\Sigma_{\ast}(t,R)$ denote the surface densities of the disk’s gaseous
and stellar components, respectively. We express the SK law as follows

${\dot{\Sigma}}_{\ast}=A\Sigma_{g}^{n}\;,$ (8)

where ${\dot{\Sigma}}_{\ast}$ is the SFR per unit area of the disk. We adopt
the parameters, $n=1.54$ and $A=10^{-3.95}$ (Kennicutt & De Los Reyes, 2021)
where $\Sigma_{\textrm{SFR}}$ and $\Sigma_{g}$ are, respectively in units of
$\textrm{M}_{\odot}\,\textrm{pc}^{-2}$ and
$\textrm{M}_{\odot}\,\textrm{pc}^{-2}$.

The basic picture is as follows. The gas already present in the disk undergoes
partial conversion into stars following the SK law, accompanied by subsequent
stellar feedback (SN and winds) that expels a fraction of the gas mass.
Additionally, the reservoir of gas in the disk is replenished and expanded as
the halo accretes additional matter.

The application of the SK law requires the gas surface density
$\Sigma_{g}(t,R)$ as a function of time and radius. This is expressed through
the mass balance equation,

${\dot{\Sigma}}_{g}={\dot{\Sigma}}_{\textrm{acc}}-{\dot{\Sigma}}_{\ast}-{\dot{\Sigma}}_{\textrm{ej}}\;,$
(9)

where ${\dot{\Sigma}}_{\textrm{acc}}$ denotes the added baryonic mass to the
disk surface density due to accretion onto the halo, ${\dot{\Sigma}}_{\ast}$
is the rate of gas mass converted into star, and
${\dot{\Sigma}}_{\textrm{ej}}$ represents the gas ejected from disk by stellar
feedback. We restrict the model to halos of mass $~{}\hbox
to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}10^{12}M_{\odot}$, where AGN
feedback is considered subdominant (e.g. Croton et al., 2006; Bower et al.,
2006; Puchwein & Springel, 2013).

The ejected mass rate is expressed in terms of the circular velocity of the
host halo as (e.g. White & Frenk, 1991b; Kauffmann et al., 1993; Mitchell et
al., 2018; Yung et al., 2019)

${\dot{\Sigma}}_{\textrm{ej}}=\left(\frac{V_{c}}{V_{SN}}\right)^{-\gamma}{\dot{\Sigma}}_{\ast}\;.$
(10)

This relation implies that stellar feedback is more effective in halos with a
shallow gravitational potential (Larson, 1974; Dekel & Silk, 1986). Mitchell
et al. (2018) adopt the $V_{SN}=380\,{\rm km}\,{\rm s}^{-1}$ and $\gamma=3.2$
while Yung et al. (2019) use the parameters $V_{SN}=240\,{\rm km}\,{\rm
s}^{-1}$ and $\gamma=2.8$, resulting in a comparatively weaker feedback.

We now detail our approach for estimating $\Sigma_{\textrm{acc}}(t,R)$. Given
the halo mass $M_{\textrm{h}}(t)$ as a function of time, and consequently,
through section 2, the halo virial radius $r_{\textrm{h}}(t)$ and circular
velocity ${V}_{c}(t)$, we determine the gas mass accreted within a short time
interval $\delta t$ as,

$\delta M_{\textrm{acc}}=f_{\textrm{b}}{\dot{M}}_{\textrm{h}}\delta t\;.$ (11)

Due to the short crossing time in high-redshift halos (e.g.,
$t_{\textrm{cross}}\approx r_{\textrm{h}}/{V}_{c}=0.1/H=70,\textrm{Myr}$ at
$z=10$) and efficient gas cooling (Dekel et al., 2023; Yung et al., 2024), we
assume that $\delta{M}_{\textrm{acc}}$ quickly settles into the disk, with a
surface density distribution given by,

$\delta\Sigma_{\textrm{acc}}(t,R)=\delta\Sigma_{0}e^{-R/R_{d}(t)}\;,$ (12)

where $R_{\textrm{d}}(t)$ could be a general function of $t$, but we rely here
on the approximation (Mo et al., 1998; Yang et al., 2023),

$R_{\textrm{d}}=0.7\lambda_{\textrm{B}}r_{\textrm{h}}\;,$ (13)

where $\lambda_{\textrm{B}}$ is the halo spin parameter $\lambda_{\textrm{B}}$
(Bullock et al., 2001). The constant $\delta\Sigma_{0}$ is determined by
equating $\delta M_{\textrm{acc}}$ to the total added mass to the disk as
derived by integrating over eq. 12. This yields,

$\delta M_{\textrm{acc}}=2\pi\delta\Sigma_{0}R_{\textrm{d}}^{2}(t)\;.$ (14)

Thus, from eqs. 11 to 14, we obtain,

${\dot{\Sigma}}_{\textrm{acc}}(t,R)=\frac{f_{\textrm{b}}\dot{M}_{\textrm{h}}}{2\pi
R_{\textrm{d}}^{2}(t)}e^{-R/R_{\textrm{d}}(t)}\;.$ (15)

For the halo growth rate, we assume
${\dot{M}}_{\textrm{h}}/M_{\textrm{h}}\propto t^{-4/3}$ as represented by the
smooth black curve in fig. 6, yielding the form,

$M_{\textrm{h}}(t)=M_{\textrm{h}}(t_{z})\mathrm{e}^{A(t_{z}^{-\beta}-t^{-\beta})}\;,$
(16)

where $A=-8.1$, $\beta=1/3$ and $t$ is in units of Gyr. Although we only
present results for this particular choice of $M_{\textrm{h}}(t)$, we have
also tested ${\dot{M}}_{\textrm{h}}/M_{\textrm{h}}\propto t^{-5/3}$, similar
to (Fakhouri et al., 2010) , and found no significant differences in the SHMR.

We perform a numerical integration of eqs. 8 and 15 from an initial time
$t=t_{i}\ll t_{z}$ to $t=t_{z}$. We set the initial conditions as
$\Sigma_{\ast}(t_{i},R)=10^{-3}\Sigma_{g}(t_{i},R)$, where
$\Sigma_{g}(t_{i},R)$ follows an exponential disk profile. The scale length of
the disk at $t_{i}$, $R_{\textrm{d}}(t_{i})$, is determined according to eq.
13 with $r_{\textrm{h}}=r_{\textrm{h}}(t_{i})$. The total initial gas mass in
the disk is equal to $f_{\textrm{b}}M_{\textrm{h}}(t_{i})$ minus the initial
mass in stars. In all numerical calculations, we maintain
$\lambda_{\textrm{B}}=0.035$ in eq. 13.

In fig. 7, the SHMR is plotted against halo mass for three redshift values,
$z$. The model incorporating the feedback parameters of Yung et al. (2019),
depicted with solid curves, agrees reasonably well with the observed local
SHMR as estimated through abundance matching techniques, according to studies
such as (e.g. Moster et al., 2013; Rodríguez-Puebla et al., 2017; Moster et
al., 2020; Shuntov et al., 2022; Girelli et al., 2020). The dashed curve for
$z=0$, which represents the stronger feedback with parameters from (Mitchell
et al., 2018), shows that our model underestimates the local SHMR by nearly an
order of magnitude at a halo mass $\approx 10^{11}M_{\odot}$.

The model SHMR acquires larger values at higher redshifts for a given mass.
This is due to the relation $V_{c}\sim M_{\textrm{h}}/t$, indicating that a
higher circular velocity $V_{c}$ occurs at earlier times for a fixed
$M_{\textrm{h}}$, leading to less efficient SN feedback and, consequently, a
larger gas reservoir for star formation. For the solid curves, the SHMR for
$M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$ increases by approximately an
order of magnitude at high redshifts compared to $z=0$, while for
$10^{10}\,{\rm M}_{\odot}$, the increase is about a factor of 25. In contrast,
predictions from the UniverseMachine (Behroozi et al., 2020) suggest that the
SHMR increases by about a factor of 10 from $z=0$ to $z=12$ for halos with
$10^{10}\,{\rm M}_{\odot}$ (their figure 12).

Figure 7: The SHMR from the simplified model at three different redshifts, as
indicated in the figure. Dashed and solid curves correspond to strong and weak
feedback recipes, respectively.

##  5 Sizes

The compactness of high-redshift galaxies have been identified as a potential
challenge for theoretical models. Observational estimates for the half-light
radius, $r_{\textrm{e}}$, are available for some galaxies at high redshifts.
Baggen et al. (2023) measure $r_{\textrm{e}}\approx 80-300\;\textrm{pc}$ for
nine galaxies at $z\approx 6-9$, while Robertson et al. (2023) provide the
estimates $r_{\textrm{e}}\approx 50-165\;\textrm{pc}$ for four galaxies at
$z\approx 10-13$. As highlighted by Baggen et al. (2023), these value are
below the prediction of the numerical simulation by (Marshall et al., 2022),
across the redshift range $z=7$ to $z=11$. However note that the gravitational
softening of these simulations is substantial, about $240\ {\rm pc}$ at $z=8$.

As a rough estimate of the effective radius of the stellar disk
$r_{\textrm{e}}$, we use $r_{\textrm{e}}\simeq 1.7R_{\textrm{d}}$, where
$R_{\textrm{d}}$ is taken from eq. 13. This yields

$r_{\textrm{e}}=550\ {\rm
pc}\frac{\lambda_{\textrm{B}}}{0.035}\left(\frac{11}{1+z}\right)\left(\frac{M}{10^{11}\,{\rm
M}_{\odot}}\right)^{1/3}\;.$ (17)

Adopting a normal distribution for ${\rm ln}\lambda_{\textrm{B}}$, where a
mean expressed as ${\rm ln}\lambda_{0}$ with $\lambda_{0}=0.035$ and an rms of
$\sigma=0.55$, we find that the observational estimates are within $\approx
2-3\sigma$ of eq. 17.

The effective radii derived from eq. 17, including the scatter in
$\lambda_{\textrm{B}}$, are broadly consistent the findings of (Marshall et
al., 2022). However, these authors do not present results as a function of
halo mass, preventing a direct comparison.

A clumpy disk at high redshifts may actually be subject to dynamical friction,
leading to smaller disk sizes. The role of dynamical friction in the simulated
galaxies of Marshall et al. (2022) remains uncertain, due to potential
limitations in gravitational force resolution that may not sufficiently
capture the effects of dynamical friction. The dynamical timescale of a clump
of mass $M_{1}$ at a distance $r$ from the center of the halo is
$t_{\textrm{df}}\approx({\rm ln}\Lambda)^{-1}M(r)/M_{1}r/V_{c}$, where $M(r)$
is the mass of the halo contained within $r$ (Binney & Tremaine, 2008).
Consider a clump of mass $M_{1}$ at a distance $r=R_{\textrm{d}}$. For an
isothermal sphere
$M(R_{\textrm{d}})=M_{\textrm{h}}(R_{\textrm{d}}/r_{\textrm{h}})=\lambda_{\textrm{B}}M_{\textrm{h}}$
and $V_{\textrm{c}}=\sqrt{GM_{\textrm{h}}/r_{\textrm{h}}}$, we find

$t_{\textrm{df}}=\left(\frac{9}{{\rm
ln}\Lambda}\right)\left(\frac{M_{\textrm{h}}/10^{11}}{M_{1}/10^{7}}\right)\left(\frac{\lambda_{\textrm{B}}}{0.035}\right)^{2}\left(\frac{11}{1+z}\right)^{3/2}\,95\,{\rm
Myr}\;,$ (18)

where we have used section 2, ${\rm ln}\Lambda=9$ corresponds to
$M_{\textrm{1}}=10^{7}\,{\rm M}_{\odot}$ and $M_{\textrm{h}}=10^{11}\,{\rm
M}_{\odot}$ (yielding a maximum impact parameter
$b_{\textrm{max}}=r_{\textrm{h}}=13\,{\rm kpc}$ at $z=10$). Therefore,
dynamical friction could play an important role over the relevant dynamical
timescales of high redshift galaxies.

##  6 summary and discussion

The approach of constrained simulations is highly beneficial for studying rare
cosmological structures. Here we have only conducted simulations with constant
resolution across the entire simulation box. However, a combination of zoom-in
techniques and constrained initial conditions is most appropriate for
resolving rare structures as well as capturing the gravitational influence of
the large scale environment.

Our simulations indicate the the average accretion history of rare halos is
similar to form found for average halos in the literature. Furthermore, the
accretion curves of halos show both long-term fluctuations on a $100\,{\rm
Myr}$ scale and short-term variations, potentially leading to variability
(stochasticity) in the associated star formation rate.

In fig. 6 we have seen that a dividing the observed SFR by a factor
$f_{t}=0.02$, leads to ${\dot{M}}_{\textrm{h}}$ that are consistent with the
simulations constrained to include a halo of mass $M_{\textrm{h}}\approx
10^{11}\,{\rm M}_{\odot}$. The agreement spans the entire considered redshift
range, $z>9$. Interestingly, dividing the observed stellar masses by the same
$f_{t}$ yields a good match with the halo masses in the simulations, but this
is only true for $z\gtrsim 10$. For galaxies at $z\simeq 9$, the factor
required is $\approx 0.004-0.02$, closer to what is seen in low redshift
galaxies.

This peculiar behaviour can stem from challenges in accurately estimating the
stellar masses. Indeed, SFRs estimated from UV magnitudes are generally more
reliable than stellar mass estimates, which require assumptions about the
entire star formation history (e.g. Whitler et al., 2023). Nonetheless, in the
estimating the SFRs, the impact of potential uncertainties due to potential
dust attenuation needs to be assessed (e.g. Kennicutt et al., 2012). However,
dust attenuation is expected to be small in these high redshift galaxies (e.g.
Bouwens et al., 2023), suggesting that the observed discrepancies could
primarily be due to the complexities of stellar mass estimation.

Another possibility for this behaviour is an abrupt change in the conditions
for star formation at $z\approx 9-10$, similar to the suggestion of Silk et
al. (2024) although their model refers to transition at $z\approx 6$.

There is a wealth of studies investigating the SHMR,
$M_{\ast}/M_{\textrm{h}}$, against the halo mass, $M_{\textrm{h}}$, in the
literature, at redshifts $z~{}\hbox to0.0pt{$<$\hss}{\lower
4.30554pt\hbox{$\sim$}}10$ using pre-JWST data (e.g. Guo et al., 2010; Yang et
al., 2012; Leauthaud et al., 2012; Moster et al., 2013; Rodríguez-Puebla et
al., 2017; Behroozi et al., 2013; Girelli et al., 2020). These studies rely on
mappings between observed stellar mass function and a halo mass function, with
the latter being motivated and calibrated through N-body simulations. The two-
point correlation function could also be invoked to further constrain the SHMR
(e.g Rodríguez-Puebla et al., 2013). Despite important differences among these
studies, they all agree on basic features. The SHMR shows a peak at
$M_{\textrm{h}}\approx 10^{12}\,{\rm M}_{\odot}$ at $z\approx 0$. The peak
tends to shift towards lower halo masses at high redshifts
($M_{\textrm{h}}\approx 2\times 10^{11}\,{\rm M}_{\odot}$ at $z\approx 7$).
The peak SHMR does not exceeds $0.05$ across $z~{}\hbox
to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}10$ (Girelli et al., 2020).
Moving from the peak to lower masses, the SHMR declined rapidly from a few
percents at the peak to a few times $10^{-3}$ at $M_{\textrm{h}}\approx
10^{11}\,{\rm M}_{\odot}$ at $z=0$. Within the observational uncertainties,
most studies suggest little evolution of the SHMR for redshifts below
$z\approx 2$ (Rodríguez-Puebla et al., 2017), although Behroozi et al. (2013)
find a decrease in the SHMR, from $z\approx 0$ to $z\approx 2$, for halo
masses below. Moreover, for halos with $M_{\textrm{h}}\approx 10^{11}\,{\rm
M}_{\odot}$, these authors find a significant increase of the SHMR from low to
high redshifts, reaching approximately $10^{-2}$ at $z\approx 8$.

The default implementations of semi-analytic models of galaxy formation (SAMs)
imply modest changes in the SHMR with redshift for halo masses below the peak
of the SHMR, where stellar dominates AGN feedback. The Durham SAMs predict a
nearly non-evolving SHMR, for $M_{\textrm{h}}~{}\hbox to0.0pt{$<$\hss}{\lower
4.30554pt\hbox{$\sim$}}10^{12}\,{\rm M}_{\odot}$ out to $z\approx 4$, with
$M_{\ast}/M_{\textrm{h}}\approx 3\times 10^{-3}$ at $M_{\textrm{h}}\approx
10^{11}\,{\rm M}_{\odot}$ (Mitchell et al., 2018). Furthermore, the Santa Cruz
models, predict very little evolution of the SHMR for $M_{\textrm{h}}~{}\hbox
to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}10^{11}\,{\rm M}_{\odot}$
between $z=8$ and $z=17$, with $M_{\ast}/M_{\textrm{h}}\approx 10^{-2}$ at
$M_{\textrm{h}}\approx 10^{11}\,{\rm M}_{\odot}$. This trend of minimal
evolution continues across the lower redshift range ($z~{}\hbox
to0.0pt{$<$\hss}{\lower 4.30554pt\hbox{$\sim$}}6$), where the SHMR remains
close to its value at higher redshifts according to the default implementation
of the Santa Cruz SAMs (Somerville et al., 2015). However, comparing between
the SHMR in figure 8 of (Somerville et al., 2015) to that in figure B2 in
(Yung et al., 2024), reveals differences at the low mass end, specifically for
$M_{\textrm{h}}~{}\hbox to0.0pt{$<$\hss}{\lower
4.30554pt\hbox{$\sim$}}10^{10}\,{\rm M}_{\odot}$.

The question of why SAMs of galaxy formation yield a nearly invariant SHMR for
halo masses below $\approx 10^{11}M_{\odot}$, even when SN feedback is
considered less efficient at $z\gtrsim 10$ due to higher values of $V_{c}$ at
those epochs, is intriguing. In these models, gas expelled from the disk by
stellar feedback is typically assumed to join the hot component of the
circumgalactic medium (CGM). This gas eventually cools on a few dynamical
timescales (e.g. Faerman et al., 2020) and settles into the disk to be
reincorporated into the star formation process. In contrast, the simplified
model presented here assumes that the gas ejected by SN feedback is not
returned to the halo and never becomes available for star formation again.
However, this process is sensitive to the circular velocity which for there
same halo mass, is more

In the model outlined in §4, we have considered the circular velocity of the
halo, $V_{c}$, as the parameter governing SN feedback. However, it can be
argued that the maximum circular velocity (the peak of the rotation curve),
$V_{\textrm{max}}$, is more relevant since it reflects the depth of the
gravitational potential of the halo.

In the regime of stable clustering, $V_{\textrm{max}}$ is expected to remain
constant over long cosmic epochs. Thus, employing $V_{\textrm{max}}$ instead
of $V_{c}$ in the model should result in a constant SHMR over extended
periods. However, we are examining halos of the same mass at different epochs.
Hence, the question is whether halos within the relevant mass range,
approximately $10^{11}\,{\rm M}_{\odot}$, also experience moderate mass growth
over cosmic time scales.

The study by (Rodríguez-Puebla et al., 2016) on the median growth of halos in
MultiDark N-body simulations indicates that indeed, halos with a mass of
$10^{11}\textrm{ M}{\odot}$ at redshift $z=0$ have not significantly increased
in mass since $z=2$. For these halos, the median $V_{\textrm{max}}$ remains
notably constant up to redshifts of $z=3-4$ (refer to figures 5&6 in their
paper). However, the scenario shifts as we approach $z=10$. The mass scale of
$10^{11}\textrm{ M}{\odot}$ at $z=10$ corresponds to the median mass of a halo
of $10^{15}\textrm{ M}{\odot}$ at $z=0$. At $z=10$, this median halo possesses
a $V_{\textrm{max}}$ of approximately $200\textrm{ km/s}$.

##  7 Data availability

No new data were generated or analysed in support of this research.

##  8 acknowledgements

The author has benefited from fruitful conversations with Andrew Benson, Enzo
Branchini, Stephane Charlot, Avishai Dekel and Joe Silk. This research is
supported by a grant from the Israeli Science Foundation and a grant from the
Asher Space Research Institute. The research in this paper made use of the
SWIFT open-source simulation code (http://www.swiftsim.com, Schaller et al.
(2018)) version 1.0.0.

## References

  * Angulo et al. (2012) Angulo, R. E., Springel, V., White, S. D., et al. 2012, MNRAS, 426, 2046, doi: 10.1111/j.1365-2966.2012.21830.x
  * Asgari et al. (2021) Asgari, M., Lin, C. A., Joachimi, B., et al. 2021, A&A, 645, A104, doi: 10.1051/0004-6361/202039070
  * Baggen et al. (2023) Baggen, J. F. W., van Dokkum, P., Labbé, I., et al. 2023, ApJL, 955, L12, doi: 10.3847/2041-8213/acf5ef
  * Behroozi et al. (2020) Behroozi, P., Conroy, C., Wechsler, R. H., et al. 2020, MNRAS, 499, 5702, doi: 10.1093/mnras/staa3164
  * Behroozi & Silk (2015) Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32, doi: 10.1088/0004-637X/799/1/32
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, The Astrophysical Journal, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20, doi: 10.1088/0067-0049/208/2/20
  * Benson (2010) Benson, A. J. 2010, New Astronomy, 17, 175, doi: 10.1016/j.newast.2011.07.004
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic dynamics (Princeton University Press), 885
  * Bouwens et al. (2023) Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036, doi: 10.1093/mnras/stad1145
  * 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
  * Boylan-Kolchin (2023) Boylan-Kolchin, M. 2023, Nature Astronomy, 7, 731, doi: 10.1038/S41550-023-01937-7
  * Bullock et al. (2001) Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240, doi: 10.1086/321477
  * Carrick et al. (2015) Carrick, J., Turnbull, S. J., Lavaux, G., & Hudson, M. J. 2015, MNRAS, 450, 317, doi: 10.1093/mnras/stv547
  * Chandrasekhar (1949) Chandrasekhar, S. 1949, Reviews of Modern Physics, 21, 383, doi: 10.1103/RevModPhys.21.383
  * Chen et al. (2023) Chen, Y., Mo, H. J., & Wang, K. 2023, MNRAS, 526, 2542, doi: 10.1093/mnras/stad2866
  * Chevallard & Charlot (2016) Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415, doi: 10.1093/mnras/stw1756
  * Cole et al. (2023) Cole, J. W., Papovich, C., Finkelstein, S. L., et al. 2023, eprint arXiv:2312.10152, arXiv:2312.10152, doi: 10.48550/ARXIV.2312.10152
  * Cole et al. (2005) Cole, S., Percival, W. J., & Peacock, J. A. 2005, MNRAS, 362, 505, doi: 10.1111/j.1365-2966.2005.09318.x
  * Cora (2006) Cora, S. A. 2006, MNRAS, 368, 1540, doi: 10.1111/j.1365-2966.2006.10271.x
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, Monthly Notices of the Royal Astronomical Society, 365, 11, doi: 10.1111/j.1365-2966.2005.09675.x
  * Davis et al. (2011) Davis, M., Nusser, A., Masters, K. L., et al. 2011, MNRAS, 413, 2906, doi: 10.1111/j.1365-2966.2011.18362.x
  * Dekel et al. (2009a) Dekel, A., Sari, R., & Ceverino, D. 2009a, ApJ, 703, 785, doi: 10.1088/0004-637X/703/1/785
  * Dekel et al. (2023) Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, Monthly Notices of the Royal Astronomical Society, 523, 3201, doi: 10.1093/mnras/stad1557
  * Dekel & Silk (1986) Dekel, A., & Silk, J. 1986, ApJ, 303, 39, doi: 10.1086/164050
  * Dekel et al. (2009b) Dekel, A., Birnboim, Y., Engel, G., et al. 2009b, Nature, 457, 451, doi: 10.1038/nature07648
  * Despali et al. (2016) Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486, doi: 10.1093/mnras/stv2842
  * D’Eugenio et al. (2023) D’Eugenio, F., Maiolino, R., Carniani, S., et al. 2023, eprint arXiv:2311.09908, arXiv:2311.09908, doi: 10.48550/ARXIV.2311.09908
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35, doi: 10.3847/1538-4365/aaee8c
  * Donnan et al. (2023) Donnan, C. T., Mcleod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011, doi: 10.1093/mnras/stac3472
  * Dressler et al. (2024) Dressler, A., Rieke, M., Eisenstein, D., et al. 2024
  * Eisenstein et al. (2007) Eisenstein, D. J., Seo, H.-J., Sirko, E., & Spergel, D. N. 2007, The Astrophysical Journal, 664, 675, doi: 10.1086/518712
  * El-Zant et al. (2001) El-Zant, A., Shlosman, I., & Hoffman, Y. 2001, ApJ, 560, 636, doi: 10.1086/322516
  * Elahi et al. (2019) Elahi, P. J., Cañas, R., Poulton, R. J., et al. 2019, Publications of the Astronomical Society of Australia, 36, e021, doi: 10.1017/pasa.2019.12
  * Faerman et al. (2020) Faerman, Y., Sternberg, A., & Mckee, C. F. 2020, ApJ, 893, 82, doi: 10.3847/1538-4357/ab7ffc
  * Fakhouri et al. (2010) Fakhouri, O., Ma, C. P., & Boylan-Kolchin, M. 2010, Monthly Notices of the Royal Astronomical Society, 406, 2267, doi: 10.1111/j.1365-2966.2010.16859.x
  * Ferrara et al. (2023) Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986, doi: 10.1093/MNRAS/STAD1095
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M., Song, M., et al. 2022, ApJ, 928, 52, doi: 10.3847/1538-4357/ac3aed
  * Fu et al. (2022) Fu, H., Shankar, F., Ayromlou, M., et al. 2022, MNRAS, 516, 3206, doi: 10.1093/mnras/stac2205
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, doi: 10.1088/1538-3873/acd1b5
  * Girelli et al. (2020) Girelli, G., Pozzetti, L., Bolzonella, M., et al. 2020, A&A, 634, A135, doi: 10.1051/0004-6361/201936329
  * Grand et al. (2021) Grand, R. J., Marinacci, F., Pakmor, R., et al. 2021, MNRAS, 507, 4953, doi: 10.1093/mnras/stab2492
  * 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
  * Harikane et al. (2024) Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56, doi: 10.3847/1538-4357/ad0b7e
  * Harikane et al. (2023) Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5, doi: 10.3847/1538-4365/acaaa9
  * Hellwing et al. (2017) Hellwing, W. A., Nusser, A., Feix, M., & Bilicki, M. 2017, MNRAS, 467, 2787, doi: 10.1093/mnras/stx213
  * Henriques et al. (2015) Henriques, B. M., White, S. D., Thomas, P. A., et al. 2015, MNRAS, 451, 2663, doi: 10.1093/MNRAS/STV705
  * Hoffman et al. (2015) Hoffman, Y., Courtois, H. M., & Tully, R. B. 2015, MNRAS, 449, 4494, doi: 10.1093/mnras/stv615
  * Hoffman & Ribak (1991) Hoffman, Y., & Ribak, E. 1991, ApJL, 380, L5, doi: 10.1086/186160
  * Jenkins et al. (2001) Jenkins, a. R., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372, doi: 10.1046/j.1365-8711.2001.04029.x
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Katz et al. (1993) Katz, N., White, S. D. M., Katz, N., & White, S. D. M. 1993, ApJ, 412, 455, doi: 10.1086/172935
  * Kauffmann et al. (1993) Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201, doi: 10.1093/mnras/264.1.201
  * Kennicutt (1998) Kennicutt, R. C. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Kennicutt & De Los Reyes (2021) Kennicutt, R. C., & De Los Reyes, M. A. C. 2021, ApJ, 908, 61, doi: 10.3847/1538-4357/abd3a2
  * Kennicutt et al. (2012) Kennicutt, R. C., Evans, N. J., Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531, doi: 10.1146/ANNUREV-ASTRO-081811-125610
  * Krumholz et al. (2018) Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716, doi: 10.1093/mnras/sty852
  * Larson (1974) Larson, R. B. 1974, MNRAS, 169, 229, doi: 10.1093/mnras/169.2.229
  * Leauthaud et al. (2012) Leauthaud, A., Tinker, J., Bundy, K., et al. 2012, ApJ, 744, 159, doi: 10.1088/0004-637X/744/2/159
  * Libeskind et al. (2020) Libeskind, N. I., Carlesi, E., Grand, R. J., et al. 2020, MNRAS, 498, 2968, doi: 10.1093/mnras/staa2541
  * Lilow et al. (2021) Lilow, R., Nusser, A., Lilow, R., & Nusser, A. 2021, MNRAS, 507, 1557, doi: 10.1093/MNRAS/STAB2009
  * Ludlow et al. (2014) Ludlow, A. D., Borzyszkowski, M., Porciani, C., et al. 2014, MNRAS, 445, 4110, doi: 10.1093/MNRAS/STU2021
  * Marshall et al. (2022) Marshall, M. A., Wilkins, S., Di Matteo, T., et al. 2022, MNRAS, 511, 5475, doi: 10.1093/mnras/stac380
  * Mason et al. (2023) Mason, C. A., Trenti, M., & Treu, T. 2023, MNRAS, 521, 497, doi: 10.1093/mnras/stad035
  * Mitchell et al. (2018) Mitchell, P. D., Lacey, C. G., Lagos, C. D., et al. 2018, MNRAS, 474, 492, doi: 10.1093/mnras/stx2770
  * Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. 1998, MNRAS, 295, 319, doi: 10.1046/j.1365-8711.1998.01227.x
  * Moster et al. (2020) Moster, B. P., Naab, T., & White, S. D. 2020, MNRAS, 499, 4748, doi: 10.1093/MNRAS/STAA3019
  * Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121, doi: 10.1093/mnras/sts261
  * Nusser & Silk (2022) Nusser, A., & Silk, J. 2022, MNRAS, 509, 2979, doi: 10.1093/mnras/stab3121
  * Okalidis et al. (2021) Okalidis, P., Grand, R. J., Yates, R. M., & Kauffmann, G. 2021, MNRAS, 504, 4400, doi: 10.1093/MNRAS/STAB1142
  * Pacifici et al. (2023) Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141, doi: 10.3847/1538-4357/ACACFF
  * Pallottini & Ferrara (2023) Pallottini, A., & Ferrara, A. 2023, A&A, 677, doi: 10.1051/0004-6361/202347384
  * Peebles (1980) Peebles, P. J. E. 1980, Research supported by the National Science Foundation. Princeton, N.J., Princeton University Press, 1980. 435 p.
  * Pérez-González et al. (2023) Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJL, 951, L1, doi: 10.3847/2041-8213/acd9d0
  * Planck Collaboration et al. (2018) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, eprint arXiv, 1807.06209.  http://arxiv.org/abs/1807.06209
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425, doi: 10.1086/152650
  * Press et al. (1974) Press, W. H., Schechter, P., Press, W. H., & Schechter, P. 1974, ApJ, 187, 425, doi: 10.1086/152650
  * Puchwein & Springel (2013) Puchwein, E., & Springel, V. 2013, MNRAS, 428, 2966, doi: 10.1093/mnras/sts243
  * Riess et al. (2023) Riess, A. G., Anand, G. S., Yuan, W., et al. 2023, ApJL, 956, L18, doi: 10.3847/2041-8213/acf769
  * Robertson et al. (2009) Robertson, B. E., Kravtsov, A. V., Tinker, J., & Zentner, A. R. 2009, ApJ, 696, 636, doi: 10.1088/0004-637X/696/1/636
  * Robertson et al. (2023) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, Nature Astronomy, 7, 611, doi: 10.1038/S41550-023-01921-1
  * Rodríguez-Puebla et al. (2013) Rodríguez-Puebla, A., Avila-Reese, V., & Drory, N. 2013, ApJ, 767, 92, doi: 10.1088/0004-637X/767/1/92
  * 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
  * Rodríguez-Puebla et al. (2017) Rodríguez-Puebla, A., Primack, J. R., Avila-Reese, V., & Faber, S. M. 2017, MNRAS, 470, 651, doi: 10.1093/mnras/stx1172
  * Schaller et al. (2018) Schaller, M., Gonnet, P., Draper, P. W., et al. 2018, Astrophysics Source Code Library, ascl:1805.020.  https://ui.adsabs.harvard.edu/abs/2018ascl.soft05020S/abstract
  * Schaller et al. (2023) Schaller, M., Borrow, J., Draper, P. W., et al. 2023, MNRAS, 000, 1.  www.swiftsim.com.
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243, doi: 10.1086/146614
  * Seppi et al. (2021) Seppi, R., Comparat, J., Nandra, K., et al. 2021, A&A, 652, A155, doi: 10.1051/0004-6361/202039123
  * Sheth & Tormen (1999) Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119, doi: 10.1046/j.1365-8711.1999.02692.x
  * Shuntov et al. (2022) Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61, doi: 10.1051/0004-6361/202243136
  * Silk et al. (2024) Silk, J., Begelman, M. C., Norman, C., Nusser, A., & Wyse, R. F. G. 2024, ApJL, 961, L39, doi: 10.3847/2041-8213/AD1BF0
  * Silk & Rees (1998) Silk, J., & Rees, M. J. 1998, A&A, 331, L1.  https://ui.adsabs.harvard.edu/abs/1998A%26A...331L...1S/abstract
  * Somerville & Davé (2015) Somerville, R. S., & Davé, R. 2015, ARAA, 53, 51, doi: 10.1146/annurev-astro-082812-140951
  * Somerville et al. (2015) Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337, doi: 10.1093/MNRAS/STV1877
  * Sun et al. (2023a) Sun, G., Faucher-Giguère, C. A., Hayward, C. C., & Shen, X. 2023a, MNRAS, 526, 2665, doi: 10.1093/mnras/stad2902
  * Sun et al. (2023b) Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023b, ApJL, 955, L35, doi: 10.3847/2041-8213/acf85a
  * Tegmark et al. (2004) Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, The Astrophysical Journal, 606, 702, doi: 10.1086/382125
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, The Astrophysical Journal, 688, 709, doi: 10.1086/591439
  * Wang et al. (2023) Wang, B., Leja, J., Labbé, I., et al. 2023, ApJS, 270, 12, doi: 10.3847/1538-4365/ad0846
  * Wang et al. (2016) Wang, H., Mo, H. J., Yang, X., et al. 2016, ApJ, 831, 164, doi: 10.3847/0004-637x/831/2/164
  * White & Frenk (1991a) White, S. D. M., & Frenk, C. S. 1991a, ApJ, 379, 52, doi: 10.1086/170483
  * White & Frenk (1991b) —. 1991b, ApJ, 379, 52, doi: 10.1086/170483
  * Whitler et al. (2023) Whitler, L., Stark, D. P., Endsley, R., et al. 2023, MNRAS, 519, 5859, doi: 10.1093/mnras/stad004
  * Yang et al. (2023) Yang, H., Gao, L., Frenk, C. S., et al. 2023, MNRAS, 518, doi: 10.1093/mnras/stac3335
  * Yang et al. (2012) Yang, X., Mo, H. J., Van Den Bosch, F. C., Zhang, Y., & Han, J. 2012, ApJ, 752, 41, doi: 10.1088/0004-637X/752/1/41
  * Yung et al. (2019) Yung, L. Y., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019, MNRAS, 483, 2983, doi: 10.1093/mnras/sty3241
  * Yung et al. (2024) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2024, MNRAS, 527, 5929, doi: 10.1093/mnras/stad3484