# Public Release of A-SLOTH: Ancient Stars and Local Observables by Tracing Halos

Tilman Hartwig, Mattis Magg, Li-Hsin Chen, Yuta Tarumi, Volker Bromm, Simon C.O. Glover, Alexander P. Ji, Ralf S. Klessen, Muhammad A. Latif, Marta Volonteri, et al.

► To cite this version:

Tilman Hartwig, Mattis Magg, Li-Hsin Chen, Yuta Tarumi, Volker Bromm, et al. Public Release of A-SLOTH: Ancient Stars and Local Observables by Tracing Halos. Astrophys.J., 2022, 936 (1), pp.45. 10.3847/1538-4357/ac7150. hal-03697143

HAL Id: hal-03697143

https://hal.science/hal-03697143

Submitted on 11 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Public Release of A-SLOTH: Ancient Stars and Local Observables by Tracing Halos

Tilman Hartwig<sup>1,2,3</sup> <sup> DID </sup>, Mattis Magg<sup>4,5</sup> <sup> DID </sup>, Li-Hsin Chen<sup>4,5</sup> <sup> DID </sup>, Yuta Tarumi<sup>2</sup> <sup> DID </sup>, Volker Bromm<sup>6</sup> <sup> DID </sup>, Simon C. O. Glover<sup>4</sup> <sup> DID </sup>, Alexander P. Ji<sup>7,8</sup> <sup> DID </sup>, Ralf S. Klessen<sup>4,9</sup> <sup> DID </sup>, Muhammad A. Latif<sup>10</sup> <sup> DID </sup>, Marta Volonteri<sup>11</sup> <sup> DID </sup>, and Naoki Yoshida<sup>1,2,3</sup> <sup> DID </sup>

<sup>1</sup> Institute for Physics of Intelligence, School of Science, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan; hartwig@phys.s.u-tokyo.ac.jp

<sup>2</sup>Department of Physics, School of Science, The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan

<sup>3</sup> Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

<sup>4</sup> Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany

<sup>5</sup> International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD), Königstuhl 17, D-69117 Heidelberg, Germany

<sup>6</sup> Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA

<sup>7</sup> Department of Astronomy & Astrophysics, University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637, USA

<sup>8</sup> Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

<sup>9</sup> Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 225, D-69120 Heidelberg, Germany

<sup>10</sup> Physics Department, College of Science, United Arab Emirates University, PO Box 15551, Al-Ain, UAE

<sup>11</sup> Institut d'Astrophysique de Paris, CNRS & Sorbonne Université, UMR 7095, 98 bis bd Arago, F-75014 Paris, France

Received 2022 March 17; revised 2022 May 8; accepted 2022 May 17; published 2022 August 29

Abstract

