# LIMFAST. I. A Semi-Numerical Tool for Line Intensity Mapping

Lluís Mas-Ribas  Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA  California
Institute of Technology, 1200 E. California Blvd, Pasadena, CA 91125, USA
Guochao Sun  California Institute of Technology, 1200 E. California Blvd,
Pasadena, CA 91125, USA
Tzu-Ching Chang  Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA  California
Institute of Technology, 1200 E. California Blvd, Pasadena, CA 91125, USA
Michael O. Gonzalez  California Institute of Technology, 1200 E. California
Blvd, Pasadena, CA 91125, USA
Richard H. Mebane  Department of Astronomy and Astrophysics, University of
California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA
©2023. All rights reserved.

###### Abstract

We present LIMFAST, a semi-numerical code for simulating high-redshift galaxy
formation and cosmic reionization as revealed by multi-tracer line intensity
mapping (LIM) signals. LIMFAST builds upon and extends the 21cmFAST code
widely used for 21 cm cosmology by implementing state-of-the-art models of
galaxy formation and evolution. The metagalactic radiation background,
including the production of various star-formation lines, together with the 21
cm line signal tracing the neutral intergalactic medium (IGM), are self-
consistently described by photoionization modeling and stellar population
synthesis coupled to the galaxy formation model. We introduce basic structure
and functionalities of the code, and demonstrate its validity and capabilities
by showing broad agreements between the predicted and observed evolution of
cosmic star formation, IGM neutral fraction, and metal enrichment. We also
present the LIM signals of 21 cm, Ly$\alpha$, H$\alpha$, H$\beta$, [O II], and
[O III] lines simulated by LIMFAST, and compare them with results from the
literature. We elaborate on how several major aspects of our modeling
framework, including models of star formation, chemical enrichment, and
photoionization, may impact different LIM observables and thus become testable
once applied to observational data. LIMFAST aims at being an efficient and
resourceful tool for intensity mapping studies in general, exploring a wide
range of scenarios of galaxy evolution and reionization and frequencies over
which useful cosmological signals can be measured.

††journal: ApJ††software: ARES (Mirocha et al., 2017), BPASS (Eldridge et al.,
2017), cloudy (Ferland et al., 2017), 21cmFAST (Mesinger et al., 2011)

##  1 Introduction

Line intensity mapping (LIM) provides a statistical approach to the study of
the formation and evolution of galaxies and large-scale structure (LSS) in the
universe. Compared to more traditional observational techniques that are
limited to the individually detectable bright sources, LIM measurements of
galaxies take into account the emission produced by the entire galaxy
population present in large areas of the sky (Madau et al., 1997; Suginohara
et al., 1999; Visbal & Loeb, 2010; Kovetz et al., 2019). This characteristic
is especially relevant for high redshift studies, including the epochs of
reionization and cosmic dawn. The LIM approach can prove beneficial here
because the faint end of the galaxy population may have played a major role
during these early times, but this is difficult to explore directly (e.g.,
Fontanot et al., 2012; Choudhury & Ferrara, 2007; Robertson et al., 2015; Yue
et al., 2018).

A number of emission lines resulting from different radiative processes and
phases of the interstellar medium (ISM) are taken into account for LIM
studies, the usual ones being the [C II] line at 158 $\mu$m (e.g., Gong et
al., 2012; Silva et al., 2015; Yue et al., 2015; Dumitru et al., 2019; Yue &
Ferrara, 2019; Sun et al., 2021b), those of the CO molecule (e.g., Righi et
al., 2008; Gong et al., 2011; Lidz et al., 2011; Pullen et al., 2013; Li et
al., 2016; Chung et al., 2019; Ihle et al., 2019), the hydrogen 21 cm spin-
flip transition (e.g., Scott & Rees, 1990; Madau et al., 1997; Furlanetto et
al., 2006; Chang et al., 2008; Visbal et al., 2009; Chang et al., 2010;
Pritchard & Loeb, 2012; Switzer et al., 2013; Liu & Shaw, 2020), and the
potentially bright rest-frame optical/ultraviolet (UV) lines such as
H$\alpha$, H$\beta$, Ly$\alpha$, He II, [O II], and [O III] among others
(e.g., Silva et al., 2013; Pullen et al., 2014; Visbal et al., 2015; Comaschi
& Ferrara, 2016; Heneka et al., 2017; Gong et al., 2017; Visbal & McQuinn,
2018; Mas-Ribas & Chang, 2020; Heneka & Cooray, 2021; Kannan et al., 2022b;
Padmanabhan et al., 2021; Parsons et al., 2022). Distinct from past
theoretical studies focusing on individual line tracers, an increasing number
of recent modeling efforts aim at building a unified and self-consistent
framework for multi-tracer investigations, which may provide a much more
detailed and complete picture of high-redshift universe of interest (e.g., Sun
et al., 2019; Yang et al., 2021; Béthermin et al., 2022; Kannan et al.,
2022b).

Numerically modeling LIM data is of major importance to guide future missions
and experiments, but it is computationally challenging because the statistical
power of LIM resides in the analysis of line emission signals over large areas
of the sky and at multiple frequencies that are sensitive to small-scale
physics. Simulations need to include both detailed processes related to star
formation and the emission and transport of radiation in large cosmological
volumes. This combination of a broad range of dynamical scales is demanding:
numerical simulations accounting for resolved ISM-scale galaxy physics
typically only exist for a small number of galaxies (e.g., Hopkins et al.,
2018; Katz et al., 2019; Pallottini et al., 2019; Kannan et al., 2020). On the
other hand, simulations covering large volumes tend to lack numerical
precision at the smallest scales, and the processes connected to star
formation are often modeled by means of sub-grid prescriptions (e.g.,
Vogelsberger et al., 2014; Eide et al., 2018, 2020; Shen et al., 2020, 2022;
Kannan et al., 2022a, b; Lewis et al., 2022, and the review by Vogelsberger et
al. 2020). Overall, in all cases these simulations typically require
substantial computational resources and are thus not well-suited for parameter
space exploration and model inference.

With these limitations and constraints in mind, we present here LIMFAST, a
semi-numerical tool designed for flexible modeling of high-redshift LIM
signals. LIMFAST aims at self-consistently simulating line emission from
galaxies and the intergalactic medium (IGM), over scales of several hundreds
of Mpc and spanning the epoch of cosmic reionization, in a matter of hours
with a current personal computer. LIMFAST builds upon and uses the 21cmFAST
code (Mesinger & Furlanetto, 2007; Mesinger et al., 2011) to compute the
underlying large-scale structure in large volumes of the universe. This
computational step is rapidly achieved because 21cmFAST uses perturbation
theory and analytical approaches to approximate the evolution of the density
and velocity fields, as well as the formation of collapsed objects. LIMFAST
inherits these calculations and applies analytic galaxy formation models
adopted from Furlanetto (2021) to it to compute the radiation fields for a
number of emission lines, in addition to the original 21 cm line from
21cmFAST. The progress of reionization is simultaneously computed also
following the approach in 21cmFAST, but with extensions that self-consistently
include the emission from galaxy populations that co-evolve with redshift.

This paper is the first of a series that introduces the main structure of the
LIMFAST code, and demonstrates the validity and capabilities of it to be used
to address science questions related to the epoch of cosmic reionization. In
the second paper, Sun et al. (2022, Paper II hereafter) presents the
computations to include the [C II] 158 $\mu$m and CO line emission and
explores the effects and observational implications of different feedback and
star formation prescriptions beyond the fiducial case presented in this work.
Furthermore, a progenitor version of LIMFAST was presented and used in Parsons
et al. (2022) to address the application of LIM to measure the average He
II/H$\alpha$ line ratio for inferring the initial mass function (IMF) of
Population III (Pop III) stars. The structure of this present paper is as
follows: the models and key calculations implemented in LIMFAST are detailed
in Section 2, and the results from our fiducial simulation runs are shown and
compared with the literature in Section 3. We discuss some key considerations
about these results, focusing on assumptions that give rise to qualitative
differences in observables that may be tested by future LIM data in Section 4,
before concluding in Section 5. A flat, $\Lambda$CDM cosmology consistent with
recent measurements by Planck (Planck Collaboration et al., 2016) is assumed
throughout.

##  2 LIMFAST: the Code

We detail the LIMFAST code below, after a brief description of 21cmFAST. In
Section 2.1, we present the galaxy model by Furlanetto (2021) and its
implementation in LIMFAST. We next describe in Section 2.2 the usage of BPASS
to create the stellar spectra, and the photoionization calculations with
cloudy to obtain the nebular radiation. The details of the ionization
computations in the IGM are presented in Section 2.3, and the calculations of
the intergalactic and background Ly$\alpha$, as well as the 21 cm emission,
are presented in Sections 2.4, 2.5 and 2.6, respectively. Finally, in Section
2.7, we describe the inclusion of redshift-space distortions (RSD) in the
calculations.

LIMFAST builds upon and extends 21cmFAST (Mesinger & Furlanetto, 2007;
Mesinger et al., 2011) after inheriting the large-scale density and velocity
fields computed by the latter at each simulated cell. In detail, 21cmFAST
computes evolving density and velocity fields from a set of initial conditions
and Lagrangian perturbation theory (Zel’Dovich, 1970; Scoccimarro, 1998).
Then, the code makes use of the extended Press-Schechter formalism (Lacey &
Cole, 1993; Somerville & Kolatt, 1999) to obtain the collapsed mass field. As
will be described in the next section, LIMFAST makes use of the Sheth-Tormen
halo mass function (HMF) to connect and derive the properties of galaxies from
halos at each cell, following the prescriptions in Furlanetto (2021). We
stress that in 21cmFAST, this approach corresponds to the “matter density
field” case, where individual halos are not resolved nor identified in the
simulation; at each cell, the halo mass distribution follows the Sheth-Tormen
HMF, and the total number of halos depends on the matter overdensity and the
collapsed fraction in that region.

After the collapsed mass field is obtained, 21cmFAST computes the ionization
state of the IGM by comparing the cumulative number of ionizing photons from
sources and the number of neutral hydrogen atoms and recombinations in
spherical regions, from large to small volumes. This method, first proposed by
Furlanetto et al. (2004), is analogous to the excursion-set formalism and it
allows for scenarios of inhomogeneous reionization. Finally, the code computes
the X-ray and Ly$\alpha$ radiation backgrounds and uses them to derive the
spin and brightness temperature of the 21 cm emission. In the following
sections, we will discuss extensions and variations of these calculations in
LIMFAST, and refer the interested reader to Mesinger & Furlanetto (2007) and
Mesinger et al. (2011) for more details on 21cmFAST.

Figure 1 summarizes the basic ingredients of the LIMFAST code and its
connection to 21cmFAST. Starting from the stellar, gas, and metal masses of a
given halo derived from the Furlanetto (2021) galaxy models, the halo star
formation rate and the metallicity are derived. These last two parameters,
combined with the tabulated stellar and nebular spectra computed with BPASS
and cloudy, then yield the amount of line and ionizing emission from the halo.
The halo luminosities in a given cell are integrated over the (conditional)
HMF computed by 21cmFAST from the density field to generate the intensity
field, with the RSDs being accounted for by the evolved velocity field.
Meanwhile, the ionizing radiation is used to calculate the ionization state of
the IGM.

Figure 1: The basic, modular structure of LIMFAST. Inheriting the baseline
method in 21cmFAST for approximating the formation of LSS and the partitioning
of mass into dark matter halos, LIMFAST starts from assigning the SFR, as well
as stellar, gas and metal masses as a function of halo mass and redshift using
the Furlanetto (2021) galaxy models. These galaxy properties are combined with
tabulated stellar and nebular SEDs and luminosities computed with BPASS and
cloudy, and then yield the amount of line and ionizing continuum emission from
the halos. The halo luminosities are integrated over the HMF computed from the
density field (and offset by the velocity field in redshift-space
calculations) for individual cells to generate the line intensity fields of
interest, whereas the ionizing radiation is used to compute the ionization of
the IGM as done in 21cmFAST.

###  2.1 Galaxy Formation and Evolution Model

We summarize next how the properties of galaxy-hosting halos are derived from
the galaxy formation model introduced by Furlanetto (2021). In Furlanetto
(2021), galaxies evolve due to the interplay between star formation and
feedback. Star formation is fueled by a smooth accretion of mass onto the dark
matter halo, and it is at the same time regulated by the feedback that ejects
gas outside the galaxy via outflows. In this scenario, the co-evolution of
gas, stellar, and metal masses can be described by

$\dot{M}_{\rm g}=f_{\mathrm{b}}\dot{M}-(\mathcal{R}+\eta)\dot{M}_{*}~{},$ (1)
$\dot{M}_{*}=M_{\rm g}/t_{\rm sf}~{},$ (2)

and

$\dot{M}_{Z}=-(1+\eta)Z\dot{M}_{*}+y_{Z}\dot{M}_{*}~{}.$ (3)

In Equation (1) above,
$f_{\mathrm{b}}=\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}\approx 1/6$ is the
baryon fraction, $\dot{M}$ is the mass accretion rate of the halo,
$\mathcal{R}\approx 0.25$ denotes the fraction of mass available for star
formation that resides in stars, and $\eta\gg 1$ accounts for the amount of
mass ejected out of the galaxy by stellar feedback. In Equation (2), $t_{\rm
sf}=t_{\rm orb}/\epsilon$ denotes the star-formation timescale, where $t_{\rm
orb}\sim 18\left[7/(1+z)\right]^{3/2}$ Myr is the orbital timescale, and
$\epsilon$ characterizes the temporal star formation efficiency per orbital
timescale. In Equation (3), $Z\equiv M_{Z}/M_{\rm g}$ denotes the metallicity,
and the metal yield factor $y_{Z}=0.03$ (Benson, 2010) is the fraction of
stellar mass returned to the ISM in the form of metals.

Now writing $M_{\rm g}$ as a fraction of the total accreted mass
$M_{\mathrm{a}}=f_{\mathrm{b}}M$, one can define $X_{\rm g}\equiv M_{\rm
g}/M_{\mathrm{a}}$ as the gas retention factor (Dekel & Mandelker, 2014), and
similarly $X_{*}\equiv M_{*}/M_{\mathrm{a}}$. The previous equations can then
be rewritten in terms of the dimensionless quantity $\tilde{M}\equiv M/M_{0}$,
where $M_{0}$ denotes the halo mass at some initial redshift $z_{0}$, and
taking derivatives with respect to redshift instead of time, namely
$\tilde{M}^{\prime}=d\tilde{M}/dz$. With these substitutions, the system of
differential equations describing the halo, gas, and stellar mass evolution
with redshift finally equates

$\frac{\tilde{M}^{\prime}}{\tilde{M}}=-|\tilde{M}_{0}^{\prime}|~{},$ (4)

$\frac{\tilde{M}_{\rm g}^{\prime}}{\tilde{M}_{\rm
g}}=-|\tilde{M}_{0}^{\prime}|\left[X_{\rm
g}^{-1}-\frac{\epsilon(\mathcal{R}+\eta)}{|\tilde{M}_{0}^{\prime}|C_{\mathrm{orb}}}\left(\frac{1+z_{0}}{1+z}\right)\right]~{},$
(5)

