# The impact of mass-dependent stochasticity at cosmic dawn

Viola Gelli  Charlotte Mason  Cosmic Dawn Center (DAWN)  Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, 2200 København N, Denmark  Christopher
C. Hayward  Center for Computational Astrophysics, Flatiron Institute, 162
Fifth Avenue, New York, NY 10010, USA

###### Abstract

JWST is unveiling a surprising lack of evolution in the number densities of
ultraviolet-selected (UV) galaxies at redshift $z\gtrsim 10$. At the same
time, observations and simulations are providing evidence for highly bursty
star formation in high-$z$ galaxies, resulting in significant scatter in their
UV luminosities. Galaxies in low-mass dark matter halos are expected to
experience most stochasticity due to their shallow potential wells. Here, we
explore the impact of a mass-dependent stochasticity using a simple analytical
model. We assume that scatter in the $M_{\rm UV}-M_{h}$ relation increases
towards lower halo masses, following the decrease in halo escape velocity,
$\sigma_{\mathrm{UV}}\sim M_{h}^{-1/3}$, independent of redshift. Since low-
mass halos are more dominant in the early universe, this model naturally
predicts an increase in UV luminosity functions (LFs) at high redshifts
compared to models without scatter. We make predictions for additional
observables which would be affected by stochasticity and could be used to
constrain its amplitude, finding: (i) galaxies are less clustered compared to
the no-scatter scenario, with the difference increasing at higher-$z$; (ii)
assuming star-bursting galaxies dominate the ionizing photon budget implies
reionization starts earlier and is more gradual compared to the no-scatter
case, (iii) at fixed UV magnitude galaxies should exhibit wide ranges of UV
slopes, nebular emission line strengths and Balmer breaks. Comparing to
observations, the mass-dependent stochasticity model successfully reproduces
the observed LFs up to $z\sim 12$. However, the model cannot match the
observed $z\sim 14$ LFs, implying additional physical processes enhance star
formation efficiency in the earliest galaxies.

High-redshift galaxies; Galaxy evolution; Cosmology

##  1 Introduction

During its first two years of observations, the James Webb Space Telescope
(JWST) has revealed a puzzling abundance of luminous galaxies at high
redshifts ($z>10$) that are challenging our understanding of early star
formation. Numerous observations of ultraviolet (UV) bright galaxies have been
reported, and the photometrically derived UV luminosity functions (LF) show a
surprising lack of evolution, across a range of luminosities, during the first
$\sim 500\leavevmode\nobreak\ \rm Myr$ after the Big Bang (e.g. Adams et al.,
2023; Atek et al., 2023a; Castellano et al., 2022; Donnan et al., 2024;
Finkelstein et al., 2023; Harikane et al., 2023; McLeod et al., 2024; Pérez-
González et al., 2023; Bouwens et al., 2023). Most pre-JWST models struggle to
reproduce the observed trends, prompting the emergence of numerous new
theories in an effort to solve this discrepancy.

One possibility is a systematic enhancement of the median UV flux of high-$z$
galaxies. This enhancement may be achieved through different mechanisms, such
as: higher star formation efficiencies (Mason et al., 2023; Dekel et al.,
2023; Inayoshi et al., 2022; Li et al., 2023), top-heavier initial mass
functions (although likely not sufficient to reproduce the observations up to
the highest $z$, e.g., Rasmussen Cueto et al., 2023; Trinca et al., 2024), or
efficient dust removal (Ferrara, 2023; Fiore et al., 2023).

An alternative explanation is the presence of a stochasticity in the galaxies
luminosities: the relation between the UV magnitude of a galaxy ($M_{\rm UV}$)
and the mass of the dark matter halo in which it resides ($M_{h}$) may not be
univocal, but rather characterized by a dispersion $\sigma_{UV}$ (e.g. Mason
et al., 2023; Mirocha & Furlanetto, 2023; Shen et al., 2023; Kravtsov &
Belokurov, 2024). In this way, the increased normalization of the LFs can be
explained by the up-scattered luminosities in the numerous low-mass halos
dominating the high-$z$ population.

A physical mechanism that can naturally produce a scatter in the UV
luminosities of galaxies is a bursty star formation history (SFH). The complex
interplay between multiple mechanisms, such as internal feedback processes,
the continuous inflow of gas from the cosmic web and mergers, can cause a
stochastic behavior in the star formation rate (SFR) of high-$z$ galaxies,
making their UV luminosity fluctuate over time (e.g. Sparre et al., 2017;
Furlanetto & Mirocha, 2022; Sun et al., 2023a). Both JWST observations and
cosmological simulations suggest that star formation in the first billion
years after the Big Bang may have occurred in a particularly bursty fashion.

From the observational side, deep JWST surveys are revealing a broad range of
SFRs in high-$z$ star forming galaxies (e.g. Curti et al., 2023; Endsley et
al., 2023b). In particular, Endsley et al. (2023b) recently found that UV-
faint galaxies ($M_{\rm UV}\gtrsim-17$) often have weaker nebular emission
lines and can be experiencing a down-turn in star formation compared to UV-
bright galaxies ($M_{\rm UV}\lesssim-20$) which are more likely to have strong
emission lines and to be in an up-turn of star formation. This is consistent
with a picture whereby the bright-end of the UV LF is dominated by galaxies
which have been up-scattered, whereas fainter galaxies show a range of star
formation history stages. Moreover, the first detections of quenched low-mass
galaxies ($M_{\star}<10^{9}{\rm M}_{\odot}$) at $z>5$ (Looser et al., 2023a,
b; Strait et al., 2023) can be interpreted as post-starburst galaxies
undergoing temporary periods of quiescence due to a highly bursty and
feedback-regulated star formation (Gelli et al., 2023, 2024; Dome et al.,
2023). The timescales on which these SFR variations occur tightly depend on
the ongoing feedback processes within the galaxy (e.g., Iyer et al., 2020).
Thus, constraining the level of stochasticity and burstiness of high-$z$
systems is essential for understanding the physical mechanisms driving their
evolution.

State-of-the-art cosmological simulations also predict a highly time variable
star formation in the early Universe (e.g. Ma et al., 2018; Pallottini et al.,
2022). However, whether the stochasticity in the SFRs of early galaxies is
sufficient to explain the high number densities observed is still debated:
Pallottini & Ferrara (2023) argue that an additional source of luminosity is
needed already at $z\sim 8$, while Sun et al. (2023b) finds that $z\sim 8-12$
LFs can be reproduced with fire-2 simulations, which show strongly time-
variable SFHs. Meanwhile, analytical studies have shown that, to match the
observed LFs, the level of scatter in the UV luminosities $\sigma_{UV}$ needs
to increase significantly above $z\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10$ (Muñoz et al., 2023; Shen et
al., 2023; Mason et al., 2023; Kravtsov & Belokurov, 2024). In these works,
the UV scatter is assumed to be independent of the halo mass. However, low-
mass halos are intrinsically characterized by shallower potential wells and
the galaxies residing in them are expected to be more sensitive to all those
feedback and environmental processes that drive bursty star formation (e.g.
Gelli et al., 2020; Furlanetto & Mirocha, 2022; Legrand et al., 2022; Hopkins
et al., 2023). For this reason, we expect stochasticity to be more relevant
for lower mass halos, which, in a $\Lambda$CDM Universe, dominate the galaxy
population at progressively higher-$z$.

In this paper, we develop a simple analytical model to explore a mass-
dependent UV scatter, i.e. increasing towards lower halo masses, independent
of the redshift, described in Section 2. We analyze the impact of this mass-
dependent UV scatter on the LF, but also on other key observables: the neutral
fraction in the intergalactic medium (IGM), the galaxy bias and spectral
features, providing unique independent tools to help us constrain the
relevance of bursty star formation in high-$z$ galaxies. These results are
presented in Section 3. We discuss our results in the context of the $z>12$
LFs in Section 4 and present our conclusions in Section 5. We assume a flat
$\Lambda$CDM cosmology with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$ and
$h=0.7$.

##  2 Modeling the mass-dependent UV scatter

Figure 1: Top panel: Probability distribution of halo masses at given UV
magnitude $M_{\rm UV}$, with corresponding scatter $\sigma$ determined by the
model in the bottom panel. The distributions are shown for three different
redshifts, highlighting the broadening of the distributions towards earlier
epochs. Bottom panel: Mass-dependent UV scatter $\sigma_{\rm UV}(M_{h})$ as a
function of the halo mass $M_{h}$. The black line shows the fiducial model
used in this work, for all redshifts.

We model the stochasticity by assuming a mass-dependent but redshift-
independent scatter $\sigma_{\rm UV}(M_{h})$, describing the dispersion in the
$M_{\rm UV}-M_{h}$ relation, which decreases with halo mass. This is motivated
by cosmological simulations such as FIRE-2 which find a decreasing trend with
the halo mass, without a significant evolution between $z\sim 8-12$ (Sun et
al., 2023b). This scatter can be well approximated by a power-law relation
$\sigma_{\rm UV}(M_{h})=a\log M_{h}/{\rm M}_{\odot}+b$. We assume a fiducial
model with $a=-0.34$ and $b=4.5$, shown in the bottom panel of Fig. 1,
obtained by fitting the results of the FIRE-2 simulations at $z\sim 8$
(derived from Fig. 1 of Sun et al. 2023b). Interestingly, this scaling with
halo mass is very similar to the $M_{h}^{-1/3}$ scaling expected if the
stochasticity scales with the inverse of halos’ escape velocities – i.e. more
tightly bound halos can retain more gas and have lower stochasticity.

Figure 2: UV LF for the mass-dependent UV scatter model (solid lines) and the
no-scatter model (dashed lines) at redshifts from $z=5$ to 14. Also shown are
observations by Bouwens et al. (2022a); Finkelstein et al. (2023); Pérez-
González et al. (2023); Donnan et al. (2024); McLeod et al. (2024).

To include the scatter in our modeling, we adopt the conditional luminosity
function (CLF) approach, following Ren et al. (2019) (see also Whitler et al.,
2020). The probability of a dark matter halo with mass $M_{h}$ to host a
galaxy with magnitude $M_{\rm UV}$ is given by:

$p(M_{\rm UV}\mid M_{h})=\\\ \frac{1}{\sqrt{2\pi}\sigma_{\rm
UV}(M_{h})}\exp\left(\frac{-[M_{\rm UV}-M_{{\rm
UV},c}(M_{h},z)]^{2}}{2\sigma_{\rm UV}^{2}(M_{h})}\right),$ (1)

where $M_{{\rm UV},c}$ is the median magnitude at a given redshift $z$ and
halo mass $M_{h}$. We derive it by following Mason et al. (2015, 2023),
assuming $\langle{\rm SFR}\rangle\propto\epsilon_{\star}\,f_{b}M_{h}$, where
$f_{b}$ is the cosmic baryon fraction and $\epsilon_{\star}$ the star
formation efficiency, and then using starburst99 (Leitherer et al., 1999)111We
use padova stellar tracks (Bertelli et al., 1994) and assume a Kroupa (2001)
initial mass function between 0.1 and 100$M_{\odot}$. However we note that the
following results are not sensitive to these specific choices. to derive
$M_{\rm UV}$. We assume a mass-dependent star formation efficiency
$\epsilon_{\star}(M_{h})\sim M_{\star}/M_{h}$, and we calibrate222We note that
when introducing mass-independent scatter, Ren et al. (2019) required a break
in $\epsilon_{\star}(M_{h})$ above $M_{h}\gtrsim 10^{11.5}$ to match the
bright-end of the LF, however, this is not necessary in our model as
$\sigma_{\rm UV}(M_{h})$ is so low at these high halo masses. it at a single
redshift $z\sim 5$ to reproduce the observed UVLFs by Bouwens et al. (2023)
(see also Mason et al., 2015; Ren et al., 2019; Mirocha & Furlanetto, 2023).
Like the UV scatter, once calibrated, the derived star formation efficiency is
assumed to be redshift independent.