The semianalytical model A-SLOTH (Ancient Stars and Local Observables by Tracing Halos) is the first public code that connects the formation of the first stars and galaxies to observables. After several successful projects with this model, we publish the source code (https://gitlab.com/thartwig/asloth) and describe the public version in this paper. The model is based on dark matter merger trees that can either be generated based on Extended Press–Schechter theory or be imported from dark matter simulations. On top of these merger trees, A-SLOTH applies analytical recipes for baryonic physics to model the formation of both metal-free and metal-poor stars and the transition between them with unprecedented precision and fidelity. A-SLOTH samples individual stars and includes radiative, chemical, and mechanical feedback. It is calibrated based on six observables, such as the optical depth to Thomson scattering, the stellar mass of the Milky Way and its satellite galaxies, the number of extremely metal-poor stars, and the cosmic star formation rate density at high redshift. A-SLOTH has versatile applications with moderate computational requirements. It can be used to constrain the properties of the first stars and high-z galaxies based on local observables, predicts properties of the oldest and most metal-poor stars in the Milky Way, can serve as a subgrid model for larger cosmological simulations, and predicts next-generation observables of the early universe, such as supernova rates or gravitational wave events.

Unified Astronomy Thesaurus concepts: Population III stars (1285); Population II stars (1284); High-redshift galaxies (734); Astronomical simulations (1857); Milky Way formation (1053); Open source software (1866)

Software reviewed by the Journal of Open Source Software JSS

<figure><img src="image_3.png" /></figure>

1. Introduction

One of the main goals of astrophysics is to understand the formation of stars and galaxies in order to explain the evolution of the universe from an almost homogeneous state shortly after the big bang until the present day. There are many open questions on small and large scales about the formation of stars and evolution of galaxies. One possible approach to answer these questions is to model the temporal evolution of the universe from high redshift to z = 0 (e.g., Peebles 1993).

Structures in the universe formed hierarchically under the influence of gravity. Small overdensities in the early universe evolved into the first gravitationally bound structures, which then hosted the first stars. These small early structures then merged with other halos to form bigger structures, which eventually became the galaxies that we observe in the local universe. In general, the evolution of the universe proceeds from simple to complex, from small to big, from pristine to chemically enriched.

While we have many observations of phenomena in the local universe, our understanding of the cosmos decreases as we move further back in time. For example, we understand star formation in the Milky Way (MW) sufficiently well to model the mass distribution of newly formed stars or predict the final fates of nearby stars. However, we have only a vague understanding of the formation of the first generations of stars owing to a lack of direct observables (Glover 2005; Bromm & Yoshida 2011; Greif 2015; Klessen 2019).

Numerical simulations and semianalytical models (SAMs) can bridge this gap. Starting from ab initio initial conditions, one can model the successive formation of stars and galaxies in time and eventually compare the results of such models with observations of the real universe. However, numerical simulations have to compromise between mass resolution and effective volume or use a zoom-in approach in which only some part of the volume is simulated with higher resolution. For example, TNG50 (Nelson et al. 2019; Pillepich et al. 2019), the highest-resolution cosmological simulation in the IllustrisTNG project, has a mass resolution of $m_{\text{DM}} = 4.5 \times 10^5 M_{\odot}$ and $m_{\text{baryon}} = 8.5 \times 10^4 M_{\odot}$ for dark matter (DM) and baryons, respectively. Therefore, the smallest galaxies reliably identified in TNG50 (Engler et al. 2021) have stellar masses of $\sim 6.5 \times 10^6 M_{\odot}$ or more, far larger than the values we expect to find in the first minihalos at high redshift. On

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
the other hand, numerical simulations with high enough resolution to reliably model individual minihalos (e.g., Schauer et al. 2019) typically involve volumes that are too small to allow them to follow the formation of massive galaxies at low redshift.

SAMs are computationally more efficient (see below). Depending on their specific flavor, they approximate galaxies with a few characteristic parameters and components and follow their evolution through time. This makes SAMs computationally very efficient and allows one to focus on a few specific aspects of the model. For example, one can trace certain types of stars from high redshift to the present-day MW.

Moreover, SAMs allow us to explore a broad parameter space due to their modest runtime. While large cosmological simulations have runtimes of several million CPU hours on modern supercomputers (Pillepich et al. 2019; Yoshikawa et al. 2021), an SAM that covers a similar effective volume can run on an ordinary computer within a few minutes.

This computational efficiency makes SAMs ideally suited to study poorly constrained phenomena in the early universe that have an observable effect in the local universe. For example, the formation of the first stars in the universe (Population III or Pop III) is not well understood. However, it affects the epoch of reionization (Bromm & Yoshida 2011; Dayal & Ferrara 2018), the formation of the first supermassive black holes (SMBHs; Latif & Ferrara 2016; Woods et al. 2019; Inayoshi et al. 2020), and the population of old and metal-poor stars in the MW and its satellites (Frebel & Norris 2015; Ishigaki et al. 2021). An SAM allows us to model the formation of the first stars under different assumptions (e.g., regarding their mass distribution or star formation efficiency (SFE)) and determine which assumption agrees best with observations, and the models that reproduce successfully known observables can be used to predict new diagnostics. Therefore, SAMs provide a fruitful connection between the high-redshift universe and local observables.

In the next section, we first present our previous models that lead to the development of A-SLOTH, then discuss and compare various other SAMs in the literature, and eventually present the benefits of A-SLOTH.

## 1.1. A-SLOTH Ancestors

A-SLOTH (Figure 1) was originally developed to model the formation and distribution of possible low-mass Pop III survivors in the MW (Hartwig et al. 2015) and is based on GALFORM (Parkinson et al. 2008). We then extended A-SLOTH to be applicable for a variety of questions relating to the high-redshift universe, as well as the oldest stars in the MW. Among the topics we examined were whether the $Ly\alpha$ emitter CR7 (Sobral et al. 2015) could be a Pop III galaxy (Hartwig et al. 2016a), surviving metal-free stars in MW satellites (Magg et al. 2018), constraints on the pristine initial mass function (IMF; Tarumi et al. 2020), and the low-mass end of the satellite mass function of the MW (Chen et al. 2022). Each of these works built on the previous models and contributed to the development of A-SLOTH.

We also implemented a weighting of merger trees based on their mass and predicted the rates of supernovae (SNe) from the first stars (Magg et al. 2016). This weighting allows us to sample halo masses from the halo mass function in order to create a cosmologically representative effective volume from merger trees generated using the Extended Press–Schechter (EPS) formalism (see Section 2.1.1 below). We also applied this updated model to predict the gravitational wave (GW) signals from the remnants of the first stars and showed that the

detection rate of Pop III GW events is lower than that of events from other channels, but that Pop III remnants could be identified in the future by their higher masses (Hartwig et al. 2016b). In an independent study, Liu & Bromm (2021) adapted A-SLOTH to study GWs from the remnants of the first stars in nuclear star clusters and found promising detection rates.

Our improved treatment of metallicity allowed us to model the transition from Pop III to Pop II star formation more realistically, and we have implemented a new model for inhomogeneous metal mixing in the first galaxies (Hartwig et al. 2018b; Tarumi et al. 2020; Ishigaki et al. 2021). This allows us to predict the properties of second-generation and extremely metal-poor (EMP) stars.

An updated model of star formation allowed us to sample individual Pop II stars and to simulate the formation of MW satellites (Chen et al. 2022). A-SLOTH can also be used as a subgrid model to larger cosmological simulations, e.g., to predict the impact of the first generations of stars on the global 21 cm signal (Magg et al. 2022a).

## 1.2. Other Semianalytical Models

SAMs are popular in astrophysics owing to their moderate computational costs. They approximate and simplify physical processes to model the universe, and each has been developed for a specific purpose, which limits their applicability. Here we summarize SAMs in the literature, discuss their strengths, and later show how A-SLOTH improves on several aspects of these existing models. We limit this discussion to only SAMs that model Pop III star formation explicitly. There are further notable SAMs that model the pre-reionization universe but that do not focus on the formation of the first generations of stars (Volonteri et al. 2003; Barausse 2012; Starkenburg et al. 2013; Agarwal et al. 2016; Pacucci et al. 2018). A comparison of further SAMs is presented in Knebe et al. (2015).

GAMETE (Salvadori et al. 2007) and its extension GAMETE/QSTDUST (Valiante et al. 2016, 2018) model early star and galaxy formation based on EPS-generated merger trees. The

<figure><img src="image_5.png" /><figcaption>Figure 1. A-SLOTH conquering a merger tree.</figcaption></figure>
model is suited for all redshifts. It has an arbitrarily fine mass resolution and accounts for radiative, chemical, and mechanical feedback. To distinguish between Pop III and Pop II star formation, the authors use a criterion based on a critical metallicity $Z_{crit}$ (Bromm & Loeb 2003), i.e., on the total metal and dust content of a halo. They calibrate their model based on EMP stars, the carbon-enhanced metal-poor (CEMP) star fraction, and the baryonic content of the MW at $z=0$ (Salvadori et al. 2007; de Bennassuti et al. 2017). They are also able to reproduce a sample of quasars at high redshift (Valiante et al. 2018) with a subgrid model for the seeding, growth, and feedback of SMBHs (Valiante et al. 2016, 2018; Sassano et al. 2021; Trinca et al. 2022). The GAMETE code has also been combined with N-body simulations to gain spatial information (Salvadori et al. 2010; Pacucci et al. 2017). Graziani et al. (2017) apply the model GAMESH (Graziani et al. 2015) to simulate the formation of a MW-like halo. They follow the baryonic evolution and radiative transfer in a semi-analytic way. This model reproduces the observed galaxy main sequence, the mass-metallicity and the fundamental plane of metallicity relations in the redshift range $z = 0-4$.

The SAM presented in Visbal et al. (2020) simulates the formation of the first stars based on merger trees from N-body simulations. It has a mass resolution of ~$10^5 M_\odot$ and operates in the redshift range $z = 6-30$. The authors include feedback in the form of ionizing radiation, Lyman-Werner (LW) photons, and metals based on a computationally efficient fast Fourier transform. They also take into account baryonic streaming and reproduce the galaxy UV luminosity function at high redshift.

DELPHI (Dayal et al. 2014, 2017, 2019; Chatterjee et al. 2020; Dayal et al. 2020) is a versatile model used to simulate the formation of galaxies and SMBHs at high redshift. It operates down to $z=4$ and allows for different cosmologies (other than ΛCDM). The authors include external radiative feedback and can reproduce various independent observables, such as the UV luminosity function of galaxies and quasars at $z=6$, star formation rate density (SFRD), the black hole mass function, and the black hole mass-stellar mass relation.

GAMMA is an SAM used to probe the impact of nuclear astrophysics uncertainties on galactic chemical evolution and to address the chemical evolution of dwarf galaxies in the early universe (Côté et al. 2017, 2018). It uses spatially resolved DM merger trees from N-body simulations and operates down to $z \sim 7$. It includes inhomogeneous metal mixing in post-processing to reproduce the metallicity distribution function (MDF) of the Renaissance simulation (O'Shea et al. 2015). GAMMA traces individual elements and even s- and r-process elements.

STARTREE models the formation of the first and second generation of stars to constrain properties of the first stars based on EMP stars in the MW. It is based on EPS trees and follows individual Pop III stars (Komiya et al. 2009, 2010). The model includes LW and ionizing feedback, a stochastic model for the inhomogeneous distribution of metal enrichment inside a protogalaxy, and faint Pop III SNe and allows the user to trace individual elements. The authors also take into account the possible surface pollution of low-mass Pop III stars from accreted metal-enriched interstellar medium (ISM; Komiya et al. 2010, 2016). The model reproduces the observed mass-metallicity relation for dwarf galaxies and various abundance ratio distributions. It operates in the redshift range $z = 0-30$ and can also be used to model the formation of EMP (Komiya 2011) and CEMP stars (Komiya et al. 2020).

The SAM by Tumlinson (2006) models the chemical evolution of MW-like galaxies to constrain the properties of the first and second generation of stars. The original version is based on EPS-generated merger trees, and an updated model uses merger trees from N-body simulations as input (Tumlinson 2010; Gómez et al. 2014). The authors include chemical, radiative, and mechanical feedback. The model is calibrated to reproduce the MDF, ionization history, relative elemental abundances of EMP stars, MW stellar mass, and satellite luminosity function.

Trenti & Stiavelli (2009) developed an SAM based on EPS merger trees to model the formation of the first and second generation of stars. The model includes radiative feedback and stochastic metal enrichment. The authors do not constrain their model based on observables but instead vary various input parameters to study the effect of such parameters on the model.

Finally, the SAM by Crosby et al. (2013) is based on N-body simulations and allows the user to study cosmic variance. It includes radiative, chemical, and mechanical feedback. The authors constrain their model based on the SFRD at $z = 10.3$.

### 1.3. Advantages of A-SLOTH

A-SLOTH is the first public SAM that models the formation of the first stars and galaxies and connects those to local observables. It is calibrated transparently, which guarantees reproducibility and has several advantages over existing SAMs.

A-SLOTH uses both EPS-based merger trees and merger trees extracted from N-body simulations as input. We provide instructions to obtain spatially resolved merger trees, or users can input their own data. Currently, A-SLOTH supports merger trees in the format of CONSISTENT-TREES (Behroozi et al. 2013b), but the input routine can be modified easily. Other SAMs that use only EPS-generated merger trees (e.g., Salvadori et al. 2007; Komiya et al. 2009; Dayal et al. 2014) lack spatial information for the feedback model.

A-SLOTH is highly efficient and parallelized. We have implemented a tree-based lookup for external feedback and minimized the memory requirements as much as practically possible. An MW-like galaxy can be simulated in just over 1 minute on four cores, and a cosmologically representative box with side length $8 Mpc h^{-1}$ can be simulated down to redshift $z = 4.5$ on 32 cores in only 14 minutes. The SAM by Visbal et al. (2020) evaluates feedback with a computationally efficient fast Fourier transform. Other SAMs do not comment on their runtime efficiency. For more details on the numerical implementation of A-SLOTH, see Magg et al. (2022c).

A-SLOTH is the only SAM that samples and traces individual Pop III and Pop II stars. The stochastic sampling of individual stars and their SNe is essential to correctly model feedback in the first galaxies (Chen et al. 2022), and the tracking of individual Pop II stars allows us to directly compare our results to metal-poor stars in the MW. Only the EPS-based models by Tumlinson (2006) and Komiya et al. (2009) sample individual Pop III stars. In the nomenclature of A-SLOTH, we distinguish between Pop III (metal-free) and Pop II (metal-poor) stars. The group of Pop II stars in A-SLOTH also includes Pop I (about Solar metallicity) stars, but the exact distinction between Pop II and Pop I is not relevant for all current applications of A-SLOTH.

A-SLOTH tracks individual elements. In the current version, we track carbon and iron explicitly, which allows us to apply a more sophisticated, observationally motivated criterion
to distinguish between Pop III and Pop II star formation (Chiaki et al. 2018). Some SAMs also track individual elements (Côté et al. 2018; Komiya et al. 2020) or even include dust (Valiante et al. 2018), while others use a simple metallicity cut to discriminate between Pop III and Pop II stars (Tumlinson 2006; de Bennassuti et al. 2017).

A-SLOTH is calibrated based on six observables, and our calibration procedure is robust and transparent. SAMs usually have various free parameters that need to be chosen. Depending on the application, one guesstimates or calibrates these input parameters to reproduce specific data. Other SAMs aim to reproduce the MDF (Tumlinson 2006; Komiya et al. 2010; de Bennassuti et al. 2017), observations of quasars or the UV luminosity function at high redshift (Dayal et al. 2014; Valiante et al. 2018; Visbal et al. 2018), or results from simulations (Côté et al. 2018). A-SLOTH is calibrated to reproduce observables from both the local and high-redshift universe simultaneously, and we show that this calibration also reproduces additional constraints that we did not actively aim for, such as the metal-enriched volume fraction, the recovery time distribution from simulations, or the observed stellar mass-halo mass (SMHM) relation.

A-SLOTH models the recovery time (Jeon et al. 2014) between Pop III and Pop II star formation self-consistently by following the baryonic contents over time. Most models assume a constant recovery time, which defines the time between the last Pop III SN explosion and the first Pop II star formation in a halo. We solve the differential equations that govern the heating and cooling of gas inside minihalos in order to calculate when there is enough gas reaccreted after the first SNe to trigger second-generation star formation.

A-SLOTH includes a subgrid recipe for inhomogeneous metal mixing inside the first galaxies, which is crucial to correctly model the metallicity of second-generation stars (Hartwig & Yoshida 2019; Tarumi et al. 2020). Only the SAM of Komiya et al. (2020) applies a similar model that takes into account inhomogeneous mixing of metals. To mimic this effect of inhomogeneous mixing, Côté et al. (2018) apply a smoothing kernel to their MDF in post-processing in order to match results from simulations.

## 2. Methodology

### 2.1. Input Data

Throughout the code, we assume a flat $\Lambda$CDM universe and use cosmological parameters from Planck Collaboration et al. (2016). Specifically, we use the Hubble constant $H_0 = 68 \text{ km s}^{-1} \text{Mpc}^{-1} = 100 \text{ h km s}^{-1} \text{Mpc}^{-1}$, a matter density of $\Omega_m = 0.31$, and a baryon mass fraction of $\Omega_b = 0.049$.

The primary input to A-SLOTH is a DM halo merger tree describing the merger history of the DM halos whose stellar and gas content we wish to model. These merger trees can be generated using the EPS formalism or extracted from a cosmological simulation, as described in more detail below.

#### 2.1.1. Extended Press-Schechter Merger Trees

The “Press-Schechter” formalism (Press & Schechter 1974) provides an analytical estimate for the number density of DM halos as a function of their mass and redshift. The comoving number density of halos of mass between $M$ and $M+dM$ is given by

$$ \frac{dn}{dM} = \sqrt{\frac{2}{\pi}} \frac{\rho_m}{M} \frac{-d(\ln \sigma(M))}{dM} \nu_c \exp(-\nu_c^2/2), \quad (1) $$

where $\rho_m$ is the matter density and $\sigma(M)$ is the standard deviation of the matter power spectrum smoothed on a mass scale $M$. The parameter $\nu_c = \delta_c(z)/\sigma(M)$ is the critical threshold for collapse, where $\delta_c$ is the overdensity for nonlinear collapse with $\delta_c(z=0) = 1.686$.

The original derivation of Equation (1) by Press & Schechter (1974) suffered from a normalization problem: the mass function they derived accounted for only half the mass in the universe, with Equation (1) following only after they multiplied their original expression by a factor of two, a procedure with little justification in their study. Later, Bond et al. (1991) showed that this normalization problem is a consequence of the cloud-in-cloud problem, i.e., the difficulty of correctly counting clouds in hierarchical models where smaller clouds can be embedded within larger ones. They showed that this problem could be avoided by an approach based on excursion sets of the density field. In this approach, the DM density field is smoothed on a series of successively smaller scales $R$, each of which is associated with a corresponding mass scale $M(R)$. Points in the density field that first exceed the overdensity $\delta_c$ when smoothed on scale $R$ are then associated with halos of mass $M$. In the simple case where the smoothing function is a sharp “top-hat” filter in wavenumber space, the mass function that one obtains from this approach is identical to Equation (1), i.e., the method avoids the original normalization problem.

Lacey & Cole (1993) subsequently extended this idea by showing that the excursion set approach allows one to derive not only the DM halo mass function but also the halo merger rate, as well as specific merger histories for individual halos. In A-SLOTH, we use the publicly available implementation of this EPS formalism presented by Parkinson et al. (2008), which is based on the code by Cole et al. (2000).

The EPS formalism as used in the original GALFORM code (Cole et al. 2000) systematically underpredicts the mass of the most massive progenitors at higher redshifts. Hence, we use the updated version of the code by Parkinson et al. (2008), which modifies the progenitor mass function with a perturbing function

$$ \frac{dN}{dM_1} \rightarrow \frac{dN}{dM_1} G(\sigma_1/\sigma_2, \delta_2/\sigma_2) \quad (2) $$

to reproduce the halo merger histories of the Millennium simulation (Springel et al. 2005). The best-fitting perturbing function is given by

$$ G(\sigma_1/\sigma_2, \delta_2/\sigma_2) = 0.57 \left(\frac{\sigma_1}{\sigma_2}\right)^{0.38} \left(\frac{\delta_2}{\sigma_2}\right)^{-0.01}. \quad (3) $$

We have chosen this specific implementation of the merger tree because it performs best compared to other EPS implementations (Jiang & van den Bosch 2014).

#### 2.1.2. Merger Trees from N-body Simulations

In addition to statistically generated merger trees, we use merger trees from cosmological simulations. In this study we use merger trees from Ishiyama et al. (2016) and from the Caterpillar project (Griffen et al. 2016, 2018). These N-body DM-only simulations both have sufficient resolution to resolve minihalos. For both simulations, halos were obtained with
ROCKSTAR (Behroozi et al. 2013a) and merger trees were generated with CONSISTENT-TREES (Behroozi et al. 2013b). From the Caterpillar simulations we use merger trees of 30 zoom-in simulations of MW-mass galaxies ($M_{\text{halo}} \sim 10^{12} M_{\odot}$), selected from a $(100 h^{-1} \text{Mpc})^3$ parent box. The high-resolution DM particle mass is $3 \times 10^4 M_{\odot}$, and the simulations are run to $z=0$ with a high-resolution region that extends at least to the virial radius of each halo. This is the largest suite of high-resolution zoom-in simulations of MW mass that can resolve the formation of low-mass minihalos, providing a good estimate of cosmic variance of $z=0$ properties related to the first stars and galaxies. They are ideal for applications related to stellar archeology. The data from Ishiyama et al. (2016) include merger trees for every halo within an $(8 h^{-1} \text{Mpc})^3$ comoving volume down to a redshift of $z=4$ with a mass resolution of $5000 h^{-1} M_{\odot}$. These merger trees allow us to model a large and well-defined cosmological volume for applications such as ionization histories, cosmic star formation histories, or the large-scale enrichment of the intergalactic medium (IGM).

Both the Caterpillar trees and the merger trees from Ishiyama et al. (2016) are not one individual tree but rather contain a larger number of separate merger trees. In order to allow for easier processing, we resort the merger trees into a strictly time-ordered manner and save them as a binary file.

For brevity, we refer later to the cosmologically representative merger trees from the box with side length $8 h^{-1} \text{Mpc}$ as the “8 Mpc box.” For comparison, each set of Caterpillar merger trees is extracted from a volume with an effective side length of $3-6 h^{-1} \text{Mpc}$.

In contrast to the Press–Schechter merger trees, trees from N-body simulations generally contain subhalos. We prevent accretion of hot gas onto subhalos because their masses tend to fluctuate significantly during mergers. Otherwise, subhalos are treated, followed, and modeled just as normal halos.

### 2.1.3. Auxiliary Input Data

SAMs do not simulate star formation from first principles, as hydrodynamic simulations do. Instead, A-SLOTH relies on analytical formulae and precomputed input quantities. In this section, we explain and justify our input data and explain how users can modify these assumptions for their own models. In the following equations, $x = \log_{10}(M_{\text{star}}/M_{\odot})$, where $M_{\text{star}}$ is the zero-age main-sequence mass of an individual star.

1. **Initial mass function:** For Pop III stars, we assume a power-law IMF of the form

$$ \frac{dN}{d \log M_{\text{star}}} \propto M_{\text{star}}^{-\alpha_{\text{III}}} \quad (4) $$

in the mass range $M_{\text{min}} \le M_{\text{star}}/M_{\odot} \le M_{\text{max}}$, where $\alpha_{\text{III}}$ is an adjustable parameter. For Pop II stars, we follow the Kroupa (2001) IMF in the mass range (0.01, 100) $M_{\odot}$ with the functional form

$$ \begin{aligned} \frac{dN}{d \log M_{\text{star}}} &\propto M_{\text{star}}^{0.7} \quad \text{for} \quad \frac{M_{\text{star}}}{M_{\odot}} < 0.08 \\ \frac{dN}{d \log M_{\text{star}}} &\propto M_{\text{star}}^{-0.3} \quad \text{for} \quad 0.08 \le \frac{M_{\text{star}}}{M_{\odot}} < 0.5 \\ \frac{dN}{d \log M_{\text{star}}} &\propto M_{\text{star}}^{-1.3} \quad \text{for} \quad 0.5 \le \frac{M_{\text{star}}}{M_{\odot}}. \end{aligned} \quad (5) $$

2. **Stellar lifetime:** We compute the stellar lifetimes of Pop III stars with the fitting function (Schaerer 2002)

$$ \log_{10}(T_{\text{III}}/\text{yr}) = 9.785 - 3.759x + 1.413x^2 - 0.186x^3 \quad (6) $$

and the stellar lifetimes of Pop II (1/50 $Z_{\odot}$) stars with the fitting function (Stahler & Palla 2005)

$$ \log_{10}(T_{\text{II}}/\text{yr}) = 10 - 3.68x + 1.17x^2 - 0.12x^3. \quad (7) $$

These fitting functions were constructed to express the lifetimes of stars in the mass ranges 5–500 $M_{\odot}$ for Pop III stars and 7–150 $M_{\odot}$ for Pop II stars. However, the values predicted by Equation (6) agree well with the individual stellar lifetimes of Pop III stars calculated by Marigo et al. (2001) down to subsolar masses. We therefore use this fitting function down to Pop III stellar masses of 0.7 $M_{\odot}$. For lower stellar masses, the lifetimes exceed the current age of the universe, and so the precise values are unimportant. For Pop II stars, our adopted fitting function agrees well with data from Ekström et al. (2008) and Stahler & Palla (2005) for stellar masses ranging from subsolar values up to 200 $M_{\odot}$. Therefore, we also use the fit of stellar lifetimes in this extended range.

3. **Ionizing photon production:** The time-averaged ionizing photon production rates of Pop III and Pop II stars are computed with the formulae$^{12}$ (Schaerer 2002)

$$ \log_{10}(Q_{\text{III}}/\text{s}^{-1}) = 43.61 + 4.9x - 0.83x^2 \quad (8) $$

and

$$ \log_{10}(Q_{\text{II}}/\text{s}^{-1}) = 27.8 + 30.68x - 14.8x^2 + 2.5x^3, \quad (9) $$

respectively. These fitting functions are valid in the mass ranges 9–500 $M_{\odot}$ for Pop III and 7–150 $M_{\odot}$ for Pop II stars. We consider all photons with energies above $E_{\gamma} = 13.6$ eV to be ionizing and neglect the ionization of helium.

4. **Metal yields ejected by SNe:** For the metal yields ejected by SNe, we use tabulated values from Nomoto et al. (2013) and Kobayashi et al. (2006) for Pop III and Pop II stars, respectively. In its current version, A-SLOTH follows carbon, iron, and the sum of all elements explicitly. The code can be modified easily to track more or less elements.

The implemented SN energy, metal yields, ionizing flux, and stellar lifetimes are illustrated as a function of stellar mass in Figure 2. We illustrate the quantities in the stellar mass ranges for which they are valid. However, please note that we sample Pop II stars only up to 100 $M_{\odot}$ and Pop III stars only in the mass range $M_{\text{min}} = 5 M_{\odot}$ to $M_{\text{max}} = 210 M_{\odot}$ in our fiducial model. Users of A-SLOTH can adapt these mass ranges as desired.

We interpolate values onto our grid of stellar masses for lifetime and ionizing flux. For the metal yields, we use the nearest-neighbor values, i.e., the nearest stellar mass available in the tabulated grid of SN yields, so that we do not create artificially mixed combinations by interpolating SN yields from

<footnote n="12">12 These formulae are for single stars. We do not currently account for the impact of binaries on the ionizing photon production rate. However, we expect that their impact will be relatively small. At early times, when photoionization is the most important form of feedback, the ionizing photon production rate is insensitive to the binary fraction, while at later times, when the sensitivity is much larger (Secunda et al. 2020), SNe will typically be the most important form of feedback.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 2. Illustration of SN energy, ionizing photon production rate, carbon and iron ejecta mass, and stellar lifetimes that are available as input data in A-SLOTH. Pop III stars are plotted in red, and Pop II stars are plotted in blue. Users can define their own specific mass ranges or use other tabulated input data.</figcaption></figure>

different progenitor masses. We assume that Pop II and Pop III stars in the mass range 10–40 $M_\odot$ explode as core-collapse SNe (CCSNe), that Pop III stars in the mass range 140–260 $M_\odot$ explode as pair-instability SNe (PISNe; Heger & Woosley 2002, 2010), and that there are no SN explosions for stellar masses between these two ranges.

## 2.2. Baryon Cycle

Each halo is initially assigned a baryonic mass fraction equivalent to the mean cosmic baryon fraction, i.e., the initial baryonic mass of a halo with virial mass $M_{vir}$ is assumed to be $(\Omega_b/\Omega_m)M_{vir}$ prior to any star formation. While baryonic streaming might lower this baryon mass in low-mass halos at high redshift, this effect is not important for star-forming halos in our fiducial model (Naoz et al. 2012). To model the evolution of this baryonic mass, we adopt a model in which the baryonic content associated with a halo is divided into four components: cold gas ($M_{cold}$), hot gas ($M_{hot}$), stars ($M_*$), and outflows ($M_{out}$). We summarize the nomenclature in Table 1. During the evolution, we add hot gas to a halo until

$$M_{\text{vir}} \frac{\Omega_b}{\Omega_m} = M_{\text{cold}} + M_{\text{hot}} + M_* + M_{\text{out}}, \quad (10)$$

which implies that halos accrete additional baryons whenever their DM mass increases by smooth accretion. Smooth accretion here refers to increases in DM mass and the associated increase in baryonic mass, which is not caused by mergers or caused by mergers with unresolved halos. During mergers, all baryonic mass components from the merging halos are added. If DM is accreted via smooth accretion, the baryonic mass budget of the halo is increased accordingly. The additional mass is gradually added to the hot component of the halo. Note that in our model hot gas is any gas that is not cold, i.e., it corresponds to the sum of the warm and hot phases of the ISM in the usual three-phase description (McKee & Ostriker 1977; Klessen & Glover 2016). The outflow mass is

the cumulative mass that is unbound from the halo and is never reaccreted.

In addition, we allow star formation to have multiple epochs instead of one starburst in each global time step $\Delta t_z$. To accurately follow star formation within individual halos, we introduce subcycling. This is necessary because the global time step of most available merger trees is larger than the typical lifetime of massive stars, and so multiple epochs of stellar birth and feedback are expected within one global model time step. At each global time step, adopted from the time evolution of the merger trees, the halo inherits the baryonic properties from all of its progenitors, which are $M_{\text{cold}}^0, M_{\text{hot}}^0, M_{*\text{,II}}^0$, and $M_{\text{out}}^0$. Note that halos at the very first global time step are initialized with hot gas only owing to the accretion shock. Then, in each sub-time step $i$, the baryonic content is updated using the following equations:

$$ 
\begin{aligned}
M_{\text{cold}}^{i+1} &= M_{\text{cold}}^i + \frac{\delta t_i M_{\text{hot}}^i}{t_{\text{dyn}}} - \delta M_{\text{out,cold}}^i - \delta M_{\text{heat}}^i - \delta M_*^i \\
M_{\text{hot}}^{i+1} &= M_{\text{hot}}^i - \frac{\delta t_i M_{\text{hot}}^i}{t_{\text{dyn}}} - \delta M_{\text{out,hot}}^i + \delta M_{\text{heat}}^i + \delta M_{\text{acc,hot}}^i \\
M_{\text{out}}^{i+1} &= M_{\text{out}}^i + \delta M_{\text{out,cold}}^i + \delta M_{\text{out,hot}}^i,
\end{aligned}
\quad (11) $$

where we assume that the timescale of the cooling of hot gas is comparable to the dynamical timescale of the disk, which is composed of stars and cold gas in a central region with a radius of $R_s$. Therefore, $t_{dyn} = R_s/v_{dyn}$, $R_s = R_{vir}/c_{dm}$ is the scale radius of a halo (Navarro et al. 1996), and the velocity

$$v_{\text{dyn}} = \sqrt{G(M_* + M_{\text{cold}})/R_s} \quad (12)$$

is the characteristic velocity of the central part of the halo. For the halo's concentration, $c_{dm}$, we use the fitting functions provided by Correa et al. (2015), which are provided in Appendix A.1. We describe how we determine each baryonic component in the following sections.
<table><caption>**Table 1**
Definitions of the Variables Used in A-SLOTH</caption><thead><tr><th></th><th>Halo-specific Masses</th></tr></thead><tbody><tr><td>M<sub>vir</sub></td><td>Virial mass of the halo</td></tr><tr><td>M<sub>vir,peak</sub></td><td>Peak virial mass of the halo</td></tr><tr><td>M<sub>cold</sub></td><td>Cold gas mass</td></tr><tr><td>M<sub>hot</sub></td><td>Hot gas mass</td></tr><tr><td>M<sub>out</sub></td><td>Cumulative outflow mass</td></tr><tr><td>M<sub>*</sub></td><td>Total stellar mass</td></tr><tr><td>M<sub>disk</sub></td><td>Disk mass, including gas and stars</td></tr><tr><td>δM<sub>out</sub></td><td>Outflow mass</td></tr><tr><td>δM<sub>out,cold</sub></td><td>Cold gas mass that enters outflow</td></tr><tr><td>δM<sub>out,hot</sub></td><td>Hot gas mass that enters outflow</td></tr><tr><td>δM<sub>heat</sub></td><td>Mass that transfers from cold to hot gas</td></tr><tr><td>δM<sub>*</sub></td><td>Stellar mass that forms per substep</td></tr><tr><td>δM<sub>acc,hot</sub></td><td>Hot gas mass that is accreted from IGM</td></tr><tr><td>δM<sub>metals</sub></td><td>Ejected metals by SNe per substep</td></tr><tr><td></td><td>Parameters of the Star Formation Model</td></tr><tr><td>M<sub>min</sub>, M<sub>max</sub></td><td>Limiting masses of the IMF</td></tr><tr><td>α</td><td>Slope of the IMF</td></tr><tr><td>η</td><td>Star formation efficiency</td></tr><tr><td>f<sub>esc</sub></td><td>Ionizing escape fraction</td></tr><tr><td>n<sub>den</sub></td><td>Number density of dense gas</td></tr><tr><td>v<sub>BC</sub></td><td>Baryonic streaming velocity</td></tr><tr><td>γ<sub>out</sub></td><td>Outflow efficiency (Equation (28))</td></tr><tr><td>α<sub>out</sub></td><td>Exponent in γ<sub>out</sub></td></tr><tr><td>M<sub>out,norm</sub></td><td>Normalization mass in γ<sub>out</sub></td></tr><tr><td></td><td>Times and Timescales</td></tr><tr><td>Δt<sub>z</sub></td><td>Time between redshift steps</td></tr><tr><td>δt<sub>i</sub></td><td>Adaptive time step i in the SF subcycle</td></tr><tr><td>t<sub>dyn</sub></td><td>Dynamical timescale of the halo</td></tr><tr><td>t<sub>cold,ff</sub></td><td>Freefall timescale of cold gas</td></tr><tr><td>t<sub>star</sub></td><td>Stellar formation timescale</td></tr><tr><td>t<sub>cool</sub></td><td>Cooling timescale</td></tr><tr><td>t<sub>acc</sub></td><td>Hot gas accretion timescale</td></tr></tbody></table>

Note. The subscripts II and III in the main text refer to Pop II-specific and Pop III-specific quantities, respectively.

## 2.3. Star Formation and Feedback

### 2.3.1. Critical Halo Mass for Cooling and Star Formation

To form stars, gas needs to cool down efficiently. Since we cannot follow the chemical evolution of the gas without information on the spatial distribution of the baryons inside a halo, we instead use a simple mass criterion: rapid gas cooling and star formation are assumed to only occur in halos that have masses greater than a critical mass $M_{crit}$. The appropriate value to take for $M_{crit}$ is somewhat uncertain, with different studies in the literature yielding different results (Yoshida et al. 2003; Latif & Khochfar 2019; Kulkarni et al. 2021; Schauer et al. 2021). A-SLOTH therefore offers three different models for this critical mass: one based on Schauer et al. (2021) (hereafter S21; our fiducial model); a second one based on O’Shea & Norman (2008), Hummel et al. (2012), and Stacy et al. (2011) (OHS); and a model based on Fialkov et al. (2013) (F13). We describe the details of our fiducial model here. The other two models are described in Appendix A.2.

In our fiducial model, we follow Equations (9) and (10) from Schauer et al. (2021), which quantify the dependence of $M_{crit}$ on the strength of the LW background and the velocity associated

with the large-scale streaming of the baryons relative to the DM, $v_{BC}$. Combining these equations yields the following expression for $M_{crit}$:

$$ \log_{10}(M_{\text{crit,S21}} / M_{\odot}) = 6.0174 (1.0 + 0.166\sqrt{J_{21}}) + 0.4159 v_{\text{BC}}. \quad (13) $$

Here $J_{21}$ is the strength of the LW background in units of $10^{-21}$ erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$ sr$^{-1}$, and $v_{BC}$ is the streaming velocity in units of the rms value, $\sigma_{rms}$. We adopt $v_{BC} = 0.8$, as this is the most likely strength of the streaming velocity (Schauer et al. 2021), and we assume a redshift-dependent form for the LW background following Greif & Bromm (2006),

$$ J_{21} = 10^{2-z/5}. \quad (14) $$

We assume a global LW background because the effective volume that we simulate with A-SLOTH is smaller than the mean free path of LW photons (Ahn et al. 2009). In this public release, we do not include near-field LW feedback, which could self-regulate Pop III star formation. However, users can easily implement this feature, for example, following the model in Hartwig et al. (2016b).

We note that the expression for $M_{crit}$ we use here corresponds to the mass scale at which an average halo of that mass can cool and form stars. The minimum halo mass for cooling, $M_{crit,min}$, is typically a factor of a few smaller (Schauer et al. 2021), but as cooling and star formation only occur in a few halos in the mass range $M_{crit,min} < M < M_{crit}$, we consider that $M_{crit}$ is the more appropriate value for our purposes.

If a halo is above the atomic-cooling limit ($T_{vir} \ge 10^4$ K), the gas can cool efficiently as long as the ionizing background is not strong (Visbal et al. 2017; see also Section 2.5.2). We account for this via another critical mass scale (Hummel et al. 2012)

$$ M_{\text{ACH}} / M_{\odot} = 10^{7.5} \left( \frac{1+z}{10} \right)^{-1.5}. \quad (15) $$

The final critical mass scale is then simply given by

$$ M_{\text{crit,final}} = \min(M_{\text{crit}}, M_{\text{ACH}}). \quad (16) $$

### 2.3.2. Stochastic Sampling of Stars

Earlier versions of A-SLOTH considered star formation only in clusters with IMF-averaged and continuous feedback (Hartwig et al. 2015; Magg et al. 2018; Tarumi et al. 2020). We have improved the star formation model by adding a new implementation such that we keep track of individual stars (Chen et al. 2022). We store the numbers, masses, and birth times of stars that formed in each time step. This novel feature of A-SLOTH allows us to account for stellar feedback at specific times and to correctly model the stochastic feedback. In addition, we are also able to keep track of those stars that survive until $z=0$ and their stellar properties.

First, we determine whether we form Pop III stars or Pop II stars based on the empirically derived critical metallicity of Chiaki et al. (2017). We form Pop II stars if

$$ 10^{[\text{C/H}]-2.30} + 10^{[\text{Fe/H}]} > 10^{-5.07}, \quad (17) $$

where [Fe/H] and [C/H] are the iron and carbon abundances of the star-forming gas. Then, we estimate the newly formed
stellar mass given an SFE $\eta_*$:

$$M_{*,est}^i = \eta_* M_{\text{cold}}^i \frac{\delta t_i}{t_{\text{cold,ff}}^i}, \qquad (18)$$

where $\eta_*$ is the SFE per freefall time, $t_{\text{cold,ff}}^i = (G\rho_{\text{cold}}^i)^{-1/2}$ is the freefall time of cold gas, $\rho_{\text{cold}}^i = M_{\text{cold}}^i/V_{\text{cold}}$ is the mean cold gas density of the halo, and $M_{\text{cold}}^i$ and $V_{\text{cold}}$ are the mass and volume of cold gas, respectively.

Next, we compute the average number of stars in the logarithmically spaced IMF bins. For Pop II stars, we assume a Kroupa IMF, and for Pop III stars, we assume an IMF with slope $\alpha_{\text{III}} = 1$. If $n_{j,\text{avg}} \le 10$, we determine the number of stars in bin j with Poisson sampling. We deactivate Poisson sampling when $n_{j,\text{avg}} > 10$ to save computational time.

### 2.3.3. Photoheating

Massive stars with $M_{\text{star}} \ge 8 M_{\odot}$ dominate the production of energetic photons in a stellar population. These photons can ionize hydrogen, leading to the heating of gas. We compute a mass conversion rate for stars at different masses by estimating the mass that is enclosed in the H II region of a massive star. Our calculation assumes that the expansion of the H II region is governed by the Spitzer (1978) solution

$$R_{\text{I}}(t) = R_{\text{D}} \left[ 1 + \frac{7}{4} \frac{c_s (t - t_{\text{D}})}{R_{\text{D}}} \right]^{4/7}, \qquad (19)$$

where $R_{\text{I}}(t)$ is the radius of the region at time t, $R_{\text{D}}$ is the radius at which the ionization front first undergoes the transition from R-type to D-type expansion, $t_{\text{D}}$ is the time at which the front reaches $R_{\text{D}}$, and $c_s$ is the sound speed of the ionized gas. The latter is given by $c_s = 11.4(T_{\text{ion}}/10^4 \text{ K})^{1/2} \text{ km s}^{-1}$, and we adopt a value $T_{\text{ion}} = 10^4 \text{ K}$ for the temperature of the ionized gas.<fnref n="13" /> Using the Spitzer solution, we derive<fnref n="14" /> the following expression for the instantaneous mass conversion rate:

$$\dot{M}_{\text{heat}} = m_{\text{H}} n_{\text{cold,den}} R_{\text{D}}^2 c_s \left[ 1 + \frac{7}{4} \frac{c_s (t - t_{\text{D}})}{R_{\text{D}}} \right]^{-1/7}, \qquad (20)$$

where $m_{\text{H}}$ is the mass of a hydrogen atom, and we assume $n_{\text{cold,den}} = 10^3 \text{ cm}^{-3}$ for the number density of the neutral gas surrounding a newly formed star. As star formation correlates with gas density (Gao & Solomon 2004), we expect this value to be larger than the mean density of the cold gas. Our choice of value here is motivated by Roman-Duval et al. (2010), but in practice $\dot{M}_{\text{heat}}$ is only weakly sensitive to this choice ($\dot{M}_{\text{heat}} \propto n_{\text{cold,den}}^{\beta}$ with $-3/7 \le \beta \le -1/3$). For other physical quantities related to cold gas mass such as the gas binding energy and the cold gas freefall time, we use the mean cold gas density, which is based on the current cold gas mass. We find that the mass conversion rate stays almost constant throughout the stellar lifetime (see Figure B1 in Chen et al. 2022). To save

computational time, we therefore use a time-averaged mass conversion rate, $\langle \dot{M}_{\text{heat}} \rangle$.

The radius of the ionization front at the time that it undergoes the transition from R type to D type is simply the initial value for the Strömgren radius:

$$R_{\text{D}} = \left( \frac{3Q}{4\pi n_{\text{cold,den}}^2 \alpha_{\text{B}}} \right)^{1/3}, \qquad (21)$$

where $\alpha_{\text{B}}$ is the case B recombination coefficient, which for our adopted temperature of $10^4 \text{ K}$ is given by $\alpha_{\text{B}} = 2.6 \times 10^{-13} \text{ cm}^3 \text{s}^{-1}$ (Ferland et al. 1992). The time required for the ionization front to reach $R_{\text{D}}$ is given by

$$t_{\text{D}} = \ln \left[ \frac{842}{23} \left( \frac{Q}{10^{48} \text{s}^{-1}} \right)^{1/3} \left( \frac{n_{\text{cold,den}}}{10^3 \text{cm}^{-3}} \right)^{1/3} \right] \tau_{\text{recomb}}, \qquad (22)$$

where $\tau_{\text{recomb}} = (\alpha_{\text{B}}n_{\text{cold,den}})^{-1}$ is the recombination time.

In dense molecular gas, stars form in clusters. Star clusters are less efficient in heating up the surrounding regions than isolated stars because their ionizing regions overlap. Therefore, we assume that 90% of massive stars in a halo form in one big cluster at the galactic center, whereas the other 10% of massive stars form in isolation (Chen et al. 2022). We sum up the ionizing photons from 90% of massive stars and compute one mass conversion rate $\dot{M}_{\text{heat,cl}}$ for the star cluster and compute individual $\langle \dot{M}_{\text{heat,iso}} \rangle$ for isolated stars. Therefore, in time step i, the total mass conversion rate from massive stars is given by

$$\delta M_{\text{heat,i}} = \left( \dot{M}_{\text{heat,cl}} + \sum_{j=1}^{0.1N} \langle \dot{M}_{\text{heat,iso}}^j \rangle \right) \delta t_i, \qquad (23)$$

where N is the total number of massive stars. For the full derivation and the two extreme cases, see the Appendix of Chen et al. (2022).

### 2.3.4. Outflows Driven by Supernovae

We compute the gas mass ejected from the halo by SNe by comparing the binding energy of gas with the total energy transferred to the gas by the SNe. DM is the dominant source of the gravitational potential, which determines the gas binding energy. The full 3D DM distribution within the halos is not provided in the merger trees, and we therefore assume that all halos follow the Navarro–Frenk–White profile (Navarro et al. 1996). We compute the gas binding energy by considering contributions from DM, cold and hot gas, and stellar components. The gravitational potential of a halo is dominated by DM; however, we find that the contribution of baryonic content is nonnegligible. We assume that the gas surrounding massive stars that die as SNe has been heated up during their stellar lifetimes. Therefore, SNe in A-SLOTH eject the hot gas first. The binding energy of hot gas is computed with the formula

$$E_{\text{bind,hot}}^i = \frac{3 \text{GM}_{\text{vir,peak}} M_{\text{hot},i}}{R_{\text{vir}} \left[ \frac{-R_{\text{vir}}}{R_s + R_{\text{vir}}} + \ln \frac{R_s + R_{\text{vir}}}{R_s} \right]} \times \left[ -\frac{1}{4} + \frac{1}{2} \left( 1 - \frac{R_s^2}{R_{\text{vir}}^2} \right) \ln \frac{R_s + R_{\text{vir}}}{R_s} + \frac{1}{2} \frac{R_s}{R_{\text{vir}}} \right] + \left( \frac{3R_s}{2R_{\text{vir}}} - \frac{13R_s^3}{30R_{\text{vir}}^3} \right) \frac{\text{GM}_{\text{disk},i} M_{\text{hot},i}}{R_s} + \frac{3}{5} \frac{\text{GM}_{\text{hot},i}^2}{R_{\text{vir}}}, \qquad (24)$$

<footnote n="13">¹³ In reality, we expect gas in H II regions to have a range of different temperatures, depending on the metallicity of the gas and of the massive stars, the age and evolutionary state of the region, etc. However, in general these temperatures will differ from $10^4 \text{ K}$ by no more than a few tens of percent (see, e.g., Ho et al. 2019), and so the error we make by adopting the same fixed value for all H II regions will be small.</footnote>

<footnote n="14">¹⁴ Details of this derivation can be found in Appendix B of Chen et al. (2022).</footnote>
where $M_{vir,peak}$ is the peak virial halo mass, $R_s$ is the scale radius of the DM halo (Navarro et al. 1996), $R_{vir}$ is the virial radius of the halo, and $M_{disk}^i$ is the combined mass of stars and cold gas since we assume that they reside in the same central region in the halo (Chen et al. 2022). We use the peak virial halo mass until the current redshift. This definition is more robust against the fluctuating results of the halo finder.

Similarly, the binding energy of cold gas can be described by

$$E_{\text{bind,cold}}^i = \frac{3 GM_{\text{vir,peak}} M_{\text{cold},i}}{4R_s \left[ \frac{-R_{\text{vir}}}{R_s + R_{\text{vir}}} + \ln \frac{R_s + R_{\text{vir}}}{R_s} \right]} + \frac{6 GM_{\text{cold},i} M_*}{5R_s} + \frac{GM_{\text{hot},i} M_{\text{cold},i}}{R_{\text{vir}}} \left( \frac{3}{2} - \frac{3R_s^2}{10R_{\text{vir}}^2} \right) + \frac{3 GM_{\text{cold},i}^2}{5R_s} \quad (25)$$

Complete derivations of both of these binding energies can be found in the Appendix of Chen et al. (2022).

We distribute the sum of SN energies in this step, $E_{\text{SNe}}^i$, to hot gas and cold gas based on their binding energies and masses. The fractions of $E_{\text{SNe}}^i$ that affect hot gas and cold gas are

$$\begin{align}
f_{\text{hot}} &= \frac{E_{\text{bind,hot}}^i M_{\text{hot}}^i}{E_{\text{bind,hot}}^i M_{\text{hot}}^i + E_{\text{bind,cold}}^i M_{\text{cold}}^i} \nonumber \\
f_{\text{cold}} &= \frac{E_{\text{bind,cold}}^i M_{\text{cold}}^i}{E_{\text{bind,hot}}^i M_{\text{hot}}^i + E_{\text{bind,cold}}^i M_{\text{cold}}^i} \quad (26)
\end{align}$$

The amount of hot gas ejected by SNe is calculated using the expression

$$\delta M_{\text{out,hot}}^i = \min \left( \frac{E_{\text{SNe}}^i f_{\text{hot}} / \gamma_{\text{out}}}{E_{\text{bind,hot}}} M_{\text{hot}}^i, M_{\text{hot}}^i \right) \quad (27)$$

Here $\gamma_{\text{out}}$ is the outflow efficiency introduced by Chen et al. (2022) and defined as

$$\gamma_{\text{out}} = \left( \frac{M_{\text{vir,peak}}}{M_{\text{out,norm}}} \right)^{\alpha_{\text{out}}}, \quad (28)$$

where the normalization mass $M_{\text{out,norm}}$ and $\alpha_{\text{out}}$ are free parameters whose values are calibrated in Section 3.

Similarly, the amount of cold gas ejected by SNe is determined with

$$\delta M_{\text{out,cold}}^i = \min \left( \frac{E_{\text{SNe}}^i f_{\text{cold}} / \gamma_{\text{out}}}{E_{\text{bind,cold}}^i} M_{\text{cold}}^i, M_{\text{cold}}^i \right) \quad (29)$$

## 2.4. Adaptive Time Step

To evolve the baryonic physics in time, we numerically integrate the corresponding differential equations with an adaptive time step. If chosen correctly, such an adaptive time step saves computational time, while still allowing for accurate integration of dynamical effects that require a fine time resolution. The adaptive time step is calculated for each halo at every time. The general idea is that we consider the main characteristic timescales of the system and then choose a time step that is smaller than all these times. We consider the star

formation timescale

$$t_{\text{SF}} = \frac{M_* + 10^2 M_{\odot}}{\dot{M}_*}, \quad (30)$$

where $M_*$ is the number of stars already present in this halo, $\dot{M}_*$ is the current SFR, and we add a small offset mass to prevent this time step from getting too small, i.e., to gain computational speed. These small offsets here and below are chosen empirically as the largest values that still guarantee a stable time integration under all conditions. We expect roughly one SN per 100 $M_{\odot}$. So this offset corresponds to the smallest unit of stellar mass whose formation we want to resolve in time. A smaller time step would mean that we follow the formation of SN progenitor stars with unnecessarily high time resolution.

The accretion timescale is given by

$$t_{\text{acc}} = \frac{M_{\text{hot}} + 10^4 M_{\odot}}{\dot{M}_{\text{acc}}}, \quad (31)$$

where $M_{\text{hot}}$ is the amount of hot gas already in the halo, $\dot{M}_{\text{acc}}$ is the accretion rate of new hot gas from the IGM, and we add a small offset mass to prevent this time step from getting too small. Here the offset mass is chosen as ~1% of the baryonic mass of a typical minihalo. This means that we allow A-SLOTH to approximate the accretion of hot gas at the percent level in order to gain computational efficiency.

We also use the dynamical time as defined in Equation (12). We found that the time integration of cold gas can be improved when we use $0.1 t_{\text{dyn}}$ as a constraint for the dynamical time step.

The merger tree is defined based on redshift steps, which is the timescale on which halos merge. To account for this, we also include the time between redshift steps, $\Delta t_z$, in our calculation of the adaptive time step.

The adaptive time step is then given by

$$t_{\text{adapt}} = 0.025 \min(t_{\text{SF}}, t_{\text{acc}}, 0.1 t_{\text{dyn}}, \Delta t_z) + 1 \text{ yr.} \quad (32)$$

The addition of 1 yr prevents extremely small time steps in some pathological cases, and the prefactor of 2.5% is the largest value that still guarantees a stable and converged evolution of the baryon budget.

The smallest of these timescales, which then defines the adaptive time step, depends on redshift, halo mass, and other factors. A-SLOTH provides routines to visualize the relative importance of these characteristic timescales in a specific setup. Figure 3 shows how the different gas phases of our model evolve over time with an adaptive time step.

This figure is an example for the main branch of a merger tree leading to the formation of an MW-like galaxy and is therefore not representative of the behavior of all minihalos. It illustrates how different baryonic components are converted into each other in a minihalo in the redshift range $z = 27-18$.

## 2.5. Spatial Feedback

One of the primary advantages of using merger trees from N-body simulations is that the positions of halos are known and subhalos are resolved. This opens up the possibility of modeling the spatially resolved reionization and enrichment history of the IGM. The method for this was introduced in Magg et al. (2018), and the algorithm was refined in Magg et al. (2022a). For completeness, we describe the implementation in this section. We
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Illustration of the baryonic subcycling and what physical quantities are available within A-SLOTH. Top: various mass contents along one merger tree branch as a function of time after the big bang. Bottom: gas metallicity of the IGM (blue) and ISM (black). The illustrated time corresponds to the redshift range $z = 27 - 18$. We witness the formation of the first Pop III stars around 110 Myr (orange), the suppression of further star formation by SNe around 120 Myr (olive), the accumulation of further cold gas around 130 Myr (blue), and the formation of second-generation stars after ~140 Myr (cyan). The ISM metallicity experiences strong variations during the first SN explosions. This is due to the small amount of residual gas and metals that remain in the halo. The new stellar mass $\delta M_*$ (purple) and the ejected metal mass $\delta M_{\text{metals}}$ (olive) are the instantaneous values at the specific time step during subcycling. The halo mass $M_{\text{halo}}$ (green) increases in steps, which illustrates the original redshift steps from the merger tree. The stellar masses (orange, cyan) and outflow mass (pink) are cumulative and therefore monotonically increasing.</figcaption></figure>

provide further technical details on the implementation of spatial feedback in Appendix A.3.

### 2.5.1. External Enrichment and IGM Metallicity

Each star-forming halo $i$ is assigned an enriched region of radius $R_{\text{en},i}$. Outflows are launched at the virial radius, i.e., $R_{\text{en},i}$ is initialized to be the virial radius. Here outflows refer to baryons that are pushed outside the virial radius of a halo.

Following Agarwal et al. (2012), we assume an outflow velocity of $v_{\text{out}} = 110 \text{ km s}^{-1}$. We chose this fixed value because, due to the calibrated efficiency factor in Equation (28), we cannot self-consistently compute the outflow velocity from momentum conservation or energetic consideration. This efficiency factor was required to be able to match the SMHM relation. The mass of these outflows is modeled according to Equations (27) and (29). The expansion of the outflow is modeled as a pressure-driven snowplow. In this approximation the enriched region expands with constant momentum and slows down as it sweeps up the IGM. The expansion velocity is

$$v_{\text{en}} = v_{\text{out}} \frac{M_{\text{out}}}{M_{\text{out}} + \frac{4}{3} \pi \rho_{b} (R_{\text{en}}^{3} - R_{\text{vir}}^{3})}, \quad (33)$$

where $\rho_b$ is the mean baryon density of the IGM, i.e.,

$$\rho_b = \Omega_b \frac{3H_0^2}{8\pi G} (1+z)^3. \quad (34)$$

The enriched regions also increase in volume owing to cosmic expansion. When two halos with enriched bubbles around them merge, we add their volumes and their matter and metal contents.

A halo $k$ is subject to external enrichment if its center is within the enriched bubble of any halo. We compute the mean IGM metallicity at the position $\vec{x}$ of a halo as

$$Z_{\text{IGM},k} = \frac{\sum_{j \in \mathcal{E}} \rho_{\text{met,out},j}}{\rho_b + \langle 3\Omega_b \Omega_m^{-1} M_{\text{vir},k} (4\pi R_{\text{en},j}^3)^{-1} \rangle_{j \in \mathcal{E}} + \sum_{j \in \mathcal{E}} \rho_{\text{out},j}}, \quad (35)$$

where $\mathcal{E}$ is the set of enriching halos, i.e., the set of halos $j$ with a distance to $\mathbf{x}$ of less than $R_{\text{en},j}$. For each enriching bubble the mass density is

$$\rho_{\text{out},j} = 3M_{\text{out},j}(4\pi R_{\text{en},j}^3)^{-1}, \quad (36)$$

and the metal density is

$$\rho_{\text{met,out},j} = \rho_{\text{out},j} Z_{\text{out},j}. \quad (37)$$

The term $\langle 3\Omega_b \Omega_m^{-1} M_{\text{vir},k} (4\pi R_{\text{en},j}^3)^{-1} \rangle_{j \in \mathcal{E}}$ in Equation (35) takes into account that the enriching bubbles have to encompass the halo that is being enriched. As we base the accretion of metals on the metallicity of the IGM, it is not limited by the mass of metals contained within the bubbles. Thus, without this term, a very small but high-metallicity bubble could lead to the accretion of more metals than it contains. With this correction term, the accretion shows the correct behavior: in the limit that the accretion is dominated by very massive bubbles and in the limit that the mass budget is dominated by the accreting halo, the accreted metal mass is limited to the available amount of metals. This method of correcting is necessary, as our algorithm currently does not allow us to remove metals from the bubbles.
### 2.5.2. Photoionization of the IGM

In addition to enriching the IGM with metals, massive stars will also photoionize the surrounding IGM. We assign each halo an ionized region with radius $R_{ion}$. Again, this radius is initialized to be the virial radius. For convenience, we phrase the expansion of each ionized region in terms of its volume $V$ and afterward compute the corresponding radius as $R_{ion} = (3V/4\pi)^{1/3}$. The change in ionized volume can be computed by modeling the expansion as an R-type ionization front, i.e., modeling the expansion assuming that the difference between the number of recombinations and the number of produced ionizing photons leads to expansion or shrinking of the ionized region. If we assume that the gas inside the ionized region is fully singly ionized, we find

$$ \dot{N}_{ion} = Vn^2 C \alpha_B + \dot{V}n \quad (38) $$

for a nucleon number density of $n$ and an ionizing photon emission rate $\dot{N}_{ion}$. Helium reionization is neglected in this model. We assume a clumping factor of $C = 3$ (Robertson et al. 2013). To ensure numerical stability, we integrate this equation with an implicit Euler method, and therefore we update the volume of ionized regions during a time step $\Delta t_z$ from $k$ to $k+1$ as

$$ V_{k+1} = \left( V_k + \frac{\dot{N}_{ion} \Delta t}{n} \right) (1 + \Delta t_z n \alpha_B C)^{-1}. \quad (39) $$

Since the temperature of ionized gas in the IGM is around $T \approx 10^4$ K, halos with a virial temperature lower than this cannot gravitationally bind the ionized gas. Thus, halos below the atomic-cooling limit (i.e., $T_{vir} = 10^4$ K) are not allowed to accrete gas from the IGM. We extend this to masses up to a factor of 10 above the atomic-cooling limit (Visbal et al. 2017) for halos that are exposed to large fluxes of ionizing photons, $F_{ion} > 6.7 \times 10^5$ photons s$^{-1}$ cm$^{-2}$. If the center of one bubble lies within another, they are joined in a new bubble, where the center of the new bubble is the center of mass of the previous two bubbles and the volumes are added. This method prevents a cluster of strong ionizing sources from ionizing the same volume multiple times, and instead creates one big joined ionized region. Similar to the enriched volumes, when two halos with ionized bubbles merge, we add the volume and increase it according to cosmic expansion.

We do not use a global, external ionizing background but calculate the reionization self-consistently within our merger trees. Hence, we implicitly assume that the ionizing radiation feedback in our simulations is dominated by the local sources. The mean free path of ionizing photons in a neutral IGM at the redshift of reionization ($z \sim 7$) is only 3 comoving kpc, so this is a good approximation at early times, while the size of most ionized bubbles remains small. It breaks down once connected ionized regions with sizes comparable to the size of our simulation volume start to become common, as at this point there is a substantial likelihood of such a region expanding into our volume from the outside. Hydrodynamical simulations of reionization in large volumes (e.g., Gnedin & Kaurov 2014; Rahmati & Schaye 2018; Kannan et al. 2022) suggest that this happens once the ionized gas volume filling fraction exceeds ~50%, and so it is likely that A-SLOTH slightly underestimates the time taken to transition from this ionized fraction to a fully

ionized state. In order to do a better job of modeling this late stage of reionization, it would be necessary to use merger trees drawn from a simulation of a much larger volume—for instance, Iliev et al. (2014) argue that a box size of at least $100 h^{-1}$ Mpc is needed to obtain a fully converged reionization history. Unfortunately, we are not aware of a publicly available N-body simulation that is larger than the one that we use from Ishiyama et al. (2016) but that also has the mass resolution and redshift coverage that we require for A-SLOTH.

Finally, we note that we do not yet include the contribution from AGNs to the reionization model. This should be a good approximation, as most studies find their contribution to the ionizing budget prior to reionization to be negligible (see, e.g., Trebitsch et al. 2021).

### 2.5.3. Metallicities and Mixing

Stars synthesize metals and contribute to chemical enrichment. We assume progenitor-mass-specific SN yields for Pop II and Pop III stars separately and calculate the mass of each element by multiplying the yields by the number of exploding SN. The synthesized metals can escape from the halo together with the outflow bubbles. We assume that the fraction of metals ejected from the galaxy is

$$ f_{eject} = \frac{\delta M_{out}}{M_{hot} + M_{cold}}. \quad (40) $$

This implementation might underestimate the metallicity of the outflow, and hence of the IGM, because SN-driven outflows are usually more metal-rich than the average gas in the galaxy (Mac Low & Ferrara 1999; Emerick et al. 2018; Pandya et al. 2021). However, our simple approximation is able to reproduce the relative number of EMP stars observed in the MW, as we show below.

It is also possible that the metallicity of star-forming gas in a galaxy could be significantly different from the mean gas-phase metallicity of the galaxy (Xu et al. 2016b; Magg et al. 2022b). Tarumi et al. (2020) investigate this scenario by analyzing 3D hydrodynamic simulations and find that star-forming gas in halos that have not experienced star formation (purely externally enriched halos) has significantly lower metallicity than the mean gas metallicity. Metals that are accreted onto the halo can enrich the gas at the outskirts but do not enrich the star-forming gas close to the center of the halo. To account for this, in A-SLOTH we define a metallicity shift factor $dZ$ as

$$ dZ = Z_{sf} - Z_{ave}, \quad (41) $$

where $Z_{sf}$ is the metallicity of the star-forming gas and $Z_{ave}$ is the average metallicity of the halo. The distribution of $dZ$ is different between internal enrichment and external enrichment (Tarumi et al. 2020). For computational efficiency, we define a metal-enriched halo as externally enriched if the total metal mass from SNe inside this halo is $<10^{-15} M_{\odot}$ (effectively zero). For internal enrichment, inhomogeneity is not significant. The distribution of $dZ$ in this case is a Gaussian function with $(\mu, \sigma) = (-0.03, 0.15)$. In the case of external enrichment, $dZ$ is typically significantly negative. We use a fitting function that we have derived in our previous study (Tarumi et al. 2020) and that is provided in Appendix A.4. We apply this metallicity correction to the metallicity of the newly formed stars.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Illustration of the spatial extent of external feedback at $z=6$. We show the fraction of ionized volume along lines of sight in orange, the average IGM metallicity along lines of sight in blue, and halos with $M_{vir} > 4 \times 10^8 M_{\odot}$ as black circles.</figcaption></figure>

To illustrate the spatial feedback, we show a projection of the 8 Mpc box in Figure 4. This figure visualizes the spatial extent of radiative and chemical feedback in our cosmologically representative 8 Mpc box. While this figure is a qualitative visualization to demonstrate the available information inside A-SLOTH, we quantify the evolution of volume filling fractions over time below.

## 2.6. Stochastic Feedback

In some situations, it is not possible or desirable to follow the detailed spatial feedback. Either one wants to save memory and runtime, or one is using EPS-generated merger trees that do not have the information on the spatial positions of the halos that our spatial feedback implementation requires. For such cases, we have also implemented and tested an option for stochastic feedback that does not rely on spatial information but that still provides reasonable results.

The main idea is that we follow the volume filling fractions of metal enrichment and of ionizing radiation and then determine stochastically whether a halo is currently in an ionized or metal-enriched region. Once a halo is in such a region, we assume that all its children are so as well. One could construct rare pathological cases in which this assumption does not hold, but we have found that this simplifying assumption adequately reproduces the desired behavior.

The bubbles around each halo of the ionized and metal-enriched volume expand as described above. At each check, the probability that a halo is affected by external feedback is equal to the respective volume filling fraction at that time. Therefore, we only check for external feedback once, when the halo first surpasses $M_{crit}$. Otherwise, if we were to check for external feedback every time step and make the time steps arbitrarily short, the effective probability for external feedback would artificially increase because we would check more often. We show below that checking for external enrichment once provides the correct results. Moreover, we do not add metal masses from external metal enrichment in the case of stochastic feedback; we only flag the halo as externally enriched in such a case.

For any application for which metal enrichment and radiative feedback are important, we recommend users of A-SLOTH to use spatially resolved merger trees from N-body simulations

instead of EPS-generated merger trees. We still provide this mode of A-SLOTH because it is a very efficient way to explore trends with input parameters, and it allows new users to immediately run simulations, without the need to download several GB of merger tree data.

## 3. Calibration

In this section, we first introduce the six observables, based on which we calibrate A-SLOTH. Then, we discuss our calibration strategy and how we found an optimal fit for the input parameters. We show how this calibrated model reproduces other constraints that A-SLOTH was not incentivized to fit. Finally, we discuss how the model performs with EPS-generated merger trees.

### 3.1. Six Observables

We use six different observables to calibrate the input parameters of A-SLOTH: the optical depth to Thomson scattering, the stellar mass of the MW, the cumulative distribution of stellar masses of the MW satellite galaxies, the fraction of EMP stars in the MW halo, the ratio of ultra-metal-poor (UMP, [Fe/H] < -4) to EMP ([Fe/H] < -3) stars, and the cosmic SFRD at high redshift. For each of these observables, we calculate a goodness-of-fit parameter, $p_i$, which quantifies how consistent the output of A-SLOTH is with observations. These parameters can have values in the range $0 \le p_i \le 1$, and higher values represent better fits. Some of these values are inspired by $p$-values (see below), for which a value of $p_i > 0.05$ conventionally means that a model is accepted (i.e., cannot be rejected).

Eventually, we multiply these six individual parameters to obtain our final goodness-of-fit measure, $p_{all}$. Maximizing $p_{all}$ during the calibration guarantees that all six observables are reproduced simultaneously in the final model. In Appendix B.1, we present how the six observables depend on the input parameters. In Table 2 we present the best-fit values, which will be discussed in more detail below.

#### 3.1.1. Ionization History

A key observable from the early universe is the Thomson scattering optical depth $\tau_e$. Free electrons can scatter cosmic microwave background (CMB) photons, leading to a slight polarization of the CMB. As only free electrons take part in this process, $\tau_e$ is influenced by the ionization fraction of the IGM as a function of time, i.e., the ionization history of the universe. The $\tau_e$ parameter is computed as

$$ \tau_e(z) = c\sigma_T n_H \int_0^z dz' f_e Q_{\text{ion}}(z')(1+z')^3 \left| \frac{dt}{dz'} \right|, \quad (42) $$

where $z$ is the redshift, $\sigma_T = 0.665 \times 10^{-24} \text{cm}^2$ is the Thomson scattering cross section, $n_H$ is the cosmological mean density of hydrogen nuclei at $z=0$, $Q_{\text{ion}}$ is the volume filling fraction of ionized regions, and $f_e$ is the number of free electrons per hydrogen nucleus in the ionized IGM. The presence of helium makes this number slightly larger than 1, and we assume for simplicity that

$$ f_e = \begin{cases} 1 + Y_p/2X_p & \text{at } z \le 4 \\ 1 + Y_p/4X_p & \text{at } z > 4, \end{cases} \quad (43) $$
<table><caption>**Table 2**
Six Observables Based on Which We Calibrate A-SLOTH, Together with Their Observed Values, Modeled Values, and Their Resulting Goodness-of-fit Parameters</caption><thead><tr><th>Observable<br>(1)</th><th>Observed<br>(2)</th><th>A-SLOTH<br>(3)</th><th>Goodness-of-fit Parameter<br>(4)</th></tr></thead><tbody><tr><td>Optical depth</td><td>0.0544 ± 0.0073</td><td>0.0539</td><td>0.989</td></tr><tr><td>M<sub>*</sub>MW</td><td>(5.43 ± 0.57) × 10<sup>10</sup> M<sub>☉</sub></td><td>(7.3 ± 3.0) × 10<sup>10</sup> M<sub>☉</sub></td><td>0.529</td></tr><tr><td>CSMF</td><td>...</td><td>...</td><td>0.136</td></tr><tr><td>log<sub>10</sub>(EMP/All)</td><td>-4.7 ± 0.5</td><td>-4.4 ± 0.6</td><td>0.697</td></tr><tr><td>log<sub>10</sub>(EMP/UMP)</td><td>(1.4 - 2) ± 0.3</td><td>1.2 ± 0.2</td><td>0.563</td></tr><tr><td>SFRD</td><td>...</td><td>...</td><td>0.479</td></tr><tr><td>...</td><td>...</td><td>...</td><td>p<sub>all</sub> = 0.0134</td></tr></tbody></table>

Note. The statistical uncertainty for A-SLOTH-based values is the standard deviation from 30 independent Caterpillar trees.

where $Y_p$ and $X_p$ are the primordial abundances of He and H, respectively (Robertson et al. 2013). This prescription assumes that $He^+$ is produced at the same time as $H^+$, primarily via radiation from stellar sources, and that the helium is later ionized further to $He^{++}$, at $z=4$, as a consequence of the steadily increasing contribution of high-redshift AGNs to the extragalactic UV background. In practice, $\tau_e$ is only weakly sensitive to the redshift at which this occurs. The time evolution of the volume filling fraction is based on Equation (39). Planck Collaboration (2016) determined<fnref n="15" /> the value of the Thomson scattering optical depth to be $\tau_e = 0.0544 \pm 0.0073$. Notably, this value is much lower than what was reported earlier (e.g., Komatsu et al. 2011), which implies much less star formation at the highest redshifts (Visbal et al. 2015). The instantaneous ionization redshift corresponding to this optical depth is $z=8.3$. We can compute $\tau_e$ from the ionization fraction of the universe that we determine based on the total volumes of the ionized regions, as described in Section 2.5.2.

In Figure 5, we show the ionization history derived from our 8 Mpc box as compared to various observations and a selection of representative models from the literature (Graziani et al. 2015; Robertson et al. 2015; de Bennassuti et al. 2017). The observational constraints are based on various different techniques (McGreer et al. 2015; Mesinger et al. 2015; Planck Collaboration 2016; Greig & Mesinger 2017; Zheng et al. 2017; Bañados et al. 2018; Konno et al. 2018; Mason et al. 2018; Ouchi et al. 2018). The models by Graziani et al. (2015) and de Bennassuti et al. (2017) are derived for a MW-like volume, and are hence not directly comparable to our cosmologically representative box. We calibrate A-SLOTH only based on the optical depth to Thomson scattering, i.e., the integrated reionization history. The ionization history that our model predicts (the blue line in the figure) agrees well with the observational constraints on the volume filling factor of ionized gas at high redshift (the data points in the figure). This is an important test of the consistency of the model.

### 3.1.2. Stellar Mass of MW and Its Satellites

We calibrate A-SLOTH to reproduce the observed MW stellar mass and the number counts of satellites at $z=0$. To do this, we use the merger trees from the Caterpillar project

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Ionization history simulated with A-SLOTH (blue), compared to various models (dashed) and observations (points).</figcaption></figure>

simulations, each of which was designed to follow the buildup of an MW-like galaxy. The mean stellar mass we obtain with our fiducial A-SLOTH model for these MW analogs is $\langle M_{\text{star,MW}} \rangle = (7.3 \pm 3.0) \times 10^{10} M_{\odot}$. For comparison, the observed stellar mass of the real MW is $M_{\text{star,MW}}^{\text{obs}} = (5.43 \pm 0.57) \times 10^{10} M_{\odot}$ (McMillan 2017).

In Figure 6, we show the cumulative stellar mass function (CSMF) of satellites from the 30 Caterpillar trees, computed using our fiducial model. Each gray curve represents the CSMF from one tree, and for comparison the blue curve shows the observed CSMF from McConnachie (2012) and Muñoz et al. (2018). Due to the limit of current surveys, the search for MW satellites is incomplete below $M_* \sim 10^5 M_{\odot}$. We shade the region below this mass in gray and only calibrate our results to the observed CSMF above this observational completeness limit. We note that these compilations of satellite galaxies do not yet include Ant II and Cra II, which both have stellar masses of $\sim 10^6 M_{\odot}$ (Ji et al. 2021), but including these two galaxies would not significantly change our results.

To quantify the difference between our predicted CSMFs and the observed one, we use a two-sample Kolmogorov–Smirnov test. We perform the tests on each CSMF produced by A-SLOTH and the observed CSMF. From each test, we retrieve

<footnote n="15">¹⁵ We adopt the value of $\tau_e$ from Planck Collaboration (2016) for consistency with earlier work with A-SLOTH, but note that the difference between this value and the most recently published value in Planck Collaboration et al. (2020) is negligible.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** CSMFs at $z=0$ extracted from our set of 30 Caterpillar trees, processed using our fiducial model parameters. The CSMFs from our Caterpillar trees are plotted in gray, and the observed CSMF for the MW is plotted in blue. The gray shaded area represents the region in which the observations are known to be incomplete. The observed CSMF in this region is plotted using a dashed line and is not used to calibrate the model.</figcaption></figure>

a p-value, which we use as a goodness-of-fit parameter. We then compute the average of these parameters from the 30 tests.

We note that, by construction, we neglect the effects of tidal stripping. While this effect is important for at least some of the nearby satellites (Simon 2018), there are several reasons why it cannot be accurately modeled in our current framework. First, both the shape of the potential well of the MW and that of the satellites are inaccurate without accounting for the baryonic component and without the presence of the MW disk. Thus, the tidal stripping cannot be accurately modeled in the DM-only simulations that we base our model on (Bullock & Boylan-Kolchin 2017). Second, halo finders often struggle to correctly identify subhalos and their masses. Indeed, in many cases we see a substantial increase in the mass of the subhalo during tidal stripping events, which is caused by the halo finders associating some of the mass of the main halo with the subhalo. Furthermore, the effect of tidal stripping is more significant on satellites that are in proximity to the MW. Thus, we neglect tidal stripping in the current work and leave it to future versions, ideally built on merger trees with disk potentials to improve on this aspect.

### 3.1.3. Extremely Metal-poor Stars

We use two metallicity-related values for the calibration. These are the fraction of stars in the MW that are EMP stars, $f_{EMP}$, and the ratio between the fraction of UMP stars and the fraction of EMP stars, $f_{UMP}/f_{EMP}$. We focus on these two individual values rather than comparing with the full MDF because it is easier to combine the values obtained from several different surveys (see below) and because it gives us a straightforward way of quantifying our agreement with the observations. We assume that EMP stars in A-SLOTH are predominantly found in the halo (Kobayashi et al. 2011). In addition, Arentsen et al. (2020) did not find any EMP stars among 8000 spectroscopically observed bulge stars. However, Sestito et al. (2019) found several EMP stars with disk-like orbits.

For the EMP fraction, we assume that (i) the stellar mass in the MW is $6 \times 10^{10} M_{\odot}$ (McMillan 2017), (ii) the stellar mass

in the halo of the MW is $1 \times 10^9 M_{\odot}$ (Bullock & Johnston 2005), and (iii) 1 out of 800 halo stars is an EMP star (Youakim et al. 2020). For the UMP-to-EMP ratio, we assume three possible values based on the literature. Youakim et al. (2020) find in the Pristine survey that $10^{-2}$ of EMP stars are UMP stars. Chiti et al. (2021) argue that the MDF $dN/dZ$ can be well fitted with a power-law function $Z^{-1.5}$, implying that $f_{UMP}/f_{EMP} = 10^{-1.5}$. Finally, Bonifacio et al. (2021) use Gaia and Sloan Digital Sky Survey photometry in the TOPoS survey to derive a UMP-to-EMP ratio of 8/217 $\simeq 10^{-1.4}$.

To evaluate the goodness of fit, we calculate how significantly the two metrics ($f_{EMP}$ and $f_{UMP}/f_{EMP}$) are away from these observations. We outline here the procedure we use in the case of $f_{EMP}$ but note that the treatment is similar for $f_{UMP}/f_{EMP}$. For the 30 Caterpillar merger trees, we first calculate the mean value of the logarithm of $f_{EMP}$, $\mu_{CTP} = \text{mean}[\log_{10}(f_{EMP})]$, and its standard deviation $\sigma_{CTP} = \text{stddev}[\log_{10}(f_{EMP})]$. Then, we calculate the “error” from an observation $i$ as $f_i = |\mu_{CTP} - \log_{10}[f_{EMP,i}]|/\sigma$. Here $i$ stands for the three observations that we compare to, and $\sigma = \sqrt{(0.5)^2 + \sigma_{CTP}^2}$ is the error estimate that includes 0.5 dex observational uncertainty. The main source of this large uncertainty is the stellar mass estimate in the MW halo that ranges from $4 \times 10^8 M_{\odot}$ (Bell et al. 2008) to $10^9 M_{\odot}$ (Bullock & Johnston 2005). Then, we use $\min(f_{\text{Pristine}}, f_{\text{Chiti}}, f_{\text{TOPoS}})$ as the goodness-of-fit evaluator.

Our approach predicts a mean EMP fraction of $\log_{10}[f_{EMP}] = -4.39$, together with a mean EMP-to-UMP fraction of $\log_{10}[f_{EMP}/f_{UMP}] = -1.23$. The value of $\log_{10}(f_{EMP})$ that we obtain lies $0.40\sigma$ away from the observed value. The value of $\log_{10}(f_{UMP}/f_{EMP})$ lies $0.58\sigma$ away from the Bonifacio et al. (2021) value and around $1\sigma$ away from the Chiti et al. (2021) value. We consider these to be good fits, which provide the second- and third-highest goodness-of-fit parameters out of the six observables.

### 3.1.4. Cosmic Star Formation Rate Density

The cosmic SFRD quantifies how many stars form per unit volume per unit time in the universe. A-SLOTH can predict the SFRD at high redshifts, and surveys can estimate the SFRD out to $z \sim 10$. Since the SFRD is a cosmic average, we use the 8 Mpc box as a cosmologically representative sample of the universe. We only have DM halo merger trees for the 8 Mpc box down to redshift $z \sim 5$. So for this observable, we compare the SFRD from A-SLOTH to observations in the redshift range $5 \le z \le 10$.

The SFRD cannot be observed directly. Instead, galaxy counts at high redshift need to be debiased, and the galaxy mass function needs to be extrapolated to the faint end at different redshifts. Hence, the observation-based SFRD is not well constrained at high redshift. Therefore, we will compare our model with two extreme cases. Specifically, we compare our modeled SFRD to the predictions by Madau & Dickinson (2014) and Behroozi et al. (2019). These two estimates differ by up to 1 dex at high redshift. We require that A-SLOTH lie between these observations or penalize the final goodness-of-fit parameter if our prediction is above both or below both observations.

We compare the maximum difference between our model and observations by using the cumulative stellar mass density that has formed up to a certain redshift. To simplify the interpretation, we calculate the cumulative stellar mass density
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Best A-SLOTH model (black) compared to two independent observation-based constraints. We incentivize the result to lie between these two constraints and penalize the maximum distance that our modeled cumulative SFRD is outside the green band.</figcaption></figure>

as a function of cosmic time (and not redshift), so that the resulting unit is stellar mass per unit volume. We compare these cumulative stellar mass densities in the range 300–1300 Myr after the big bang, which is given by the highest redshift for which we have observational constraints and the lowest redshift of the DM merger tree.

Figure 7 shows the two observation-based constraints and our modeled prediction. At higher redshift, the SFRD of A-SLOTH is below the prediction of both Madau & Dickinson (2014) and Behroozi et al. (2019). At each time, we measure how far we are outside the band that is allowed by observations, and we normalize it by the width of this band. Then, we quantify the discrepancy by choosing the time at which this relative discrepancy is largest. In this case, we would then calculate the relative distance as

$$d_{SFR} = \frac{SFR_{B19} - SFR_{asloth}}{SFR_{MD14} - SFR_{B19}}. \qquad (44)$$

Then, we convert this maximum relative discrepancy into a goodness-of-fit parameter via

$$p_{SFR} = \exp(-4d_{SFR}). \qquad (45)$$

This formula provides the correct limits of $p_{SFR} = 1$ for a good agreement between model and observation and an exponential penalty if the prediction of A-SLOTH is outside the observed values. The factor of 4 is an empirical weighting to give this observational constraint similar importance compared to the five observables.

Dubois et al. (2021) show with a 3D hydrodynamical simulation that the largest uncertainty to the SFRD is the cosmic variance in their zoomed-in simulated volume of approximately 4187 cMpc³. Our effective volume is smaller, which might even worsen this effect. Crosby et al. (2013) find with their SAM that cosmic variance between different random realizations of the same cosmic volume can result in SFRD that differs by up to a factor of two. This is of the same order as the observational uncertainty of the SFRD at high redshift. We plan to investigate this variance in detail in a future study in a larger effective volume.

## 3.2. Calibration Strategy

Testing each new combination of parameters requires us to run A-SLOTH using the 30 Caterpillar trees and the 8 Mpc box and takes about 1 hr of wall-clock time. Exploring the entire high-dimensional parameter space to obtain the full posterior distribution is prohibitively expensive for this optimization problem.

Therefore, we first identify three input parameters that are degenerate with other parameters and to which the observables are not very sensitive: $M_{min}$, $\alpha_{III}$, and $f_{esc,III}$. The minimum mass of Pop III stars is degenerate with the SFE, as a smaller value of $M_{min}$ can be compensated with a larger value of $\eta_{III}$. The escape fraction for Pop III stars is degenerate with that for Pop II stars, and the slope of the Pop III IMF is also degenerate with the maximum mass of Pop III stars and their SFE. Furthermore, we have verified that these three parameters have only a small effect on the fit of our six observables. We therefore keep them fixed during the calibration and will study their influence in more detail in a follow-up study.

The values of these three fixed parameters are chosen and motivated as follows. There are no observed metal-free stars in the MW, which constrains $M_{min} \gtrsim 0.8 M_{\odot}$ (Hartwig et al. 2015; Ishiyama et al. 2016; Magg et al. 2019; Rossi et al. 2021). On the other side, we have seen evidence of the chemical fingerprint of SNe from metal-free progenitor stars with masses of $\sim 10 M_{\odot}$ (Ishigaki et al. 2018). We have chosen our fiducial $M_{min} = 5 M_{\odot}$ between those two limits.

A logarithmically flat IMF slope ($\alpha_{III} = 1.0$) puts equal importance on each logarithmically spaced mass bin. This avoids any prior assumptions about dominant mass ranges. Moreover, several hydrodynamical simulations of Pop III star formation yield results that are consistent with this IMF slope (Greif et al. 2011; Stacy et al. 2016; Wollenberg et al. 2020; Sharda et al. 2021; Latif et al. 2022).

The Pop III SFRD in A-SLOTH is rather high compared to other simulations at $z > 20$ as we show below. In order to avoid that the Pop III contribution to reionization requires an even higher SFRD, we set $f_{esc,III} = 60\%$, the highest value that is allowed by theory and observations (Dayal & Ferrara 2018).

In Table 3, we list the six remaining free parameters of our model and the parameter ranges that we have explored while searching for an optimal model fit. We aim to fit six observables with six free parameters. Two of these observables (CSMF and SFRD) are functions (and not just single values), so that the final calibration problem is overconstrained.

We expect further degeneracies of input parameters, and we do not know whether there is one unique global optimum or several local optima instead. The goal of this optimization strategy is therefore not to prove that we found the global optimum. Instead, we try to find a model configuration that reproduces all observables and that is robust to small variations in the input parameters.

We calibrate our model in two steps. First, we run a random parameter exploration in which we sample all input parameters randomly. During this step, we sample 128 different combinations of parameters. From these 128 parameter combinations, we choose the top four models that provide the highest value of $p_{all}$. These four models are then used as the starting point of “Quadient Descent” (QD), a combination of quadratic fits and gradient descent that we have developed to calibrate A-SLOTH.

Instead of optimizing all six input parameters at once, we calibrate A-SLOTH with a sequence of 1D optimizations. At
<table><caption>**Table 3**
List of Free Parameters in A-SLOTH</caption><thead><tr><th>Parameter<br>(1)</th><th>Description<br>(2)</th><th>Explored Range<br>(3)</th><th>Best Fit<br>(4)</th></tr></thead><tbody><tr><td><i>M</i><sub>max</sub></td><td>Maximum mass of Pop III stars</td><td>100 – 330 <i>M</i><sub>☉</sub></td><td>210 <i>M</i><sub>☉</sub></td></tr><tr><td>η<sub>III</sub></td><td>Pop III star formation efficiency</td><td>0.03–1.0</td><td>0.38</td></tr><tr><td>η<sub>II</sub></td><td>Pop II star formation efficiency</td><td>0.01–1.0</td><td>0.19</td></tr><tr><td>α<sub>out</sub></td><td>Slope of outflow efficiency</td><td>0–2.5</td><td>0.86</td></tr><tr><td><i>M</i><sub>out,norm</sub></td><td>Normalization mass of outflow efficiency</td><td>10<sup>9</sup>–10<sup>11</sup> <i>M</i><sub>☉</sub></td><td>7.5 × 10<sup>9</sup> <i>M</i><sub>☉</sub></td></tr><tr><td><i>f</i><sub>esc,II</sub></td><td>Pop II ionizing photon escape fraction</td><td>0–0.4</td><td>0.37</td></tr><tr><td><i>M</i><sub>min</sub></td><td>Minimum mass of Pop III stars</td><td>5 <i>M</i><sub>☉</sub></td><td>5 <i>M</i><sub>☉</sub></td></tr><tr><td>α<sub>III</sub></td><td>Slope of the Pop III IMF</td><td>1</td><td>1</td></tr><tr><td><i>f</i><sub>esc,III</sub></td><td>Pop III ionizing photon escape fraction</td><td>0.6</td><td>0.6</td></tr></tbody></table>

**Note.** Top: six parameters that we calibrate based on observables. Bottom: three further parameters that we do not vary as we explain in the text. The explored range of $M_{max}$ should enclose the mass range of PISNe, and the upper limit for the escape fraction is based on Dayal & Ferrara (2018). Other explored ranges are based on previous experience (Chen et al. 2022) and extensive experiments.

each QD step, we choose a random input parameter that should be explored. Ideally, we want to calculate the gradient $dp_{all}/dx_i$, where $x_i$ is the ith input parameter that is currently under investigation. However, A-SLOTH includes random sampling of variables in various places, so that the resulting gradient $dp_{all}/dx_i$ is not a continuous function of $x_i$.

Therefore, we cannot just sample one additional point to calculate the gradient. Instead, we calculate four additional models and fit a quadratic function to the obtained $p_{all}(x)$. A quadratic fit is the lowest-order polynomial fit that is able to match local optima in the p-landscape, and four sampled points are an efficient trade-off between sample density and computational time. The four explored values, $x_i$, of one parameter are chosen randomly based on a normal distribution centered at the previous value. Eventually, we choose the value of $x_i$ for which the fit predicts the highest $p_{all}$ (inside the sampled parameter range). An example of one such QD step can be seen in Figure 8. After several loops over all input parameters, the QD chains should converge to a similar solution. We stop these QD chains after around 350 iterations and analyze their explored and preferred parameter ranges in Figure 18 in Appendix B.

All chains converge to a local optimum, and we cannot find any combination of input parameters that provides a significantly better fit. Three out of four chains converge to the same optimal parameters. One chain prefers a slightly higher $α_{out}$. This higher value results in a better fit to the CSMF but a worse fit to the MW stellar mass. The values for $p_{all}$ for this chain are similar to the values for the other three chains.

To obtain one set of optimal parameters, we proceed in the following way. We cut the first 60 iterations, which guarantees that we only take sampled parameters into account for which the chains are already exploring near the optimum. Then, we take all the best-fit parameters from all chains together and use the median value as our final best-fit value. Using the median is less sensitive to outlying chains. The resulting best-fit values and observables of our best-fit model are summarized in Table 2. The best-fitting parameters are reported in the last column of Table 3. While the Pop III SFE of $η_{III} = 0.38$ might appear high, we note that it is defined as the fraction of cold gas that forms stars per freefall time. It is therefore not directly comparable to numerical simulations, which usually define the SFE as mass fraction and not per characteristic timescale (e.g.,

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Example of one QD step. We start at the previous value (vertical dotted line), sample four additional values (blue points), fit a quadratic function (blue line), and choose the maximum value as a new parameter (orange star).</figcaption></figure>

Latif & Khochfar 2019). We plan to explore the full posterior distributions of these parameters in a future study.

### 3.3. Further Constraints

Besides these six observables that we actively try to reproduce, we also can compare our results to further results from the literature and from independent simulations. For example, we find a similar metal enrichment history to that in other simulations (Figure 9). Our results are in excellent agreement with Visbal et al. (2020) and in good agreement with other models.

We show the SMHM relations from 30 Caterpillar trees produced by our fiducial model with gray curves in Figure 10. The circles show the mean stellar mass of galaxies at given $M_{vir,peak}$ along with $1\sigma$ among them. The SMHM relation derived by Garrison-Kimmel et al. (2014) is plotted as a red curve, where the curve ends at $M_* = 2.9 \times 10^5 M_{☉}$ owing to the observational completeness. We plot the SMHM relation of Nadler et al. (2020) in gray contours. We find that the difference among SMHM relations from 30 Caterpillar trees is small. The larger scatter at $>10^{11} M_{☉}$ comes from the small statistics. The SMHM relations are consistent with the relations presented by both Garrison-Kimmel et al. (2014) and Nadler et al. (2020) above the observational completeness and consistent with the relation of Nadler et al. (2020) below the observational completeness. We find a steepening of the SMHM relation below $M_{vir,peak} \approx 10^8 M_{☉}$. This is contrary to
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Metal-enriched volume fraction simulated with A-SLOTH in the 8 Mpc box (blue solid) compared to other results from the literature (Pallottini et al. 2014; Salvadori et al. 2014; Jaacks et al. 2018; Visbal et al. 2020).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** We show the SMHMs at $z=0$ of 30 Caterpillar trees from our best model. At given $M_{vir,peak}$, we show the mean $M_*$ (gray circles) and the standard deviation among the galaxies from each tree. We also plot the SMHM derived by Nadler et al. (2020) as a gray contour, where the dark gray shows the $1\sigma$ region and the light gray shows the $2\sigma$ region. The SMHM derived with the AM technique by Garrison-Kimmel et al. (2014) is plotted with a red curve until the observational completeness.</figcaption></figure>

the finding in an earlier study (Chen et al. 2022, which utilized an earlier version of A-SLOTH), where they found a flattening of SMHM relation below $M_{vir,peak} \approx 10^9 M_\odot$. The flattening presented by Chen et al. (2022) comes from the high Pop II SFE (they adopted a fiducial value of $\eta_{II}=2$). This high $\eta_{II}$ introduces a minimal stellar mass that forms before the stellar feedback kicks in and regulates further star formation. In our fiducial model, we have $\eta_{II}=0.19$, which avoids such a threshold. Note that the physical models have also been improved, compared to the one used by Chen et al. (2022), where the Pop III star formation and stellar feedback model and the general chemical model were more simplified than our fiducial model.

In Figure 11, we show MDFs of the 30 Caterpillar trees and three observational estimates (Youakim et al. 2020; Bonifacio et al. 2021; Chiti et al. 2021). For low-metallicity stars, the distribution is reproduced well. For metal-rich stars ([Fe/H] $\gtrsim$-2.5), the contribution from Type Ia SNe is expected to

<figure><img src="image_15.png" /><figcaption>**Figure 11.** We compare the MDF produced with the calibrated model and observations of Youakim et al. (2020), Bonifacio et al. (2021), and Chiti et al. (2021). The black lines show model results for all the Caterpillar trees. Other lines show three different observations. The histogram shows the number of stars in each bin. The sum of the histograms at [Fe/H] < -2 is normalized to be 1 for both model results and observations.</figcaption></figure>

be significant; therefore, we do not claim that A-SLOTH can reproduce the MDF also at higher metallicities. Around [Fe/H] = -3.5, A-SLOTH seems to systematically underproduce the observed MDF. This dip is related to the rather high Pop III SFE in our calibrated model. However, a low SFE would not reproduce the fractions of EMP to UMP stars and of EMP to all stars. So this shape of the MDF is a compromise between other observables. We show below that the model with stochastic feedback is able to reproduce the shape of the MDF better.

We illustrate the cosmic SFRDs for Pop III stars (top panel) and Pop II stars (bottom panel) in Figure 12. Our Pop III SFRD is higher at high redshift. However, when we compare the cumulative stellar mass density of Pop III stars of A-SLOTH with other models, our prediction is in the middle of other literature values. With A-SLOTH we obtain a cumulative Pop III stellar mass density of $9.5 \times 10^4 M_\odot \text{ Mpc}^{-3}$. This is within the range of other models that provide predictions from $1.8 \times 10^3 M_\odot \text{ Mpc}^{-3}$ (Mebane et al. 2018) to $2.9 \times 10^5 M_\odot \text{ Mpc}^{-3}$ (Sarmento et al. 2019).

We also analyze the recovery time between Pop III and Pop II star formation. For this comparison, we define the recovery time as the time between the first Pop III SN and the first Pop II SN in a halo. This timescale is set by the radiative and mechanical feedback of the Pop III star formation and the ability of the halo to retain or accrete new gas. We find recovery times in the range of ~1 Myr – 1 Gyr. This spread is in good agreement with numerical simulations (Jeon et al. 2014; Chiaki et al. 2018; Latif & Schleicher 2020; Magg et al. 2022b; Hicks et al. 2021). Although such a large range is seen in several different simulations, its cause is still poorly understood. Chiaki et al. (2018) relate it to the ratio of the energy injected by feedback (from both photoionization and SNe) to the binding energy of the halo. In this picture, if the total energy output by a star is larger than the binding energy, the halo will undergo a long recovery process of at least several tens of millions of years. Alternatively, Magg et al. (2022b) suggest that it could be caused by variations in the number of substructures present in the minihalos.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Cosmic SFRDs for Pop III stars (top panel) and Pop II stars (bottom panel). The blue lines show our model prediction, the dashed lines show other models from the literature (Trenti & Stiavelli 2009; Johnson et al. 2013; de Souza et al. 2014; Xu et al. 2016a; Jaacks et al. 2018; Mebane et al. 2018; Sarmento et al. 2019; Visbal et al. 2020), the dotted lines in the bottom panel show observations (Madau & Dickinson 2014; Behroozi & Silk 2015; Finkelstein 2016; Behroozi et al. 2019), and the gray dashed-dotted lines in both panels show the total SFRD (Pop II+Pop III).</figcaption></figure>

Altogether, these additional comparisons are further confirmation that our model for baryonic physics in minihalos works well.

### 3.4. Calibration of EPS Mode

To test EPS trees and our implementation of stochastic feedback, we compare three different cases. First, we run a Caterpillar tree with spatial feedback (our fiducial model). Then, we run exactly the same merger tree but ignore the spatial information and use the stochastic feedback instead. Those two runs should produce similar results. Then, we run a third case in which we use an EPS-generated merger tree of a halo with the same mass as the MW with stochastic feedback. Ideally, the essential physics can be implemented appropriately, despite the difference in merger tree details, so that also this last case should produce similar results. However, with EPS, we follow only the formation history of the MW main halo and ignore subhalos and environmental effects. Therefore, we do not expect that this EPS-based run produces identical results to the Caterpillar-based runs. Previously, we have confirmed that external feedback affects the observables and

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Comparison of SFRD (top) and MDF (bottom) between different MW-like simulations. The gray lines are models from the literature (compare Figures 11 and 12). We focus on the difference between our three models, and the SFRD in an MW-like volume should not reproduce the cosmologically representative SFRD. Hence, we do not focus on the comparison to literature here, but rather provide them to guide the eye. The two models that are based on Caterpillar merger trees (blue, orange) provide similar results, whereas the EPS-based model differs at high redshift and low metallicity.</figcaption></figure>

that about 10% of halos are enriched externally before they are enriched internally (Tarumi et al. 2020).

The comparison between different models of external feedback can be seen in Figures 13 and 14. The two runs that are based on the same Caterpillar merger tree result in very similar star formation rates. The MDFs are slightly different because we do not accrete externally enriched metals in the run with stochastic feedback. The run based on EPS merger trees deviates from the spatially resolved merger trees at high redshift and low metallicities.

Figure 14 shows how the volume filling fractions of ionized (top) and metal-enriched (bottom) regions evolve over time. Again, the two runs based on Caterpillar merger trees are very similar. To obtain an estimate of the total volume, we divide
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Volume filling fraction of ionized regions (top) and metal-enriched regions (bottom). We compare three different approaches to model an MW-like galaxy with and without spatial information. The gray lines and symbols are models from the literature (compare Figures 5 and 9). These literature models should only guide the eye since an MW-like volume is not representative and these simulations do not intend to reproduce these literature values. The early peak of ionizing radiation in the EPS-based model is due to a smaller effective volume and due to a higher Pop III SFRD at early times.</figcaption></figure>

the total mass of a merger tree at $z=0$ by the mean cosmic density. Since the MW is not cosmologically representative, we do not expect to reproduce the ionization history or cosmic evolution of metals in these test scenarios.

Overall, we see that our implementation of stochastic feedback is working. However, EPS trees are in general not very well suited to reproduce observables because of their simplified structure that cannot reproduce an MW-like galaxy with its satellites and cosmic environment. A-SLOTH is calibrated and intended for use with spatially resolved DM merger trees. We recommend the use of EPS merger trees in A-SLOTH only for initial tests and the use of spatially resolved merger trees for all scientific applications.

## 4. Summary

A-SLOTH is the first public SAM that was designed to simulate the formation of the first stars and galaxies. In this paper, we have shown how A-SLOTH can reproduce six observables with our calibration strategy: the ionization history

of the universe, the stellar mass of the MW and its satellite galaxies, the relative fractions of EMP stars, and the SFRD at high redshift.

We have shown that A-SLOTH is computationally efficient and has versatile applications (Section 1.1). Users can start experimenting with the memory-friendly EPS mode and then run scientific experiments on publicly available merger trees.

A-SLOTH includes sophisticated models for external feedback (Section 2.5), inhomogeneous metal mixing (Section 2.5.3), and baryonic physics (Section 2.2) inside the first galaxies. All modules and input parameters can be modified by the user in order to study new physics. It is the first SAM that samples the formation of individual Pop III and Pop II stars, which allows us to model feedback by individual SNe in dwarf galaxies.

In addition to this code release paper, we also plan to publish the main data products of A-SLOTH that can directly be used by the community, such as SFRD, recovery time distribution, and MW satellite properties with convenient fitting functions.

We also expect that A-SLOTH can contribute to new scientific discoveries in the next decade, by optimizing surveys with the James Webb Space Telescope or interpreting LISA and third-generation terrestrial detector observations. For example, A-SLOTH can predict the optimal redshift to search for direct detections of Pop III SNe (Magg et al. 2016; Hartwig et al. 2018a) or the contribution of Pop III remnant binary black holes to the data stream of GW detectors (Hartwig et al. 2016b). With these next-generation telescopes and surveys, we will greatly extend the frontier of observational cosmology, which is bound to result in serendipitous discoveries of unknown objects and events. Given its inherent adaptability and computational speed, A-SLOTH is ideal for quickly interpreting such frontier observations and helping to define follow-up campaigns.

We thank Boyuan Liu, Britton Smith, Anna Schauer, and Betül Uysal for valuable discussions. We also thank the anonymous referee for constructive feedback. We acknowledge the work of Lacey & Cole (1993) and Parkinson et al. (2008), who developed GALFORM, based on which we started the implementation of A-SLOTH. We acknowledge funding from JSPS KAKENHI grant Nos. 19K23437 and 20K14464. The team in Heidelberg is thankful for support from the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center (SFB 881, Project-ID 138713538) “The Milky Way System” (subprojects A1, B1, B2, and B8) and from the Heidelberg Cluster of Excellence (EXC 2181-390900948) “STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data,” funded by the German Excellence Strategy. They also acknowledge funding from the European Research Council in the ERC Synergy Grant “ECOGAL—Understanding our Galactic ecosystem: From the disk of the Milky Way to the formation sites of stars and planets” (project ID 855130). M.V. acknowledges funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 grant agreement No. 614199, project “BLACK”).