$\tilde{M}_{*}^{\prime}=-\mathcal{R}\tilde{M}_{\rm
g}\left[\frac{\epsilon}{C_{\mathrm{orb}}}\left(\frac{1+z_{0}}{1+z}\right)\right]~{},$
(6)

$\frac{\tilde{M}_{Z}^{\prime}}{\tilde{M}_{Z}}=\mathcal{R}\left(-1-\eta+y_{Z}Z^{-1}\right)\left[\frac{\epsilon}{C_{\mathrm{orb}}}\left(\frac{1+z_{0}}{1+z}\right)\right]~{}.$
(7)

Here, $C_{\rm orb}=(1+z_{0})t_{\rm orb}H(z)$ characterizes the parameter
dependence of the orbital timescale, and $H(z)$ is the Hubble parameter at
redshift $z$.

In order to solve the above equations, the parameters denoting the feedback
model, $\eta$ and $\epsilon$, also need to be specified. For the fiducial
LIMFAST case introduced here, we adopt the momentum-driven feedback model by
Furlanetto et al. (2017), and present detailed dependencies on feedback
prescriptions, as well as other star formation recipes, in Paper II. In the
momentum-driven feedback case, we set $\epsilon=0.015$, and compute the term
denoting the relation between the rate of gas expelled from the galaxy and the
star formation rate as

$\eta(M,z)=C\left(\frac{10^{11.5}}{M}\right)^{\xi}\left(\frac{9}{1+z}\right)^{\sigma}~{},$
(8)

where we have adopted the parameter values $C=5$, $\xi=1/3$, and $\sigma=1/2$,
consistent qualitatively with findings from abundance matching to the observed
galaxy UV luminosity functions (UVLF) at $z\gtrsim 5$ (e.g., Mason et al.,
2015; Sun & Furlanetto, 2016; Mirocha et al., 2017).

For the implementation of this galaxy model in LIMFAST, we have solved the
above equations considering an initial redshift of $z=30$ and tabulated the
results as a function of halo mass and redshift. These tables cover the
redshift range between $z=5$ and $z=30$ in steps of $\Delta z=0.1$, and the
halo masses range from $M=10^{7}\,M_{\odot}$ to $M=10^{16}\,M_{\odot}$ in 900
evenly-distributed logarithmic bins. LIMFAST then interpolates the tables to
obtain the result for any combination of mass and redshift within these
ranges.

Once the above quantities are calculated, we can derive two other important
observables, the star formation rate density (SFRD) and the mass-averaged
metallicity. The SFRD in a simulation cell is obtained by integrating the star
formation rate per halo over the (conditional) HMF of the cell at position
$\boldsymbol{x}$ as

$\dot{\rho}_{*}(\boldsymbol{x},z)=\int{d}n/{d}M(\boldsymbol{x},M,z)\,\dot{M}_{*}(M,z)\,{d}M~{},$
(9)

where ${\rm d}n/{\rm d}M$ denotes the HMF in 21cmFAST, consistent with the
Sheth-Tormen (Sheth & Tormen, 1999) formalism with the correction by Jenkins
et al. (2001). For an individual halo, the metallicity at a given redshift is
simply defined as $Z(M,z)\equiv M_{Z}(M,z)/M_{\rm g}(M,z)$, whereas for a cell
the metallicity is computed as the ratio of total metal mass to total gas mass
in such a cell, namely

$Z(\boldsymbol{x},z)=\frac{\int{d}n/{d}M(\boldsymbol{x},M,z)\,M_{Z}(M,z)\,{d}M}{\int{d}n/{d}M(\boldsymbol{x},M,z)\,M_{\rm
g}(M,z)\,{d}M}~{}.$ (10)

This calculation differs from that of the star formation rate because the
metallicity is not an additive quantity and, therefore, one cannot integrate
the metallicity per halo over the HMF. However, one may want to perform the
metallicity per halo integration, and then divide it by the total number
density of halos in the respective cell, to obtain the average metallicity per
halo in such a volume. The integrals in the above two equations are performed
from a minimum halo mass, $M_{\mathrm{min}}$, set to the atomic cooling halo
mass at a virial temperature of $T_{\rm vir}=10^{4}\,$K, to a maximum halo
mass of $M_{\mathrm{max}}=10^{16}\,M_{\odot}$. The value of $M_{\mathrm{min}}$
varies with redshift from $\sim 10^{7}\,M_{\odot}$ at $z=20$ to $\sim
10^{8}\,M_{\odot}$ at $z=5$. Considering the atomic cooling threshold here is
valid because we do not account for Pop III star formation in this work; one
may want to include smaller halo masses corresponding to molecular cooling
when accounting for that stellar population (see, e.g., Mebane et al., 2018;
Muñoz et al., 2022; Parsons et al., 2022, and references therein for further
discussions). Meanwhile, it should also be considered as an optimistic case
for the formation of Pop II stars, since their formation in low-mass (but
still above this limit) systems may be regulated by radiative feedback due to
reionization (Yue et al., 2014, 2018), although the most recent constraints on
the UVLF faint-end slope from deep lensed galaxies have yet to show the
presence of a turnover (Bouwens et al., 2022). Moreover, while so far only
kinetic feedback from supernova explosions are considered in LIMFAST, we note
that reionization feedback can also leave detectable imprints on intensity
mapping signals (see e.g., Mirocha et al., 2022).

Finally, we define the comoving luminosity density from halos in a cell as

$\rho_{l}(\boldsymbol{x},z)=\int{d}n/{d}M(\boldsymbol{x},M,z)\,l(M,z)\,{d}M~{},$
(11)

where $l(M,z)$ is the luminosity of a halo at a given redshift derived from
the SFR of the halo and the line luminosity per unit SFR (to be detailed in
the next section), and the integration limits are the same as before111We do
not account for scatter in the relation between luminosity and halo mass in
this work. This could be incorporated by considering a distribution of
luminosity values instead of a single value in Equation (11) and a random
sampling of halos from the HMF.. Then, the observed specific intensity is
derived from the luminosity density as

$I_{\nu}(\boldsymbol{x},z)=\frac{c}{4\pi}\frac{\rho_{l}(\boldsymbol{x},z)}{\nu_{0}\,H(z)}~{},$
(12)

where $c$ is the speed of light and $\nu_{0}$ is the rest-frame frequency of
the emission line of interest. In Section 3 and beyond, we will refer to
values of these quantities averaged over the entire co-eval simulation box in
order to assess their redshift evolution, e.g., $\bar{I}_{\nu}(z)\equiv
N^{-1}_{\mathrm{cell}}\sum_{i=0}^{N_{\mathrm{cell}}}I_{\nu}(\boldsymbol{x_{i}},z)$.

###  2.2 Stellar and Nebular SEDs

In this work, we account only for normal, Population II (Pop II) stellar
populations, and leave the inclusion of Pop III stars to future
implementations (see, e.g., Mebane et al., 2018; Qin et al., 2020, 2021;
Tanaka & Hasegawa, 2021; Muñoz et al., 2022; Parsons et al., 2022, for
previous work on implementations of Pop III stars in 21cmFAST).

The original 21cmFAST code considers one single, redshift-independent stellar
spectral energy distribution (SED) describing the Pop II stellar population,
and uses its properties for the calculations of the ionization of the IGM and
the Ly$\alpha$ radiation background. The SED used by 21cmFAST assumes a Scalo
IMF (Scalo, 1998), a metallicity of 0.05$Z_{\odot}$, and continuous star
formation for 100 Myr (Barkana & Loeb, 2005). LIMFAST replaces this SED by a
set of 13 metallicity-dependent stellar SEDs, and uses them to compute the
ionization state of the IGM, the Ly$\alpha$ background, and the nebular line
emission. This parameterization allows us to connect the processes of galaxy
evolution and reionization consistently by using different SEDs that trace the
redshift evolution of metallicity.

The stellar SEDs used by LIMFAST are computed by using BPASS v2.1 (Eldridge et
al., 2017), assuming a single-star and constant star formation mode with an
age of 100 Myr, and a Salpeter IMF (Salpeter, 1955) with stellar masses within
the range $0.5–100\,M_{\odot}$. The 13 SEDs differ from each other by their
metallicity value, and they respectively account for the default absolute
BPASS metallicity values of $Z_{\star}=$ [$10^{-5}$, $10^{-4}$, $10^{-3}$,
$2\times 10^{-3}$, $3\times 10^{-3}$, $4\times 10^{-3}$, $6\times 10^{-3}$,
$8\times 10^{-3}$, $10^{-2}$, $1.4\times 10^{-2}$, $2\times 10^{-2}$, $3\times
10^{-2}$, $4\times 10^{-2}$]. BPASS provides this range of metallicities that
is well suited to reach the low metallicity values that may occur at the
redshifts of reionization, as well as the possible higher values in massive
objects. The number of ionizing photons per stellar baryon in these SEDs spans
the range $N_{\mathrm{ion}}\sim 2500$–6000, where higher numbers are for more
metal-poor SEDs. For comparison, the single value for the number of ionizing
photons used in 21cmFAST is 4361.

For the calculation of the nebular line emission, we use our stellar SEDs as
the incident spectrum in the photoionization code cloudy (version 17.02,
Ferland et al., 2017), and the following quantities describing the nebular
medium; we consider a gas density of $n_{\mathrm{H}}=10^{2}$
$\mathrm{cm}^{-3}$ typical for H II regions (Byler et al., 2017), a distance
between the radiation source and the medium of $r=10^{19}$ cm, and the
ionization parameter values $\log U=[-4,\,-3.5,\,-3,\,-2.5,\,-2,\,-1.5,\,-1]$.
Changing the value of $U$ with the other quantities fixed would imply that the
number of ionizing photons also changes. In practice, however, the number of
photons is fixed by the incident SED, so we renormalize the resulting nebular
emission by the corresponding default number of photons in the SEDs in all
cases. In other words, this approach would be equivalent to vary the value of
$r$ or $n_{\rm H}$ to obtain different ionization parameter values while
keeping a constant number of photons (see similar procedures in Byler et al.,
2017, and Xiao et al. 2018). For each SED case, we assume the gas and the SED
to have the same metallicity (i.e., equal stellar and gas-phase
metallicities). We then perform photoionization calculations combining these
metallicity and ionization parameter values and tabulate the emission results
in units of luminosity per unit of star formation rate. For a given pair of
metallicity and ionization parameter values describing a halo, LIMFAST then
linearly interpolates the tabulated results to derive the luminosity per unit
star formation rate in that halo, $l(M,z)$. In Appendix A, we show the
luminosity of star-formation lines of interest as a function of metallicity
and ionization parameter as the key input to LIMFAST from the spectral
synthesis and photoionization modeling using BPASS and cloudy (see Figure 11).

We do not consider dust in these sub-grid calculations and leave its
implementation to future versions of the code, because the effect of dust
attenuation remains highly uncertain for high-redshift galaxies (Casey et al.,
2014, 2018; Capak et al., 2015; Popping et al., 2017). Not including dust
implies that the default escape fraction of nebular radiation in LIMFAST is
100%, where we have presently also neglected the effect of neutral hydrogen on
the escape and transfer of Ly$\alpha$ emission. This is also the case for the
damping wing absorption of Ly$\alpha$ photons by intervening patches of
neutral IGM along the line of sight (LOS). In other words, the line emission
presented here is the intrinsic one. However, the calculation of the intrinsic
emission allows the user to apply desired custom attenuation effects a
posteriori both at a cell level and along any specific LOS, e.g., applying
analytical extinction prescriptions that depend on redshift or other
parameters (Gong et al., 2017), or treating the absorption and scattering of
Ly$\alpha$ photons by the intergalactic H I (Mesinger & Furlanetto, 2008;
Heneka et al., 2017). Alternatively, one could implement the effects of dust
directly on the SEDs (e.g., through a wavelength-dependent attenuation curve)
resulting in attenuated nebular spectra (Byler et al., 2017), or by linking
the dust properties to the halo metallicities resulting from the LIMFAST
simulations.

Finally, the default calculations assume a constant escape fraction of
ionizing photons of 10% throughout. This value is chosen such that the
resulting reionization history is broadly consistent with current constraints
as further discussed in Section 3. In practice, one can obtain different
reionization histories in LIMFAST by changing the value of the escape
fraction, as well as by varying the star formation parameters, such as
feedback mode, star formation law, or $M_{\mathrm{min}}$ in the above galaxy
model. We explore how these variations may impact the reionization history and
morphology, as well as the corresponding multi-tracer LIM signals in Paper II.

###  2.3 Ionization Calculation

As mentioned in the previous section, 21cmFAST uses a fixed value of 4361
ionizing photons per stellar baryon in the ionization calculations at all
redshifts. Following our approach of metallicity-dependent SEDs, LIMFAST
instead computes a number of ionizing photons that depends and evolves with
metallicity. In detail, LIMFAST computes the number of ionizing photons for a
halo of a given metallicity by linearly interpolating the number of ionizing
photons from the two SEDs with metallicities closer to that of the halo.
Therefore, because the metallicity of the halos changes with time, the number
of ionizing photons also evolves with redshift. Finally, we have not
considered radiation transfer effects that may produce further spatial
variations in the ionization calculations (see, e.g., Davies & Furlanetto 2022
for a recent discussion on the propagation of ionizing radiation and its
effects on 21cmFAST, and Lewis et al. 2022).

###  2.4 IGM Ly$\alpha$ Emission

The IGM Ly$\alpha$ emission denotes here the Ly$\alpha$ radiation produced in
situ in the IGM due to the recombination of ionized gas leading to the
formation of H I. We ignore collisional effects that may also lead to the
production of Ly$\alpha$ because these are expected to be subdominant compared
to recombination (see further discussions in Silva et al. 2013 and Comaschi &
Ferrara 2016).

The comoving Ly$\alpha$ luminosity density from the recombination of IGM gas
in a cell can be expressed as

$\rho_{\rm Ly\alpha,\,rec}^{\rm IGM}(z)=f_{\rm rec}\,E_{\rm
Ly\alpha}\,\dot{N}_{\rm rec}(z)\,n_{\rm H\textsc{ii}}(z)~{},$ (13)

where $f_{\rm rec}=0.66$ denotes the fraction of recombinations producing
Ly$\alpha$ photons, $E_{\rm Ly\alpha}=1.637\times 10^{-11}\,{\rm erg}$ is the
energy of the Ly$\alpha$ transition, and $\dot{N}_{\rm rec}(z)$ represents the
recombination rate per baryon computed by the original 21cmFAST code. The last
term above equates

$n_{\rm H\textsc{ii}}(z)=[1-x_{\rm H\textsc{i}}(z)]\,n_{\rm
b}\,[1+\delta(z)]~{},$ (14)

and it describes the comoving number density of ions in the IGM. Here, $x_{\rm
H\textsc{i}}(z)$ denotes the neutral hydrogen fraction, $\delta(z)$
corresponds to the matter overdensity, and $n_{\rm b}$ is the present day
comoving number density of baryons. The comoving luminosity density can be
finally converted to the observed intensity by means of Equation (12) as
above.