We can derive halo mass distributions $p(M_{h}\mid M_{\rm UV})$ using a
numerical approach, evaluating the probability in Equation 1 for an array of
halo masses at fixed $M_{\rm UV}$ and $z$ (see also Whitler et al., 2020). The
results are shown in the top panel of Fig. 1 for two different magnitudes and
redshifts. The position of the peak of the distribution is determined by the
median $M_{{\rm UV},c}(M_{h},z)$, while the width of the distribution is
determined by the scatter $\sigma_{\rm UV}(M_{h})$ and is independent of the
redshift. At fixed redshift, brighter galaxies are more likely to be hosted by
larger dark matter halos. Note however that here we are not weighting by the
halo mass function. As lower mass halos are in reality more numerous, galaxies
we observe at fixed UV magnitudes will be more likely to be hosted in lower
mass halos, as it will be evident in the next section. The dispersion around
the median $M_{h}$ is larger for lower luminosities. This is a direct
consequence of the increasing $\sigma_{\rm UV}$ towards lower masses, making
lower mass halos have high probability to appear with a wider range of
luminosities.

Comparing the different redshifts, we can see that there is a broadening of
the distribution towards high-$z$, meaning that galaxies at a given $M_{\rm
UV}$ can be hosted at higher $z$ by a broader range of halos with respect to
their lower $z$ analogues. This behavior is due to the combined effect of: (i)
the increasing median luminosity $M_{{\rm UV},c}(M_{h},z)$ with the redshift –
i.e. at fixed halo mass, galaxies are on-average more luminous at higher
redshift, due to higher SFRs, and (ii) the increasing scatter $\sigma_{\rm
UV}$ in lower mass halos.

In the following sections we determine and analyze the expected impact of such
mass-dependent UV scatter on different observables (UVLFs, SFR density,
neutral hydrogen fractions, spectral features).

##  3 Results

We present our model for the UV luminosity function in Section 3.1 and the
luminosity density in Section 3.2. In Section 3.3 we make predictions for
other independent observables: galaxy clustering, the reionization timeline,
and galaxy SEDs; which could be combined with the LF to constrain the level of
stochasticity.

###  3.1 Luminosity function

Following the CLF approach, we derive the luminosity functions by integrating
the number density of dark matter halos, weighted by the probability that a
halo $M_{h}$ hosts a galaxy with magnitude $M_{\rm UV}$:

$\Phi(M_{\rm UV})=\int_{M_{h,min}}^{\infty}dM_{h}\,p(M_{\rm UV}\mid
M_{h})\,\frac{dn}{dM_{h}};$ (2)

where $p(M_{\rm UV}\mid M_{h})$ is given by Equation 1, and for $dn/dM_{h}$ we
use the halo mass function (HMF) from Reed et al. (2007). We set a lower limit
to the integral, $M_{h,min}$, to prevent the unphysical effect of very low-
mass halos forming stars with extremely large UV scatter. Its value is set at
each redshift as the minimum halo mass for which atomic cooling can occur,
i.e. corresponding to a virial temperature of $T_{\rm vir}\sim 10^{4}\rm K$
(e.g. Bromm & Yoshida, 2011). Following Mason et al. (2015, 2023) we add dust
attenuation using the Meurer et al. (1999) attenuation law
$A_{\mathrm{UV}}=4.43+1.99\beta$, including the UV slope $\beta(M_{\rm UV},z)$
empirically from observations by Bouwens et al. (2014). We remove dust
attenuation at $z>10$ motivated by recent results suggesting dust attenuation
is negligible in $z>10$ galaxies (Topping et al., 2024a; Cullen et al., 2024;
Morales et al., 2024), though we note that beyond $z\gtrsim 10$ dust has only
a small impact on the UV LFs (Mason et al., 2023).

The resulting LFs for our mass-dependent UV scatter model are shown in Fig. 2
from $z=5$ to 14. We compare them with the model without UV scatter (dashed
lines), where the evolution of the LF directly follows the evolution of the
underlying halo mass function. The solid curves are the LF obtained including
the mass-dependent scatter $\sigma_{\rm UV}(M_{h})$ described in the previous
section.

At progressively higher redshifts, the scatter in the UV luminosities has a
stronger impact on the LF, which deviates more and more from the no-scatter
case. The presence of a UV scatter that decreases with halo mass naturally
produces larger number densities of $M_{\rm UV}\lesssim-14$ galaxies at higher
redshifts with respect to the no-scatter case. This is due to low mass halos
(characterized by higher UV scatter) being progressively more dominant at
earlier epochs.

As evident from the bottom right panel of Fig. 2 where we show the LFs up to
fainter luminosities, the model also predicts a faint end turnover of the LF
at magnitudes $M_{\rm UV}\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}-13$. This is because we have
imposed a cut-off in galaxy formation at the atomic cooling limit $M\sim
10^{8}\,M_{\odot}$, i.e. lower mass halos cannot host a galaxy and thus will
not be scattered up to populate the LF. While galaxies in higher mass halos
are scattered down to lower luminosities during lulls in star formation, their
number densities are lower than low mass halos, and they cannot boost the
faint-end of the LF. If ever probed with deep observations (e.g. in lensed
fields with JWST Atek et al., 2023b), this could provide an additional
constraint on the dominant driver of star formation at high redshifts.

We find the mass-dependent UV scatter model successfully reproduces the
observed HST and JWST up to $z\sim 12$ (Bouwens et al., 2022b; Finkelstein et
al., 2023; Pérez-González et al., 2023; Donnan et al., 2024; McLeod et al.,
2024), producing a $\sim 1$ dex increase in number densities of $M_{\rm
UV}<-18$ galaxies at $z\gtrsim 9$ compared to the no-scatter model. The model
slightly overpredicts the bright-end of the LF at $z\sim 8$, which is likely
due to a combination of the dust modelling and the precise form of
$\sigma_{\rm UV}(M_{h})$. However, at redshifts greater than $>12$, even the
additional contribution from an enhanced UV scatter at low masses is not
sufficient to reproduce the observed high number densities of galaxies. This
suggests that additional physical processes may be at play at these redshifts
to enhance the UV luminosities, as we will discuss in Sec. 4.

Figure 3: Luminosity density $\rho_{\rm UV}$ and cosmic star formation rate
density $\rm SFRD$ as a function of redshift, derived integrating the UV LFs
down to $M_{\rm lim}=-17$. The solid and dashed curves show the mass-dependent
UV scatter and no-scatter models respectively. Also shown are observations
from Bouwens et al. (2022a); Finkelstein et al. (2023); Pérez-González et al.
(2023); Donnan et al. (2024); McLeod et al. (2024). We note that the apparent
match to the observations at $z\sim 14$ is due to the steep faint-end slope of
our LFs relative to that assumed by Donnan et al. (2024) (see Section 3.2).

###  3.2 Luminosity density evolution

By integrating the LFs we derive the luminosity density $\rho_{\rm UV}$ for
the mass-dependent UV scatter model, shown in Fig. 3 as a function of
redshift. To be consistent with observations, we obtain $\rho_{\rm UV}$ by
integrating the LFs down to a magnitude limit of $M_{\rm lim}=-17$. We also
calculate the corresponding cosmic star formation rate density (SFRD) using
the empirical relation $\rm SFR/({\rm M}_{\odot}\,{\rm yr}^{-1})=8\times
10^{27}L/(erg\,s^{-1}Hz^{-1})$ from Madau et al. (1999).

The mass-dependent UV scatter and the no scatter models, shown with solid and
dashed lines respectively, predict similar luminosity density at $z\sim 7$,
but progressively diverge towards higher redshifts. While the no-scatter model
underpredicts the observed luminosity density already from $z>9$, the mass-
dependent UV scatter model is in good agreement with the latest JWST
observations up to $z\sim 11$, but then diverges from the observations at
progressively earlier epochs.

We note that the luminosity density of the mass-dependent UV scatter model is
consistent with the data from Donnan et al. (2024) up to $z\sim 14$, despite
their LF exhibiting higher values at this redshift. This happens because their
double power-law fit yields a shallower LF slope than our model towards
fainter magnitudes, resulting in lower luminosity density when integrating up
to $M_{\rm lim}=-17$. It is important to note that the luminosity density only
reflects the integrated LF and caution should be taken when drawing
conclusions based solely on its value. Therefore, while our luminosity density
only slightly diverges from the data towards higher-$z$, this still suggests
the need for extra luminosity production beyond that provided by stochastic
star formation.

###  3.3 Constraining the level of stochasticity

Even if it cannot fully explain the number density of $z\gtrsim 12$ galaxies,
bursty star formation still plays a role at cosmic dawn and it results from
the complex interplay of multiple feedback processes acting on different
timescales, which are not yet fully understood. Constraining the level of
stochasticity in galaxies would offer key insights into the physical and
feedback processes driving the evolution of galaxies of different masses
across cosmic time, and their associated timescales (e.g., Iyer et al., 2020).

However, the LF and luminosity density are intrinsically degenerate with
respect to multiple physical effects, such as an enhancement of the star
formation efficiency or of the UV scatter. As discussed by e.g., Mirocha
(2020); Muñoz et al. (2023), this makes it impossible to use them alone as a
tool to disentangle these physical processes. For this reason, it is important
to identify additional key observables that can help us to independently
constrain the importance of stochastic star formation at high-$z$.

To this aim, we discuss the expected implications of mass-dependent UV scatter
on multiple quantities: galaxy clustering (Section 3.3.1), the evolution of
the intergalactic medium (IGM) neutral fraction (Section 3.3.2), and spectral
energy distributions (Section 3.3.3), providing predictions and tools to test
the model against future observations.

####  3.3.1 Clustering

Figure 4: Redshift evolution of the predicted galaxy bias for the mass-
dependent UV scatter model (solid lines) and the no-scatter model (dashed
lines). Observations by Qiu et al. (2018) and Dalmasso et al. (2024a) are
shown.

A key observational test of the importance of stochasticity in high-$z$
galaxies is to connect galaxies to their dark matter halos. The clustering of
galaxies is driven by the underlying distribution of the dark matter halos
that they inhabit (see e.g., Wechsler & Tinker, 2018, for a recent review). As
discussed by Muñoz et al. (2023), constraints on galaxy clustering can thus
break the degeneracies of the UVLF and distinguish the presence of UV
stochasticity from the effect of an increased average UV emission.

In a $\Lambda$CDM Universe, the galaxy bias (quantifying the halo-galaxy
connection) is expected to increase with the halo mass: galaxies hosted by
massive halos cluster more strongly than those residing in lower mass ones.
Estimates of Lyman-break galaxy clustering at $z\sim 4-7$ have revealed the
brightest galaxies tend to reside in the most massive halos (e.g., Barone-
Nugent et al., 2014; Harikane et al., 2016, 2022; Qiu et al., 2018; Dalmasso
et al., 2024b). As discussed by Muñoz et al. (2023), scatter in the $M_{\rm
UV}-M_{h}$ relation will lower the bias of galaxies, since UV luminous
galaxies can be found in a broad range of halo masses.

It is therefore important to provide predictions for the galaxy bias in our
mass-dependent UV-scatter model. We estimate the luminosity-weighted effective
bias of galaxies as (see e.g., Muñoz et al., 2023):