The project benefited from computing resources provided by the State of Baden-Württemberg through bwHPC and DFG through grant INST 35/1134-1 FUGG, and from the data storage facility SDS@hd supported through grant INST 35/1314-1 FUGG. We also acknowledge the Texas Advanced Computing Center (TACC) for providing HPC resources under
XSEDE allocation TG-AST120024. And we acknowledge the Leibniz Computing Center (LRZ) for providing HPC resources in the project pr74nu. N.Y. acknowledges financial support by a Grant-in-Aid from Japan Society for the Promotion of Science 20H05845.

*Software:* GALFORM (Lacey & Cole 1993; Parkinson et al. 2008), python (Van Rossum & Drake 2009), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), pandas (McKinney 2010; Reback 2022).

# Appendix A
## Supplementary Information on Method

### A.1. Fitting Function of Halo Concentration

We follow the fitting functions in Correa et al. (2015) to compute the halo's concentration at $0 \le z < 10$ and adopt a constant halo concentration of 3.5 at $z > 10$. At $0 \le z < 4$, we have

$$ \begin{aligned} \log_{10} c_{\text{dm}} &= a + b \log_{10}(M_{\text{vir}})(1 + c M_{\text{vir}}^2) \\ a &= 1.7543 - 0.2766(1 + z) + 0.02039(1 + z)^2 \\ b &= 0.2753 + 0.00351(1 + z) - 0.3038(1 + z)^{0.0269} \\ c &= -0.01537 + 0.02102(1 + z)^{-0.1475}, \end{aligned} \quad \text{(A1)} $$