Although the IGM recombination process may produce diffuse emission of other
hydrogen lines, we ignore their contribution. Besides Ly$\alpha$, H$\alpha$
would be the brightest of these hydrogen lines, but given its recombination
coefficient and the transition probabilities between the atomic energy level
of the hydrogen atom, the expected luminosities for H$\alpha$ are expected to
be around one order of magnitude fainter than those of Ly$\alpha$.

###  2.5 Ly$\alpha$ Background

Another component contributing to the Ly$\alpha$ radiation field is that
produced by the scattering of high energy UV photons that, while traveling
through the IGM, redshift into the frequency of the Lyman series lines in the
rest frame of the IGM gas. When these initially high energy photons reach the
frequencies of the Lyman lines, they are susceptible to be absorbed by the
neutral H I and they can subsequently lead to the processes of resonant
scattering or a down-cascade by the electron in the atom that may ultimately
produce Ly$\alpha$ radiation.

We follow the calculation of the Ly$\alpha$ background in 21cmFAST, but we use
our set of SEDs consistently instead of a single SED. In detail, each of the
13 stellar SEDs is tabulated to account for the number of photons in between
the first 23 energy levels of the hydrogen atom as in 21cmFAST (see details of
this calculation in Barkana & Loeb 2005 and Mesinger et al. 2011). Then,
interpolation is used to find the photon number corresponding to metallicities
within those of the two nearest SEDs as done for the line and ionizing
emission.

With the Ly$\alpha$ background in the frame of the gas in units of photon rate
per unit frequency, area and steradian, $J_{\alpha}$, the observed intensity
can be obtained as (Silva et al., 2013)

$I^{\mathrm{BG}}_{\mathrm{Ly\alpha}}=\frac{6\,E_{\rm
Ly\alpha}}{(1+z)^{4}}J_{\alpha}~{},$ (15)

where the term in the denominator accounts for the redshift dimming of the
surface brightness.

###  2.6 21 cm Signal

The 21 cm signal sourced by the neutral hydrogen in the IGM is often expressed
as a differential brightness temperature that can take a positive or negative
sign depending on whether the gas is emitting or absorbing 21 cm radiation.
This differential brightness temperature can be written as (Furlanetto et al.,
2006)

$\displaystyle\delta T_{b}(z)$ $\displaystyle\approx 27\,x_{\rm
H\textsc{i}}(z)\,[1+\delta(z)]\,{\left|1+\frac{1+z}{H(z)}\frac{dv_{\parallel}}{dr_{\parallel}}\right|}^{-1}$
$\displaystyle\times\left[1-\frac{T_{\gamma}(z)}{T_{S}(z)}\right]\left[\frac{1+z}{10}\frac{0.15}{\Omega_{m}h^{2}}\right]\left[\frac{\Omega_{b}h^{2}}{0.023}\right]\,{\rm
mK}~{}.$ (16)

Here, $T_{\gamma}$ and $T_{S}$ denote the temperature of the cosmic microwave
background (CMB) radiation and the spin temperature of the intergalactic
neutral hydrogen, respectively. The term $dv_{\parallel}/dr_{\parallel}$
represents the peculiar velocity gradient of the gas along the LOS, and it is
responsible for introducing the RSD effect as will be detailed in the next
section.

The spin temperature term in Equation (2.6) traces the CMB temperature and the
thermal history of the IGM via different coupling mechanisms. These include
the collisional (thermal) coupling of the intergalactic gas, as well as
radiative coupling effects through CMB and Ly$\alpha$ photons. The coupling
through resonant scattering of Ly$\alpha$ photons is known as the Wouthuysen-
Field (WF) effect, and it is sourced by the cosmic soft-UV Ly$\alpha$
background described in the previous section. Taking into account these
processes, we can express the offset $1-T_{\gamma}/T_{S}$ as

$1-\frac{T_{\gamma}}{T_{S}}=\frac{x_{c}+x_{\alpha}}{1+x_{c}+x_{\alpha}}\left(1-\frac{T_{\gamma}}{T_{K}}\right)~{},$
(17)

where $T_{K}$ is the kinetic temperature of the intergalactic gas, and $x_{c}$
and $x_{\alpha}$ are the collisional and WF radiative coupling coefficients,
respectively. While the Ly$\alpha$ background is the responsible for the WF
coupling as mentioned above, the heating of the intergalactic gas is mostly
dominated by the X-ray background.

For the calculation of the 21 cm signal, LIMFAST follows the methodology
developed in 21cmFAST, except for the Ly$\alpha$ background derivation. As
mentioned in the previous section, for the Ly$\alpha$ background we adopt our
metallicity-dependent SEDs, which results in a varying Ly$\alpha$ background
due to the evolution of metallicity, in addition to that of star formation,
with redshift. For the X-ray background, we use the prescriptions relating
luminosity and star formation adopted by Park et al. (2019). Finally, for both
radiation background calculations, we have modified the original code to take
into account our star formation formalism.

###  2.7 Redshift-Space Distortions

The implementation of RSD calculations in LIMFAST generally follows that of
21cmFAST, as described in Mesinger et al. (2011) and Greig & Mesinger (2018),
which includes the numerical approach to transform the data from real to
redshift space outlined by Jensen et al. (2013). This method divides each
simulation cell into sub-cells of equal intensity and computes their new
position according to the peculiar velocities. Then, the sub-cells are re-grid
to the original resolution (see Section 3.1.2 in Jensen et al., 2013).
Finally, for the calculation of intensities with RSD, we use the default
21cmFAST cut-off velocity gradient value of
${dv_{\parallel}}/{dr_{\parallel}}=0.2\,H(z)$. This limit is set to avoid
extremely large intensities when the velocity gradient presents large values.
The exact value for this cut-off is somewhat arbitrary, but Mao et al. (2012)
found that values around the default one used here do not produce extreme
departures from the true results (their Section 9), although this depends on
the specific simulation and case. We have verified that using
${dv_{\parallel}}/{dr_{\parallel}}=0.1\,H(z)$ or
${dv_{\parallel}}/{dr_{\parallel}}=0.3\,H(z)$ does not significantly change
our results. This cut-off, however, is not required and, therefore, not used
in the 21 cm calculations. For the 21 cm line, the full radiative transfer
derivation ignoring the optically thin assumption results in terms that vanish
instead of going to infinity (see Section 5.1.1 in Mao et al., 2012). In
principle, a similar full derivation would be possible for other lines, but
this requires accounting for the specific radiative processes affecting each
line of interest, which goes beyond the scope of LIMFAST in its current form.
We leave these calculations to future work and here simply use the cut-off
method for the star-formation lines. For interested readers, a more detailed
summary of the RSD-related calculations in LIMFAST and effects on the LIM
signals of interest is provided in Appendix B.

##  3 Results

We present in the remainder of this section the results from LIMFAST runs
considering simulation boxes of 0.5 Gpc on a side, with a cell size of
$2.5^{3}\,\mathrm{Mpc}^{3}$, from $z=15$ to $z=5$. In this paper, we focus on
the 21 cm line tracing the neutral IGM and various optical/UV lines generally
emitted from H II regions created by the formation of massive stars in
galaxies. These star-formation lines are some of the strongest emission
features in the optical/UV spectrum of a typical star-forming galaxy,
including Ly$\alpha$, H$\alpha$, H$\beta$, [O II] $3727\,{\rm\AA}$, and [O
III] $5007\,{\rm\AA}$. It is noteworthy that the intensity field of Ly$\alpha$
(and in principle other hydrogen recombination lines like H$\alpha$, though to
a much less extent) is not only sourced by recombinations of ionized gas
inside galaxies; it also contains emission from the diffused, ionized IGM and
the redshifted background radiation associated with sources within the
“horizons” (see e.g., Furlanetto et al., 2006; Holzbauer & Furlanetto, 2012).
These additional components are self-consistently calculated and generated
separately. In Paper II, we present extensions that model emission lines from
the neutral ISM, such as the [C II] $158$ $\mu$m and CO rotational lines. For
predictions displayed in this section, we present the fiducial case where we
adopt a fixed value of the ionization parameter $\log U=-2$ consistent with
the enhanced $U$ values typically observed for high-$z$ galaxies (Sanders et
al., 2016; Strom et al., 2017; Katz et al., 2023). A plausible variation that
assumes metallicity-dependent $U$ is also considered. We will revisit these
assumptions later again in Section 4.3 and discuss their implications for
future LIM observations. In Section 3.1, we present the results for the
redshift evolution of the variety of quantities of interest, such as the SFRD,
metallicity, IGM neutral fraction, and mean line intensities. In Section 3.2,
we show power spectra derived from the simulated line intensity fields.
Results from the literature for are also displayed for comparison, and we
discuss the origin of several key differences between our results and previous
studies in Section 4.

###  3.1 Redshift Evolution of the Cosmic Means

####  3.1.1 Star Formation Rate Density

Figure 2: The evolution of the cosmic star formation rate density computed by
LIMFAST with the fiducial model described in previous sections. For
comparison, we also plot similar results from the literature (Silva et al.,
2013; Park et al., 2019; Muñoz et al., 2022; Kannan et al., 2022a) as colored
lines, together with observational data from Bouwens et al. (2021) after
extrapolating the luminosity functions to a limiting magnitude of $M_{\rm
UV}=-10$. The error bars represent the 16th and 84th percentiles derived from
the Schechter parameters reported. The dashed lines show the fractions of
cosmic star formation associated with massive halos with
$M>10^{10}\,M_{\odot}$.

By averaging the SFRDs of individual cells computed as explained in Section
2.1 over the coeval box at different redshifts, we can determine the mean
cosmic SFRD evolution. In Figure 2, we show results from LIMFAST, assuming the
fiducial model and default setups specified in Section 2 (black solid line).
The colored lines illustrate a compilation of predictions from the literature
(Silva et al., 2013; Park et al., 2019; Muñoz et al., 2022; Kannan et al.,
2022a), which typically involve different assumptions of $M_{\mathrm{min}}$
that can host star-forming galaxies. Unless otherwise labeled, the atomic
cooling limit is assumed, as is the case for LIMFAST (see Section 2.1). These
theoretical predictions are compared with the data points that represent SFRDs
implied by the observed UVLFs (Bouwens et al., 2021). Note that to make fair
comparisons to our results, we use the Schechter parameters (including their
variances but not covariances) reported in Bouwens et al. (2021) to
extrapolate and integrate the UVLF measurements down to a magnitude of $M_{\rm
UV}=-10$, which roughly corresponds to a halo mass at the atomic cooling limit
in the redshift range of interest.

Overall, the SFRD evolution from LIMFAST agrees well with the latest
observations after plausible extrapolations, although these constraints become
rather weak beyond $z\sim 7$. Our SFRD is generally comparable (within a
factor of a few) to models from the literature at $5\lesssim z\lesssim 8$,
with the difference in the shape and amplitude growing larger towards higher
redshift. As will be further discussed in Section 4.1, the overall higher
amplitude and shallower slope are mainly because our model allows star
formation to occur in low-mass halos (i.e., faint galaxies) down to the atomic
cooling limit with a non-negligible efficiency. To demonstrate this, we also
include in Figure 2 cases showing the SFRD from only the massive halos with
$M>10^{10}\,M_{\odot}$.

####  3.1.2 Metallicity

Similar to the SFRD, we can compute the mean cosmic metal enrichment history
by averaging the metallicity per cell over coeval boxes. The solid black line
in Figure 3 shows the mean metallicity evolution of the collapsed structures,
following Equation (10). Specifically, these metallicity values represent the
mean metallicities of the gas (and stars) residing in galaxies in the
simulation cells at different redshifts, but they should not be confused with
the mean halo metallicity. For comparison, a number of metallicity estimates
of $z>5$ star-forming galaxies and gamma-ray burst (GRB) hosts available to
date, together with the metal enrichment history and its extrapolation to
$z>7$ implied by the cosmic star formation history (Madau & Fragos, 2017), are
also shown. The galaxy metallicities are determined from recent analysis of
JWST/NIRSpec (Taylor et al., 2022) using the direct $T_{e}$ method, whereas
the GRB metallicities are measured from absorption spectra of their afterglows
(Thöne et al., 2013; Bolmer et al., 2019). Because GRBs may preferentially
reside and thus represent the ISM of low-mass (i.e., metal-poor) galaxies
(e.g., Graham & Fruchter, 2017; Palmerio et al., 2019, and references
therein), we also show the metallicity evolution considering only halos with
$M<10^{10}\,M_{\odot}$. Given the small sample size and large uncertainties
associated with the metallicity estimates, especially for GRBs where the
analysis is highly sensitive to properties of their host halos (see, e.g.,
Metha et al., 2021), we find broad agreement between these observational
constraints/implications and our model predictions.

Figure 3: Metallicity evolution of the collapsed objects in LIMFAST, as
derived from Equation (10) considering all halos more massive than
$M_{\mathrm{min}}$ (black solid curve) or those with
$M_{\mathrm{min}}<M<10^{10}\,M_{\odot}$ (black dashed curve). The data points
represent the metallicity estimates of five $z>5$ galaxies measured from
JWST/NIRSpec data using the direct method (Taylor et al., 2022), as well as
those of the two highest redshift GRB hosts to date with confident,
absorption-based metal detections (Thöne et al., 2013; Bolmer et al., 2019).
The blue curve shows a best-fit metallicity evolution of the observed galaxy
population as a whole (extrapolated to $z>7$) adopted from Madau & Fragos
(2017).

####  3.1.3 Neutral Fraction

The redshift evolution of the IGM neutral fraction is a key diagnostic for our
simulations of the reionization process. Figure 4 shows the evolution of the
gas neutral fraction as computed by LIMFAST, assuming a LyC escape fraction of
10% that yields a CMB optical depth $\tau_{\mathrm{CMB}}\approx 0.06$
consistent with recent estimates (see also Paper II for a comparison of the
reionization scenarios and corresponding optical depths from plausible
variations of the stellar feedback). The colored data points represent
observational constraints from quasar absorption spectra (McGreer et al.,
2015; Greig et al., 2017; Davies et al., 2018; Greig et al., 2019; Zhu et al.,
2022) and Ly$\alpha$ transmission (Hoag et al., 2019; Mason et al., 2019;
Whitler et al., 2020), whereas the colored lines denote some alternative
models from previous 21cmFAST-based simulations (Park et al., 2019; Muñoz et
al., 2022). We note that some recent observations have provided strong
evidence for a late completion of reionization by $z\lesssim 5.5$ (e.g., Zhu
et al., 2022), which poses challenges on existing reionization models. As our
fiducial model still agrees with the updated upper limits at $z<6$, we do not
attempt to match them and consider scenarios producing $\tau_{\mathrm{CMB}}$
values consistent with CMB observations to be generally valid.