$b_{\rm eff}(M_{\rm UV})=\frac{1}{\Phi(M_{\rm UV})}\int dM_{h}\,p(M_{\rm
UV}\mid M_{h})\,\frac{dn}{dM_{h}}b(M_{h})$ (3)

where we use the linear bias $b(M_{h})$ fits from Tinker et al. (2010). Our
predictions for the redshift evolution of the bias at fixed $M_{\rm UV}$ are
shown in Fig. 4.

As expected, biases from the UV scatter model always have lower values than
the no-scatter case, and the difference is progressively more pronounced
towards higher redshifts as galaxies of fixed $M_{\rm UV}$ occupy lower mass
halos. As discussed by Dalmasso et al. (2024b), measurements of galaxy bias
are extremely challenging at high-$z$ as estimates are limited by both Poisson
noise due to low number counts and cosmic variance in small fields. Pre-JWST
studies have measured the bias up to $z\sim 7$ for $M_{\rm UV}\lesssim-19$
galaxies, but the errors are too large to draw definite conclusions and
constrain the models, as we can see when comparing our model with data from
Qiu et al. (2018).

However, the sensitivity of JWST makes more accurate estimates of the bias
feasible by enabling us to measure clustering for sub-$L^{\star}$ galaxies at
$z>6$ for the first time (e.g., Endsley et al., 2020; Dalmasso et al., 2024a).
We compare our results with those of recent estimates from JWST/NIRCam
observations of $M_{\rm UV}\lesssim-15$ galaxies at $z\sim 5-10$ in the
GOODS-S field by Dalmasso et al. (2024a). The bias estimated by Dalmasso et
al. (2024a) shows a very steep increase with redshift. At $z<9$ the difference
between our two models is too small to be probed by these observations.
However, $z>9$, the Dalmasso et al. (2024a) constraints are more consistent
with the no-scatter case, and even slightly above both models at $z\sim 10$.

This could suggest galaxies experience a lower level of stochasticity than
what we have assumed and/or are hosted in higher mass halos than expected at
higher redshifts. This is somewhat surprising, as at fixed UV magnitude,
models predict and clustering results imply that galaxies reside in lower mass
halos at higher redshifts as star formation rates are higher (e.g., Mason et
al., 2015; Tacchella et al., 2018; Harikane et al., 2022). This may
potentially imply a change in the star formation efficiency with halo mass at
high redshift.

However, we note that Dalmasso et al. (2024a) assumed a different cosmology
than we have, and full interpretation of JWST observations likely requires a
Halo Occupation Distribution (HOD) model approach (e.g., Harikane et al.,
2022) to separate non-linear clustering on small scales due to satellite
galaxies, as noted by Dalmasso et al. (2024a). In addition, large area,
unbiased, surveys with JWST will be essential to constrain models by probing
the number density and clustering of UV bright galaxies ($M_{\rm
UV}\lesssim-20$) at $z\gtrsim 7$.

Figure 5:  Left: $M_{\rm UV}-M_{h}$ relation at $z=9$, populated with galaxies
following the halo mass function and coloured with the product $f_{\rm
esc}\xi_{\rm ion}$ (assuming an exponentially increasing $f_{\rm esc}$ and
linearly increasing $\xi_{\rm ion}$ with $\Delta M_{\rm UV}$), showing how up-
scattered galaxies will contribute the most to the number density of ionizing
photons. Centre: escape fraction $f_{\rm esc}$ and ionising photon production
efficiency $\xi_{\rm ion}$ as a function of $\Delta M_{\rm UV}$. For the
escape fraction we show three different cases of constant, linear and
exponential increase towards up-scattered galaxies. Right: Redshift evolution
of the fraction of neutral hydrogen in the IGM for the same three scenarios,
as well as for the case without scatter, assuming fixed $f_{\rm esc}=0.08$ and
$\log_{10}\xi_{\mathrm{ion}}/(\rm Hz\,erg^{-1})=25.7$. We also plot
measurements from observations of: the evolution of the $\rm Ly\alpha$
equivalent width (EW, Mason et al. 2018, 2019; Bolan et al. 2022), the
clustering of $\rm Ly\alpha$ emitters (LAE, Ouchi et al. 2010; Sobacchi &
Mesinger 2015), $\rm Ly\alpha$ forest dark pixel fraction (McGreer et al.,
2015), and quasar damping wings (Davies et al., 2018; Greig et al., 2019).

####  3.3.2 Reionization history

Constraints on the IGM reionization process offer vital information about star
formation in all sources, even those too faint to detect directly with JWST.
Muñoz et al. (2024) recently described how early JWST estimates of high
ionizing photon production efficiency in high-$z$ galaxies, taken at face-
value, may imply a potential ‘overproduction’ of ionizing photons, which is at
odds with current constraints on the timeline of reionization. Bursty star
formation may alleviate this tension. In particular, Muñoz et al. (2024)
assumed constant values of the ionizing photon escape fraction $f_{\rm esc}$
and ionizing photon production rate, $\xi_{\rm ion}$ at fixed $M_{\rm UV}$,
which may not be realistic assumptions in a stochastic star formation scenario
(see also Section 3.3.3).

Here, first, we discuss results from observations and simulations implying
that ionizing photon production and escape may be linked to bursty star
formation. Then, we combine our mass-dependent UV scatter model with simple
models for $f_{\rm esc}$ and $\xi_{\rm ion}$ as a function to offset from the
median $M_{\mathrm{UV,c}}(M_{h})$ relation to explore how bursty star
formation impacts the reionization history.

Observationally, high ionizing photon production rate ($\xi_{\rm ion}$) is
linked to high equivalent width nebular emission lines, and thus young,
massive stars (e.g., Chevallard et al., 2018; Tang et al., 2019, 2021).
Therefore, if galaxies undergo bursty star formation, the number of ionizing
photons they produce will vary significantly throughout their star formation
history. Indeed, as discussed in the introduction, Endsley et al. (2023b)
recently found that UV-faint galaxies ($M_{\rm UV}\gtrsim-17$) are more likely
to have weaker nebular emission lines and be experiencing a down-turn in star
formation compared to UV-bright galaxies ($M_{\rm UV}\lesssim-20$). Endsley et
al. (2023b) discuss that one interpretation of the weaker emission lines in
UV-faint galaxies is due to their star formation histories (rather than an
increase in $f_{\mathrm{esc}}$ in faint galaxies, which produces a degenerate
effect in photometric data). In this case, they find $\xi_{\rm ion}$ is a
broad distribution that can be strongly dependent on UV magnitude.

Furthermore, hydrodynamical simulations predict that the escape fraction of
ionizing photons may also coincide with (or shortly after) star formation
bursts (Kimm & Cen, 2014; Paardekooper et al., 2015; Ma et al., 2020a; Barrow
et al., 2020). In particular, hard radiation from massive stars may allow
ionizing photons to escape by creating ‘porous’ density-bounded HII regions
and/or stellar winds or supernovae feedback can open low-density channels
(e.g., Zackrisson et al., 2013; Steidel et al., 2018). From the observational
side, this picture is not completely clear: while many Lyman-continuum leakers
detected at $z\sim 0.3-0.4$ are in a period of intense star formation, with
specific star formation rates $\gtrsim 10$ Gyr-1 and star formation rate
surface densities $\Sigma_{\mathrm{SFR}}\gtrsim 10\,M_{\odot}$ yr-1 kpc-2
(e.g., Izotov et al., 2016, 2018; Flury et al., 2022), observations at $z\sim
3$ show no clear trend of $f_{\mathrm{esc}}$ with specific star formation
rate, age or $\Sigma_{\mathrm{SFR}}$ (Pahl et al., 2022, 2023).

In the context of our mass-dependent UV scatter model, we combine these
observational and theoretical insights with our UV LFs to assess the impact of
stochasticity on the timeline of reionization. In particular, as discussed by
Tang et al. (2019), Naidu et al. (2022) and Endsley et al. (2023b), among
others, we consider that galaxies experiencing a burst of star formation, i.e.
that will be up-scattered in the $M_{\rm UV}-M_{h}$ distribution, could have
an outsized contribution to the ionizing photon budget.

This effect can be translated into $\xi_{\mathrm{ion}}$ and $f_{\rm esc}$
being described by a distribution at fixed UV magnitudes rather than fixed
values for every galaxy. The production rate of ionizing photons available for
reionization can therefore be expressed as:

$\dot{n}_{\rm ion}(z)=\int_{-\infty}^{M_{lim}}dM_{\rm UV}L(M_{\rm
UV})\Phi(M_{\rm UV})\\\ \times f_{\rm esc}\,\xi_{\rm ion}\,p(f_{\rm
esc},\xi_{\rm ion}\mid M_{\rm UV});$ (4)

which can be solved when an empirical relation between $\xi_{\rm ion}$ and
$M_{\rm UV}$ is available.

Within our modeling framework, we model $\xi_{\mathrm{ion}}$ and $f_{\rm esc}$
as functions of $\Delta M_{\rm UV}=M_{\rm UV,c}(M_{h})-M_{\rm UV}$. As
illustrated in the middle panel of Fig. 5, we assume $\log_{10}\xi_{\rm ion}$
increases linearly with $\Delta M_{\rm UV}$, which naturally produces a log-
normal distribution of $p(\xi_{\rm ion}\mid M_{\rm UV})$, similar to
observations (e.g., Prieto-Lyon et al., 2023; Endsley et al., 2023b). For the
escape fraction we explore three different scenarios: (i) exponentially
increasing $f_{\rm esc}(\Delta M_{\rm UV})$, which produces an exponential
distribution of $p(f_{\rm esc}\mid M_{\rm UV})$ (similar to that derived for
$z\sim 3.5$ star forming galaxies by Kreilgaard et al., 2024); (ii) linearly
increasing $f_{\rm esc}(\Delta M_{\rm UV})$; and (iii) constant $f_{\rm
esc}(\Delta M_{\rm UV})$ for all galaxies.

In Fig. 5, left, we show the $M_{\rm UV}-M_{h}$ relation at $z=9$, populated
following the halo mass function and our mass-dependent
$\sigma_{\mathrm{UV}}(M_{h})$ model. Each galaxy is coloured with the
corresponding product of escape fraction and ionizing efficiency for the
exponentially increasing $f_{\rm esc}$ case.

Using these assumptions, we can derive $\dot{n}_{\rm ion}(z)$ as:

$\dot{n}_{\rm ion}(z)=\int_{-\infty}^{M_{lim}}dM_{\rm UV}L(M_{\rm
UV})\Phi(M_{\rm UV})\\\ \times\int_{M_{h,min}}^{\infty}dM_{h}\,f_{\rm
esc}[\Delta M_{\rm UV}(M_{h})]\,\xi_{\rm ion}[(\Delta M_{\rm UV}(M_{h})]\\\
\times p(M_{\rm UV}\mid M_{h})\,\frac{dn}{dM_{h}};$ (5)

where we set $M_{\mathrm{lim}}=-13$, but note that our results are relatively
insensitive to the exact value of this choice due to the turn-over of the UV
scatter LFs at $M_{\rm UV}\sim-13$ (Fig. 2). The ionizing photon rate
production can be then used to calculate the cosmic evolution of the ionized
hydrogen fraction by solving Eq. (6) in Mason et al. (2015).

The resulting reionization histories for the three scenarios are shown on the
right of Fig. 5. We also compare to the no-scatter case, for which we assume
fixed $f_{\rm esc}$ and $\xi_{\rm ion}=10^{25.7}\rm Hz\,erg^{-1}$ (consistent
with the latest JWST results; e.g. Endsley et al. 2023a; Atek et al. 2024).
For each model, $f_{\rm esc}$ is tuned to obtain the end of reionization by
$z\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}6$, as inferred from
Lyman-alpha forest observations (Qin et al., 2021; Bosman et al., 2022; Jin et
al., 2023). In the no-scatter scenario this leads to a value of $f_{\rm
esc}=0.08$. As discussed by Muñoz et al. (2024) this may be in tension with
estimates from low-$z$ analogs which would predict higher $f_{\rm esc}$.