and when $4 \le z < 10$, we have

$$ \begin{aligned} \log_{10} c_{\text{dm}} &= a + b \log_{10}(M_{\text{vir}}) \\ a &= 1.3081 - 0.1078(1 + z) + 0.00398(1 + z)^2 \\ b &= 0.0223 - 0.0944(1 + z)^{-0.3907}. \end{aligned} \quad \text{(A2)} $$

At $z \ge 10$, the halo's concentration in the relevant virial mass range ($M_{\text{vir}} > 10^5 M_{\odot}$) has little variation (Figure 7 in Correa et al. 2015). Therefore, we adopt a constant concentration value, and our results are insensitive to the adopted value.

### A.2. Alternative Models for the Critical Halo Mass

The following two models for the critical halo mass of star formation are implemented as alternatives to our fiducial model (Schauer et al. 2021). In model OHS, we compute $M_{\text{crit}}$ by considering prescriptions in O'Shea & Norman (2008), Hummel et al. (2012), and Stacy et al. (2011). The idea is to combine models that take into account the LW background (O'Shea & Norman 2008) and baryonic streaming (Stacy et al. 2011) with our baseline model (Hummel et al. 2012) that we had used in previous versions of A-SLOTH. O'Shea & Norman (2008) assumed that the critical mass of a halo depends solely on $J_{21}$. We denote this critical mass as $M_{\text{crit,O}}$, given by