Figure 4: Evolution of the volume-averaged IGM neutral fraction
$\bar{x}_{\mathrm{HI}}$ as computed by LIMFAST with a 10% escape fraction of
ionizing photons (black solid curve). For comparison, the colored data points
show constraints from the literature based on different methods, including
quasar absorption spectra McGreer et al. (2015); Greig et al. (2017, 2019);
Davies et al. (2018); Zhu et al. (2022) and Ly$\alpha$ emission (Hoag et al.,
2019; Mason et al., 2019; Whitler et al., 2020). Furthermore, the dashed
curves represent $\bar{x}_{\mathrm{HI}}$ evolution predicted by other previous
21cmFAST-based simulations (Park et al., 2019; Muñoz et al., 2022) that also
broadly match observational constraints available to date. Figure 5: LIMFAST
light cones over redshift $5\leq z\leq 15$. The 8 panels shown correspond to
the hydrogen neutral fraction, the brightness temperature of the 21 cm line,
the SFRD, the metallicity of collapsed structures, and the intensities of
Ly$\alpha$ (from star formation, the ionized IGM, and the scattering
background, respectively) and [O III] $5007\,{\rm\AA}$ emission. The color
schemes are identically normalized (to the maximum value) for [O III] and
Ly$\alpha$ from star formation to highlight their difference.

####  3.1.4 Line Emission

To showcase the output of LIMFAST, we show in Figure 5 the LIMFAST light cones
for various quantities. As in 21cmFAST, these light cones are created by first
generating coeval boxes at discrete redshift steps over the redshift range of
interest and then interpolating (after being rotated to avoid repeating
structures) between redshifts along the LOS to describe the continuous
evolution of each quantity. From the top to the bottom row, the signals
displayed are the neutral fraction of hydrogen, the 21 cm brightness
temperature, the SFRD, the metallicity of collapsed structures, Ly$\alpha$
intensities from star formation, recombinations in the diffuse ionized IGM,
and the scattering of UV photons in the IGM, as well as the [O III]
$5007\,{\rm\AA}$ line intensity, respectively.

Figure 5 visualizes qualitatively the correlation between different
quantities. For instance, beyond the well-studied correlation of the 21 cm
signal with the hydrogen neutral fraction, it is easy to see the increase of
metallicity (metal enrichment) that corresponds to the SFRD evolution, which
is directly traced by the evolving Ly$\alpha$ emission from star-forming
galaxies. The Ly$\alpha$ emission from the diffused ionized IGM, on the other
hand, shows a stronger dependence on the ionization state of the IGM than the
amount of cosmic star formation, and its small-scale spatial fluctuations
appear to be weaker compared to the star-formation component of Ly$\alpha$
emission. Small-scale fluctuations are even weaker for the Ly$\alpha$
background radiation from the scattering of UV photons, due to the large
Ly$\alpha$ “horizon distance” (a few hundred comoving Mpc at $z\sim 10$)
associated with the distance beyond which UV photons can no longer redshift
into the Ly$n$ series and produce Ly$\alpha$ radiation (Holzbauer &
Furlanetto, 2012). Finally, compared with hydrogen lines, metal lines like [O
III] $5007\,{\rm\AA}$ are much more sensitive to the metallicity, which
controls the hardness of the stellar ionizing spectrum. As a result, the non-
trivial metallicity dependence leads to a stronger redshift evolution for the
[O III] line compared with Ly$\alpha$, as shown by the steeper color gradient.

A more quantitative view of the redshift evolution of these LIM signals is
provided in Figure 6, where we also show results from the literature for
comparison. In addition to our fiducial model, for the Balmer and oxygen lines
we consider two alternative cases that exemplify how variations of the
photoionization modeling affect the predicted mean intensity evolution of
different lines. As denoted by the gray dashed lines, one of such alternative
cases is to adopt a local, redshift-independent relation between the line
luminosity and the SFR, which otherwise evolves with redshift through the
metallicity dependence in our fiducial model. Here, we take scaling relations
similar to those in Gong et al. (2017), which are calibrated against galaxies
observed in the local universe. Specifically, considering the scaling relation
$L_{\mathrm{line}}=C\,\dot{M}_{*}$, we have $C=1.3\times
10^{41}\,\mathrm{erg\,s^{-1}}\,\dot{M}_{\star}^{-1}$ for H$\alpha$ (Kennicutt,
1998), $C=4.4\times 10^{40}\,\mathrm{erg\,s^{-1}}\,\dot{M}_{\star}^{-1}$ for
H$\beta$ from the relation ${\rm H}\alpha/{\rm H}\beta=2.86$, $C=7.2\times
10^{40}\,\mathrm{erg\,s^{-1}}\,\dot{M}_{\star}^{-1}$ for [O II] from the
relation [O II]$/{\rm H}\alpha=0.57$ (Kennicutt, 1998), and $C=1.3\times
10^{41}\,\mathrm{erg\,s^{-1}}\,\dot{M}_{\star}^{-1}$ for [O III] (Ly et al.,
2007). Besides the clear impact on the overall normalization of the redshift
evolution, assuming locally-calibrated scaling relations also results in a
slightly shallower slope, especially for the metal lines, because in this case
the ionizing spectrum does not harden as the metallicity decreases (towards
higher redshift).

The other alternative case, as denoted by the gray dotted lines, assumes that
the ionization parameter $U$ anti-correlates with metallicity following
$U\propto Z^{-0.8}$, instead of assigning a fixed value like $\log U=-2$ as in
the fiducial case. Such an anti-correlation is broadly in line with the ISM
physical condition of high-$z$ galaxies probed by nebular emission line ratios
(Sanders et al., 2016) and the expectation of simple analytic models of wind-
driven bubble (Dopita et al., 2006). Due to the way the oxygen line ratio
$\mathrm{O32=[O\,III]/[O\,II]}$ depends on $U$, the assumption of $U\propto
Z^{-0.8}$ leads to a steeper (shallower) redshift evolution of
$\bar{I}_{\mathrm{O\,II}}$ ($\bar{I}_{\mathrm{O\,III}}$), though the
difference from the fiducial case with $\log U=-2$ is modest. While neither a
universal, fixed $U$ nor the exact anti-correlation considered may be true in
reality, they illustrate how much LIM predictions might be affected by the
detailed photoionization modeling of H II regions associated with the
formation of massive stars. Without a good understanding of the chemical
abundance and ionization conditions of H II regions in high-$z$ galaxies, it
would therefore be difficult to gauge how reliably LIM measurements of lines
like [O II] and [O III] can be applied to constraining the cosmic SFRD at high
redshift (Gong et al., 2017; Sun, 2022). We will further elaborate on this and
its implications for future LIM observations in Section 4.3.

Qualitatively, the redshift evolution of the various lines extracted from our
simulations share a few similarities with results in the literature. For line
intensities that directly trace the ongoing cosmic SFRD (including the
accumulated Ly$\alpha$ background that encodes the past star formation
history), the steepness of the redshift evolution is comparable to other works
(e.g., Comaschi & Ferrara 2016 and THESAN-1) that allow stars to form
efficiently in low-mass halos, in such a way that consistently predicts the
faint-end slope of the UVLFs observed. The calibration against the UVLF,
typically done at at $z\sim 5$–6, also leads to reasonable agreement in the
amplitude at these redshifts. Quantitatively, discrepancies ranging from
orders of a few at $z\sim 5$ to several hundred at $z>10$ exist, as a result
of different assumptions made for galaxy formation and additional effects like
dust attenuation. Impacts of systematic differences caused by the higher and
slowly evolving SFRD and the exclusion of attenuation corrections, as well as
complications due to the photoionization modeling that lead to deviations from
the local relations (e.g., Silva et al., 2017; Kewley et al., 2019), will be
further discussed in Section 4.2.

Figure 6: Redshift evolution of the brightness of emission lines computed by
LIMFAST (solid black curves), compared with results from the literature shown
by the colored curves (Pullen et al., 2014; Comaschi & Ferrara, 2016; Mesinger
et al., 2016; Heneka et al., 2017; Kannan et al., 2022b). Two alternative
cases assuming local scaling relations and the anti-correlation $U\propto
Z^{-0.8}$ for the Balmer and oxygen lines are shown by the dashed and dotted
curves in gray, respectively.

###  3.2 Power Spectra

Figure 7: Power spectra of emission lines at $z=7$, 10, and 12, derived from
the fiducial LIMFAST simulations (solid curves), together with two alternative
cases assuming local scaling relations (dashed curves) and the anti-
correlation $U\propto Z^{-0.8}$ (dotted curves) for the Balmer and oxygen
lines. The RSD effect on the redshift-space power spectrum is displayed for
the 21 cm line only by the dash-dotted curves.

Figure 7 shows the spherically-averaged power spectra of intrinsic, dust-
unattenuated line emission at redshifts $z=7$, $z=10$, $z=12$, in the form of
$kP(k)$. Results assuming local luminosity–SFR relations ($U\propto Z^{-0.8}$)
are again shown by the dashed (dotted) curves for comparison. At all
redshifts, clear differences exist between the shapes of power spectra for
lines tracing star formation and the IGM. In particular, Ly$\alpha$ emission
from the IGM shows a clear deficit of fluctuation power on small scales
($k\gtrsim 0.5\,h/\mathrm{Mpc}$), especially for the scattering background
radiation which encodes the large Ly$\alpha$ “horizon distance” that smooths
out small-scale fluctuations (note the many orders of magnitude shown in the
$y$-axis of the Ly$\alpha_{\mathrm{BG}}$ panel; see also Figure 5). The
emission of Ly$\alpha$ (as well as other star-formation lines) from star-
forming galaxies, on the other hand, has an overall flatter power spectrum,
whose shape reflects how biased the line is as a LSS tracer. Compared with
hydrogen lines, oxygen lines show stronger shape and amplitude evolution over
$7<z<12$ due to the strong metallicity dependence of their luminosities. The
(different) ways [O II] and [O III] depend on $U$ can also impact their power
spectra differently, which may be used to probe ionization conditions of the
ISM (see Section 4.3). On small scales ($k\gtrsim 1\,h/\mathrm{Mpc}$) and high
redshifts ($z\gtrsim 10$), in particular, the shapes of [O II] and [O III]
power spectra deviate, with [O II] being steeper and [O III] being more
flattened out, which can be seen from the comparison with the case assuming
the local scaling relation where the line luminosity is simply proportional to
the SFR (dashed lines). Such a systematic difference in the power spectrum
shape can be attributed to the fact that [O II] has a steeper luminosity–halo
mass relation than [O III], which gives rise to the (relatively) higher large-
scale clustering power for [O II].

Insofar as the model variations are considered, we see that orders of
magnitude difference in the power spectrum amplitude can easily arise from
seemingly modest changes of model assumptions like applying locally-calibrated
scaling relations to high redshifts. Understanding such model dependence has
important implications for making sensitivity forecasts and survey strategies
of future LIM experiments, especially pathfinder programs that are usually
instrument-noise-limited. A few additional examples of model variations and
observational implications will be discussed in Section 4.3.

Finally, we note that considering the RSDs generally makes the spherically-
averaged power spectra in real space and in redshift space only modestly
different. We therefore only show a comparison for the 21 cm signals in Figure
7 (dash-dotted curves) for visual clarity. The net effect is a small increase
in the fluctuation power on all scales, with the only exception of the 21 cm
signal which instead shows a small decrease in power at $z=7$ due to the
growing anti-correlation between the density and ionization fields. To better
elucidate effects of the RSDs, we examine in Appendix B the 2D power spectra
of a few representative line signals at $z=7$ and 11 (see Figure 12).
Difference between real-space and redshift-space 2D power spectra is only
discernible for the 21 cm and Ly$\alpha$ emission from the ionized IGM
(without considering Ly$\alpha$ scattering effects), and is largely negligible
for line signals tracing star formation that occurs in more overdense and thus
more strongly biased regions. Interested readers are referred to Appendix B
for a pedagogical overview of the implementation and impact of the RSDs in
LIMFAST.

##  4 Discussion

Below, we highlight the main features that distinguish LIMFAST from, and thus
put our simulations into the context of, previous work from the literature,
including the parent code 21cmFAST. Our goal here is to focus on model
assumptions giving rise to qualitative differences that can be observationally
tested (at least in principle) by various LIM signals, rather than
scrutinizing quantitatively the likely causes of specific discrepancies
between LIMFAST and the literature. Along these lines, we also discuss a few
planned applications of LIMFAST in studies of the high-redshift universe
exploiting the intensity mapping concept.

###  4.1 Improvements with Respect to 21cmFAST

Besides an extension of the original 21cmFAST code to the more general regime
of LIM, LIMFAST also has a few noteworthy improvements on simulations of the
metagalactic and 21 cm radiation backgrounds, including the implementation of
(1) a more physically-grounded galaxy formation model, and (2) redshift-
dependent galaxy SEDs consistent with the cosmic chemical enrichment history
implied by the galaxy formation model. While these changes made do not hugely
impact the predicted 21 cm signals, understanding possible systematic effects
caused by them is critical for both high-accuracy forward modeling and
constraining astrophysics and cosmology with the simulated observables at
cosmic dawn and reionization eras.

In terms of the galaxy formation model, we note that in 21cmFAST, for fast
simulation of the ionization and spin temperature fields as key inputs to the
21 cm signal, the rate of star formation in dark matter halos is simply
described by an efficiency parameter $f_{*}$ as a power-law of the halo mass
and a star formation timescale given by $t_{*}/H(z)$, where $t_{*}$ takes a
fixed value between 0 and 1. While these assumptions give a reasonable
description of the co-evolution of galaxies and the IGM during the EoR, in
detail some non-trivial discrepancies may occur, especially when the mass
growth histories of individual halos are considered.

To illustrate this, in Figure 8 we compare the SFR as a function of halo mass
and redshift given by the galaxy formation models adopted in LIMFAST and
21cmFAST, respectively. Note that for the model parameters in 21cmFAST, we
take the best-fit values in Table 2 of Park et al. (2019) derived from
observations of the UVLF, CMB optical depth, and quasar absorption spectra,
instead of fiducial values reported there, which overall predict 10% lower
SFRs because they are not calibrated against observations. From the
comparison, it is clear that even though both ways of parameterizing the SFR
fit the observed UVLFs at $5\lesssim z\lesssim 8$ reasonably well, the
fiducial model adopted in LIMFAST implies a weaker mass dependence and a
stronger redshift dependence of the SFR–halo mass relation, thus resulting in
a more gradual SFRD evolution when integrated over the halo mass function as
shown in Figure 2. The difference in the mass dependence is mainly attributed
to the different ways feedback on star formation in low-mass halos is
considered in LIMFAST and 21cmFAST. In LIMFAST, the SFR–halo mass relation is
derived physically from star formation sustained by continuous accretion of
baryonic matter onto dark matter halos and regulated by stellar feedback, of
which the implied observational consequences are then checked against existing
data. In 21cmFAST, the mass dependence and timescale of star formation is
tuned to match the observational data without a clear physical interpretation
in the context of feedback. Consequently, discrepant SFR–halo mass relations
arise and lead to the difference in the SFRD that grows larger when
extrapolated to higher redshifts. Different halo mass dependence of the SFR
also leads to testable differences in the history and morphology of the
reionization process tracked by statistics such as the bubble size
distribution (BSD), as demonstrated in Paper II. These differences underscore
the impact of model dependence in simulating the onset of galaxy formation and
its impact on the early phase of the EoR, especially in semi-numerical
simulations where simple galaxy formation models are employed.