When adding stochastic star formation, we obtain an earlier onset and more
gradual evolution of reionization, as low mass halos can start producing
significant ionizing photons at early times (see also Furlanetto & Mirocha,
2022). We find the shape of the reionization timeline is strongly dependent on
our $f_{\mathrm{esc}}$ model. The steeper the dependency of $f_{\mathrm{esc}}$
with $\Delta M_{\rm UV}(M_{h})$ (i.e. going from the constant to the
exponential scenario), the earlier the start and the shallower the evolution
of the reionization process.

At face-value, stochastic star formation exacerbates the potential ionizing
photon crisis: due to the higher luminosity densities in our scatter model, we
require an even lower value of $f_{\rm esc}\sim 0.05$ than the no-scatter case
($f_{\rm esc}\sim 0.08$) to obtain a $z\sim 6$ end of reionization. We note
here that our linear $\xi_{\mathrm{ion}}$ model also assumes
$\xi_{\mathrm{ion}}$ is higher for starbursting galaxies, but even when using
a constant $\xi_{\mathrm{ion}}$ model the higher luminosity density of the
scatter model implies a higher injection rate of ionizing photons at $z\gtrsim
8$ compared to the no-scatter case. This motivates the need to consider a
distribution of escape fractions and $\xi_{\mathrm{ion}}$, different for
galaxies in different stages of their bursty SFH, to properly model early
galaxies’ contribution to reionization.

We note however that while there are observational constraints on
$\xi_{\mathrm{ion}}$ at high-$z$, the escape fraction is not directly
constrained. Thus our predictions here are merely illustrative, but
demonstrate that stochasticity can have a marked impact on the timing of
reionization and is an important additional constraint to understand early
star formation, and can be constrained now to $z\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10$ with JWST observations (e.g.,
Nakane et al., 2023; Umeda et al., 2023).

As stochasticity implies lower mass halos (which are less biased, see above)
are more important for reionization this would also impact the morphology of
reionization. We would expect more numerous, smaller bubbles at fixed neutral
fraction compared to a case with no stochasticity (e.g. Hutter et al., 2021,
2023; Lu et al., 2024). This could potentially be measured via Lyman-alpha
transmission and 21-cm intensity mapping.

####  3.3.3 Spectral Energy Distributions

Figure 6: Examples of star formation histories (left) and continuum spectral
energy distributions (right) for three galaxies at $z=10$ with same UV
magnitude $M_{\rm UV}=-20.$

Highly time variable, bursty star formation histories that can produce drastic
variability of the luminosity of a galaxy across its evolution, should also
cause large scatter in galaxies’ spectral features when we observe galaxies at
fixed $M_{\rm UV}$ or stellar mass, as we are observing galaxies at different
phases in their SFHs (e.g. Iyer et al., 2024; Sun et al., 2023a; Endsley et
al., 2023b). Different spectral features are sensitive to different timescales
and galaxy properties, and are commonly used to probe galaxy star formation
histories (e.g. Weisz et al., 2012; Sparre et al., 2017; Ciesla et al., 2023;
Cole et al., 2023; Topping et al., 2024a; Langeroodi & Hjorth, 2024). In the
context of understanding the role of stochasticity at high-$z$ it is important
to identify which are the best observable features to help us constrain
stochasticity

JWST makes it possible to measure variations across galaxy populations at
$z\sim 5-15$. Several recent works have quantified stochasticity at $z\gtrsim
5$ via the scatter in the star-formation main sequence using JWST photometry
(SFR - stellar mass, e.g., Ciesla et al., 2023; Cole et al., 2023; Caputi et
al., 2023), finding evidence for a potential increase in scatter with
increasing redshift. However, it becomes increasingly challenging to measure
accurate stellar masses from photometry for bursty galaxies, due to the
‘outshining’ of rest-optical stellar emission by nebular emission from younger
stars (e.g., Stark et al., 2013; Whitler et al., 2023; Giménez-Arteaga et al.,
2023; Narayanan et al., 2024). Analyses at fixed UV magnitude may provide more
reliable comparisons with models during this bursty epoch of galaxy formation
(e.g., Endsley et al., 2023b).

To this aim, we explore the impact of simple bursty star formation histories
on galaxy SEDs at fixed $M_{\rm UV}$, to avoid the issues of stellar mass
uncertainties.

Following Mason et al. (2015), we set the average star formation rate of a
halo of mass $M_{h}$ as $\langle{\rm
SFR}\rangle\propto\epsilon_{\star}(M_{h})\,f_{b}M_{h}$. On top of this
constant evolution we artificially add a “burstiness” by modeling the SFH as a
periodic pattern in which burst-phases and lull-phases alternates.
Specifically the SFR varies between $\rm\log\langle SFR\rangle+\Delta\log
SFR$, during bursts, and $\rm\log\langle SFR\rangle-\Delta\log SFR$, during
lulls. The value of $\rm\Delta SFR$ and the duration of burst and lull phases
is calibrated so that the stellar mass assembled over the total period remains
the same as for the constant evolution.

To check how this burstiness affects the emission of galaxies, we produce
synthetic continuum SEDs. We sample the SFH as a series of single bursts every
1 Myr and use starburst99 (Leitherer et al., 1999) with padova stellar tracks
(Bertelli et al., 1994) and assuming a Kroupa (2001) initial mass function
between 0.1 and 100 $M_{\odot}$ to derive the rest-frame spectra and $M_{\rm
UV}$ (see also Gelli et al., 2021).

We then select galaxies at the same magnitude $M_{\rm UV}$, finding that they
can show a wide range of different SFH and properties due to their stochastic
behaviour.

In Fig. 6, left panel, we show as an example the SFHs of three galaxies with
same magnitude $M_{\rm UV}=-20$ at $z=10$, but hosted in halos of different
mass. The one in the lowest mass halos is experiencing a burst of star
formation, and is therefore up-scattered in the $M_{\rm UV}-M_{h}$ relation,
whereas the ones in higher mass halo are in a lull phase.

In the right panel of Fig. 6 we show with the same color coding their
redshifted observed flux $z=10$. By definition, their rest-frame 1500A∘
luminosity is the same, but we can see that significant differences are
present in their UV slopes $\beta_{UV}$ and Balmer breaks at
$\lambda_{rest}\sim
3645\textrm{A\kern-5.59721pt\raisebox{2.58334pt}{${}^{\circ}$}}$. Galaxies
experiencing a recent rise in the SFR have steeper UV slopes, since their
emission is dominated by young massive stars. On the other hand, galaxies
undergoing a constant or recently decreasing period of SFR have shallower UV
slopes and much more prominent Balmer breaks.

While not modelled here, we also expect the strength of nebular emission lines
that trace star formation on short timescales ($<10$ Myr), such as $H\alpha$,
to be enhanced during bursty phases (see e.g., Weisz et al., 2012; Sparre et
al., 2017). As a consequence, the $\beta_{UV}$, Balmer breaks and line
equivalent widths represent quantities that can give us precious direct
insight on bursty star formation (see also e.g. Endsley et al., 2023b;
Trussler et al., 2024; Langeroodi & Hjorth, 2024; Topping et al., 2024a).
Larger samples of galaxies are expected to show wide distributions of these
quantities at fixed magnitude. In a forthcoming work, we will extend our
predictions, analyzing in detail the expected distributions of $\beta_{UV}$
slopes and Balmer breaks, as well as of emission line equivalent widths and
other spectral features, all essential to constrain the bursty behavior of
high-$z$ galaxies and to infer the timescales and physical processes
regulating them.

##  4 discussion: The abundance of galaxies at $z>12$

Despite the success of the mass-dependent UV scatter model in naturally
producing higher UVLFs at progressively higher redshifts, it still cannot
reproduce the observed high number densities of galaxies at $z>12$.

Our analysis was based on some fundamental assumptions regarding: (i) the
cosmological parameters, (ii) the halo mass function, and (iii) the form used
for the matter transfer function. We tested how these assumptions impact the
predicted LFs in Appendix A, and find that their effect is mostly negligible
(with some exception for the transfer function).

Thus, the persisting discrepancy at $z>12$ suggests that stochastic star
formation alone is not sufficient to fill the gap and that additional physical
processes may be at play.

As mentioned in the introduction, another way to boost the abundance of
high-$z$ galaxies is through an enhancement of the average $M_{\rm UV}-M_{h}$
relation. This could be achieved if galaxies are forming stars in a more
efficient manner at high-$z$. Our model currently employs a redshift
independent star formation efficiency, a strong assumption that likely needs
to be relaxed to match the $z>12$ observations.

The fact that our model reflects a stochasticity proportional to the inverse
of the halo escape velocity offers an interesting insight. As
$v_{\mathrm{esc}}\sim M_{h}^{1/3}(1+z)^{1/2}$ we should expect a decrease in
stochasticity at fixed halo mass at higher redshifts, as halos are more
compact and thus better able to retain their gas. This implies higher redshift
halos should host higher gas masses at fixed halo mass, which may provide fuel
for more efficient star formation.

As proposed by Dekel et al. (2023), a possible cause for higher star formation
efficiency at high redshift is an increase in the density of star forming
regions. Cloud-scale simulations show that when a surface density threshold is
reached ($\Sigma\gtrsim 10^{3}\,M_{\odot}$ pc-2) stellar feedback may not be
fast enough to counteract gravitational collapse and the star formation
efficiency can be extremely enhanced from usual values of $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10\%$ to $\sim 100\%$ (e.g. Grudić
et al., 2018; Lancaster et al., 2021; Menon et al., 2023, 2024). Galaxy-scale
hydrodynamic simulations show that to create such high density clouds likely
requires a turbulent gas-rich environment where gas is compressed by feedback-
driven winds and shocks from inflowing gas (Ma et al., 2020b). Dekel et al.
(2023) and Li et al. (2023) discuss that we may expect more galaxies at high
redshift to satisfy these high star formation efficiency conditions as halos
become increasingly compact and accretion rates increase.

Indeed, observations show galaxies appear to be increasingly compact at $z>4$
(e.g. Langeroodi & Hjorth, 2023; Morishita et al., 2024). Furthermore,
extremely high electron densities ($n_{e}\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10^{4}$ cm-3) have been detected in
two UV-bright, compact $z\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}6$
sources (GNz11 and RXCJ2248-ID) which show UV emission line signatures of
massive stars (Senchyna et al., 2023; Topping et al., 2024b, though c.f.
Maiolino et al. 2023 where the high electron density of GNz11 is attributed to
an AGN broad-line region). Topping et al. (2024a) suggest such high electron
densities could result from the high density birth clouds described above,
which can fuel efficient star formation. Larger samples of $z>6$ sources with
deep, high-resolution rest-UV spectroscopy trace both massive star signatures
and ionized gas densities, combined with morphological information, should
enable a better understanding of the links between efficient star formation
and high densities.

##  5 Conclusions

We have explored the impact of a mass-dependent scatter of the $M_{\rm
UV}-M_{h}$ relation on the UV luminosity functions of high-$z$ galaxies, and
other observables, to understand the relevance of bursty star formation in the
early Universe. We assumed a decreasing trend $\sigma_{\rm UV}(M_{h})\sim
M_{h}^{-1/3}$ with increasing $M_{h}$, reflecting the higher sensitivity of
low-mass halos to stochastic star formation processes due to their lower halo
escape velocities. Our main conclusions are the following:

  * •