$$ M_{\text{crit,O}}/M_{\odot} = 4(1.25 \times 10^5 + 8.7 \times 10^5(4\pi J_{21})^{0.47}). \quad \text{(A3)} $$

Hummel et al. (2012) assumed that gas reaches high density and collapses to form Pop III stars when a halo reaches a critical virial temperature, $T_{\text{crit}} = 2200$ K. The critical mass $M_{\text{crit,H}}$ is then dependent on redshift

$$ M_{\text{crit,H}}/M_{\odot} = 10^6 \left( \frac{T_{\text{crit}}}{1000 \text{ K}} \right)^{1.5} \left( \frac{1+z}{10} \right)^{-1.5}. \quad \text{(A4)} $$

Finally, we follow Stacy et al. (2011) and compute a critical mass $M_{\text{crit,S}}$,

$$ M_{\text{crit,S}}/M_{\odot} = \frac{\pi v_{\text{eff}}^3}{6G^{3/2} \rho^{1/2}}, \quad \text{(A5)} $$

where $\rho$ is the mean DM density of the halo and $v_{\text{eff}} = \sqrt{v_{\text{BC}}^2 \sigma_{\text{rms}}(z)^2 + c_s^2}$ is the effective velocity of the gas. This depends on the sound speed in the IGM, $c_s = \sqrt{k_B T / \mu m_{\text{H}}}$, where $T = 0.017(1+z)^2$ K prior to the onset of Pop III star formation (Schneider 2015), and on the redshift-dependent streaming velocity