Figure 8: Top: comparison of the halo mass and redshift dependence of the SFR
in LIMFAST (solid) and 21cmFAST (dashed), with the latter assuming best-fit
parameters from Park et al. (2019) constrained by the UVLF, CMB, and quasar
absorption spectra. Bottom: the fractional difference in the SFR from galaxy
formation models adopted by LIMFAST and 21cmFAST. The horizontal line marks
where the difference is 50% of the SFR predicted by the Park et al. (2019)
model.

The updated galaxy formation model implemented in LIMFAST also makes possible
simulations of the metagalactic background consistent with histories of star
formation and metal enrichment. Particularly, we can trace how SEDs of the
radiation sources evolve over time as galaxies build up their stellar mass and
become more metal-enriched. As described in Section 2, the fixed-value SED
template used in 21cmFAST is replaced by the metallicity-dependent SED
templates coupled with the metal production and enrichment history from our
galaxy formation model. The impact of such redshift evolution of the source
SEDs on the 21 cm signal can then be examined.

Figure 9 compares the 21 cm sky-averaged signal and the power spectrum as a
function of redshift at a fixed scale ($k=0.2\,h/\mathrm{Mpc}$), predicted by
our simulations with and without invoking metallicity-dependent SEDs. While
the overall impact from metallicity-driven SED evolution is modest, clear
differences in the timings of Ly$\alpha$ coupling ($z\sim 20$) and
reionization ($z\sim 8$) can be seen, whereas the timing of the epoch of
heating (EoH) driven by X-ray binaries ($z\sim 15$) is little affected because
of the metallicity-independent X-ray luminosity assumed. As a result, the
asymmetries of the 21 cm sky-averaged signal is modified. These observations
are consistent with how we expect the SED evolution with metallicity to
modulate the 21 cm signals — the overall harder spectrum resulted from an
evolving, metallicity-dependent SED shifts the Ly$\alpha$ coupling and
reionization epochs to higher redshifts. It also slightly reduces the contrast
of fluctuations on large versus small scales and speeds up its evolution
during the EoR. These effects can potentially introduce significant systematic
uncertainties to efforts on using the statistics of the 21 cm signal to probe
astrophysics and cosmology, such as the presence of massive Pop III stars
(e.g., Mirocha et al., 2018; Mebane et al., 2020), effects of cosmic rays
(e.g., Bera et al., 2022), and alternative dark matter models (e.g., Jones et
al., 2021).

Finally, we note that any metallicity dependence of the luminosity of X-ray
binaries can similarly modulate the 21 cm signal during the EoH, whose
implementation in LIMFAST is left for future work. In a recent study, Kaur et
al. (2022) find that different plausible prescriptions for the metallicity
dependence of the X-ray luminosity correspond to order of unity changes in the
21 cm signals during the EoH. It is therefore interesting to understand, as a
future step, how exactly the combined effect of metal enrichment on galaxy UV
and X-ray spectra modifies the thermal and ionization histories of the IGM to
be revealed by 21 cm and other complementary observations.

Figure 9: Effects of evolving galaxy SEDs on the 21 cm sky-averaged signal
(top) and fluctuations (bottom). The same underlying galaxy formation model is
assumed for both cases with and without invoking metallicity-dependent SEDs.
Also shown in both panels is the relative difference defined as the fractional
deviation of the red curve from the gray curve.

###  4.2 Comparison to Other Previous Work

In Section 3, we presented results from LIMFAST fiducial model and some simple
variations, juxtaposed with a number of forecasts from previous LIM studies in
the literature, without providing a careful explanation for the differences
observed. Here, we outline a few factors that distinguish our results from the
others, and identify the specific differences that they are responsible for.

####  4.2.1 High-Redshift Galaxy Formation

As already pointed out in Sections 3.1.1 and 4.1, the most significant
difference in terms of high-$z$ galaxy formation between LIMFAST and other
studies is that we predict overall higher and less steep SFRD as a function of
redshift, especially at $z>8$ (see Figure 2). Two major factors embedded in
our galaxy formation model together make the fractional contribution from
faint galaxies to the total SFRD higher and thus are responsible for this
discrepancy — a small minimum halo mass, $M_{\mathrm{min}}$, for star
formation that equals the atomic cooling threshold (assuming
$T_{\mathrm{vir}}=10^{4}\,$K) and a relatively shallow, power-law relation
between the amount of star formation and halo mass above that. Both factors
modulate the SFRD evolution through the integral over the HMF, though
$M_{\mathrm{min}}$ is relatively less important at the low-$z$ end considered,
when a considerable fraction of halos have grown massive enough, e.g.,
$M>10^{10}\,M_{\odot}\gg M_{\mathrm{min}}$.

A useful example to consider is the SFRD from Silva et al. (2013), which has a
much steeper redshift dependence such that it predicts about 3 times more star
formation at $z\sim 6$ and about 5 times less at $z\sim 10$. This is mainly
because the Silva et al. (2013) model predates and thus is not calibrated
against the much improved UVLFs measurements at $z>6$ (e.g., Bouwens et al.,
2015), which put stringent constraints on the SFR–halo mass relation
especially for faint galaxies. Even though star formation is allowed down to a
similar halo mass cutoff of $10^{8}\,M_{\odot}$, Silva et al. (2013) infer a
steeper SFR–halo mass relation that implies a UVLF faint-end slope now
disfavored by observations. Indeed, for more recent, UVLF-calibrated models
(e.g., Muñoz et al., 2022; Kannan et al., 2022b), the predicted steepness of
SFRD evolution is more comparable to that in LIMFAST, with the remaining
discrepancy in amplitude largely due to the varying minimum mass cutoff
assumed.

Finally, we stress that the difference in the SFRD evolution, combined with
the evolution it drives such as the chemical enrichment, serves as a major
source of the discrepancies observed for the LIM signals, as will be discussed
next.

####  4.2.2 LIM Signals at Cosmic Dawn and Reionization

A fundamental distinction between the LIM signals predicted by LIMFAST and
from the literature is that we have presented the intrinsic line emission,
without accounting for any attenuation. The attenuation effect by intervening
dust or neutral IGM can be treated on-the-fly or post-processed, analytically
or numerically. Because of these different approaches and complications, and
that the main goal of LIMFAST at this stage is to establish the connection
between galaxy formation physics and the production of different emission
lines, we leave the implementation of such attenuation effects to future work.
Yet, it is still useful to compare the line signals at face value in terms of
both their redshift evolution and spatial fluctuations measured by the power
spectrum.

For the 21 cm line, the fiducial case in LIMFAST predicts a global signal that
reaches a minimum of depth $\sim-110\,$mK centered at $z\approx 17$, which
appears somewhat earlier in time compared with some recent studies involving
similar assumptions about high-$z$ galaxy formation (see e.g., Mirocha et al.,
2017; Mirocha & Furlanetto, 2019; Park et al., 2019, where the trough locates
at $z\sim 12$). The difference can be best understood by considering the fact
that our fiducial case predicts an earlier onset of star formation, which
presents a non-negligible amount of Ly$\alpha$ and X-ray backgrounds at $z>15$
that can modulate the 21 cm spin temperature. Indeed, models that produce
similarly high SFRD by resorting to efficient star formation in low-mass halos
or the presence of Pop III stars (which we neglect in our fiducial
simulations; but see Parsons et al. 2022) predict 21 cm signals centered at
$z\gtrsim 15$ as well (e.g., Mesinger et al., 2016; Muñoz et al., 2022). At
lower redshifts when reionization began, the 21 cm signal is increasingly
driven by the ionization field, which is in turn sourced by the metagalactic
ionizing radiation background that rapidly builds up as the SFRD rises.

For Ly$\alpha$, Figure 6 shows good agreement between the shapes of redshift
evolution from LIMFAST and Comaschi & Ferrara (2016) for all three sources of
Ly$\alpha$ emission considered, even though we purposely choose to leave out
the dust correction in this work. This is likely because (1) both calculations
assume the atomic cooling threshold for $M_{\mathrm{min}}$ and are calibrated
against the observed UVLFs and CMB optical depth constraints, and (2) sources
that contribute the majority of Ly$\alpha$ emission at these high redshifts
are not strongly attenuated by dust, as is the case in Comaschi & Ferrara
(2016) at least. Generally, other studies like Pullen et al. (2014) and Heneka
et al. (2017) show a steeper decrease towards high redshift with an overall
lower amplitude because of the significant lower contribution from faint
galaxies as discussed above, and the IGM attenuation in the Ly$\alpha$ damping
wing may also play a role, especially when the neutral fraction is high
(Heneka et al., 2017). On the other hand, for Ly$\alpha$ emission from
recombinations in the diffuse ionized IGM, LIMFAST and Comaschi & Ferrara
(2016) both predict intensities a few to a hundred times lower than the
previous findings by Pullen et al. (2014) and Heneka et al. (2017). Because
the global reionization history is roughly the same for all these studies, the
offset in normalization might be caused by different model assumptions made
for the recombination rate, such as the choice between Case A (Heneka et al.,
2017) vs. Case B (LIMFAST) recombinations and homogeneous vs. inhomogeneous
recombinations (Park et al., 2019). We also note that the Ly$\alpha$ emission
due to IGM recombinations in Pullen et al. (2014) is derived from an SFRD that
seems exceedingly high ($\dot{\rho}_{*}\sim
10\,M_{\odot}\,\mathrm{yr^{-1}\,Mpc^{-3}}$ at $z=6$).

Similarly, for other star-formation lines like H$\alpha$ and [O II], the
redshift evolution in LIMFAST matches well that from THESAN 1 simulations,
despite that small discrepancies exist for both the slope and the
normalization as a result of differences in the cosmic star formation and
photoionization modeling. There is not such a match for THESAN 2 because these
simulations do not resolve the small halos that drive the signal at high
redshift, resulting in a steeper evolution (Kannan et al., 2022b). Finally, we
emphasize again the subtle yet non-trivial differences in the redshift
evolution and power spectrum arising from the usage or not of locally-
calibrated scaling relations, and from different photoionization modeling
assumptions (e.g., $\log U=-2$ vs. $U\propto Z^{-0.8}$) — a useful feature
that can be exploited to retrieve information about physical properties of the
ISM (Silva et al., 2017), which we will discuss in more detail next in Section
4.3.

####  4.2.3 Other Tools for Simulating Multi-tracer LIM

Over the course of recent years, several other numerical tools have been
developed for self-consistently simulating LIM signals of multiple tracers
across the electromagnetic spectrum. These tools make use of either semi-
analytic/empirical models of galaxy evolution to “paint on” galaxy properties
(including emission lines) to the distribution of dark matter halos extracted
from cosmological $N$-body simulations (see e.g., Yang et al., 2021; Béthermin
et al., 2022; Sato-Polito et al., 2022), or fully cosmological
magnetohydrodynamic (MHD) simulations to directly generate different line
intensity fields of interest (e.g., Kannan et al., 2022a). It is therefore
interesting to compare them with LIMFAST, which has distinctive features of
being semi-numerical, in order to understand the pros and cons of different
tools available for the community.

Compared to fully numerical simulations like THESAN (Kannan et al., 2022a),
the most notable difference of semi-analytic/empirical/numerical tools is the
computational cost. Even though the former captures much more physical
information about the evolution of radiation sources (e.g., details of gas
dynamics and radiative transfer) in a highly integrated manner, compromises
still need to be made due to the computational cost, which prevents the
simulations from resolving the ISM-scale physics or reaching large, Gpc-scale
volumes needed for probing astrophysics on large scales and cosmology. On the
other hand, by running only $N$-body simulations without including baryonic
physics, the “paint-on” method provides the access to larger cosmological
volumes while still achieving the mass resolution necessary for properly
capturing the statistics of halos and sub-halos. Depending on the type of
questions that need to be addressed, prescriptions motivated by semi-analytic
(e.g., Yang et al., 2021; Béthermin et al., 2022) or semi-empirical (e.g.,
Sato-Polito et al., 2022) galaxy evolution models may be adopted to paint
light onto halos.

The semi-numerical formalism LIMFAST inherits from 21cmFAST (Section 2) makes
it even more computationally efficient to approximate the LSS formation and
evolution as revealed by various LIM signals. The density-field-based method
allows fast realizations (in a matter of hours) of multiple line intensity
fields on a personal computer. While certain (potentially important) aspects
of simulating the intensity fields and their synergies with other LSS tracers
that are characteristic of “paint-on” simulations have not been implemented in
LIMFAST yet, such as the inclusion of scatter in various relations (Reis et
al., 2022; Murmu et al., 2023) and point sources like Lyman alpha emitters
(LAEs), they can be included through relatively simple extensions of the
current formalism using halo finders or subgrid sampling.

Another important and related distinction to be considered is how complete and
efficient the multi-scale and multi-phase cosmic gas may be simulated. None of
the aforementioned “paint-on” simulations has so far incorporated calculations
of the 21 cm signal during the EoR, although it is straightforward to post-
process the density field realizations with methods like the excursion set
formalism used in 21cmFAST/LIMFAST to describe the phase transition of the IGM
and therefore the 21 cm signal. As demonstrated in Paper II, capturing the
complementarity between the 21 cm line and star-formation lines tracing
opposite phases of the IGM during reionization plays a pivotal role in
constraining the astrophysics of high-$z$ galaxy formation. The highly
efficient (though approximate) simulations of structure formation and the
production of multiple LIM signals that LIMFAST provides pave the way for
fully Bayesian model inference studies in the near future, similar to what has
been established for analyzing the 21 cm signal (Greig & Mesinger, 2015) to
study both astrophysics and cosmology (Park et al., 2019; Jones et al., 2021;
Rudakovskyi et al., 2021; Muñoz et al., 2022).

###  4.3 Photoionization Modeling: Effects and Diagnosis

As introduced in Section 2.2, while we have explored and generated grids of
photoionization parameters like the metallicity $Z$ and ionization parameter
$U$ that allow flexible modeling of the star-formation lines from H II
regions, it is unclear how exactly $U$ behaves among the entire source
population and whether it connects to other galaxy properties like $Z$ at all.
Therefore, potentially large uncertainties exist in our photoionization
modeling, which substantially affect the predicted LIM signals, as has been
demonstrated explicitly in Sections 3.1.4 and 3.2. In order to further
illustrate this, we compare our fiducial case that assumes a fixed,
arbitrarily-chosen value of $U$ with two ad hoc model variations that adopt
physically-grounded but qualitatively different correlations between $U$ and
$Z$. We further compare the predicted LIM signals of these models against the
detectability of SPHEREx, in order to demonstrate how future LIM experiments
might help improve the photoionization modeling of high-$z$ galaxies.