The mass-dependent UV scatter model naturally produces higher UV LFs at
progressively higher redshifts with respect to the no-stochasticity case, due
to low-mass halos being more dominant at earlier epochs;

  * •

The predicted LF and luminosity density are in agreement with the latest JWST
observations up to $z\sim 12$, but start to diverge at progressively earlier
epochs, suggesting the need for additional UV luminosity production on top of
that provided by stochastic star formation. This may be linked to the
increasing compactness, and thus density, of halos at fixed mass with
redshift;

  * •

The model predicts significantly lower galaxy bias values than the no-scatter
case as a function of redshift, with the models diverging at $z\gtrsim 8$ for
$M_{\rm UV}\sim-21$ galaxies;

  * •

Assuming that up-scattered galaxies are the dominant source of ionizing
photons yields an earlier start to reionization and a more gradual
reionization timeline than a model without scatter;

  * •

Balmer breaks, nebular emission line strengths and UV-slopes are also
sensitive to bursty star formation histories and can have significantly
different values for galaxies with same $M_{\rm UV}$.

Our results show that bursty star formation leading to a mass-dependent
stochasticity in the UV luminosity - halo mass relation is expected to have an
impact on multiple observable quantities. JWST observations are probing the
high-$z$ UV luminosity functions with increasing depth and precision, but, as
discussed by Muñoz et al. (2023), the contribution of bursty star formation
cannot be constrained by the LF alone due to degeneracies with other physical
effects. For this reason, it is essential to use multiple independent
observables simultaneously to definitively understand the role of bursty star
formation in high-$z$ galaxies’ evolution (see also e.g., Mirocha &
Furlanetto, 2023; Endsley et al., 2023b; Topping et al., 2024a; Langeroodi &
Hjorth, 2024; Trussler et al., 2024).

Since our model assumes a redshift-independent star formation efficiency and
still struggles to reproduce the extremely high redshift observations at
$z>12$, there is a clear need for further investigation of the physics of the
earliest galaxies to understand how they formed the first stellar populations.

##

We thank Andrea Ferrara, Anne Hutter, Julian Muñoz and Stefania Salvadori for
useful discussions, and Dan Stark for providing comments on a draft of the
manuscript. VG and CAM acknowledge support by the Carlsberg Foundation under
grant CF22-1322. The Cosmic Dawn Center (DAWN) is funded by the Danish
National Research Foundation under grant DNRF140. The model LFs can be
obtained upon request to the authors and will be available online upon
publication.

##  Appendix A The impact of cosmological assumptions on the UVLF

To make our predictions for the UVLFs, we had to make some assumptions
regarding: (i) the cosmological parameters, (ii) the halo mass function, and
(iii) the form used for the matter transfer function. Previous work modelling
high-$z$ LFs have assumed a variety of choices for the above. To assess the
impact of these choices on our results, we have tested our UV-scatter model
results varying these set of assumptions one at a time.

By changing the cosmological parameters from those of a flat $\Lambda$CDM with
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$ and $h=0.7$, to the Planck
Collaboration et al. (2020) measurements through the cosmic microwave
background (CMB) anisotropies ($\Omega_{m}=0.315$, $\Omega_{\Lambda}=0.685$
and $h=0.67$), we find a completely negligible change in the LF.

For the HMF, we use the hmf python package by Murray et al. (2013) and assumed
a Reed et al. (2007) function, which had been simulated specifically for halos
at $z\sim 10-30$. Using instead the functions from Sheth et al. (2001) and
Behroozi et al. (2013) results in slightly higher and lower number densities
respectively, with the maximum change produced in the LFs of $\sim 0.5$ dex
(see Figure 7, left panel). We can conclude that the choice of the HMF
produces a minor effect on the LFs, but that our overall conclusions do not
depend on its choice.

Finally, we test the impact of the assumed linear transfer function by
changing it from the Eisenstein & Hu (1998) form to the one calculated with
the Code for Anisotropies in the Microwave Background (CAMB, Lewis et al.,
2000), both included in the hmf package. This time, the change in the LFs is
significant, especially towards the faint end, with the LF being higher by up
to $\sim 1$ dex when using CAMB. This could be linked to an excess of the
power spectrum at very small scales in the CAMB function implementation,
and/or possibly due to the default extrapolation of the CAMB transfer function
at high wavenumber in the hmf code333https://github.com/halomod/hmf/issues/90.
We therefore caution that care must be taken when using the default CAMB
option for the transfer function in hmf as it can lead to an overestimation of
the LF at high-$z$.

Figure 7: UVLF at $z=7,10,14$ for different halo mass functions (left) and
transfer functions (right).