$$ \frac{v_{\text{BC}} \sigma_{\text{rms}}(z)}{\text{cm s}^{-1}} = v_{\text{BC}} \times \frac{6 \times 10^5}{201} \times (1 + z), \quad \text{(A6)} $$

where $v_{\text{BC}}$ is the streaming velocity in units of the rms value. In model OHS, the combined critical mass is determined by taking the maximum of $M_{\text{crit,O}}$, $M_{\text{crit,H}}$, and $M_{\text{crit,S}}$, i.e.,

$$ M_{\text{crit,OHS}} = \max(M_{\text{crit,O}}, M_{\text{crit,H}}, M_{\text{crit,S}}). \quad \text{(A7)} $$

The third $M_{\text{crit}}$ model we implement in A-SLOTH is based on the formula given in Fialkov et al. (2013), and we denote this model as F13. This critical mass $M_{\text{crit,F13}}$ is computed by

$$ M_{\text{crit,F13}}/M_{\odot} = M_0(1 + 6.96(4\pi J_{21})^{0.47}), \quad \text{(A8)} $$

where $J_{21}$ depends on redshift via $J_{21} = 10^{2-z/5}$, and $M_0$ is the critical mass when there is no LW background. We compute $M_0$ based on Fialkov et al. (2012):

$$ M_0/M_{\odot} = \left( \frac{v_{\text{cool}}(z)}{146.6 \text{km s}^{-1}} \right)^3 \times \Omega_{m,0}^{-0.5} \times (1+z)^{-1.5} \times 10^{12} h^{-1}, \quad \text{(A9)} $$