The first case, which has already been briefly discussed and shown in Figures
6 and 7, builds on the wind-driven bubble model proposed by Dopita et al.
(2006), which predicts an anti-correlation $U\sim 3U_{c}(Z/Z_{\odot})^{-0.8}$,
where $U_{c}$ is a characteristic ionization parameter corresponding to
$Z=0.1Z_{\odot}$. Qualitatively, the negative correlation in this case can be
attributed to two effects — stellar winds that become both more opaque and
stronger at higher metallicities lead to more absorption of ionizing photons
and more diluted ionizing flux due to larger shocked regions. The second case
follows empirically the findings from photoionization models of nearby, MaNGA
galaxies by Ji & Yan (2022), which suggests $U\sim
0.15U_{c}(Z/Z_{\odot})^{0.5}$. Such a positive correlation contradicts the
expectation from the bubble model and might be a result of geometric effects
of H II regions in photoionization modeling (e.g., spherical vs. plane-
parallel), although the exact cause of the correlation is still an open
question that needs to be addressed by better data and models of the
photoionization and gas dynamics in H II regions. To be consistent with our
fiducial choice of $\log U=-2$, we normalize both $U$–$Z$ relations such that
they return $U=U_{c}=0.01$ for $Z=0.1\,Z_{\odot}$.

Effects of combining $U$ and $Z$ differently in photoionization models can be
conveniently projected into the observable space using, e.g., the O32–R23
diagram. The ionization-/excitation-sensitive line index $\mathrm{O32}=$ [O
III]/[O II] and the oxygen-abundance-sensitive line index $\mathrm{R23}=$ ([O
III]+[O II])/$\mathrm{H\beta}$ are often combined to jointly constrain the
photoionization model in the analysis of individual H II regions or galaxies
(e.g., Sanders et al., 2016; Strom et al., 2017). Even though LIM observations
do not resolve individual H II regions or galaxies, it is helpful to
understand how counterparts of these line indices, in terms of the mean line
intensities directly probed by LIM data, compare under different assumptions
of $Z$ and $U$. Both forward modeling and the interpretation of the mean line
ratios as key observables in LIM data analysis benefit from such close-up
views of the photoionization model.

Figure 10: The relationship between O32 and R23 as a diagnostic of the
photoionization model parameters $Z$ and $U$. Colored sequences of dots show
O32 and R23 values of individual emitters with different combinations of $Z$
and $U$. The sequences with fixed marker size in blue, green, and red
represent cases assuming $\log U=-2$ (fiducial), $-2.5$, and $-3$, with
metallicities in range $-1.5<\log(Z/Z_{\odot})<0.5$ being indicated by the
slowly evolving shades. For comparison, model variations assuming anti-
correlated (correlated) $U$ and $Z$ are shown by the sequences in gray
(purple), where $U$ decreases (increases) with $Z$ as indicated by the
changing marker size (Dopita et al., 2006; Ji & Yan, 2022), with bigger
circles representing larger $U$. Dashed lines show O32 and R23 values computed
from the mean intensities of individual lines at $z=5$, as simulated by
LIMFAST in three different cases ($U=0.01$, $U\propto Z^{-0.8}$, and $U\propto
Z^{0.5}$). The horizontal and vertical shaded bands shown for the fiducial
case ($U=0.01$) indicate the projected constraints from the SPHEREx deep-field
survey of the lines at $z=5$.

Figure 10 shows how O32 and R23 respond to changes in $U$ and its possible
dependence on $Z$. Predictions from five cases of the $U$–$Z$ relation are
shown by the sets of points in different colors, including the two
aforementioned model variations based on findings of Dopita et al. (2006) and
Ji & Yan (2022). Cases with fixed $U$ show similar trends for O32 versus R23,
with the normalization primarily set by the exact $U$ value. The non-monotonic
and double-valued behavior at high R23 and O32 is caused by the $Z$ dependence
of the oxygen line strengths, which increase and decrease with $Z$ at the
low-$Z$ and high-$Z$ ends, respectively (see also Figure 11). The same general
behavior persists for cases with correlated $U$ and $Z$, despite that when $U$
and $Z$ are anti-correlated the metallicity effect is partially counteracted
by the corresponding $U$ evolution, making the turnover at intermediate $Z$
much more gradual.

On top of pairs of O32 and R23 shown by the scattered points, we also overplot
the corresponding line indices formed by the derived mean line intensities
$\bar{I}_{\mathrm{[O\,II]}}$, $\bar{I}_{\mathrm{[O\,III]}}$, and
$\bar{I}_{\mathrm{H\beta}}$ at $z=5$ (the lowest redshift covered by current
LIMFAST simulations), as shown by the horizontal and vertical dashed lines for
cases with $U=0.01$, $U\propto Z^{-0.8}$, and $U\propto Z^{0.5}$. While these
quantities cannot be directly mapped to the mean physical conditions of
photoionization in $z=5$ galaxies due to its highly heterogeneous nature among
the galaxy population, measuring and comparing them with predictions from the
photoionization model still provides useful insights into the permitted
parameter space. With the shaded bands, we illustrate the levels of
uncertainty that can be achieved by the 200 deg2 SPHEREx deep-field survey for
these line indices. Assuming the nominal, deep-field survey depth222See the
publicly-available SPHEREx sensitivity data product
https://github.com/SPHEREx/Public-
products/blob/master/Surface_Brightness_v28_base_cbe.txt., we expect SPHEREx
to measure both O32 and R23 with $\mathrm{S/N}\approx 3$ at $z=5$. Clearly,
extreme scenarios that deviate significantly from our fiducial case, e.g.,
$\log U=-3$ (red), can be readily tested with the LIM-derived O32 and R23, but
it is challenging for LIM measurements to distinguish more similar cases like
those differing only in how $U$ scales with $Z$ at high redshift. Given how
complicated the detailed picture of photoionization might become and vary for
different galaxies, reliably forward modeling or inferring from the LIM
signals of lines such as [O II] and [O III] can be quite challenging,
especially at high redshifts where prior knowledge on the physical conditions
of H II regions is limited. Imminent JWST observations of these lines in
individual high-$z$ galaxies will therefore serve as a resourceful guide for
how to better model and interpret their LIM signals in the near future (e.g.,
Katz et al., 2023).

###  4.4 Applications of LIMFAST

The density field-based nature of LIMFAST simulations makes them suitable for
efficiently generating many realizations in the forward modeling and analysis
pipeline, which is pivotal to the success of any intensity mapping experiment.
As has been demonstrated by 21CMMC (Greig & Mesinger, 2015) for analyzing the
21 cm signal, it is also possible to build up a full Bayesian inference
framework for parameter studies, provided that a flexible, generative model
can be constructed, which is what we aim to demonstrate with LIMFAST for
investigations of the multi-tracer LIM concept. A few plans have been made so
far to implement LIMFAST into the simulation pipeline of forthcoming intensity
mapping surveys, including the Tomographic Ionized-carbon Mapping Experiment
(TIME; Crites et al. 2014; Sun et al. 2021b) and the Spectro-Photometer for
the History of the Universe, Epoch of Reionization, and Ices Explorer
(SPHEREx; Doré et al. 2014). While LIMFAST only runs down to $z=5$, the lack
of correlation between large-scale structure at distinct redshifts allows us
to combine with simulations dedicated to the post-reionization universe
($0<z<5$) and create complete light cones that include both the target signal
and different sources of contamination.

Meanwhile, the modular nature of LIMFAST makes it convenient to implement
model variations and theoretically study their effects on various observables.
Examples of this include the different forms of stellar feedback and the star
formation law investigated in Paper II, and the impact of Pop III star
formation with different IMF assumptions on the LIM signals of He II and
H$\alpha$ at cosmic dawn explored in Parsons et al. (2022). We have also
extended our metagalactic radiation background model to simulate intensity
mapping observations of the continuum emission at near-infrared, where the
cosmic near-infrared background reveals rich information about the collective
properties of the first stars and galaxies as well as their impact on
reionization (Sun et al., 2021a), especially when jointly analyzed with the 21
cm signal (Sun et al. in prep).

##  5 Conclusions

In this paper, we have introduced LIMFAST, a semi-numerical code to physically
model the production of emission lines in galaxies and the IGM during the
cosmic dawn and reionization eras and simulate their LIM signals. LIMFAST
implements the analytic galaxy formation and evolution models by Furlanetto
(2021) into the framework of structure formation and reionization simulations
in 21cmFAST, and utilizes pre-computed look-up tables from stellar population
synthesis (BPASS) and photoionization modeling (cloudy) to obtain emissivities
of the lines of interest for LIM studies.

We have detailed the galaxy formation model, the calculations of stellar and
nebular SEDs, and the various kinds of line emission that can be simulated
(Section 2). We have showcased physical evolution of the galaxy population and
the IGM derived from LIMFAST runs, followed by a collection of predicted LIM
signals in terms of their light cones, mean line intensities, and power
spectra. Emission line signals presented include the 21 cm line, the
Ly$\alpha$ IGM and background radiation fields, as well as the Ly$\alpha$,
H$\alpha$, H$\beta$, [O II] $3727\,{\rm\AA}$, and [O III] $5007\,{\rm\AA}$
lines arising from star formation (Section 3). We have compared our results
with a number of previous studies from the literature including the parent
code 21cmFAST, discussed key aspects that distinguish LIMFAST from them, and
introduced some case studies enabled by LIMFAST (Section 4). Several aspects
of our simulations are noteworthy, which we summarize and reiterate as
follows:

  1. 1.

With LIMFAST, we are able to generate histories of cosmic star formation,
reionization, and metal enrichment in a fully self-consistent manner that is
in good agreement with the observational constraints available so far. Due to
some assumptions we make — that star formation in galaxies is regulated by
momentum-driven feedback from supernovae and stars can form in halos of mass
as low as the atomic cooling limit with $T_{\mathrm{vir}}=10^{4}\,$K, we
predict overall a higher rate of cosmic star formation at $z>6$ that evolves
less strongly with redshift compared with most previous studies in the
literature. We note, though, that this picture of high-$z$ galaxy formation
still agrees well with observations and implies a plausible reionization
history. Neglecting other complications that also modulate the line strengths,
we argue that higher SFRD in general leads to stronger emission in the various
star-formation lines considered in this work.

  2. 2.

Tracing how chemical enrichment drives the SED evolution of the galaxy
population, as sources of both the ionizing radiation and the line intensity
fields, is important for accurately modeling the LIM signals of interest. Even
for the 21 cm line whose emissivity is not directly related to the galaxy
SEDs, order of unity differences in both sky-averaged and fluctuation signals
may result from ignoring the SED evolution. As can be seen from the comparison
with cases assuming locally-calibrated relations, this has an even larger
impact on the star-formation lines whose emissivities are usually strong
functions of metallicity $Z$. This marks an important source of astrophysical
uncertainty to be considered in future modeling and analysis of LIM
observations, especially for studying cosmology.

  3. 3.

Since we choose to not apply corrections for the attenuation effects on line
intensities due to dust and neutral gas, which themselves are highly uncertain
and model dependent, it is somewhat difficult to make fair comparisons between
results from our simulations and those in the literature with corrections
applied. Nonetheless, qualitative features like the steepness of redshift
evolution help understand the source(s) of discrepancies. Overall, we find
good qualitative agreement among predictions based on models calibrated to
latest observational data including the UVLFs at $z>6$.

  4. 4.

Detailed photoionization (or more generally ISM) modeling is sophisticated and
usually difficult to be done consistently with simplistic models of galaxy
formation and evolution. Using simple assumptions and ad hoc variations, we
have demonstrated challenges posed by this issue for modeling the line signals
(e.g., [O II] and [O III]) and also showed that it is useful to consider
combinations of the line ratios as a way to constrain key parameters in the
photoionization modeling such as the ionization parameter $U$. By analogy to
the strong line calibration method in the study of chemical evolution of
individual galaxies, we have illustrated how $Z$ and $U$ may be jointed probed
by observational data from LIM surveys like SPHEREx.

  5. 5.

Although LIMFAST in its current form is obviously too simplistic to be applied
to all kinds of problems in the study of the high-redshift universe using the
intensity mapping technique, interesting applications of it to real-world LIM
efforts are already underway. Thanks to its modular nature (Figure 1), we
expect new features to be implemented and complement the existing framework as
new astrophysical and cosmological applications being pursued.

To conclude, following up on this methodology paper, in Paper II of this
series we take into account models of the neutral ISM for the inclusion of the
[C II] 158 $\mu$m and CO lines, which are popular target lines for sub-mm LIM.
We anticipate further implementations such as the aforementioned attenuation
effects due to neutral gas and dust, which can be self-consistently included
into the existing formalism. LIMFAST enables time-efficient multi-tracer
simulations of large cosmological volumes, which make it a useful complement
to more sophisticated numerical simulations such as the recently presented
THESAN simulations (Kannan et al., 2022a), for not only forward modeling but
also model inference from the LIM data sets.

We thank the anonymous reviewer for their helpful comments that improve this
paper. We are grateful to Steven Furlanetto, Adam Trapp and Fred Davies for
useful discussions and comments throughout this project, as well as Rahul
Kannan and Yuxiang Qin for their valuable comments on the early draft of this
paper. We acknowledge support from the JPL R&TD strategic initiative grant on
line intensity mapping. This research was carried out at the Jet Propulsion
Laboratory, California Institute of Technology, under a contract with the
National Aeronautics and Space Administration.

Figure 11: Dependence of the luminosity with metallicity and ionization
parameter for lines arising from star formation. Note the drastically
different scales in the vertical axes for different lines.

##  Appendix A Luminosity Dependence on $U$ and $Z$

Figure 11 shows the dependence of the luminosity of lines from star formation
on metallicity and ionization parameter. The luminosity of the oxygen lines
changes by orders of magnitude when varying the metallicity, while the changes
in the hydrogen line luminosities are only of a factor of a few. Furthermore,
varying the ionization parameter also has a major impact (up to an order of
magnitude) on the oxygen lines, with opposite effects for the O II and O III
ions, whereas the impact on the hydrogen lines is less than order of unity.

##  Appendix B Redshift-Space Distortions: A Pedagogical Description

We briefly revisit here the formalism of RSD, and show its application to LIM
observations targeting 21 cm and other emission lines. Our discussion is based
on that in Mao et al. (2012), to which we refer the interested reader for
details.

Observing an emission line at a specific frequency indicates, in principle,
the position along the LOS of the sources producing the line. In reality,
however, the observed redshift can be different from the real cosmological one
because not only the Hubble expansion determines its value; the LOS peculiar
velocity of the sources affects the value of the observed frequency via the
Doppler effect which, in turn, results in deviations from the true redshift.
Quantitatively, when peculiar velocities along the LOS, $v_{\parallel}$, are
present, the position $\boldsymbol{r}$ of a source in real space appears to a
position $\boldsymbol{s}$ in redshift space through the expression

$\boldsymbol{s}=\boldsymbol{r}+\frac{1+z_{\rm obs}}{H(z_{\rm
obs})}v_{\parallel}\hat{r}~{}.$ (B1)