## References

  * Adams et al. (2023) Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755, doi: 10.1093/mnras/stac3347
  * Atek et al. (2023a) Atek, H., Shuntov, M., Furtak, L. J., et al. 2023a, MNRAS, 519, 1201, doi: 10.1093/mnras/stac3144
  * Atek et al. (2023b) Atek, H., Chisholm, J., Alavi, A., et al. 2023b, JWST’s GLIMPSE: gravitational lensing & NIRCam imaging to probe early galaxy formation and sources of reionization, JWST Proposal. Cycle 2, ID. #3293
  * Atek et al. (2024) Atek, H., Labbé, I., Furtak, L. J., et al. 2024, Nature, 626, 975, doi: 10.1038/s41586-024-07043-6
  * Barone-Nugent et al. (2014) Barone-Nugent, R. L., Trenti, M., Wyithe, J. S. B., et al. 2014, ApJ, 793, 17, doi: 10.1088/0004-637X/793/1/17
  * Barrow et al. (2020) Barrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, ApJL, 902, L39, doi: 10.3847/2041-8213/abbd8e
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Bertelli et al. (1994) Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994, A&A Supp., 106, 275
  * Bolan et al. (2022) Bolan, P., Lemaux, B. C., Mason, C., et al. 2022, MNRAS, 517, 3263, doi: 10.1093/mnras/stac1963
  * Bosman et al. (2022) Bosman, S. E. I., Davies, F. B., Becker, G. D., et al. 2022, MNRAS, 514, 55, doi: 10.1093/mnras/stac1046
  * Bouwens et al. (2022a) Bouwens, R. J., Illingworth, G., Ellis, R. S., Oesch, P., & Stefanon, M. 2022a, ApJ, 940, 55, doi: 10.3847/1538-4357/ac86d1
  * Bouwens et al. (2014) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115, doi: 10.1088/0004-637X/793/2/115
  * Bouwens et al. (2022b) Bouwens, R. J., Smit, R., Schouws, S., et al. 2022b, ApJ, 931, 160, doi: 10.3847/1538-4357/ac5a4a
  * Bouwens et al. (2023) Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036, doi: 10.1093/mnras/stad1145
  * Bromm & Yoshida (2011) Bromm, V., & Yoshida, N. 2011, ARA&A, 49, 373, doi: 10.1146/annurev-astro-081710-102608
  * Caputi et al. (2023) Caputi, K. I., Rinaldi, P., Iani, E., et al. 2023, arXiv e-prints, arXiv:2311.12691, doi: 10.48550/arXiv.2311.12691
  * Castellano et al. (2022) Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJL, 938, L15, doi: 10.3847/2041-8213/ac94d0
  * Chevallard et al. (2018) Chevallard, J., Charlot, S., Senchyna, P., et al. 2018, MNRAS, 479, 3264, doi: 10.1093/mnras/sty1461
  * Ciesla et al. (2023) Ciesla, L., Elbaz, D., Ilbert, O., et al. 2023, arXiv e-prints, arXiv:2309.15720, doi: 10.48550/arXiv.2309.15720
  * Cole et al. (2023) Cole, J. W., Papovich, C., Finkelstein, S. L., et al. 2023, arXiv e-prints, arXiv:2312.10152, doi: 10.48550/arXiv.2312.10152
  * Cullen et al. (2024) Cullen, F., McLeod, D. J., McLure, R. J., et al. 2024, MNRAS, doi: 10.1093/mnras/stae1211
  * Curti et al. (2023) Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2304.08516, doi: 10.48550/arXiv.2304.08516
  * Dalmasso et al. (2024a) Dalmasso, N., Leethochawalit, N., Trenti, M., & Boyett, K. 2024a, arXiv e-prints, arXiv:2402.18052, doi: 10.48550/arXiv.2402.18052
  * Dalmasso et al. (2024b) Dalmasso, N., Trenti, M., & Leethochawalit, N. 2024b, MNRAS, 528, 898, doi: 10.1093/mnras/stad3901
  * Davies et al. (2018) Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, 864, 142, doi: 10.3847/1538-4357/aad6dc
  * Dekel et al. (2023) Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201, doi: 10.1093/mnras/stad1557
  * Dome et al. (2023) Dome, T., Tacchella, S., Fialkov, A., et al. 2023, arXiv e-prints, arXiv:2305.07066, doi: 10.48550/arXiv.2305.07066
  * Donnan et al. (2024) Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, arXiv e-prints, arXiv:2403.03171, doi: 10.48550/arXiv.2403.03171
  * Eisenstein & Hu (1998) Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605, doi: 10.1086/305424
  * Endsley et al. (2020) Endsley, R., Behroozi, P., Stark, D. P., et al. 2020, MNRAS, 493, 1178, doi: 10.1093/mnras/staa324
  * Endsley et al. (2023a) Endsley, R., Stark, D. P., Whitler, L., et al. 2023a, MNRAS, 524, 2312, doi: 10.1093/mnras/stad1919
  * Endsley et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.05295, doi: 10.48550/arXiv.2306.05295
  * Ferrara (2023) Ferrara, A. 2023, arXiv e-prints, arXiv:2310.12197, doi: 10.48550/arXiv.2310.12197
  * Finkelstein et al. (2023) Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2023, arXiv e-prints, arXiv:2311.04279, doi: 10.48550/arXiv.2311.04279
  * Fiore et al. (2023) Fiore, F., Ferrara, A., Bischetti, M., Feruglio, C., & Travascio, A. 2023, ApJL, 943, L27, doi: 10.3847/2041-8213/acb5f2
  * Flury et al. (2022) Flury, S. R., Jaskot, A. E., Ferguson, H. C., et al. 2022, ApJ, 930, 126, doi: 10.3847/1538-4357/ac61e4
  * Furlanetto & Mirocha (2022) Furlanetto, S. R., & Mirocha, J. 2022, MNRAS, 511, 3895, doi: 10.1093/mnras/stac310
  * Gelli et al. (2024) Gelli, V., Salvadori, S., Ferrara, A., & Pallottini, A. 2024, ApJ, 964, 76, doi: 10.3847/1538-4357/ad23ec
  * Gelli et al. (2021) Gelli, V., Salvadori, S., Ferrara, A., Pallottini, A., & Carniani, S. 2021, ApJL, 913, L25, doi: 10.3847/2041-8213/abfe6c
  * Gelli et al. (2023) —. 2023, ApJL, 954, L11, doi: 10.3847/2041-8213/acee80
  * Gelli et al. (2020) Gelli, V., Salvadori, S., Pallottini, A., & Ferrara, A. 2020, MNRAS, 498, 4134, doi: 10.1093/mnras/staa2410
  * Giménez-Arteaga et al. (2023) Giménez-Arteaga, C., Oesch, P. A., Brammer, G. B., et al. 2023, ApJ, 948, 126, doi: 10.3847/1538-4357/acc5ea
  * Greig et al. (2019) Greig, B., Mesinger, A., & Bañados, E. 2019, MNRAS, 484, 5094, doi: 10.1093/mnras/stz230
  * Grudić et al. (2018) Grudić, M. Y., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 475, 3511, doi: 10.1093/mnras/sty035
  * Harikane et al. (2016) Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123, doi: 10.3847/0004-637X/821/2/123
  * Harikane et al. (2022) Harikane, Y., Ono, Y., Ouchi, M., et al. 2022, ApJS, 259, 20, doi: 10.3847/1538-4365/ac3dfc
  * Harikane et al. (2023) Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5, doi: 10.3847/1538-4365/acaaa9
  * Hopkins et al. (2023) Hopkins, P. F., Gurvich, A. B., Shen, X., et al. 2023, MNRAS, 525, 2241, doi: 10.1093/mnras/stad1902
  * Hutter et al. (2021) Hutter, A., Dayal, P., Legrand, L., Gottlöber, S., & Yepes, G. 2021, MNRAS, 506, 215, doi: 10.1093/mnras/stab877
  * Hutter et al. (2023) Hutter, A., Trebitsch, M., Dayal, P., et al. 2023, MNRAS, 524, 6124, doi: 10.1093/mnras/stad2230
  * Inayoshi et al. (2022) Inayoshi, K., Harikane, Y., Inoue, A. K., Li, W., & Ho, L. C. 2022, ApJL, 938, L10, doi: 10.3847/2041-8213/ac9310
  * Iyer et al. (2024) Iyer, K. G., Speagle, J. S., Caplar, N., et al. 2024, ApJ, 961, 53, doi: 10.3847/1538-4357/acff64
  * Iyer et al. (2020) Iyer, K. G., Tacchella, S., Genel, S., et al. 2020, MNRAS, 498, 430, doi: 10.1093/mnras/staa2150
  * Izotov et al. (2016) Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016, MNRAS, 461, 3683, doi: 10.1093/mnras/stw1205
  * Izotov et al. (2018) Izotov, Y. I., Worseck, G., Schaerer, D., et al. 2018, MNRAS, 478, 4851, doi: 10.1093/mnras/sty1378
  * Jin et al. (2023) Jin, X., Yang, J., Fan, X., et al. 2023, ApJ, 942, 59, doi: 10.3847/1538-4357/aca678
  * Kimm & Cen (2014) Kimm, T., & Cen, R. 2014, ApJ, 788, 121, doi: 10.1088/0004-637X/788/2/121
  * Kravtsov & Belokurov (2024) Kravtsov, A., & Belokurov, V. 2024, arXiv e-prints, arXiv:2405.04578, doi: 10.48550/arXiv.2405.04578
  * Kreilgaard et al. (2024) Kreilgaard, K. C., Mason, C. A., Cullen, F., Begley, R., & McLure, R. J. 2024.  https://arxiv.org/abs/2405.10364
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Lancaster et al. (2021) Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021, ApJL, 922, L3, doi: 10.3847/2041-8213/ac3333
  * Langeroodi & Hjorth (2023) Langeroodi, D., & Hjorth, J. 2023, arXiv e-prints, arXiv:2307.06336, doi: 10.48550/arXiv.2307.06336
  * Langeroodi & Hjorth (2024) —. 2024, arXiv e-prints, arXiv:2404.13045, doi: 10.48550/arXiv.2404.13045
  * Legrand et al. (2022) Legrand, L., Hutter, A., Dayal, P., et al. 2022, MNRAS, 509, 595, doi: 10.1093/mnras/stab3034
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473, doi: 10.1086/309179
  * Li et al. (2023) Li, Z., Dekel, A., Sarkar, K. C., et al. 2023, arXiv e-prints, arXiv:2311.14662, doi: 10.48550/arXiv.2311.14662
  * Looser et al. (2023a) Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023a, arXiv e-prints, arXiv:2302.14155.  https://arxiv.org/abs/2302.14155
  * Looser et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.02470, doi: 10.48550/arXiv.2306.02470
  * Lu et al. (2024) Lu, T.-Y., Mason, C. A., Hutter, A., et al. 2024, MNRAS, 528, 4872, doi: 10.1093/mnras/stae266
  * Ma et al. (2020a) Ma, X., Quataert, E., Wetzel, A., et al. 2020a, MNRAS, 498, 2001, doi: 10.1093/mnras/staa2404
  * Ma et al. (2018) Ma, X., Hopkins, P. F., Garrison-Kimmel, S., et al. 2018, MNRAS, 478, 1694, doi: 10.1093/mnras/sty1024
  * Ma et al. (2020b) Ma, X., Grudić, M. Y., Quataert, E., et al. 2020b, MNRAS, 493, 4315, doi: 10.1093/mnras/staa527
  * Madau et al. (1999) Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648, doi: 10.1086/306975
  * Maiolino et al. (2023) Maiolino, R., Scholtz, J., Witstok, J., et al. 2023, arXiv e-prints, arXiv:2305.12492, doi: 10.48550/arXiv.2305.12492
  * Mason et al. (2015) Mason, C. A., Trenti, M., & Treu, T. 2015, ApJ, 813, 21, doi: 10.1088/0004-637X/813/1/21
  * Mason et al. (2023) —. 2023, MNRAS, 521, 497, doi: 10.1093/mnras/stad035
  * Mason et al. (2018) Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, 856, 2, doi: 10.3847/1538-4357/aab0a7
  * Mason et al. (2019) Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947, doi: 10.1093/mnras/stz632
  * McGreer et al. (2015) McGreer, I. D., Mesinger, A., & D’Odorico, V. 2015, MNRAS, 447, 499, doi: 10.1093/mnras/stu2449
  * McLeod et al. (2024) McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2024, MNRAS, 527, 5004, doi: 10.1093/mnras/stad3471
  * Menon et al. (2023) Menon, S. H., Federrath, C., & Krumholz, M. R. 2023, MNRAS, 521, 5160, doi: 10.1093/mnras/stad856
  * Menon et al. (2024) Menon, S. H., Lancaster, L., Burkhart, B., et al. 2024, arXiv e-prints, arXiv:2405.00813, doi: 10.48550/arXiv.2405.00813
  * Meurer et al. (1999) Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64, doi: 10.1086/307523
  * Mirocha (2020) Mirocha, J. 2020, MNRAS, 499, 4534, doi: 10.1093/mnras/staa3150
  * Mirocha & Furlanetto (2023) Mirocha, J., & Furlanetto, S. R. 2023, MNRAS, 519, 843, doi: 10.1093/mnras/stac3578
  * Morales et al. (2024) Morales, A. M., Finkelstein, S. L., Leung, G. C. K., et al. 2024, ApJL, 964, L24, doi: 10.3847/2041-8213/ad2de4
  * Morishita et al. (2024) Morishita, T., Stiavelli, M., Chary, R.-R., et al. 2024, ApJ, 963, 9, doi: 10.3847/1538-4357/ad1404
  * Muñoz et al. (2024) Muñoz, J. B., Mirocha, J., Chisholm, J., Furlanetto, S. R., & Mason, C. 2024, arXiv e-prints, arXiv:2404.07250, doi: 10.48550/arXiv.2404.07250
  * Muñoz et al. (2023) Muñoz, J. B., Mirocha, J., Furlanetto, S., & Sabti, N. 2023, MNRAS, 526, L47, doi: 10.1093/mnrasl/slad115
  * Murray et al. (2013) Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astronomy and Computing, 3, 23, doi: 10.1016/j.ascom.2013.11.001
  * Naidu et al. (2022) Naidu, R. P., Matthee, J., Oesch, P. A., et al. 2022, MNRAS, 510, 4582, doi: 10.1093/mnras/stab3601
  * Nakane et al. (2023) Nakane, M., Ouchi, M., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2312.06804, doi: 10.48550/arXiv.2312.06804
  * Narayanan et al. (2024) Narayanan, D., Lower, S., Torrey, P., et al. 2024, ApJ, 961, 73, doi: 10.3847/1538-4357/ad0966
  * Ouchi et al. (2010) Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869, doi: 10.1088/0004-637X/723/1/869
  * Paardekooper et al. (2015) Paardekooper, J.-P., Khochfar, S., & Dalla Vecchia, C. 2015, MNRAS, 451, 2544, doi: 10.1093/mnras/stv1114
  * Pahl et al. (2022) Pahl, A. J., Shapley, A., Steidel, C. C., Reddy, N. A., & Chen, Y. 2022, MNRAS, 516, 2062, doi: 10.1093/mnras/stac1767
  * Pahl et al. (2023) Pahl, A. J., Shapley, A., Steidel, C. C., et al. 2023, MNRAS, 521, 3247, doi: 10.1093/mnras/stad774
  * Pallottini & Ferrara (2023) Pallottini, A., & Ferrara, A. 2023, A&A, 677, L4, doi: 10.1051/0004-6361/202347384
  * Pallottini et al. (2022) Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621, doi: 10.1093/mnras/stac1281
  * 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. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Prieto-Lyon et al. (2023) Prieto-Lyon, G., Strait, V., Mason, C. A., et al. 2023, A&A, 672, A186, doi: 10.1051/0004-6361/202245532
  * Qin et al. (2021) Qin, Y., Mesinger, A., Bosman, S. E. I., & Viel, M. 2021, MNRAS, 506, 2390, doi: 10.1093/mnras/stab1833
  * Qiu et al. (2018) Qiu, Y., Wyithe, J. S. B., Oesch, P. A., et al. 2018, MNRAS, 481, 4885, doi: 10.1093/mnras/sty2633
  * Rasmussen Cueto et al. (2023) Rasmussen Cueto, E., Hutter, A., Dayal, P., et al. 2023, arXiv e-prints, arXiv:2312.12109, doi: 10.48550/arXiv.2312.12109
  * Reed et al. (2007) Reed, D. S., Bower, R., Frenk, C. S., Jenkins, A., & Theuns, T. 2007, MNRAS, 374, 2, doi: 10.1111/j.1365-2966.2006.11204.x
  * Ren et al. (2019) Ren, K., Trenti, M., & Mason, C. A. 2019, ApJ, 878, 114, doi: 10.3847/1538-4357/ab2117
  * Senchyna et al. (2023) Senchyna, P., Plat, A., Stark, D. P., & Rudie, G. C. 2023, arXiv e-prints, arXiv:2303.04179, doi: 10.48550/arXiv.2303.04179
  * Shen et al. (2023) Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254, doi: 10.1093/mnras/stad2508
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1, doi: 10.1046/j.1365-8711.2001.04006.x
  * Sobacchi & Mesinger (2015) Sobacchi, E., & Mesinger, A. 2015, MNRAS, 453, 1843, doi: 10.1093/mnras/stv1751
  * Sparre et al. (2017) Sparre, M., Hayward, C. C., Feldmann, R., et al. 2017, MNRAS, 466, 88, doi: 10.1093/mnras/stw3011
  * Stark et al. (2013) Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129, doi: 10.1088/0004-637X/763/2/129
  * Steidel et al. (2018) Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123, doi: 10.3847/1538-4357/aaed28
  * Strait et al. (2023) Strait, V., Brammer, G., Muzzin, A., et al. 2023, ApJL, 949, L23, doi: 10.3847/2041-8213/acd457
  * 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
  * Tacchella et al. (2018) Tacchella, S., Bose, S., Conroy, C., Eisenstein, D. J., & Johnson, B. D. 2018, ApJ, 868, 92, doi: 10.3847/1538-4357/aae8e0
  * Tang et al. (2019) Tang, M., Stark, D. P., Chevallard, J., & Charlot, S. 2019, MNRAS, 489, 2572, doi: 10.1093/mnras/stz2236
  * Tang et al. (2021) Tang, M., Stark, D. P., Chevallard, J., et al. 2021, MNRAS, 503, 4105, doi: 10.1093/mnras/stab705
  * Tinker et al. (2010) Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878, doi: 10.1088/0004-637X/724/2/878
  * Topping et al. (2024a) Topping, M. W., Stark, D. P., Endsley, R., et al. 2024a, MNRAS, 529, 4087, doi: 10.1093/mnras/stae800
  * Topping et al. (2024b) Topping, M. W., Stark, D. P., Senchyna, P., et al. 2024b, MNRAS, 529, 3301, doi: 10.1093/mnras/stae682
  * Trinca et al. (2024) Trinca, A., Schneider, R., Valiante, R., et al. 2024, MNRAS, 529, 3563, doi: 10.1093/mnras/stae651
  * Trussler et al. (2024) Trussler, J., Conselice, C., Adams, N., et al. 2024, arXiv e-prints, arXiv:2404.07163, doi: 10.48550/arXiv.2404.07163
  * Umeda et al. (2023) Umeda, H., Ouchi, M., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2306.00487, doi: 10.48550/arXiv.2306.00487
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435, doi: 10.1146/annurev-astro-081817-051756
  * Weisz et al. (2012) Weisz, D. R., Johnson, B. D., Johnson, L. C., et al. 2012, ApJ, 744, 44, doi: 10.1088/0004-637X/744/1/44
  * Whitler et al. (2023) Whitler, L., Stark, D. P., Endsley, R., et al. 2023, MNRAS, 519, 5859, doi: 10.1093/mnras/stad004
  * Whitler et al. (2020) Whitler, L. R., Mason, C. A., Ren, K., et al. 2020, MNRAS, 495, 3602, doi: 10.1093/mnras/staa1178
  * Zackrisson et al. (2013) Zackrisson, E., Inoue, A. K., & Jensen, H. 2013, ApJ, 777, 39, doi: 10.1088/0004-637X/777/1/39