where $\Omega_{m,0} = 0.3086$ is the matter density at $z=0$ and the redshift-dependent circular velocity threshold $v_{\text{cool}}(z)$ is given by

$$ \frac{v_{\text{cool}}(z)}{\text{km s}^{-1}} = \sqrt{(3.714)^2 + (4.015 \times 10^{-5} v_{\text{BC}} \sigma_{\text{rms}}(z))^2}, \quad \text{(A10)} $$

where $v_{\text{BC}}\sigma_{\text{rms}}(z)$ is given by Equation (A6) above. The full derivation of $M_0$ can be found in Appendix D in Chen et al. (2022).

In Figure 15 we show the evolution of $M_{\text{crit}}$ versus redshift for several different values of $v_{\text{BC}}$ for the three implemented models.

### A.3. Technical Implementation of Spatial Feedback

Determining which ionized and enriched bubbles a halo sits inside of in its simplest implementation requires computing the distance between each possible pair of halos. This would lead to the computational cost of the simulations scaling as $N^2$ for N halos and implies having to compute trillions of distances on each time step. To avoid this, we use a tree-based approach that for a given target halo limits the number of halos to which we need to compute the distances.

In this approach we build an oct-tree of bubbles and consecutively add all bubbles during a time step to it. The tree is started with a head-node encompassing the entire physical volume. This node is split into eight subnodes (i.e., it is split in two along each dimension), and the subnodes are split into eight smaller nodes. We limit the tree to a depth of 20 levels but only allocate nodes that actually contain at least one bubble.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Redshift evolution of $M_{\text{crit}}$ in our fiducial (S21) model and two other models, computed for several different values of the streaming velocity: $v_{\text{BC}} = 0, 0.8, 2$, and 3 times $\sigma_{\text{rms}}$. The values of $M_{\text{crit}}$ in the OHS, S21, and F13 models are plotted in green, red, and blue, respectively. We additionally plot $M_{\text{ACH}}$ in gray for reference.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Illustration of the oct-tree construction for determining which bubbles a halo lies inside of. The colors represent the levels (black: first level; blue: second level; orange: third level). The bubbles (circles) are colored according to the node they are associated with. Each bubble is associated with the smallest tree node it is completely contained within.</figcaption></figure>

Each bubble is associated with exactly one node of this oct-tree, namely, the smallest node that fully encompasses the bubble. Thus, large bubbles and bubbles that are very close to the edges of the larger nodes are located in the first few levels, and small bubbles are located in the lower levels. This is illustrated in Figure 16.

<table><caption>**Table 4**

Parameters to Be Used for the External Enrichment Fitting Function</caption><thead><tr><th>Metallicity</th><th>μ</th><th>λ</th><th>σ</th></tr></thead><tbody><tr><td>-2 &lt; Z &le; -1</td><td>-0.07</td><td>0.14</td><td>2.78</td></tr><tr><td>-3 &lt; Z &le; -2</td><td>-0.01</td><td>0.24</td><td>1.41</td></tr><tr><td>-4 &lt; Z &le; -3</td><td>0.08</td><td>0.20</td><td>0.77</td></tr><tr><td>-5 &lt; Z &le; -4</td><td>0.16</td><td>0.19</td><td>0.38</td></tr><tr><td>-8 &lt; Z &le; -5</td><td>1.63</td><td>1.01</td><td>0.35</td></tr></tbody></table>

**Note.** These parameter sets are for “-dZ,” not dZ itself.

In order to find all bubbles that overlap with a point (e.g., the center of a halo), one only has to test bubbles associated with tree nodes the point is inside of. This algorithm results in identical results to testing all pairwise combinations of halos at a substantially reduced computational cost.

### A.4. Fitting Functions for Inhomogeneous Metal Mixing

In Tarumi et al. (2020), we obtain fitting functions for internal and external enrichment separately. For internal enrichment, we use a Gaussian distribution function with $\mu = -0.03$, $\sigma = 0.15$:

$$p(x) = \frac{1}{\sqrt{2\pi}(0.15)^2} \exp\left[\frac{(x + 0.03)^2}{2(0.15)^2}\right]. \quad (A11)$$

For external enrichment, we use an exponentially modified Gaussian distribution. The distribution has three free parameters, ($K$, $\mu$, $\sigma$), and the probability distribution function $p(x; \mu, \sigma, \lambda)$ is

$$p(x; \mu, \sigma, \lambda) = \frac{\lambda}{2} \exp\left[\frac{\lambda}{2}(2\mu + \lambda\sigma^2 - 2x)\right] \times \text{erfc}\left(\frac{\mu + \lambda\sigma^2 - x}{\sqrt{2}\sigma}\right), \quad (A12)$$

where

$$\text{erfc}(x) = 1 - \text{erf}(x) \quad (A13)$$

$$= \frac{2}{\sqrt{\pi}} \int_{x}^{\infty} e^{-t^2} dt, \quad (A14)$$

and the three parameters are given in Table 4.

## Appendix B

### Supplementary Information on Calibration

#### B.1. Random Parameter Exploration

The first step of our calibration is that we run 128 different models with random input parameter combinations. In Figure 17 we show how the goodness-of-fit parameters depend on the input parameters. For example, the optical depth is sensitive to the ionizing escape fraction, and the CSMF of MW satellites is sensitive to our normalization mass of SN feedback. Furthermore, the maximum mass and SFE of Pop III stars are sensitive to EMP-based quantities. In the last row ($p_{\text{all}}$) we see which parameter ranges are preferred or disfavored. For example, models with $M_{\text{max}} > 200 M_{\odot}$ can reach higher values of $p_{\text{all}}$ than models with $M_{\text{max}} < 200 M_{\odot}$. Overall, this representation demonstrates that such a random parameter exploration in this 6D space is not effective to find a global optimum. Hence, we only use it as a first exploration of the parameter space and then optimize the calibration with QD.
<figure><img src="image_21.png" /><figcaption>**Figure 17.** Correlations between the input parameters (horizontal) and goodness-of-fit parameters (vertical) for 128 randomly chosen combinations of input parameters. This plot illustrates which observables are sensitive to which input parameters.</figcaption></figure>
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 18.** Iterative exploration and optimization of the input parameters (left) and how the goodness-of-fit parameter changes during this optimization (right).</figcaption></figure>
## B.2. Quadrant Descent

In Figure 18, the four colors show four independent chains. Each point shows one explored model, and each step consists of four different models. The four chains converge to optimal values for all input parameters. Only the orange chain prefers a slightly higher value for $\alpha_{out}$. As we can see in the bottom two panels, this higher value of the feedback efficiency exponent yields a better fit to the CSMF but a worse fit to the MW stellar mass. This is a legitimate trade-off between two individual goodness-of-fit parameters, and $p_{all}$ converges to similar values for all four chains.

## ORCID iDs

Tilman Hartwig https://orcid.org/0000-0001-6742-8843

Mattis Magg https://orcid.org/0000-0002-9022-5136

Li-Hsin Chen https://orcid.org/0000-0003-0475-1947

Yuta Tarumi https://orcid.org/0000-0002-9801-7788

Volker Bromm https://orcid.org/0000-0003-0212-2979