Here, $1+z_{\rm obs}=(1+z)(1-v_{\parallel}/c)^{-1}$, $\hat{r}$ is the unit
vector along the LOS direction, and $z_{\rm obs}$ and $z$ are the observed and
cosmological redshifts, respectively. Because the peculiar velocities of the
sources are typically unknown in observations, the redshift space information
is the only one that is observable. Therefore, to predict observed quantities
from our simulations, we need to account for the effect of peculiar velocities
on calculations in real space. We use the term redshift-space distortions, or
RSDs, here to refer to the change in apparent position induced by the peculiar
velocity of the sources (see Hamilton, 1998, for a review on RSDs).

As detailed in Section 4 in Mao et al. (2012), deriving the cosmological
radiative transfer equations for radiation emitted from sources with peculiar
velocity, and considering an optically-thin medium along the LOS, one finds
that the observed intensity can be written as (Equation 12 in Mao et al.,
2012)

$I_{\nu_{\rm obs}}=\frac{cL_{{\nu_{0}}}\,a\,\nu^{2}_{\rm
obs}}{4\pi\nu^{3}_{0}H(a)}\frac{n^{\prime}}{\left|1+(dv_{\parallel}/dr_{\parallel})/[aH(a)]\right|}~{}.$
(B2)

Here, $L_{{\nu_{0}}}$ is the rest-frame luminosity of the sources, $H(a)$
represents the Hubble parameter at scale factor $a$, and $\nu_{0}$ and
$\nu_{\rm obs}$ are the rest-frame and observed frequencies of the emission,
respectively. The rightmost term denotes the apparent change in the number
density of sources due to their peculiar velocities, where $n^{\prime}$ is the
proper number density of sources in real space, and
${dv_{\parallel}}/{dr_{\parallel}}$ is the velocity gradient along the LOS.
Expressing now the observed frequency in terms of the rest-frame one via
$\nu_{\rm obs}=\nu_{0}(1+z)^{-1}(1-v_{\parallel}/c)$ and considering the
comoving number density of sources, $n$, we can rewrite the previous equation
as

$\displaystyle I_{\nu_{\rm obs}}$
$\displaystyle=\frac{cL_{{\nu_{0}}}(1-v_{\parallel}/c)^{2}}{4\pi\nu_{0}H(z)}\frac{n}{\left|1+(1+z)(dv_{\parallel}/dr_{\parallel})/H(z)\right|}~{}$
$\displaystyle\simeq\frac{I_{\nu}(z)}{\left|1+(1+z)(dv_{\parallel}/dr_{\parallel})/H(z)\right|}~{}.$
(B3)

This highlights the fact that the observed intensity considering the RSDs is
simply roughly equivalent to the specific intensity computed by LIMFAST
(Equation 12), with a velocity gradient correction denoted by the denominator.
Here, we have omitted the correction term $(1-v_{\parallel}/c)^{2}$ in the
numerator because its effect is much smaller than that from the denominator.
The correction shown in the last term of Equation (B) is the one that we will
apply to the intensity of the optical and UV lines computed by LIMFAST. This
is the same correction as for the 21 cm radiation field shown in Equation
(2.6), but there the 21 cm signal is expressed as a differential brightness
temperature instead of as an intensity.

As mentioned before, the above derivation considers an optically-thin medium,
which is valid for all optical and UV emission lines considered in our work
except for Ly$\alpha$. The reason for this difference is that Ly$\alpha$ is
severely affected by the neutral hydrogen gas, both within halos and in the
IGM, due to its resonant nature (see a review on the physics of Ly$\alpha$ in
Dijkstra, 2017)333See a discussion about the optically thin assumption and its
validity for 21 cm emission in Mao et al. (2012).. For simplicity, we do not
address further corrections for Ly$\alpha$ here, but we refer the interested
reader to Zheng et al. (2011) for an investigation of the impact of Ly$\alpha$
radiative transfer to the RSDs.

Figure 12: From left to right, 2D power spectra of the 21 cm, Ly$\alpha$ (from
star formation and the ionized IGM), and [O III] lines simulated by LIMFAST at
$z=7$ (top row) and 11 (bottom row), where $\bar{x}_{\mathrm{HI}}\approx 0.3$
and 0.9, respectively. The overlaid contours highlight the shape difference of
the power spectra, with the solid and dashed ones denoting cases with and
without the RSD effect, respectively. The amplification/suppression and
anisotropies induced by the RSDs vary for different lines at different
redshifts.

Finally, the observable of interest that is sensitive to the RSDs is the power
spectrum. Thus, it is also useful to write the intensity (or brightness
temperature) in terms of fluctuations to visualize the impact from the
clustering bias of the sources and from the RSDs when comparing the power
spectra of different emission lines. For the intensity case and considering
linear scales,

$I_{\nu_{\rm obs}}=\bar{I}_{\nu}(z)\,\frac{1+b\,\delta_{\rm
DM}(\bf{r})}{|1+\delta_{\partial_{r}v}({\bf r})|}.$ (B4)

Here, $\bar{I}_{\nu}(z)$ is the cosmic mean intensity of a given emission
line, $\delta_{\rm DM}(\bf{r})$ denotes the fluctuations of the dark matter
density field, and $\delta_{\partial_{r}v}({\bf
r})=[(1+z)/H(z)]dv_{\parallel}/dr_{\parallel}$, with all three quantities
considered in real space. The term $b$ denotes the bias factor of the line
emission with respect to the dark matter field, where we assume this emission
bias to be the same as that of the sources since we ignore radiative transfer
effects. Overall, Equation (B4) readily shows the known result that the impact
of the RSDs will be more significant when the sources of radiation are less
biased. This case corresponds to emission sourced by the intergalactic gas,
i.e., 21 cm and IGM Ly$\alpha$444We remind here that we do not account for the
Ly$\alpha$ radiative transfer that can also affect the fluctuations of this
emission line., while emission from star formation tracing much more biased
galaxies will be less sensitive to the RSDs.

Figure 12 displays the two dimensional power spectra in redshift space for a
selection of lines simulated by LIMFAST at redshifts $z=7$ (top row) and
$z=11$ (bottom row). Each panel is color coded by the power, i.e.,
$I_{\nu}^{2}/V$, in units of $\mathrm{mK^{2}\,Mpc^{3}}$ for the 21 cm line or
$\mathrm{(Jy/sr)^{2}\,Mpc^{3}}$ otherwise. The overlaid contours contrast the
shape difference between real-space (dashed) and redshift-space (solid) power
spectra, to highlight the anisotropies in redshift space induced by the RSDs.
Overall, the RSDs are most visible for lines that trace and are sourced by the
intergalactic gas, i.e., the 21 cm line and the intergalactic Ly$\alpha$
radiation. When the clustering bias of the sources is large, as is the case
for star-formation lines, the impact of the RSDs is less important, as
indicated by Equation (B4). Also noteworthy is that the RSDs appear to affect
the 21 cm 2D power spectra in opposite ways at $z=7$ and 11, which correspond
to IGM neutral fractions $\bar{x}_{\mathrm{HI}}\approx 0.3$ and 0.9,
respectively. This is mainly due to the progress of the reionization process,
when an anti-correlation between the ionization and density fields contributes
more to the 21 cm fluctuations as the IGM becomes more ionized (see the
discussion in Jensen et al., 2013, for more details).