# Mixtape Application:
Last.fm Data Characterization

Luciana Fujii Pontello Pedro H. F. Holanda
Bruno Guilherme Joao Paulo V. Cardoso
Olga Goussevskaia Ana Paula C. da Silva

Computer Science Department
Universidade Federal de Minas Gerais (UFMG), Brazil

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
  2. 2 Related Work
  3. 3 Data collection
  4. 4 User profile
  5. 5 Songs and Artists
  6. 6 Tags
  7. 7 Conclusions

##  Chapter 1 Introduction

There is a never-ending debate about music and what happened to it after the
dawn of the streaming age. A few very quick revolutions happened between the
first appearance of MP3 files and illegal downloads and the current scenario,
where most artists have their full catalog on at least one free streaming
service, and many questions regarding the consequences of that process are
left without answers.

Back in the day, the only ways for people to listen to music were to turn on
the radio or, when choice was important, to buy large vinyl records that would
eventually be stored somewhere in the house, waiting to be chosen to leave the
comfort of their thick sleeves for a spin on the turntable. These records had
two sides, and required an inevitable pause between the first and the second
half of every record, since someone would have to go there and flip the
record. This single fact, for instance, had a strong influence on how records
were planned and how the sequence of songs was laid out.

Music albums were also a notable case of product tying, since they were a way
for record companies to sell a large number of songs at once even though many
customers only wanted a couple songs off each album. This situation started to
change with the release of singles, which were cheaper, smaller vinyl records
containing only a couple of songs, and, a little later, by the introduction of
portable cassette tape machines, which allowed people to copy songs from
different sources (radio, records or other tapes) to magnetic tapes, which
were cheap, small and practical.

This was the birth of the personal playlist. Tapes were the first time people
could piece together whatever tracks they wanted to listen to in whatever
order they wanted and then effectively listen to them, regardless of the
popularity of the artists involved, their likeliness to appear on the radio or
the price of their LPs. The CD came as a more practical alternative to that
system, but, down deep, the album mentality was the same.

Of course, these are baby steps when compared to the flexibility and speed of
the current music industry. Many people don’t even store music on their own
computers anymore and, in their cloud profiles, have access to music
collections so large that they would probably take up the space of a modern
house in terms of physical records. Entire albums can be bought, leased,
shared or “stolen” in a few minutes over broadband connection, and crafting
and delivering a mixed tape to a friend only takes a few clicks.

This abundance of material creates a very urgent need for services that
somehow make it easier to find music that is interesting and relevant. There
are many recommendation services that try to address this problem in different
ways. Some systems are able make taste-based recommendations based on usage
patterns, after a user has used them for long enough. Other systems are based
on curated lists of newly released and trending items or artists, such as the
user-made playlists on mixcloud111www.mixcloud.com or the professional dj-
curated streams on 22tracks22222tracks.com. Otherwise, the typically available
navigation functions in media players and online streaming services are mostly
based on filtering by attributes, like title, artist or genre, or return a
list of similar items, computed using collaborative filtering techniques.

Now that we’re witnessing the consolidation of these streaming services, it is
more important than ever to have a better understanding of how people listen
to music. This is crucial for the success of any music-related service,
especially the ones that intend to somehow recommend new music, or create
sequences or playlists for different people to listen to.

One of the simplest ways to do that is to look at how people behave in the
internet \- more specifically, on Online Social Networks (OSNs). These
websites are meant to extend the individual experience of society by bringing
social activities to the internet, allowing people to discuss their favorite
subjects and share their thoughts, opinions and feelings in a distributed and
independent way, making them an interesting source of many different kinds of
information.

Some of these networks carry information that is particularly valuable because
it does not require active behavior to gather itself. This way, it reflects
true habit and behavior that would otherwise have to be observed in an
intrusive way. A great example of an ingenious way to collect such information
is Last.Fm333www.last.fm:

Last.Fm is an online social network for music fans. It has a very ingenious
way of operating: Users deliberately install a lightweight crawler on their
personal computers to keep track of what they are listening to. That is called
scrobbling. Once users have their musical history up on the website, they can
interact with other users, make virtual “friends” (or find real ones that have
profiles) and add them to their own profiles, leave messages, navigate through
other users’ profiles to see what they listen to, among other things. It even
integrates with streaming services such as Spotify444www.spotify.com, allowing
users to upload their streaming histories directly to last.fm and listen to
songs.

Even though the way people listen to music has changed drastically over the
last few years, there is still little research characterizing these phenomena.
In this work, we focus our attention on analyzing how people share their music
listening habits with other people via the internet. In particular, we analyze
last.fm: We seek to characterize last.fm user behavior using our analysis to
shed light on how users interact in this Online Social Network (OSN), how
their preferences and activities may affect and be affected by content
popularity dynamics and, especially, how they listen to music.

We seek to answer the following key questions:

  * •

What is the Last.fm user profile, age, gender, location?

  * •

How is user activity distributed among content?

  * •

How is the distribution of tags associated with songs?

  * •

How can we use the knowledge of our data collection to improve the similarity
measures we intend to derive from it?

Our analysis reveals interesting details about the operation of last.fm. In
Particular, we show that users are young, and most are located in the United
States. We also learn that there is a large difference between the amount of
people who listen to popular songs and obscure songs, and that popularity
should be taken into account to build concepts of similarity between songs,
for recommendation purposes. And finally, we notice that users explore the tag
functionality of the network quite heavily, and that important information can
be derived from these tags and used for similarity and recommendation
purposes.

The report is structured as follows. We discuss related work about music
similarity and recommendation in chapter 2. In chapter 3 we present detailed
information about the dataset used for this study. In chapter 4 we describe
our detailed analysis of the user-related information in the dataset. In
chapter 5 we present our analysis of artists and tracks. In chapter 6 we
describe the analysis of tags. Finally, in chapter 7, we present our
conclusions.

##  Chapter 2 Related Work

The most relevant aspects of recommendation systems are the premises they are
based on, and how they build upon them. One of the most common premises is
content similarity, and that it would be desirable to recommend content that
is somehow similar to content that a given user enjoys. So, We present a small
compilation of related works in the area of Media Similarity, Music
Recommendation and Playlist Generation.

Media item-to-item similarity computation: Several research studies offer
technical solutions to define similarity between a given pair of items, and
they divide themselves among two schools. Some studies develop similarity
measures that are based solely on content (objective approach), being that
audio or video information — for instance, the spectral or rhythmic content of
songs. On the other hand, other studies develop similarity measures that are
based on user-generated data and tags, also known as collaborative filtering
(subjective approach). In the context of music, different approaches to define
item-to-item similarity have been studied extensively, such as content-based
measures [14, 2, 12, 7, 1, 21], that analyze spectral or rhythmic properties
of songs, and user-based measures [18, 9, 19, 8, 5, 10], that analyze user
listening habits in online social networks or user-generated tags [19]. Many
of these studies aim at building recommendation systems [14, 2, 21, 5, 20,
19].

Music playlist generation: Another related line of research to music
navigation systems is playlist generation. Several people have addressed this
problem from different perspectives. There are techniques that use statistical
analysis of radio streams [13, 20, 4, 3], are based on multidimensional metric
spaces [3, 8, 16, 15, 10], explore audio content [11], and user skipping
behavior [17].

In particular, authors in [17] create playlists based on audio music
similarity and skipping behavior, while authors in [6] use network flow
analysis to generate playlists from a friendship graph of artists on MySpace.
Maillet et al [13], in turn, present an approach to generating steerable
playlists from tags linked to songs played in professional radio station
playlists, and Chen et al [3] model playlists as Markov chains, which are
generated through the Latent Markov Embedding machine learning algorithm. The
heuristics in [17] use acoustic similarity and present linear complexity to
return the next item to the user, which is feasible to navigate in a
collection of 2,500 items, as the authors did, but does not scale to larger
collections, which is not only desirable but necessary for modern standards.

##  Chapter 3 Data collection

Last.fm provides a public API for collecting information about songs, artists,
albums and tags that have been contributed by millions of users.

The collection of Last.fm was comprised of 2 steps: First we collected 0.28%
of all songs from the Last.fm service (around 100 thousand songs) and their
top fans, and then we continued collecting the top songs of these fans and
their friends. We collected the top-25 most listened to songs of each user.

Our dataset was collected from November, 2014 to July, 2015, and contains
372,899 users (with their respective top-25 song lists), 2,060,173 songs, and
374,402 artists. Moreover, we also collected a total of 1,006,236 user-
generated tags, associated with songs. In particular, 47% of songs have had at
least one associated tag in our dataset.

From a total of 2,060,173 songs in our database, 983,010 have MusicBrainz
Identifiers (MBID)111MBID is a reliable and unambiguous form of music
identification (musicbrainz.org).. This is an important source of reliability,
since these entries are guaranteed to be unique, which eliminates duplicate
problems, and refer to the correct media items, which avoids the occasional
mismatch between the name of a song or artist in a user’s computer and the
actual correct name of the item, not to mention non-released live versions,
non-official covers and so on.

##  Chapter 4 User profile

We collected a total of 372,899 users from 237 countries. Table 4.1 shows the
countries with most users we have collected.

Country | # of users
---|---
US | 54248
BR | 31230
PL | 25882
RU | 22979
UK | 22515
DE | 18314
NL | 7237
CA | 6931
FI | 6356
FR | 6139
Table 4.1: Top countries with most users

Many users did not declare their gender. Of those who did, 201,096 said they
were male and 108,034 said they were female. Table 4.2 shows the number of
average playcounts per gender. We can see that the average number of
playcounts of male users is bigger than that of female users. Users that have
not declared their gender have the smallest average playcount number.