Simon C. O. Glover https://orcid.org/0000-0001-6708-1317

Alexander P. Ji https://orcid.org/0000-0002-4863-8842

Ralf S. Klessen https://orcid.org/0000-0002-0560-3172

Muhammad A. Latif https://orcid.org/0000-0003-2480-0988

Marta Volonteri https://orcid.org/0000-0002-3216-1322

Naoki Yoshida https://orcid.org/0000-0001-7925-238X

## References

Agarwal, B., Johnson, J. L., Zackrisson, E., et al. 2016, MNRAS, 460, 4003

Agarwal, B., Khochfar, S., Johnson, J. L., et al. 2012, MNRAS, 425, 2854

Ahn, K., Shapiro, P. R., Iliev, I. T., Mellema, G., & Pen, U.-L. 2009, ApJ, 695, 1430

Arentsen, A., Starkenburg, E., Martin, N. F., et al. 2020, MNRAS, 496, 4964

Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Natur, 553, 473

Barausse, E. 2012, MNRAS, 423, 2533
Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013b, ApJ, 763, 18

Bell, E. F., Zucker, D. B., Belokurov, V., et al. 2008, ApJ, 680, 295

Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440

Bonifacio, P., Monaco, L., Salvadori, S., et al. 2021, A&A, 651, A79

Bromm, V., & Loeb, A. 2003, Natur, 425, 812

Bromm, V., & Yoshida, N. 2011, ARA&A, 49, 373

Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343

Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931

Chatterjee, A., Dayal, P., Choudhury, T. R., & Schneider, R. 2020, MNRAS, 496, 1445

Chen, L.-H., Magg, M., Hartwig, T., et al. 2022, MNRAS, 513, 934

Chiaki, G., Susa, H., & Hirano, S. 2018, MNRAS, 475, 4378

Chiaki, G., Tominaga, N., & Nozawa, T. 2017, MNRAS, 472, L115

Chiti, A., Mardini, M. K., Frebel, A., & Daniel, T. 2021, ApJL, 911, L23

Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168

Correa, C. A., Wyithe, J. S. B., Schaye, J., & Duffy, A. R. 2015, MNRAS, 452, 1217

Côté, B., Ritter, C., Herwig, F., et al. 2017, in 14th Int. Symp. on Nuclei in the Cosmos (NIC2016), ed. S. Kubono et al. (Tokyo: JPS)

Côté, B., Silvia, D. W., O'Shea, B. W., Smith, B., & Wise, J. H. 2018, ApJ, 859, 67

Crosby, B. D., O'Shea, B. W., Smith, B. D., Turk, M. J., & Hahn, O. 2013, ApJ, 773, 108

Dayal, P., Choudhury, T. R., Pacucci, F., & Bromm, V. 2017, MNRAS, 472, 4414

Dayal, P., & Ferrara, A. 2018, PhR, 780, 1

Dayal, P., Ferrara, A., Dunlop, J. S., & Pacucci, F. 2014, MNRAS, 445, 2545

Dayal, P., Rossi, E. M., Shirailou, B., et al. 2019, MNRAS, 486, 2336

Dayal, P., Volonteri, M., Choudhury, T. R., et al. 2020, MNRAS, 495, 3065
de Bennassuti, M., Salvadori, S., Schneider, R., Valiante, R., & Omukai, K. 2017, MNRAS, 465, 926

de Souza, R. S., Ishida, E. E. O., Whalen, D. J., Johnson, J. L., & Ferrara, A. 2014, MNRAS, 442, 1640

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Ekström, S., Meynet, G., Chiappini, C., Hirschi, R., & Maeder, A. 2008, A&A, 489, 685

Emerick, A., Bryan, G. L., Mac Low, M.-M., et al. 2018, ApJ, 869, 94

Engler, C., Pillepich, A., Pasquali, A., et al. 2021, MNRAS, 507, 4211

Ferland, G. J., Peterson, B. M., Horne, K., Welsh, W. F., & Nahar, S. N. 1992, ApJ, 387, 95

Fialkov, A., Barkana, R., Tseliakhovich, D., & Hirata, C. M. 2012, MNRAS, 424, 1335

Fialkov, A., Barkana, R., Visbal, E., Tseliakhovich, D., & Hirata, C. M. 2013, MNRAS, 432, 2909

Finkelstein, S. L. 2016, PASA, 33, e037

Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631

Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271

Garrison-Kimmel, S., Boylan-Kolchin, M., Bullock, J. S., & Lee, K. 2014, MNRAS, 438, 2578

Glover, S. 2005, SSRv, 117, 445

Gnedin, N. Y., & Kaurov, A. A. 2014, ApJ, 793, 30

Gómez, F. A., Coleman-Smith, C. E., O'Shea, B. W., Tumlinson, J., & Wolpert, R. L. 2014, ApJ, 787, 20

Graziani, L., Salvadori, S., Schneider, R., et al. 2015, MNRAS, 449, 3137

Graziani, L., de Bennassuti, M., Schneider, R., Kawata, D., & Salvadori, S. 2017, MNRAS, 469, 1101

Greif, T. H. 2015, ComAC, 2, 3

Greif, T. H., & Bromm, V. 2006, MNRAS, 373, 128

Greif, T. H., Springel, V., White, S. D. M., et al. 2011, ApJ, 737, 75

Greig, B., & Mesinger, A. 2017, MNRAS, 465, 4838

Griffen, B. F., Dooley, G. A., Ji, A. P., et al. 2018, MNRAS, 474, 443

Griffen, B. F., Ji, A. P., Dooley, G. A., et al. 2016, ApJ, 818, 10

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hartwig, T., Bromm, V., Klessen, R. S., & Glover, S. C. O. 2015, MNRAS, 447, 3892

Hartwig, T., Bromm, V., & Loeb, A. 2018a, MNRAS, 479, 2202

Hartwig, T., Latif, M. A., Magg, M., et al. 2016a, MNRAS, 462, 2184

Hartwig, T., Volonteri, M., Bromm, V., et al. 2016b, MNRAS, 460, L74

Hartwig, T., & Yoshida, N. 2019, ApJL, 870, L3

Hartwig, T., Yoshida, N., Magg, M., et al. 2018b, MNRAS, 478, 1795

Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532

Heger, A., & Woosley, S. E. 2010, ApJ, 724, 341

Hicks, W. M., Wells, A., Norman, M. L., et al. 2021, ApJ, 909, 70

Ho, I. T., Kreckel, K., Meidt, S. E., et al. 2019, ApJL, 885, L31

Hummel, J. A., Pawlik, A. H., Milosavljević, M., & Bromm, V. 2012, ApJ, 755, 72

Hunter, J. D. 2007, CSE, 9, 90

Iliev, I. T., Mellema, G., Ahn, K., et al. 2014, MNRAS, 439, 725

Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27

Ishigaki, M. N., Hartwig, T., Tarumi, Y., et al. 2021, MNRAS, 506, 5410

Ishigaki, M. N., Tominaga, N., Kobayashi, C., & Nomoto, K. 2018, ApJ, 857, 46

Ishiyama, T., Sudo, K., Yokoi, S., et al. 2016, ApJ, 826, 9

Jaacks, J., Thompson, R., Finkelstein, S. L., & Bromm, V. 2018, MNRAS, 475, 4396

Jeon, M., Pawlik, A. H., Bromm, V., & Milosavljević, M. 2014, MNRAS, 444, 3288

Ji, A. P., Koposov, S. E., Li, T. S., et al. 2021, ApJ, 921, 32

Jiang, F., & van den Bosch, F. C. 2014, MNRAS, 440, 193

Johnson, J. L., Dalla Vecchia, C., & Khochfar, S. 2013, MNRAS, 428, 1857

Kannan, R., Garaldi, E., Smith, A., et al. 2022, MNRAS, 511, 4005

Klessen, R. 2019, in Formation of the First Black Holes, ed. M. Latif & D. Schleicher (Singapore: World Scientific Publishing), 67

Klessen, R. S., & Glover, S. C. O. 2016, Star Formation in Galaxy Evolution: Connecting Numerical Models to Reality, Saas-Fee Advanced Course, Vol. 43 (Berlin: Springer-Verlag), 85

Knebe, A., Pearce, F. R., Thomas, P. A., et al. 2015, MNRAS, 451, 4029

Kobayashi, C., Karakas, A. I., & Umeda, H. 2011, MNRAS, 414, 3231

Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Komiya, Y. 2011, ApJ, 736, 73

Komiya, Y., Habe, A., Suda, T., & Fujimoto, M. Y. 2009, ApJL, 696, L79

Komiya, Y., Habe, A., Suda, T., & Fujimoto, M. Y. 2010, ApJ, 717, 542

Komiya, Y., Suda, T., & Fujimoto, M. Y. 2016, ApJ, 820, 59

Komiya, Y., Suda, T., Yamada, S., & Fujimoto, M. Y. 2020, ApJ, 890, 66

Konno, A., Ouchi, M., Shibuya, T., et al. 2018, PASJ, 70, S16
Kroupa, P. 2001, MNRAS, 322, 231

Kulkarni, M., Visbal, E., & Bryan, G. L. 2021, ApJ, 917, 40

Lacey, C., & Cole, S. 1993, MNRAS, 262, 627

Latif, M. A., & Ferrara, A. 2016, PASA, 33, e051

Latif, M. A., & Khochfar, S. 2019, MNRAS, 490, 2706

Latif, M. A., & Schleicher, D. 2020, ApJL, 902, L31

Latif, M. A., Whalen, D., & Khochfar, S. 2022, ApJ, 925, 28

Liu, B., & Bromm, V. 2021, MNRAS, 506, 5451

Mac Low, M.-M., & Ferrara, A. 1999, ApJ, 513, 142

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magg, M., Hartwig, T., Agarwal, B., et al. 2018, MNRAS, 473, 5308

Magg, M., Hartwig, T., Chen, L.-H., & Tarumi, Y. 2022c, JOSS, 7, 4417

Magg, M., Hartwig, T., Glover, S. C. O., Klessen, R. S., & Whalen, D. J. 2016, MNRAS, 462, 3591

Magg, M., Klessen, R. S., Glover, S. C. O., & Li, H. 2019, MNRAS, 487, 486

Magg, M., Reis, I., Fialkov, A., et al. 2022a, MNRAS, 514, 4433

Magg, M., Schauer, A. T. P., Klessen, R. S., et al. 2022b, ApJ, 929, 119

Marigo, P., Girardi, L., Chiosi, C., & Wood, P. R. 2001, A&A, 371, 152

Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, 856, 2

McConnachie, A. W. 2012, AJ, 144, 4

McGreer, I. D., Mesinger, A., & D'Odorico, V. 2015, MNRAS, 447, 499

McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148

McKinney, W. 2010, in Proc. of the 9th Python in Science Conf., ed. S van der Walt & J. Millman (Austin, TX: SciPy), 56

McMillan, P. J. 2017, MNRAS, 465, 76

Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2018, MNRAS, 479, 4544

Mesinger, A., Aykutalp, A., Vanzella, E., et al. 2015, MNRAS, 446, 566

Muñoz, R. R., Côté, P., Santana, F. A., et al. 2018, ApJ, 860, 66

Nadler, E. O., Wechsler, R. H., Bechtol, K., et al. 2020, ApJ, 893, 48

Naoz, S., Yoshida, N., & Gnedin, N. Y. 2012, ApJ, 763, 27

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457

O'Shea, B. W., & Norman, M. L. 2008, ApJ, 673, 14

O'Shea, B. W., Wise, J. H., Xu, H., & Norman, M. L. 2015, ApJL, 807, L12

Ouchi, M., Harikane, Y., Shibuya, T., et al. 2018, PASJ, 70, S13

Pacucci, F., Loeb, A., & Salvadori, S. 2017, MNRAS Lett., 471, L72

Pacucci, F., Loeb, A., Mezcua, M., & Martín-Navarro, I. 2018, ApJL, 864, L6

Pallottini, A., Ferrara, A., Gallerani, S., Salvadori, S., & D'Odorico, V. 2014, MNRAS, 440, 2498

Pandya, V., Fielding, D. B., Anglés-Alcázar, D., et al. 2021, MNRAS, 508, 2979

Parkinson, H., Cole, S., & Helly, J. 2008, MNRAS, 383, 557

Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton, NJ: Princeton Univ. Press)

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Planck Collaboration 2016, A&A, 596, A108

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Planck Collaboration et al. 2020, A&A, 641, A6

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Rahmati, A., & Schaye, J. 2018, MNRAS, 478, 5123

Reback, J., Jbrockmendel, McKinney, W., et al. 2022, pandas-dev/pandas: Pandas v1.4.1 doi:10.5281/zenodo.6053272

Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJL, 802, L19

Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71

Roman-Duval, J., Jackson, J. M., Heyer, M., Rathborne, J., & Simon, R. 2010, ApJ, 723, 492

Rossi, M., Salvadori, S., & Skúladóttir, Á. 2021, MNRAS, 503, 6026

Sassano, F., Schneider, R., Valiante, R., et al. 2021, MNRAS, 506, 613

Salvadori, S., Ferrara, A., Schneider, R., Scannapieco, E., & Kawata, D. 2010, MNRAS Lett., 401, L5

Salvadori, S., Schneider, R., & Ferrara, A. 2007, MNRAS, 381, 647

Salvadori, S., Tolstoy, E., Ferrara, A., & Zaroubi, S. 2014, MNRAS, 437, L26

Sarmento, R., Scannapieco, E., & Côté, B. 2019, ApJ, 871, 206

Schaerer, D. 2002, A&A, 382, 28

Schauer, A. T. P., Glover, S. C. O., Klessen, R. S., & Ceverino, D. 2019, MNRAS, 484, 3510

Schauer, A. T. P., Glover, S. C. O., Klessen, R. S., & Clark, P. 2021, MNRAS, 507, 1775

Schneider, P. 2015, Extragalactic Astronomy and Cosmology: An Introduction (Berlin: Springer-Verlag Berlin Heidelberg)

Secunda, A., Cen, R., Kimm, T., Göteborg, Y., & de Mink, S. E. 2020, ApJ, 901, 72

Sestito, F., Longeard, N., Martin, N. F., et al. 2019, MNRAS, 893, 89

Sharda, P., Federrath, C., Krumholz, M. R., & Schleicher, D. R. G. 2021, MNRAS, 503, 2014

Simon, J. D. 2018, ApJ, 863, 89

Sobral, D., Matthee, J., Darvish, B., et al. 2015, ApJ, 808, 139

Spitzer, L. 1978, Physical Processes in the Interstellar Medium (New York: Wiley-Interscience)

Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Natur, 435, 629

Stacy, A., Bromm, V., & Lee, A. T. 2016, MNRAS, 462, 1307

Stacy, A., Bromm, V., & Loeb, A. 2011, ApJL, 730, L1

Stahler, S. W., & Palla, F. 2005, The Formation of Stars (New York: Wiley)

Starkenburg, E., Helmi, A., De Lucia, G., et al. 2013, MNRAS, 429, 725

Tarumi, Y., Hartwig, T., & Magg, M. 2020, ApJ, 897, 58

Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2021, A&A, 653, A154

Trenti, M., & Stiavelli, M. 2009, ApJ, 694, 879

Trinca, A., Schneider, R., Valiante, R., et al. 2022, MNRAS, 511, 616

Tumlinson, J. 2006, ApJ, 641, 1

Tumlinson, J. 2010, ApJ, 708, 1398

Valiante, R., Schneider, R., Volonteri, M., & Omukai, K. 2016, MNRAS, 457, 3356

Valiante, R., Schneider, R., Zappacosta, L., et al. 2018, MNRAS, 476, 407

Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (CreateSpace: Scotts Valley, CA)

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Visbal, E., Bryan, G. L., & Haiman, Z. 2017, MNRAS, 469, 1456

Visbal, E., Bryan, G. L., & Haiman, Z. 2020, ApJ, 897, 95

Visbal, E., Haiman, Z., & Bryan, G. L. 2015, MNRAS, 453, 4456

Visbal, E., Haiman, Z., & Bryan, G. L. 2018, MNRAS, 475, 5246

Volonteri, M., Haardt, F., & Madau, P. 2003, ApJ, 582, 559

Wollenberg, K. M. J., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2020, MNRAS, 494, 1871

Woods, T. E., Agarwal, B., Bromm, V., et al. 2019, PASA, 36, e027

Xu, H., Ahn, K., Norman, M. L., Wise, J. H., & O'Shea, B. W. 2016a, ApJL, 832, L5

Xu, H., Wise, J. H., Norman, M. L., Ahn, K., & O'Shea, B. W. 2016b, ApJ, 833, 84

Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645

Yoshikawa, K., Tanaka, S., & Yoshida, N. 2021, SC 21: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, arXiv:2110.15867

Youakim, K., Starkenburg, E., Martin, N. F., et al. 2020, MNRAS, 492, 4986

Zheng, Z.-Y., Wang, J., Rhoads, J., et al. 2017, ApJL, 842, L22