## References

  * Barkana & Loeb (2005) Barkana, R., & Loeb, A. 2005, ApJ, 626, 1, doi: 10.1086/429954
  * Benson (2010) Benson, A. J. 2010, Phys. Rep., 495, 33, doi: 10.1016/j.physrep.2010.06.001
  * Bera et al. (2022) Bera, A., Samui, S., & Datta, K. K. 2022, MNRAS, doi: 10.1093/mnras/stac3814
  * Béthermin et al. (2022) Béthermin, M., Gkogkou, A., Van Cuyck, M., et al. 2022, A&A, 667, A156, doi: 10.1051/0004-6361/202243888
  * Bolmer et al. (2019) Bolmer, J., Ledoux, C., Wiseman, P., et al. 2019, A&A, 623, A43, doi: 10.1051/0004-6361/201834422
  * Bouwens et al. (2022) Bouwens, R. J., Illingworth, G., Ellis, R. S., Oesch, P., & Stefanon, M. 2022, ApJ, 940, 55, doi: 10.3847/1538-4357/ac86d1
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47, doi: 10.3847/1538-3881/abf83e
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66
  * Capak et al. (2015) Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455, doi: 10.1038/nature14500
  * Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45, doi: 10.1016/j.physrep.2014.02.009
  * Casey et al. (2018) Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018, ApJ, 862, 77, doi: 10.3847/1538-4357/aac82d
  * Chang et al. (2010) Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Nature, 466, 463, doi: 10.1038/nature09187
  * Chang et al. (2008) Chang, T.-C., Pen, U.-L., Peterson, J. B., & McDonald, P. 2008, Phys. Rev. Lett., 100, 091303, doi: 10.1103/PhysRevLett.100.091303
  * Choudhury & Ferrara (2007) Choudhury, T. R., & Ferrara, A. 2007, Monthly Notices of the Royal Astronomical Society: Letters, 380, L6, doi: 10.1111/j.1745-3933.2007.00338.x
  * Chung et al. (2019) Chung, D. T., Viero, M. P., Church, S. E., et al. 2019, ApJ, 872, 186, doi: 10.3847/1538-4357/ab0027
  * Comaschi & Ferrara (2016) Comaschi, P., & Ferrara, A. 2016, MNRAS, 455, 725, doi: 10.1093/mnras/stv2339
  * Crites et al. (2014) Crites, A. T., Bock, J. J., Bradford, C. M., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, ed. W. S. Holland & J. Zmuidzinas, 91531W, doi: 10.1117/12.2057207
  * Davies & Furlanetto (2022) Davies, F. B., & Furlanetto, S. R. 2022, MNRAS, 514, 1302, doi: 10.1093/mnras/stac1005
  * Davies et al. (2018) Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, 864, 142, doi: 10.3847/1538-4357/aad6dc
  * Dekel & Mandelker (2014) Dekel, A., & Mandelker, N. 2014, MNRAS, 444, 2071, doi: 10.1093/mnras/stu1427
  * Dijkstra (2017) Dijkstra, M. 2017, arXiv e-prints, arXiv:1704.03416.  https://arxiv.org/abs/1704.03416
  * Dopita et al. (2006) Dopita, M. A., Fischera, J., Sutherland, R. S., et al. 2006, ApJ, 647, 244, doi: 10.1086/505418
  * Doré et al. (2014) Doré, O., Bock, J., Ashby, M., et al. 2014, arXiv e-prints, arXiv:1412.4872.  https://arxiv.org/abs/1412.4872
  * Dumitru et al. (2019) Dumitru, S., Kulkarni, G., Lagache, G., & Haehnelt, M. G. 2019, MNRAS, 485, 3486, doi: 10.1093/mnras/stz617
  * Eide et al. (2020) Eide, M. B., Ciardi, B., Graziani, L., et al. 2020, MNRAS, 498, 6083, doi: 10.1093/mnras/staa2774
  * Eide et al. (2018) Eide, M. B., Graziani, L., Ciardi, B., et al. 2018, MNRAS, 476, 1174, doi: 10.1093/mnras/sty272
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMXAA, 53, 385\.  https://arxiv.org/abs/1705.10877
  * Fontanot et al. (2012) Fontanot, F., Cristiani, S., & Vanzella, E. 2012, MNRAS, 425, 1413, doi: 10.1111/j.1365-2966.2012.21594.x
  * Furlanetto (2021) Furlanetto, S. R. 2021, MNRAS, 500, 3394, doi: 10.1093/mnras/staa3451
  * Furlanetto et al. (2017) Furlanetto, S. R., Mirocha, J., Mebane, R. H., & Sun, G. 2017, MNRAS, 472, 1576, doi: 10.1093/mnras/stx2132
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181, doi: 10.1016/j.physrep.2006.08.002
  * Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, The Astrophysical Journal, 613, 1, doi: 10.1086/423025
  * Gong et al. (2012) Gong, Y., Cooray, A., Silva, M., et al. 2012, ApJ, 745, 49, doi: 10.1088/0004-637X/745/1/49
  * Gong et al. (2011) Gong, Y., Cooray, A., Silva, M. B., Santos, M. G., & Lubin, P. 2011, ApJL, 728, L46, doi: 10.1088/2041-8205/728/2/L46
  * Gong et al. (2017) Gong, Y., Cooray, A., Silva, M. B., et al. 2017, ApJ, 835, 273, doi: 10.3847/1538-4357/835/2/273
  * Graham & Fruchter (2017) Graham, J. F., & Fruchter, A. S. 2017, ApJ, 834, 170, doi: 10.3847/1538-4357/834/2/170
  * Greig & Mesinger (2015) Greig, B., & Mesinger, A. 2015, MNRAS, 449, 4246, doi: 10.1093/mnras/stv571
  * Greig & Mesinger (2018) —. 2018, MNRAS, 477, 3217, doi: 10.1093/mnras/sty796
  * Greig et al. (2019) Greig, B., Mesinger, A., & Bañados, E. 2019, MNRAS, 484, 5094, doi: 10.1093/mnras/stz230
  * Greig et al. (2017) Greig, B., Mesinger, A., Haiman, Z., & Simcoe, R. A. 2017, MNRAS, 466, 4239, doi: 10.1093/mnras/stw3351
  * Hamilton (1998) Hamilton, A. J. S. 1998, in Astrophysics and Space Science Library, Vol. 231, The Evolving Universe, ed. D. Hamilton, 185, doi: 10.1007/978-94-011-4960-0_17
  * Heneka & Cooray (2021) Heneka, C., & Cooray, A. 2021, MNRAS, 506, 1573, doi: 10.1093/mnras/stab1842
  * Heneka et al. (2017) Heneka, C., Cooray, A., & Feng, C. 2017, ApJ, 848, 52, doi: 10.3847/1538-4357/aa8eed
  * Hoag et al. (2019) Hoag, A., Bradač, M., Huang, K., et al. 2019, ApJ, 878, 12, doi: 10.3847/1538-4357/ab1de7
  * Holzbauer & Furlanetto (2012) Holzbauer, L. N., & Furlanetto, S. R. 2012, MNRAS, 419, 718, doi: 10.1111/j.1365-2966.2011.19752.x
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Ihle et al. (2019) Ihle, H. T., Chung, D., Stein, G., et al. 2019, ApJ, 871, 75, doi: 10.3847/1538-4357/aaf4bc
  * Jenkins et al. (2001) Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372, doi: 10.1046/j.1365-8711.2001.04029.x
  * Jensen et al. (2013) Jensen, H., Datta, K. K., Mellema, G., et al. 2013, MNRAS, 435, 460, doi: 10.1093/mnras/stt1341
  * Ji & Yan (2022) Ji, X., & Yan, R. 2022, A&A, 659, A112, doi: 10.1051/0004-6361/202142312
  * Jones et al. (2021) Jones, D., Palatnick, S., Chen, R., Beane, A., & Lidz, A. 2021, ApJ, 913, 7, doi: 10.3847/1538-4357/abf0a9
  * Kannan et al. (2022a) Kannan, R., Garaldi, E., Smith, A., et al. 2022a, MNRAS, 511, 4005, doi: 10.1093/mnras/stab3710
  * Kannan et al. (2020) Kannan, R., Marinacci, F., Vogelsberger, M., et al. 2020, MNRAS, 499, 5732, doi: 10.1093/mnras/staa3249
  * Kannan et al. (2022b) Kannan, R., Smith, A., Garaldi, E., et al. 2022b, MNRAS, doi: 10.1093/mnras/stac1557
  * Katz et al. (2019) Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902, doi: 10.1093/mnras/stz1672
  * Katz et al. (2023) Katz, H., Saxena, A., Cameron, A. J., et al. 2023, MNRAS, 518, 592, doi: 10.1093/mnras/stac2657
  * Kaur et al. (2022) Kaur, H. D., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 513, 5097, doi: 10.1093/mnras/stac1226
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ARAA, 36, 189, doi: 10.1146/annurev.astro.36.1.189
  * Kewley et al. (2019) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARAA, 57, 511, doi: 10.1146/annurev-astro-081817-051832
  * Kovetz et al. (2019) Kovetz, E., Breysse, P. C., Lidz, A., et al. 2019, BAAS, 51, 101.  https://arxiv.org/abs/1903.04496
  * Lacey & Cole (1993) Lacey, C., & Cole, S. 1993, MNRAS, 262, 627, doi: 10.1093/mnras/262.3.627
  * Lewis et al. (2022) Lewis, J. S. W., Ocvirk, P., Sorce, J. G., et al. 2022, arXiv e-prints, arXiv:2202.05869.  https://arxiv.org/abs/2202.05869
  * Li et al. (2016) Li, T. Y., Wechsler, R. H., Devaraj, K., & Church, S. E. 2016, ApJ, 817, 169, doi: 10.3847/0004-637X/817/2/169
  * Lidz et al. (2011) Lidz, A., Furlanetto, S. R., Oh, S. P., et al. 2011, ApJ, 741, 70, doi: 10.1088/0004-637X/741/2/70
  * Liu & Shaw (2020) Liu, A., & Shaw, J. R. 2020, PASP, 132, 062001, doi: 10.1088/1538-3873/ab5bfd
  * Ly et al. (2007) Ly, C., Malkan, M. A., Kashikawa, N., et al. 2007, ApJ, 657, 738, doi: 10.1086/510828
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, ApJ, 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Madau et al. (1997) Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429, doi: 10.1086/303549
  * Mao et al. (2012) Mao, Y., Shapiro, P. R., Mellema, G., et al. 2012, MNRAS, 422, 926, doi: 10.1111/j.1365-2966.2012.20471.x
  * Mas-Ribas & Chang (2020) Mas-Ribas, L., & Chang, T.-C. 2020, Phys. Rev. D, 101, 083032, doi: 10.1103/PhysRevD.101.083032
  * Mason et al. (2015) Mason, C. A., Trenti, M., & Treu, T. 2015, ApJ, 813, 21, doi: 10.1088/0004-637X/813/1/21
  * Mason et al. (2019) Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947, doi: 10.1093/mnras/stz632
  * McGreer et al. (2015) McGreer, I. D., Mesinger, A., & D’Odorico, V. 2015, MNRAS, 447, 499, doi: 10.1093/mnras/stu2449
  * Mebane et al. (2018) Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2018, Monthly Notices of the Royal Astronomical Society, 479, 4544, doi: 10.1093/mnras/sty1833
  * Mebane et al. (2020) Mebane, R. H., Mirocha, J., & Furlanetto, S. R. 2020, MNRAS, 493, 1217, doi: 10.1093/mnras/staa280
  * Mesinger & Furlanetto (2007) Mesinger, A., & Furlanetto, S. 2007, ApJ, 669, 663, doi: 10.1086/521806
  * Mesinger et al. (2011) Mesinger, A., Furlanetto, S., & Cen, R. 2011, Monthly Notices of the Royal Astronomical Society, 411, 955, doi: 10.1111/j.1365-2966.2010.17731.x
  * Mesinger & Furlanetto (2008) Mesinger, A., & Furlanetto, S. R. 2008, MNRAS, 385, 1348, doi: 10.1111/j.1365-2966.2007.12836.x
  * Mesinger et al. (2016) Mesinger, A., Greig, B., & Sobacchi, E. 2016, MNRAS, 459, 2342, doi: 10.1093/mnras/stw831
  * Metha et al. (2021) Metha, B., Cameron, A. J., & Trenti, M. 2021, MNRAS, 504, 5992, doi: 10.1093/mnras/stab892
  * Mirocha & Furlanetto (2019) Mirocha, J., & Furlanetto, S. R. 2019, MNRAS, 483, 1980, doi: 10.1093/mnras/sty3260
  * Mirocha et al. (2017) Mirocha, J., Furlanetto, S. R., & Sun, G. 2017, MNRAS, 464, 1365, doi: 10.1093/mnras/stw2412
  * Mirocha et al. (2022) Mirocha, J., Liu, A., & La Plante, P. 2022, MNRAS, 516, 4123, doi: 10.1093/mnras/stac2530
  * Mirocha et al. (2018) Mirocha, J., Mebane, R. H., Furlanetto, S. R., Singal, K., & Trinh, D. 2018, MNRAS, 478, 5591, doi: 10.1093/mnras/sty1388
  * Muñoz et al. (2022) Muñoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657, doi: 10.1093/mnras/stac185
  * Murmu et al. (2023) Murmu, C. S., Olsen, K. P., Greve, T. R., et al. 2023, MNRAS, 518, 3074, doi: 10.1093/mnras/stac3304
  * Padmanabhan et al. (2021) Padmanabhan, H., Breysse, P., Lidz, A., & Switzer, E. R. 2021, arXiv e-prints, arXiv:2105.12148.  https://arxiv.org/abs/2105.12148
  * Pallottini et al. (2019) Pallottini, A., Ferrara, A., Decataldo, D., et al. 2019, MNRAS, 487, 1689, doi: 10.1093/mnras/stz1383
  * Palmerio et al. (2019) Palmerio, J. T., Vergani, S. D., Salvaterra, R., et al. 2019, A&A, 623, A26, doi: 10.1051/0004-6361/201834179
  * Park et al. (2019) Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, Monthly Notices of the Royal Astronomical Society, 484, 933, doi: 10.1093/mnras/stz032
  * Parsons et al. (2022) Parsons, J., Mas-Ribas, L., Sun, G., et al. 2022, ApJ, 933, 141, doi: 10.3847/1538-4357/ac746b
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Popping et al. (2017) Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152, doi: 10.1093/mnras/stx1545
  * Pritchard & Loeb (2012) Pritchard, J. R., & Loeb, A. 2012, Reports on Progress in Physics, 75, 086901, doi: 10.1088/0034-4885/75/8/086901
  * Pullen et al. (2013) Pullen, A. R., Chang, T.-C., Doré, O., & Lidz, A. 2013, ApJ, 768, 15, doi: 10.1088/0004-637X/768/1/15
  * Pullen et al. (2014) Pullen, A. R., Doré, O., & Bock, J. 2014, ApJ, 786, 111, doi: 10.1088/0004-637X/786/2/111
  * Qin et al. (2021) Qin, Y., Mesinger, A., Greig, B., & Park, J. 2021, MNRAS, 501, 4748, doi: 10.1093/mnras/staa3408
  * Qin et al. (2020) Qin, Y., Mesinger, A., Park, J., Greig, B., & Muñoz, J. B. 2020, MNRAS, 495, 123, doi: 10.1093/mnras/staa1131
  * Reis et al. (2022) Reis, I., Barkana, R., & Fialkov, A. 2022, MNRAS, 511, 5265, doi: 10.1093/mnras/stac411
  * Righi et al. (2008) Righi, M., Hernández-Monteagudo, C., & Sunyaev, R. A. 2008, A&A, 489, 489, doi: 10.1051/0004-6361:200810199
  * Robertson et al. (2015) Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJL, 802, L19, doi: 10.1088/2041-8205/802/2/L19
  * Rudakovskyi et al. (2021) Rudakovskyi, A., Mesinger, A., Savchenko, D., & Gillet, N. 2021, MNRAS, 507, 3046, doi: 10.1093/mnras/stab2333
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Sanders et al. (2016) Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23, doi: 10.3847/0004-637X/816/1/23
  * Sato-Polito et al. (2022) Sato-Polito, G., Kokron, N., & Bernal, J. L. 2022, arXiv e-prints, arXiv:2212.08056, doi: 10.48550/arXiv.2212.08056
  * Scalo (1998) Scalo, J. 1998, in Astronomical Society of the Pacific Conference Series, Vol. 142, The Stellar Initial Mass Function (38th Herstmonceux Conference), ed. G. Gilmore & D. Howell, 201.  https://arxiv.org/abs/astro-ph/9712317
  * Scoccimarro (1998) Scoccimarro, R. 1998, MNRAS, 299, 1097, doi: 10.1046/j.1365-8711.1998.01845.x
  * Scott & Rees (1990) Scott, D., & Rees, M. J. 1990, MNRAS, 247, 510
  * Shen et al. (2022) Shen, X., Vogelsberger, M., Nelson, D., et al. 2022, MNRAS, 510, 5560, doi: 10.1093/mnras/stab3794
  * Shen et al. (2020) Shen, X., Vogelsberger, M., Nelson, D., et al. 2020, Monthly Notices of the Royal Astronomical Society, 495, 4747, doi: 10.1093/mnras/staa1423
  * Sheth & Tormen (1999) Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119, doi: 10.1046/j.1365-8711.1999.02692.x
  * Silva et al. (2017) Silva, B. M., Zaroubi, S., Kooistra, R., & Cooray, A. 2017, Monthly Notices of the Royal Astronomical Society, 475, 1587, doi: 10.1093/mnras/stx3265
  * Silva et al. (2015) Silva, M., Santos, M. G., Cooray, A., & Gong, Y. 2015, ApJ, 806, 209, doi: 10.1088/0004-637X/806/2/209
  * Silva et al. (2013) Silva, M. B., Santos, M. G., Gong, Y., Cooray, A., & Bock, J. 2013, ApJ, 763, 132, doi: 10.1088/0004-637X/763/2/132
  * Somerville & Kolatt (1999) Somerville, R. S., & Kolatt, T. S. 1999, MNRAS, 305, 1, doi: 10.1046/j.1365-8711.1999.02154.x
  * Strom et al. (2017) Strom, A. L., Steidel, C. C., Rudie, G. C., et al. 2017, ApJ, 836, 164, doi: 10.3847/1538-4357/836/2/164
  * Suginohara et al. (1999) Suginohara, M., Suginohara, T., & Spergel, D. N. 1999, ApJ, 512, 547, doi: 10.1086/306787
  * Sun (2022) Sun, G. 2022, ApJL, 931, L29, doi: 10.3847/2041-8213/ac7138
  * Sun & Furlanetto (2016) Sun, G., & Furlanetto, S. R. 2016, MNRAS, 460, 417, doi: 10.1093/mnras/stw980
  * Sun et al. (2019) Sun, G., Hensley, B. S., Chang, T.-C., Doré, O., & Serra, P. 2019, ApJ, 887, 142, doi: 10.3847/1538-4357/ab55df
  * Sun et al. (2022) Sun, G., Mas-Ribas, L., Chang, T.-C., et al. 2022, arXiv e-prints, arXiv:2206.14186.  https://arxiv.org/abs/2206.14186
  * Sun et al. (2021a) Sun, G., Mirocha, J., Mebane, R. H., & Furlanetto, S. R. 2021a, MNRAS, 508, 1954, doi: 10.1093/mnras/stab2697
  * Sun et al. (2021b) Sun, G., Chang, T. C., Uzgil, B. D., et al. 2021b, ApJ, 915, 33, doi: 10.3847/1538-4357/abfe62
  * Switzer et al. (2013) Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, Monthly Notices of the Royal Astronomical Society: Letters, 434, L46, doi: 10.1093/mnrasl/slt074
  * Tanaka & Hasegawa (2021) Tanaka, T., & Hasegawa, K. 2021, MNRAS, 502, 463, doi: 10.1093/mnras/stab072
  * Taylor et al. (2022) Taylor, A. J., Barger, A. J., & Cowie, L. L. 2022, ApJL, 939, L3, doi: 10.3847/2041-8213/ac959d
  * Thöne et al. (2013) Thöne, C. C., Fynbo, J. P. U., Goldoni, P., et al. 2013, MNRAS, 428, 3590, doi: 10.1093/mnras/sts303
  * Visbal et al. (2015) Visbal, E., Haiman, Z., & Bryan, G. L. 2015, MNRAS, 450, 2506, doi: 10.1093/mnras/stv785
  * Visbal & Loeb (2010) Visbal, E., & Loeb, A. 2010, Journal of Cosmology and Astro-Particle Physics, 2010, 016, doi: 10.1088/1475-7516/2010/11/016
  * Visbal et al. (2009) Visbal, E., Loeb, A., & Wyithe, S. 2009, JCAP, 2009, 030, doi: 10.1088/1475-7516/2009/10/030
  * Visbal & McQuinn (2018) Visbal, E., & McQuinn, M. 2018, ApJL, 863, L6, doi: 10.3847/2041-8213/aad5e6
  * Vogelsberger et al. (2020) Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nature Reviews Physics, 2, 42, doi: 10.1038/s42254-019-0127-2
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * Whitler et al. (2020) Whitler, L. R., Mason, C. A., Ren, K., et al. 2020, MNRAS, 495, 3602, doi: 10.1093/mnras/staa1178
  * Xiao et al. (2018) Xiao, L., Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 477, 904, doi: 10.1093/mnras/sty646
  * Yang et al. (2021) Yang, S., Somerville, R. S., Pullen, A. R., et al. 2021, ApJ, 911, 132, doi: 10.3847/1538-4357/abec75
  * Yue & Ferrara (2019) Yue, B., & Ferrara, A. 2019, MNRAS, 490, 1928, doi: 10.1093/mnras/stz2728
  * Yue et al. (2015) Yue, B., Ferrara, A., Pallottini, A., Gallerani, S., & Vallini, L. 2015, MNRAS, 450, 3829, doi: 10.1093/mnras/stv933
  * Yue et al. (2014) Yue, B., Ferrara, A., Vanzella, E., & Salvaterra, R. 2014, Monthly Notices of the Royal Astronomical Society: Letters, 443, L20, doi: 10.1093/mnrasl/slu073
  * Yue et al. (2018) Yue, B., Castellano, M., Ferrara, A., et al. 2018, ApJ, 868, 115, doi: 10.3847/1538-4357/aae77f
  * Zel’Dovich (1970) Zel’Dovich, Y. B. 1970, A&A, 500, 13
  * Zheng et al. (2011) Zheng, Z., Cen, R., Trac, H., & Miralda-Escudé, J. 2011, ApJ, 726, 38, doi: 10.1088/0004-637X/726/1/38
  * Zhu et al. (2022) Zhu, Y., Becker, G. D., Bosman, S. E. I., et al. 2022, ApJ, 932, 76, doi: 10.3847/1538-4357/ac6e60