Gender | Average playcounts
---|---
Female | 82907.0770
Male | 113573.4874
Not informed | 56797.2135
Table 4.2: Average playcount per gender

The age range reported by the users ranges from 0 to 115 years. Figure 4.1
shows the distribution of users’ ages. We can see that the great majority of
Last.fm users are from 18 to 30 years old. Over 113,000 users did not report
their age.

Figure 4.1: Last.fm: Users’ age distribution. Figure 4.2: Last.fm: CDF of
users’ playcounts.

The Figure 4.2 shows that approximately 27% of the users have listened to less
than 10,000 songs, whereas 62% of the users have listened to 10,000 to 100,000
songs. In general users listen to a large amount of songs, making a vast
listening history, showing that their engagement with Last.fm service lasts
for a significant period of time.

Figure 4.3: Last.fm: CDF of users’ friends.

Figure 4.3 shows the CDF of users’ friends. Approximately 5% of the users have
zero friends in Last.fm, and 50% of the users have at most 25 friends. We have
only collected a maximum of 50 friends of each user, explaining the peak in
the 50 number of friends point in the plot. Approximately 30% of the users
have 50 friends or more.

##  Chapter 5 Songs and Artists

Sometimes, Last.fm registers different items which are actually the same song,
written in different ways. So, For this part of the analysis, we only took
into account the songs we could match with Musicbrainz IDs, in order to avoid
the risk of considering multiple issues related to track or artist name
mismatch or repetition.

Figures 5.2 and 5.2 show the Cumulative Density Function (CDF) of the
popularity of songs by listeners and playcount. We can observe that, in our
dataset, 50% of songs were listened to by at most 1,280 unique users (0.3%);
10%, in turn, attracted attention of more than 22,685 unique users (5.9%). The
most popular songs (top 1%) were listened to by more than 155,468 users (41%).
The top 10 songs by number of listeners can be seen in Table 5.1, each of them
with more than 1.4M listeners.

In Figure 5.2 we can see that 10% of the tracks have 1,000 playcounts or less
and more than 50% have more than 10,000 playcounts.

Track name | Artist | Number of listeners
---|---|---
Smells Like Teen Spirit | Nirvana | 1,806,180
Mr. Brightside | The Killers | 1,716,969
Wonderwall | Oasis | 1,685,703
Come as You Are | Nirvana | 1,597,611
Clocks | Coldplay | 1,507,981
Somebody Told Me | The Killers | 1,490,787
Take Me Out | Franz Ferdinand | 1,462,621
Karma Police | Radiohead | 1,431,055
Viva la Vida | Coldplay | 1,431,034
The Scientist | Coldplay | 1,404,877
Table 5.1: Top 10 songs

Figure 5.1: Last.fm: CDF of song popularity by listeners.

Figure 5.2: Last.fm: CDF of song popularity by playcount.

Figure 5.4 shows the Cumulative Density Function (CDF) of the popularity of
artists by their total listeners, retrieved from Last.fm. The top 1% most
popular artists ($\approx$ 3,577) were listened to by more than 764,000 users,
while approximately 20% of the artists were listened to by only 10 users.

Figure 5.4 shows the Artist CDF by the number of users that contain that
artist in their top 25, in the users’ history collected by us. More than 35%
of the artists appear in only one user’s top 25 songs, while 1% appear in more
than 600 user’s top 25 history.

Figure 5.3: Last.fm: CDF of artist popularity by number of total listeners.

Figure 5.4: Last.fm: CDF of artist popularity by users’ top 25.

Figure 5.5: Last.fm: CDF of song co-occurrence.

Figure 5.6: Last.fm: CDF of artist co-occurrence.

Figure 5.6 shows the co-occurrence CDF between songs that co-occurred at least
once. We define as a co-occurring relationship only the pair of songs that
appear in at least one user’s top 25 listened songs. We can see that more than
90% of the co-occurrences happen in only one user listening history, while
only 1% of the relationships co-occur 10 or more times. If we consider all the
NxN pairs of all songs, most of the songs don’t even co-occur in any user’s
top 25 songs. Therefore, initially, we can only infer the similarity between a
small percentage of the songs we have collected, and have to infer similarity
of the other songs from those similarities.

Figure 5.6 shows the co-occurrence CDF between artists that co-occurred at
least once. We define as a co-occurring relationship only the pair of artists
that appear in at least one user top 25 listened songs. We can see that 75% of
the artists co-occur only in one user’s top 25 songs, whereas approximately 1%
of the artists co-occur in 22 users’ top 25.

##  Chapter 6 Tags

Last.fm allows users to create and associate tags with the songs they listen
to. One tag can be associated to one song from one to 100 times, allowing the
users to show their agreement of a certain tag association by adding to it. We
collected a total of 1,006,236 user-generated tags, associated with songs. 47%
of the songs have had at least one associated tag in our dataset and
considering only the songs with MusicBrainz ID, 75% of the songs were
associated with at least one tag.

The top 5 most popular tags (rock, alternative, pop, indie and electronic)
were associated with 541,527 songs. Table 6.1 shows the most popular tags and
the number of songs associated with them.

Tag name | Tag count
---|---
rock | 167,610
alternative | 101,061
pop | 96,654
indie | 88,786
electronic | 87,416
alternative rock | 56,643
favorites | 56,508
beautiful | 51,870
love | 50,918
awesome | 42,364
Table 6.1: Top 10 tags

Figure 6.1 shows the Cumulative Density Function (CDF) of tag popularity using
the number of songs associated with the tag: 62% of tags have been associated
with only one song. In contrast, the top 5 most popular tags were associated
with more than 87,000 songs each.

Figure 6.1: Last.fm: CDF of tag popularity.

##  Chapter 7 Conclusions

In this technical report we presented a broad characterization of Last.fm, an
Online Social Network (OSN) for music fans. From that characterization we
understood more about the profile of Last.fm users and how they listen to
music and interact with the OSN. We also collected insights about the data we
have collected, which is a small fraction of the information available in
Last.fm database.

We learned that the majority of the users we have collected are young, from 18
to 30 years old, and from various countries, but the majority of users are
still in the US. This will probably impact our similarity measures, which will
be more accurate for the songs that appeal to this young audience and that are
popular in the US, because we will have more co-occurrence data about these
songs. If we want to have better similarity metrics for a specific country or
region, it is fair to assume that it will be beneficial to collect more data
from users from that region as a priority.

We discovered that there are songs that are listened by only a handful of
users, while the top songs are listened by thousands users. This means that we
will probably achieve better similarity measures for the most popular songs.
Based on that information, we could develop new methods to gather more data
about the less popular songs. For instance, we could collect those songs top
fans histories, in order to gather more co-occurrence data from these songs.
Or, in the future, we could try to ally our similarity measures based on
social data to other techniques, for example, techniques that use the analysis
of the content of the songs, to derive similarity measures for new or unknown
songs we want to recommend.

Users also heavily used tags to classify the songs they listen to. From tags
we learned that the most common genre of song in Last.fm is rock, the tag
associated with the most songs in the social network. Since tags are present
in most of the songs that contain a MusicBrainz ID, we can take advantage of
tags to improve our similarity metrics, for instance, using tags to derive
similarity of songs with zero co-occurrences to other songs. Or we can develop
navigation techniques that take advantage of tags. Finally, we can even use
tags to evaluate our similarity metrics and Euclidean space, by analyzing if
songs with a given tag tend to be close together or spread around the map.

## References

  * [1] J-J Aucouturier, François Pachet, and Mark Sandler.  “the way it sounds”: timbre models for analysis and retrieval of music signals.  Multimedia, IEEE Transactions on, 7(6):1028–1035, 2005.
  * [2] Pedro Cano, Markus Koppenberger, and Nicolas Wack.  Content-based music audio recommendation.  In Proceedings of the 13th annual ACM international conference on Multimedia, pages 211–212. ACM, 2005.
  * [3] Shuo Chen, Josh L Moore, Douglas Turnbull, and Thorsten Joachims.  Playlist prediction via metric embedding.  In 18th ACM SIGKDD, 2012.
  * [4] Shuo Chen, Jiexun Xu, and Thorsten Joachims.  Multi-space probabilistic sequence modeling.  In 19th ACM SIGKDD, 2013.
  * [5] Katayoun Farrahi, Markus Schedl, Andreu Vall, David Hauger, and Marko Tkalcic.  Impact of listening behavior on music recommendation.  ISMIR, 2014.
  * [6] Ben Fields, Kurt Jacobson, Christophe Rhodes, and Michael Casey.  Social playlists and bottleneck measurements: Exploiting musician social graphs using content-based dissimilarity and pairwise maximum flow values.  In ISMIR, 2008.
  * [7] Arthur Flexer.  On inter-rater agreement in audio music similarity.  In ISMIR, 2014.
  * [8] Olga Goussevskaia, Michael Kuhn, Michael Lorenzi, and Roger Wattenhofer.  From web to map: Exploring the world of music.  In Web Intelligence and Intelligent Agent Technology, 2008. WI-IAT’08, volume 1, 2008.
  * [9] Ioannis Konstas, Vassilios Stathopoulos, and Joemon M Jose.  On social networks and collaborative recommendation.  In Proceedings of the 32nd international ACM SIGIR Conference on Research and development in Information Retrieval, pages 195–202. ACM, 2009\.
  * [10] Michael Kuhn, Roger Wattenhofer, and Samuel Welten.  Social audio features for advanced music retrieval interfaces.  In Proceedings of the international conference on Multimedia, pages 411–420. ACM, 2010.
  * [11] Beth Logan.  Content-based playlist generation: Exploratory experiments.  In ISMIR, 2002.
  * [12] Beth Logan and Ariel Salomon.  A music similarity function based on signal analysis.  In ICME, 2001.
  * [13] François Maillet, Douglas Eck, Guillaume Desjardins, and Paul Lamere.  Steerable playlist generation by learning song similarity from radio station playlists.  In ISMIR, 2009.
  * [14] Brian McFee, Luke Barrington, and Gert Lanckriet.  Learning content similarity for music recommendation.  Audio, Speech, and Language Processing, IEEE Transactions on, 20(8):2207–2218, 2012.
  * [15] Joshua L Moore, Shuo Chen, Thorsten Joachims, and Douglas Turnbull.  Learning to embed songs and tags for playlist prediction.  In ISMIR, pages 349–354, 2012.
  * [16] Joshua L Moore, Shuo Chen, Douglas Turnbull, and Thorsten Joachims.  Taste over time: The temporal dynamics of user preferences.  In ISMIR, 2013.
  * [17] Elias Pampalk, Tim Pohle, and Gerhard Widmer.  Dynamic playlist generation based on skipping behavior.  In ISMIR, 2005.
  * [18] Minsu Park, Ingmar Weber, Mor Naaman, and Sarah Vieweg.  Understanding musical diversity via online social media.  In 9th International AAAI Conference on Weblogs and Social Media (ICWSM), 2015.
  * [19] Andriy Shepitsen, Jonathan Gemmell, Bamshad Mobasher, and Robin Burke.  Personalized recommendation in social tagging systems using hierarchical clustering.  In Proceedings of the 2008 ACM conference on Recommender systems, pages 259–266. ACM, 2008.
  * [20] Douglas R Turnbull, Justin A Zupnick, Kristofer B Stensland, Andrew R Horwitz, Alexander J Wolf, Alexander E Spirgel, Stephen P Meyerhofer, and Thorsten Joachims.  Using personalized radio to enhance local music discovery.  In CHI’14, 2014.
  * [21] Aaron Van den Oord, Sander Dieleman, and Benjamin Schrauwen.  Deep content-based music recommendation.  In Advances in Neural Information Processing Systems, 2013.