# LIMFAST. II. Line Intensity Mapping as a Probe of High-Redshift Galaxy
Formation

Guochao Sun  California Institute of Technology, 1200 E. California Blvd.,
Pasadena, CA 91125, USA  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  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  Steven R. Furlanetto  Department of Physics and Astronomy,
University of California, Los Angeles, CA 90024, USA  Richard H. Mebane
Department of Physics and Astronomy, University of California, Los Angeles, CA
90024, USA  Department of Astronomy and Astrophysics, University of
California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA  Michael
O. Gonzalez  California Institute of Technology, 1200 E. California Blvd.,
Pasadena, CA 91125, USA  Jasmine Parsons  McGill University, Department of
Physics & McGill Space Institute, 3600 Rue University, Montréal, QC, H3A 2T8
A. C. Trapp  Department of Physics and Astronomy, University of California,
Los Angeles, CA 90024, USA
©2022. All rights reserved.

###### Abstract

The epoch of reionization (EoR) offers a unique window into the dawn of galaxy
formation, through which high-redshift galaxies can be studied by observations
of both themselves and their impact on the intergalactic medium. Line
intensity mapping (LIM) promises to explore cosmic reionization and its
driving sources by measuring intensity fluctuations of emission lines tracing
the cosmic gas in varying phases. Using LIMFAST, a novel semi-numerical tool
designed to self-consistently simulate LIM signals of multiple EoR probes, we
investigate how building blocks of galaxy formation and evolution theory, such
as feedback-regulated star formation and chemical enrichment, might be studied
with multi-tracer LIM during the EoR. On galaxy scales, we show that the star
formation law and the feedback associated with star formation can be indicated
by both the shape and redshift evolution of LIM power spectra. For a baseline
model of metal production that traces star formation, we find that lines
highly sensitive to metallicity are generally better probes of galaxy
formation models. On larger scales, we demonstrate that inferring ionized
bubble sizes from cross-correlations between tracers of ionized and neutral
gas requires a detailed understanding of the astrophysics that shape the line
luminosity–halo mass relation. Despite various modeling and observational
challenges, wide-area, multi-tracer LIM surveys will provide important high-
redshift tests for the fundamentals of galaxy formation theory, especially the
interplay between star formation and feedback by accessing statistically the
entire low-mass population of galaxies as ideal laboratories, complementary to
upcoming surveys of individual sources by new-generation telescopes.

††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

The enormous amount of observational and modeling efforts over the past two
decades have revealed an increasingly detailed and precise picture of the
epoch of reionization (EoR). Following the onset of the first galaxy formation
at $z>10$ (Oesch et al., 2016; Naidu et al., 2020; Bouwens et al., 2021;
Harikane et al., 2022) and being completed by $z\approx 5$–6 (McGreer et al.,
2015; Becker et al., 2021; Cain et al., 2021), the neutral intergalactic
medium (IGM) after cosmic recombination was ionized again by an accumulating
background of energetic UV photons emerged from the evolving populations of
early star-forming galaxies (Fan et al. 2006; Stark 2016; Dayal & Ferrara
2018; Robertson 2021; but see also e.g., Qin et al. 2017 for alternative
sources like quasars).

An emerging technique in observational cosmology, line intensity mapping (LIM)
has been widely recognized as a powerful method to study the co-evolution of
galaxies and the IGM during the EoR (Kovetz et al., 2017; Chang et al., 2019).
Particularly, the prospects of synergies among LIM surveys targeting at
different (and usually complementary) tracers have attracted considerable
attention in recent years, as more and more target lines being identified and
observed at wavelengths across the electro-magnetic spectrum. Substantial
theoretical effort has been made in recent years to investigate the scientific
potentials of multi-tracer LIM. One important objective is to employ the
large-scale complementarity between tracers of ionized and neutral regions in
the IGM to tomographically measure the reionization history (e.g., Lidz et
al., 2011; Gong et al., 2012; Feng et al., 2017; Heneka et al., 2017; Dumitru
et al., 2019; Kannan et al., 2022b). Such joint analyses can trace the growth
of ionized regions and alleviate observational challenges like systematics and
foreground contamination. Another major objective is to infer physical
properties of the source population through simultaneous diagnosis of multiple
spectral lines emitted from the multi-phase interstellar medium (ISM) and/or
IGM (e.g., Heneka et al., 2017; Sun et al., 2019; Yang et al., 2021; Schaan &
White, 2021; Bethermin et al., 2022). Even though the coarse-grain averaged
nature of these statistical measurements makes the interpretation challenging
in many circumstances, these efforts have showcased the richness of
astrophysical information about the EoR to be gleaned from multi-tracer LIM
datasets.

Nevertheless, the majority of modeling efforts in the LIM context can be
broadly considered as “semi-empirical”, which leverage a small number of
simple, observationally-motivated trends to describe the source population and
create mock LIM signals. Although these models provide quantitative
expectations of LIM signals anchored to observations, clear physical
connections among properties of the source population and different
observables are often missing (but see Kannan et al., 2022a, b, who employ
fully-detailed, radiation-magneto-hydrodynamic simulations of reionization
with photoionization and radiative transfer modeling to study multi-tracer
LIM). Relatively little effort has been devoted so far to the development of
physical yet efficient modeling frameworks that capture the essential
astrophysical information, while being flexible enough to allow model testing
and inference from multi-tracer LIM datasets. For these reasons, we have
developed LIMFAST, a semi-numerical toolkit for self-consistently simulating a
multitude of LIM signals during the EoR, as introduced in detail in Mas-Ribas
et al. (2022, henceforth Paper I). LIMFAST extends the 21cmFAST code, and
implements significantly improved models of galaxy formation and line emission
in the high-$z$ universe.

In this work, we present a generalization and applications of the basic
framework of LIMFAST introduced in Paper I, by considering physically-
motivated variations of stellar feedback and star formation law prescriptions.
Given the important consequences these variations have on galaxy and IGM
evolution during the EoR, we investigate their implications for a number of
promising LIM targets for probing the EoR, including the 21 cm line of H I and
nebular lines at optical/UV (e.g., H$\alpha$, Ly$\alpha$) and far-infrared
(e.g., [C II], [O III], CO) wavelengths. Such a generalization allows us to
relate specific LIM observables to a fundamental picture of high-$z$ galaxy
formation described by a balance maintained by star formation from the ISM and
stellar feedback typically from supernovae (Furlanetto et al., 2017;
Furlanetto, 2021). By characterizing how these astrophysical processes impact
the summary statistics of different LIM signals, especially the auto- and
cross-power spectra, we investigate how the underlying astrophysics of
feedback-regulated star formation can be better understood from future LIM
observations combining different line tracers.

The remainder of this paper is structured as follows. In Section 2, we briefly
review LIMFAST, including its basic code structure and functionalities. In
Section 3, we specify key features of the galaxy formation model and its
variations considered in this work, namely prescriptions for stellar feedback
and the star formation law. We also introduce a supplement to the baseline
nebula model discussed in Paper I, which allows us to model lines emitted from
the photodissociation regions (PDRs) and molecular gas irradiated by the
interstellar radiation field sourced by star formation. In Section 4, we
present the main quantitative results of this work about how variations of the
galaxy model affect the reionization history, followed by how distinct forms
of feedback and the star formation law may be revealed by multi-tracer LIM
observations. We compare our results with previous work, discuss some
limitations and caveats of the analyses presented, and outline several
possible extensions of the current framework in Section 5, before summarizing
our main conclusions in Section 6. Throughout the paper, we assume
cosmological parameters consistent with recent measurements by Planck
Collaboration XIII (Planck Collaboration et al., 2016).

##  2 The LIMFAST Code

In Paper I, we describe in detail the general features and applications of the
LIMFAST code. Here, we only briefly review the key features of LIMFAST and
refer interested readers to the paper for further details.

Built on top of the 21cmFAST code Mesinger et al. (2011); Park et al. (2019),
LIMFAST shares with it the efficient, semi-numerical configuration adopted to
approximate the formation of the large-scale structure and the partitioning of
mass into dark matter halos. Specifically, the evolution of density and
velocity fields is calculated with the Lagrangian perturbation theory, whereas
the hierarchical formation of structures and the growth of ionized regions are
described by the excursion set formalism Lacey & Cole (1993); Mesinger &
Furlanetto (2007), without resolving individual halos. Using the overdensity
field derived, LIMFAST replaces the simplistic galaxy formation model used in
21cmFAST by a quasi-equilibrium model of high-$z$, star-forming galaxies
introduced by Furlanetto et al. (2017) and Furlanetto (2021), which predicts a
range of physical properties important for LIM studies, including the gas
mass, stellar mass, star formation rate (SFR), metallicity, and so on. The
line intensity fields of interest are then computed by integrating
emissivities predicted by the photoionization and radiative transfer
simulation code, cloudy Ferland et al. (2017), over the halo mass function
conditional on the local overdensity. Following Mesinger et al. (2011), we
normalize integrals of the subgrid conditional halo mass function in the
Press-Schechter formalism (Bond et al., 1991) to match the mean values from
the Sheth-Tormen formalism.

LIMFAST coherently simulates a variety of LIM signals that trace the
reionization and the underlying galaxy formation histories. In Paper I, the
simulated cosmic star formation rate density and the IGM neutral fraction
evolution are verified by comparing against latest observations of high-$z$
galaxies/quasars and the cosmic microwave background (CMB), whereas LIM
statistics of a suite of (mainly optical/UV) nebular lines typically from H II
regions (e.g., Ly$\alpha$, H$\alpha$, [O II], and [O III]) are compared with
other high-$z$ model predictions in the literature. Extensions and variations
of the baseline model presented in Paper I, including an extended model of
emission lines that predominantly originate from the neutral ISM (e.g., [C II]
and CO), are introduced in this work to facilitate the analysis.

##  3 Models

To understand the connection between astrophysics of galaxy formation and LIM
signals originating from different environments, especially the multi-phase
ISM, and demonstrate the astrophysical applications of multi-tracer LIM
studies, we need to consider some plausible variations of the galaxy formation
model, and preferably a large set of line signals sensitive to the variations
and different gas phases. In what follows, we will describe how the baseline
LIMFAST simulation presented in Paper I is extended for such purposes.

###  3.1 Models of Galaxy Formation and Evolution

Following the galaxy formation model described in Paper I based on Furlanetto
et al. (2017) and Furlanetto (2021), the star formation and chemical evolution
of individual halos can be described by a set of simple, coupled ordinary
differential equations.

Expressing any given mass 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 respective to redshift for the
time evolution (i.e., $\tilde{M}^{\prime}=d\tilde{M}/dz$), we can express the
evolution of halo mass, gas mass, stellar mass, and metal mass as

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

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

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

$\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]~{},$
(4)

where $C_{\rm orb}=(1+z_{0})A_{\rm dyn}\Delta^{-1/2}_{\rm vir}$ characterizes
the parameter dependence of the orbital timescale, which is set by some
normalization factor $A_{\rm dyn}$ and the virial overdensity of a collapsed
halo $\Delta_{\rm vir}=18\pi^{2}$. A constant return fraction
$\mathcal{R}=0.25$ is taken to describe the fraction of stellar mass recycled
back to the star-forming gas due to stellar evolution (Benson, 2010; Tacchella
et al., 2018). A metal yield $y_{Z}=0.03$ is adopted, which is appropriate for
metal-poor, Population II (Pop II) stars with a typical initial mass function
(Benson, 2010).

Figure 1: Top: the mass growth histories of a sample halo reaching $M\approx
10^{11}\,M_{\odot}$ at $z=5$ calculated by solving the system of differential
equations from $z_{i}=30$. The black set of curves show halo properties
calculated by the bathtub model with constant $\eta=10$ and $\epsilon=0.015$.
The sets of curves in blue, green, yellow, and red show the results of models
“momentum” (Model Ia), “momentum/KS” (Model Ib), “momentum/FQH” (Model Ic),
and “energy” (Model II), respectively. For reference, the dotted curve in grey
indicates the total accreted baryonic mass, which is more than 10 times larger
than the mass of gas and stars formed as a result of strong feedback
regulation. Bottom: the gas/stellar mass–halo mass relations at $z=7$
predicted by different choices of the feedback mode and star formation law.
Figure 2: The mass and redshift dependence of the instantaneous (top row) and
time-averaged (bottom row) star formation efficiencies calculated using
different feedback prescriptions.

To investigate what astrophysics of galaxy formation may be inferred from LIM
data, we consider a total of 6 model variations involving different
assumptions for the underlying feedback mode and star formation law, which are
described by the value (or functional form) of the mass loading factor,
$\eta$, and the temporal efficiency factor, $\epsilon$, respectively. Our
baseline model assumes that stellar feedback is momentum-driven ($\eta\propto
M^{-2/3}(1+z)^{-1}$) and sets $\epsilon$ to a fiducial value of 0.015
consistent with local observations (Krumholz et al., 2012), which is referred
to as Model Ia and is the fiducial model assumed throughout Paper I. A set of
feedback variations are considered, where we explore a range of feedback modes
leading to different star formation efficiency (SFE), especially in low-mass
halos, while fixing $\epsilon$. In Model II, we consider energy-driven
feedback ($\eta\propto M^{-1/3}(1+z)^{-1/2}$), whereas in Model III (IV) we
envisage a more extreme scenario where the coupling between stellar feedback
and the star-forming gas is weaker (stronger) than the momentum-driven
(energy-driven) case. Specifically, a weak coupling in Model III assumes
$\eta\propto M^{1/6}$, whereas a strong coupling in Model IV assumes
$\eta\propto M$, with the redshift dependence being dropped in both cases for
simplicity.

Table 1: Specifications of the baseline model (Model Ia, adopted for the fiducial simulations presented in Paper I) and its variations considered in this work. The value of the escape fraction is varied accordingly to ensure that reionization completes by $z\approx 5.5$.  Model | Feedback Mode | Star Formation Law | $f_{\rm esc}$
---|---|---|---
Ia | momentum | $\epsilon=0.015$ | 10%
Ib | momentum | KS | 10%
Ic | momentum | FQH | 10%
II | energy | $\epsilon=0.015$ | 12.5%
III | weak | $\epsilon=0.015$ | 5%
IV | strong | $\epsilon=0.015$ | 30%

On the other hand, a set of star formation law variations are explored in
Models Ib and Ic for momentum-driven feedback, where we further allow
$\epsilon$ to vary moderately with the gas mass and thus yield star formation
surface density–gas surface density relations111Throughout this work, we
define star formation law to be the relation between the star formation rate
surface density $(\dot{\Sigma}_{*})$ ans the gas surface density
$(\Sigma_{\mathrm{g}})$, for which we follow Furlanetto (2021) to compute the
galaxy size from the half-mass radius of the host halo, assuming that galaxy
disc has uniform surface density. corresponding to those implied by
observations and/or theoretical predictions. In Model Ib, we adopt
$\epsilon\propto M_{g}^{0.2}$ which reproduces the well-known Kennicutt-
Schmidt law Kennicutt (1998) with a power-law index of 1.4, whereas in Model
Ic we adopt $\epsilon\propto M_{g}^{0.4}$ to approximate the star formation
law with a power-law index of 2 as suggested by Faucher-Giguère et al. (2013),
where the gas disc is assumed to be entirely turbulence-supported by stellar
feedback.

Solving Equations (1)–(4), we obtain the growth histories of stellar and gas
masses in dark matter halos as they continuously accrete from an initial
redshift of $z_{i}=30$. Figure 1 shows the growth histories of gas, stellar,
and metal masses for a sample dark matter halo with $M\approx
10^{11}\,M_{\odot}$ at $z=5$, derived from models with different feedback and
star formation law combinations considered in this work. An averaged star
formation efficiency (SFE) can be defined consequently as the stellar
mass–halo mass ratio, namely $\tilde{f}_{*}=M_{*}/M_{a}$, which can be
interpreted as the time-integrated value of the instantaneous SFE
$f_{*}=\dot{M}_{*}/\dot{M}_{a}$ — both are often derived in the literature
with halo abundance matching (HAM), namely by matching the halo number density
described by the halo mass function to the galaxy number density described by
the galaxy UV luminosity function (e.g., Mason et al., 2015; Mashian et al.,
2016; Sun & Furlanetto, 2016; Furlanetto et al., 2017; Tacchella et al., 2018;
Behroozi et al., 2019). For simplicity, we limit the stellar population
resulting from our star formation prescriptions to Pop II stars only, and
defer a systematic introduction of Population III (Pop III) stars to future
studies (see Section 5.4). Figure 2 shows the instantaneous and time-averaged
SFEs derived in different feedback prescriptions as a function of halo mass
and redshift. The steepness of the color gradient indicates how strongly the
fraction of accreted gas turned into stars is regulated by stellar feedback.

Figure 3: Left: galaxy rest-frame UV luminosity functions under different
feedback assumptions. The predicted luminosity functions are compared against
the observational constraints from Bouwens et al. (2015) and Bouwens et al.
(2021), as represented by the filled squares and empty circles, respectively.
The $z=6$, $z=8$ and $z=9$ luminosity functions are offset vertically by a
multiplicative factor of 2, 0.5 and 0.25, respectively, for ease of
comparison. Right: comparison of the galaxy stellar-to-halo mass ratios
implied by different feedback assumptions to the latest estimates from
observations based on clustering analysis (Harikane et al., 2018) and HAM
(Finkelstein et al., 2015; Stefanon et al., 2021).

In Figure 3, we show a comparison between the observed galaxy UV luminosity
functions and stellar-to-halo mass relation (SHMR) and our model predictions
at $z\gtrsim 6$. As illustrated in the left panel, we verify that luminosity
functions implied by the four feedback prescriptions considered are all
reasonably well-consistent with constraints on the faint end from the Hubble
Space Telescope (HST) data (Bouwens et al., 2015, 2021). In the right panel,
we show that the SHMRs predicted by our model variations are roughly
consistent with observations in the low-mass regime. At the high-mass end, our
predictions only agree well with estimates based on the galaxy clustering
(Harikane et al., 2018), but not those based on HAM (Finkelstein et al., 2015;
Stefanon et al., 2021), which are a factor of 2–3 larger. We note that lots of
these uncertainties associated with EoR galaxies will be greatly reduced by
new-generation telescopes like the James Webb Space Telescope (JWST), but as
will be illustrated in what follows the information from multi-tracer LIM
observations, which cover much wider areas wherein the entire galaxy
population is accessed, will still be extremely valuable and complementary.

###  3.2 A Multi-Phase Extension of the Nebula Model

As described in Paper I, the numerical photoionization code cloudy (version
17.02, Ferland et al., 2017) is supplemented to galaxy properties predicted by
the galaxy formation model in LIMFAST to simulate the production of various
emission lines as target LIM signals. Here, we extend the baseline nebular
model introduced in Paper I, which mainly accounts for lines produced in H II
regions, to include bright emission lines of particular interest to LIM
studies from the neutral (atomic/molecular) ISM, such as [C II] 158 $\mu$m,
and CO(1–0) 2601 $\mu$m lines. We note that because any legitimate nebular
model based on cloudy can be used as the input of LIMFAST, in what follows we
do not repeat the analysis of optical/UV lines discussed in Paper I with the
new nebular model. For H$\alpha$, Ly$\alpha$, and [O III] 88 $\mu$m lines
considered in this work, we simply reuse the results from Paper I, although in
principle they can be captured together with lines originating from the
neutral ISM by a generalized nebular model.

For lines emitted from atomic or molecular gas, namely in PDRs or
$\mathrm{H}_{2}$ cores of gas clouds, because their strengths depend on the
gas content, we define an equivalent surface area $\mathcal{S}$ according to
the distribution of gas density in giant molecular clouds (GMCs) following the
prescription from Vallini et al. (2018). For simplicity, the gas mass $M_{g}$
of a given halo is assumed to be evenly distributed among the population of
GMCs with the same fixed mass $M_{\rm GMC}$, such that the total line
luminosity can be simply scaled from that of one single GMC. Specifically, to
describe the internal structure of GMCs, we first define a volumetric
distribution of gas density $\rho$ that follows a log-normal probability
distribution function (PDF), as suggested by models of isothermal, non-self-
gravitating, turbulent gas (e.g., Passot & Vázquez-Semadeni, 1998; Padoan &
Nordlund, 2002). Namely,

$P_{V}(\rho)\propto\frac{1}{\sqrt{2\pi}\sigma}\exp{\left[-\frac{\left(\ln(\rho/\rho_{0})-\langle\ln(\rho/\rho_{0})\rangle\right)^{2}}{2\sigma^{2}}\right]}~{},$
(5)

where $\rho_{0}=\mu m_{\rm p}n_{\rm H,0}$ is the mean gas density of the GMC,
with $n_{\rm H,0}=100\,\mathrm{cm^{-3}}$ and $\mu=1.36$ accounting for helium.
The logarithmic scatter $\sigma$ satisfies
$\langle\ln(\rho/\rho_{0})\rangle=-\sigma^{2}/2$, maintaining a fixed
expectation of $\ln({\rho/\rho_{0}})$ as $\sigma$ varies. The distribution of
$\sigma$ depends on the turbulence level characterized by the Mach number
$\mathcal{M}$ through

$\sigma^{2}=\ln(1+b^{2}\mathcal{M}^{2})~{},$ (6)

where $b=0.3$ describes the efficiency of turbulence production and we take
$\mathcal{M}=5$, a plausible value for high-redshift galaxies that tend to
have more supersonic structures (see discussion in e.g., Safarzadeh &
Scannapieco, 2016). Since in this work we focus on lines with low to
intermediate critical densities (rather than e.g., high-$J$ CO lines tracing
the densest regions in GMCs), we ignore self-gravity, which has the critical
effect of modifying the density distribution into a power law at high enough
densities. The PDF is normalized such that its integral over gas density gives
the total volume of the GMC,

$V_{\rm GMC}^{\rm tot}=\int dV=\int
P_{V}(\rho|\mathcal{M})d\rho=\frac{4\pi}{3}R^{3}_{\rm GMC}~{},$ (7)

where $R_{\rm GMC}=15\,$pc specifies the size of one GMC, implying a mass of
$M_{\rm GMC}=\rho_{0}V_{\rm GMC}^{\rm tot}=4.7\times 10^{4}\,M_{\odot}$. This
allows one to conveniently define a characteristic length scale corresponding
to each density $\rho$,

$r(\rho)=\left[\int_{\rho-\delta\rho}^{\rho+\delta\rho}P_{V}(\rho^{\prime}|\mathcal{M})d\rho^{\prime}\right]^{1/3}~{},$
(8)

with which an equivalent surface area $\mathcal{S}=4\pi r^{2}$ can be defined
for the cumulative line flux $f_{\rm line}=4\pi j_{\rm line}$ in units of
$\mathrm{erg\,s^{-1}\,cm^{-2}}$. The line luminosity density (i.e., $l_{\rm
line}\equiv dL_{\rm line}/dV$) can be then expressed as

$l_{\rm line}(\rho)=\frac{4\pi r^{2}(\rho)}{V_{\rm GMC}^{\rm tot}}f_{\rm
line}(\rho,Z_{g},U)~{},$ (9)

where $f_{\rm line}(\rho,Z_{g},U)$ means that the line luminosity is a
function of gas density, gas metallicity, and ionization parameter in cloudy.
Figure 4 shows the cumulative emissivities of the three FIR/mm-wave lines
considered in this work, for a gas density of
$n_{\mathrm{H,0}}=100\,\mathrm{cm^{3}}$ and an interstellar radiation field
(ISRF) of strength $G_{0}=10^{3}$ in units of the Habing flux, which we adopt
as fiducial parameters in our model, consistent with observations of the ISM
at high redshifts (Gullberg et al., 2015; Wardlow et al., 2017). The line
luminosity of the GMC is consequently given by

$\mathcal{L}^{\rm tot}_{\rm line}=\int l_{\rm
line}(\rho)P_{V}(\rho|\mathcal{M})\,d\rho~{}.$ (10)

Finally, to arrive at the total line luminosity of a halo of gas mass $M_{g}$,
we simply scale $\mathcal{L}^{\rm tot}_{\rm line}$ by

$L^{\rm tot}_{\rm line}=\mathcal{L}^{\rm tot}_{\rm line}M_{g}/M_{\rm GMC}~{}.$
(11)

Figure 4: Cumulative emissivities of [O III] 88 $\mu$m, [C II] 158 $\mu$m,
CO(1–0) 2601 $\mu$m, and CO(2–1) 1300 $\mu$m lines calculated with cloudy
assuming a gas cloud with density $n_{\mathrm{H}}=100\,\mathrm{cm^{3}}$ and
metallicity $Z=0.002$ illuminated by the interstellar radiation field
$G_{0}=10^{3}$ in the Habing flux. Note that CO lines are enlarged by 1000
times for the ease of comparison. Background colors indicate different ISM
phases where the lines predominantly originate from, with the boundaries
corresponding to where sharp changes in the gas kinetic temperature profile
occur.  Figure 5: The luminosity–star formation rate relations of [O III], [C
II], and CO(1–0) lines at $z\sim 6$ predicted by LIMFAST assuming different
stellar feedback and star formation law prescriptions. Empirical fittings to
observations of low-/intermediate-$z$ galaxies Lidz et al. (2011); De Looze et
al. (2014); Chung et al. (2022) and high-$z$ galaxies Harikane et al. (2020);
Schaerer et al. (2020), as well as predictions from physically-motivated
high-$z$ ISM models Lagache et al. (2018); Muñoz & Furlanetto (2013) from the
literature are shown by the gray lines. Because of the relatively dearth of
high-$z$ [O III] emitters observed, a selected number of recent [O III]
observations at $z>6$ are also shown for comparison (Inoue et al., 2016;
Hashimoto et al., 2019; Laporte et al., 2019; Harikane et al., 2020).

Figure 5 shows the luminosity–star formation rate relations for [O III], [C
II], and CO lines from $z=6$ (solid) and $z=8$ (dashed) galaxies as predicted
by our galaxy models. For comparison, we also plot empirical representations
of the observed/inferred luminosity–star formation rate relations for
low-/intermediate-$z$ galaxies Lidz et al. (2011); De Looze et al. (2014);
Chung et al. (2022), together with high-$z$ fitting relation (Harikane et al.,
2020) and predictions based on physically-motivated ISM models Muñoz &
Furlanetto (2013); Lagache et al. (2018). For [C II] emission, recent
observations of $4\lesssim z\lesssim 6$ galaxies (e.g., Schaerer et al., 2020)
suggest a lack of evolution in the $L_{\mathrm{[CII]}}$–SFR relation from that
in the local universe De Looze et al. (2014), which slightly disfavors some
theoretical predictions derived from a combination of semi-analytical models
of galaxies and photoionization simulations (Lagache et al., 2018). Notably
for CO(1–0) emission, the models displayed from LIMFAST and the literature
differ significantly in terms of the shape and normalization, partly because
of the different redshift regimes where these models are evaluated/calibrated.
Future CO galaxy and LIM surveys with ngVLA and COMAP (see e.g., Breysse et
al., 2022, and references therein) will greatly improve the constraints on
scaling relations like the one shown here.

###  3.3 Emission Lines From the IGM

Besides nebular emission lines produced in the multi-phase ISM, for the IGM
emission LIMFAST also inherits and improves on the detailed 21 cm calculations
from 21cmFAST, while adopting a simple prescription for recombination emission
from the diffuse, ionized IGM in Ly$\alpha$. As detailed in Sections 2.5 and
2.6 of Paper I, where interested readers are referred to for a complete
description, the modeling of these emission lines from the IGM is also fully
coupled with the galaxy formation model and its variations implemented in
LIMFAST. This allows the influence of galaxy astrophysics on the statistics of
21 cm and IGM Ly$\alpha$ emission to be studied self-consistently with other
emission lines from the ISM as tracers of the galaxy distribution.

Figure 6: Left: the redshift evolution of the cosmic SFRD (top) and the mean
neutrality of the IGM (bottom) simulated by LIMFAST models with different
feedback assumptions. Also shown in the top panel are observational
constraints on the SFRD from Oesch et al. (2018) integrated down to $M_{\rm
UV}=-17$, along with empirical fits from Robertson et al. (2015) and Harikane
et al. (2022) assuming different amounts of extrapolation for the faint
populations. Observational constraints on the mean IGM neutral fraction from
the dark fraction in the Ly$\alpha$ and Ly$\beta$ forests (McGreer et al.,
2015), IGM damping wing signatures in quasar spectra (Davies et al., 2018),
and Ly$\alpha$ emission from Lyman-break galaxies (Mason et al., 2018, 2019)
are shown in the bottom panel. Curves of different feedback mode are also
labeled by an asymmetry measure, $A_{\mathrm{s}}$, of the EoR history, and the
implied CMB electron scattering optical depth, $\tau_{e}$, which is verified
to be consistent with the latest observational constraints (Pagano et al.,
2020; Qin et al., 2020). Right: the bubble size distribution derived with the
“mean free path” approach when the IGM is approximately half-ionized.

##  4 Results

In this section, we present the main quantitative results of this paper
derived from the set of model variations specified in Table 1. We first show
the global reionization histories implied by models with varying feedback
prescriptions (Section 4.1), which supplements the model predictions presented
in Paper I based on Model Ia only, and then present how the corresponding sky-
averaged signals of various lines are sensitive to the changes in feedback.
Next, we demonstrate how variations of the feedback mode (Section 4.3) and the
star formation law (Section 4.4) in play affect summary statistics, namely the
auto- and cross-correlation power spectra of tracers of neutral and ionized
IGM. By examining the shape and amplitude evolution of power spectra, we
elaborate on how astrophysical information about ionizing sources and the IGM
may be extracted in turn from joint analyses of multi-tracer LIM observations.
For clarity, all results presented in the remainder of this section are shown
in real space, without considering observational effects such as redshift
space distortions (RSDs), whose treatment in LIMFAST is elaborated in Section
2.7 of Paper I.

###  4.1 Reionization Histories

Figure 7: The redshift evolution of the sky-averaged (global) 21 cm different
brightness temperature (left) and intensities of [C II], H$\alpha$ and CO
lines (middle) under different feedback assumptions. Also shown in the right
panel is the ratio of redshift dependence of different line intensities, which
serves as a measure of the feedback-sensitive metal enrichment history.

Effects of feedback regulation on star-forming galaxies as ionizing sources
can be revealed by both the global history and the detailed morphology of
reionization. Figure 6 shows two important measures of the reionization, the
volume-averaged neutral fraction and the ionized bubble size distribution
(BSD), simulated by LIMFAST assuming different feedback prescriptions, which
yield the different redshift evolution of the cosmic SFRDs shown in the upper
left panel. Notably, the cosmic SFRD directly relates to the strength of
stellar feedback through the SFE in low-mass halos. As illustrated in Figure
2, more efficient feedback coupling results in a steeper SFE gradient with
halo mass, implying less efficient star formation in low-mass halos and thus
an overall lower and steeper cosmic SFRD dominated by massive halos. The SFRDs
predicted by our momentum- and energy-regulated models are comparable to the
extrapolation to observations from Robertson et al. (2015) out to $z\sim 15$
(despite the opposite curvature), whereas stronger or weaker feedback can
result in SFRDs close to or substantially higher ($>1$ dex) than the
observational constraints available to date (Oesch et al., 2018; Harikane et
al., 2022). Since we tune $f_{\mathrm{esc}}$ such that the reionization
completes roughly at the same time at $z\approx 6$ in each feedback scenario
(see Table 1), a steeper SFRD corresponds to an overall more rapid and
asymmetric reionization history, as measured by the factor
$A_{\mathrm{s}}\equiv(z_{05}-z_{50})/(z_{50}-z_{95})$ which uses the
reionization midpoint $z_{50}$ and 5% (95%) completion point $z_{05}$
($z_{95}$) to characterize the asymmetry of the full extent of the EoR (Trac,
2018).

The impact of stellar feedback on the size of ionized regions is illustrated
in the right panel of Figure 6, which shows the BSD in different feedback
modes when the IGM is about half-ionized. Following the “mean free path”
method introduced by Mesinger & Furlanetto (2007), we describe the BSD with
the probability density function of the logarithmic bubble radius $R$, which
is calculated by repeatedly sampling the size of H II regions from random
ionized points and in random directions with a Monte Carlo process. At a fixed
neutral fraction $\langle x_{\rm H\textsc{i}}\rangle$, the BSD shifts towards
larger bubble radius when the feedback regulation is stronger. Although a
mass-dependent $f_{\mathrm{esc}}$ may lead to a similar effect, the degeneracy
can be reduced by direct constraints on $f_{\mathrm{esc}}$ from future
observations of individual EoR galaxies by e.g., JWST. As will be shown in
what follows, cross-correlations between the 21 cm signal and tracers of star-
forming galaxies turn out to be sensitive probes of the typical bubble size
encoded by the BSD, even though the exact correspondence relies on a good
understanding of the astrophysics.

###  4.2 Sky-Averaged Line Intensities

The sky-averaged intensity of spectral line emission, especially that of the
21 cm line (often referred to as the global 21 cm signal), as a spatial
monopole measurement is known to be a useful probe of the EoR history and
source population (Mirocha et al., 2015, 2017; Cohen et al., 2017; Mirocha &
Furlanetto, 2019; Park et al., 2019). Figure 7 shows the redshift evolution of
the sky-averaged signals of various lines and their ratios. For the 21 cm
global signal, $\delta\bar{T}_{\mathrm{b}}(z)$, as revealed by the timing and
strength of its extrema, a stronger feedback implies delayed Ly$\alpha$
coupling and heating, which lead to an overall smaller signal amplitude in
both absorption ($\delta\bar{T}_{\mathrm{b}}<0$) at cosmic dawn and emission
($\delta\bar{T}_{\mathrm{b}}>0$) during cosmic reionization. The absorption
trough varies between $z\sim 12$ and 18 in the central redshift and between
$\delta\bar{T}_{\mathrm{b}}=-80$ and $-120$ mK in the depth for the four
feedback modes considered, suggesting an intimate connection between feedback-
regulated star formation in the first galaxies and the 21 cm spin temperature
evolution during cosmic dawn. The is consistent with the overall shift of the
21 cm global signal towards lower redshift/higher frequency, as projected by
the recent literature taking into account of the observed UV luminosity
function at $z\gtrsim 6$ (Mirocha et al. 2017; Park et al. 2019; but see also
Cohen et al. 2017).

Figure 8: Left: a comparison between the cosmic SFRD attributed to galaxies
measurable by LIM experiments and a nominal JWST/NIRCAM ultra-deep survey
reaching a limiting magnitude of $m\approx 32$ with four 200 h pointings
covering a total area of $\approx 0.01\,\mathrm{deg^{2}}$. The thick and thin
sets of curves represent LIM- and JWST-detectable galaxies, respectively.
Right: the fractional SFRD deviation from the fiducial momentum-driven
feedback model (Model Ia) in other feedback models.

The middle panel of Figure 7 shows the mean intensities of line tracers of
galaxies [C II], H$\alpha$, and CO, whose redshift evolution is largely driven
by that of the SFRD. Nonetheless, the subtle difference in the steepness of
redshift evolution, caused by the different metallicity dependence of these
tracers, serves as a potential probe of feedback through the implicit metal
enrichment history. The effect is illustrated in the right panel of Figure 7,
where ratios of slopes $\bar{I}^{-1}d\bar{I}/dz$ as a function of redshift are
contrasted with each other. Clearly, the slope ratio of lines with high
contrast in metallicity dependence (e.g., [C II] and H$\alpha$) is sensitive
to feedback, with less efficient feedback producing a slope ratio with
stronger redshift evolution, whereas the slope ratio of lines with similar
metallicity dependence (e.g., [C II] and CO) is largely a constant insensitive
to the exact feedback mechanism. Such sensitivity to feedback is not
surprising though, given that the mean intensity evolution is mainly driven by
the much more abundant low-mass galaxies that are most feedback-sensitive and
least metal-enriched. In the case of very strong feedback, there are too few
low-mass galaxies to make a significant variation in the slope ratio, and thus
it remains roughly constant with redshift.

Measurements of the mean intensity evolution of spectral line tracers of
galaxies, especially SFR tracers like H$\alpha$ and [C II], can often be
translated into constraints on the SFRD, provided that the $L$–SFR relation
can be reliably determined. This in turn provides an angle to compare the
information from LIM observations to what will be available from forthcoming
surveys of individual high-$z$ sources by new-generation telescopes like the
JWST. In Figure 8, we illustrate a simple comparison between the
distinguishing power on the cosmic SFRD in different feedback modes available
from JWST and LIM observations in general. As an example, we consider an
potential ultra-deep (UD) configuration for a galaxy dropout survey with
JWST/NIRCAM that reaches a limiting magnitude of $m\approx 32$, similar to the
strategies considered in Mason et al. (2015) and Furlanetto et al.
(2017)222Even though more realistic survey plans are now available (see e.g.,
Williams et al., 2018; Robertson, 2021), the approximate configuration, as
presented, is sufficient for our purpose.. As shown in the left panel of
Figure 8, a JWST UD survey will likely still miss a considerable fraction
($\gtrsim 50\%$) of the total star formation in galaxies at $z\gtrsim 6$,
unless the SFRD declines steeply with redshift due to a shallow faint-end
slope of the galaxy luminosity function, a likely result of very efficient
stellar feedback. On the contrary, the statistical nature of LIM makes the
measurements sensitive to the collective star formation activity sourcing the
aggregate line emission, although in some cases the conversion between line
luminosity and the SFR can be sophisticated.

Figure 9: Left: slices of 21 cm differential brightness temperature $\delta
T_{b}$ fields at various stages of reionization ($\langle x_{\rm
H\textsc{i}}\rangle\approx 0.2$, 0.5, and 0.8, respectively) simulated by
LIMFAST, assuming different stellar feedback prescriptions as specified in
Table 1. Right: slices of H$\alpha$, Ly$\alpha$, and [C II] line intensity
fields at $z\approx 7$ ($\langle x_{\rm H\textsc{i}}\rangle\approx 0.5$)
simulated by LIMFAST. Note that the Ly$\alpha$ intensity fields displayed have
contributions from both Ly$\alpha$ emitters and recombinations in the diffuse
ionized IGM, though without accounting for the damping wing absorption due to
the intervening neutral IGM. Each slice is 256 Mpc on a side and 1 Mpc thick.

To further contrast the two types of measurements in the context of probing
stellar feedback, we show in the right panel of Figure 8 the fractional
deviation of the total, measurable SFRD in other feedback modes from that in
the fiducial, momentum-driven case of feedback. Namely,
$|\dot{\rho}_{*}-\dot{\rho}^{\mathrm{Ia}}_{*}|/\dot{\rho}^{\mathrm{Ia}}_{*}$
reflects how easily one might disprove a simple momentum-driven feedback model
using deviations (if any) of the observed SFRD from the expected one. Due to
the insufficient sensitivity of galaxy surveys to faint objects, for which the
effect of feedback regulation is most pronounced, a JWST UD survey tends to
have less distinguishing power than LIM observations especially towards higher
redshift. The exception, again, is when comparing a very strong feedback to
Model Ia, in which case the difference in distinguishing power decreases with
increasing redshift as galaxies to which LIM is uniquely sensitive diminish
rapidly, although LIM observations still offer more distinguishing power. We
note that the example presented here only represents an extremely-simplified,
special case of inferring stellar feedback from the SFRD evolution. In
practice, the individual source detection and LIM methods further complement
each other by the different quantities that are directly probed (e.g., the
luminosity function vs. moments of the luminosity function) and the different
sources of uncertainty involved (e.g., cosmic variance and the dust correction
vs. foreground contamination and the $L$–SFR relation), and thus are both
valuable probes of galaxy formation and evolution in the high-$z$ universe.

###  4.3 Characterizing Stellar Feedback With LIM

In the left panel of Figure 9, we show slices of $\delta T_{b}$ fluctuations
in different feedback scenarios at various stages of reionization, when the
average IGM neutral fraction is $\langle x_{\rm H\textsc{i}}\rangle=0.2$, 0.5,
and 0.8, respectively. It is obvious that at a given stage, the typical size
of ionized regions is on average larger with stronger feedback. This is
because under stronger feedback regulation, star formation tends to occur in
more massive halos, which are more clustered and have a higher ionization rate
to ionize a larger volume of gas thanks to their higher SFR. Note that
although across each row the volume filling factor of fully ionized regions
appears higher in the case of a stronger feedback, the volume-averaged neutral
fraction $\langle x_{\rm H\textsc{i}}\rangle$ of individual simulated boxes
are in fact comparable — because the product of ionization efficiency $\zeta$
and local collapse fraction $f_{\mathrm{coll}}(x,z,R)$ is more evenly
distributed when feedback is less efficient, more partially ionized cells with
ionized fraction equal to $\zeta^{-1}f_{\mathrm{coll}}(x,z,R_{\mathrm{cell}})$
are allowed to exist (Mesinger et al., 2011), which compensate for the deficit
in fully ionized regions.

In Section 4.2, we have demonstrated the impact of feedback on the history of
cosmic dawn and reionization eras as revealed by the 21 cm global signal from
the neutral IGM. Using LIMFAST, we supplement such a picture with the
complementary LIM signals of UV/optical and far-infrared nebular emission
lines tracing star-forming galaxies, which are considered to provide the
majority of ionizing photons required to complete the reionization by
$z\approx 5.5$ (Robertson et al., 2015; Naidu et al., 2020).

The right panel of Figure 9 shows slices through the boxes of H$\alpha$,
Ly$\alpha$, and [C II] intensity fluctuations simulated by LIMFAST when
$\langle x_{\rm H\textsc{i}}\rangle\approx 0.5$ in different feedback
scenarios, color-coded by the logarithmic line intensity in units of
$\mathrm{Jy/sr}$. In contrast to maps of the 21 cm signal, these line
intensity maps generally trace the ionizing sources residing in overdense
regions, whose spatial anti-correlation with the 21 cm signal is clearly
visible on scales larger than the typical size of ionized regions. On finer
scales, information about sources of line emission and the luminosity
distribution of the source population is encoded in detailed features of the
intensity fluctuations. For Ly$\alpha$, a spatially-extended component is
apparent, especially in the case of a strong feedback, which is sourced by
recombinations in the diffuse ionized IGM surrounding ionizing sources. We
note that, as discussed in Paper I, we present intrinsic line intensity fields
throughout, without including the attenuation effect due to dust grains in the
ISM or patches of neutral hydrogen in the IGM. Such effects can be readily
incorporated via post-processing the simulation boxes that LIMFAST outputs, as
have been demonstrated in the literature (Silva et al., 2013; Heneka et al.,
2017).

Figure 10: The redshift evolution at a fixed comoving scale of
$k=0.2\,h/\mathrm{Mpc}$ (top) and scale dependence (bottom) of auto- and
cross-power spectra (absolute value) between [C II] emission from galaxies and
the 21 cm line in different feedback models. The cross-power spectrum changes
sign from negative on large scales to positive on small scales at
$k_{\mathrm{trans}}\sim 2\,h/\mathrm{Mpc}$. Note the different format in which
the [C II] auto-power spectra are plotted as a function of $k$ to facilitate
visual comparison.

On the other hand, the fact that the [C II] intensity field shows a larger
spatial gradient compared with that of H$\alpha$ indicates that the former is
preferentially sourced by more massive and therefore more biased sources. This
results from the difference in the luminosity–halo mass ($L$–$M$) relation of
the two lines. As demonstrated in Appendix A, where we contrast the $L$–$M$
relation of several nebular lines under varying assumptions of the metallicity
dependence, [C II] luminosity is a steeper function of halo mass compared with
lines like H$\alpha$ due to its much stronger metallicity dependence. The
paucity of contribution from low-mass halos not only leads to a steep
$\bar{I}(z)$ evolution shown in Figure 7, but also implies that [C II]
intensity fluctuations will be more dominated by the Poisson noise from rare
[C II]-bright sources, as can be shown by the power spectrum. As we will see,
a sensitive luminosity–halo mass relation makes statistical measurements of
lines like [C II] and CO promising ways of testing models of galaxy formation
involving different feedback assumptions.

####  4.3.1 Information From Auto-Power Spectra

The statistical information about spatial fluctuations of a given LIM signal
is directly available from its auto-power spectrum. As an example, we
illustrate in Figure 10 the power spectra of 21 cm and [C II] lines, whereas
similar illustrations for the power spectra of other nebular lines considered
(and their cross-correlations with the 21 cm line) are provided in Appendix B.
The left two columns of Figure 10 show the shape and redshift evolution of [C
II] and H I 21 cm power spectra calculated from simulations boxes in the four
cases of feedback considered. Even though the power spectrum only partially
describes these potentially highly non-gaussian fields, it is encouraging to
see that useful information about the feedback mode in play can be probed by
either the shape or amplitude evolution of the auto-power spectrum.

From the redshift evolution of the power spectrum at $k=0.2\,h/\mathrm{Mpc}$
shown in the top row, it is evident that the amplitude of large-scale
fluctuations encapsulates statistics of the key drivers for the sky-averaged
signal evolution. In the case of 21 cm power spectrum amplitude, the three
characteristic peaks (from high $z$ to low $z$) corresponds to the eras of
Ly$\alpha$ coupling, X-ray heating, and reionization, when high-amplitude
fluctuations in $\delta T_{\rm b}$ are concurrent with rapid changes (i.e.,
steep slopes) in the 21 cm global signal as shown in Figure 7. Different
feedback prescriptions modulate these peaks in significantly different ways,
with strong feedback yielding peaks later and more squeezed in redshift and
less contrasted in amplitude.

The redshift evolution of [C II] power spectrum amplitude, on the other hand,
largely reflects the sky-averaged intensity $\bar{I}(z)$ evolution of the
signal, which in turn traces the cosmic SFRD evolution as discussed in Section
4.2. Thanks to the quadratic dependence $\Delta^{2}(k)\propto\bar{I}^{2}$,
different feedback modes become more distinguishable, provided that the power
spectrum amplitude can be monitored over a wide enough redshift range.

From the shape of power spectra at $z=6$ and 9 shown in the bottom row, it is
also straightforward to see the modulation effect by feedback. For 21 cm power
spectrum, the impact of feedback on the ionized bubble size is manifested by
the shift of the scale at which the power spectrum peaks, which is most
discernible at $z=6$ ($k\sim 0.1\,h/\mathrm{Mpc}$ for the “strong” model and
$k\sim 0.4\,h/\mathrm{Mpc}$ for the “weak” model) when different feedback
models predict similar $\langle x_{\rm H\textsc{i}}\rangle$ but different
BSDs. The difference appears to be a lot smaller at $z=9$ when the remains
close to high, except for the case of very strong feedback which shows a
qualitative difference from other cases. This is because the 21 cm spin
temperature field is dominated by highly-biased, massive sources in the
presence of strong feedback, thereby showing a distinctive large-scale power
excess at $k\sim 0.1\,h/\mathrm{Mpc}$ due to the source clustering. The
generally much lower amplitude on smaller scales in case, compared with other
feedback cases, is due to the delayed reionization by strongly-suppressed
cosmic star formation.

The shape evolution of [C II] power spectrum is much more subtle in the plot,
although the effect of feedback can still be inferred from the shape contrast
between two redshifts. Overall, stronger feedback leads to both a steeper
$\Delta^{2}_{\mathrm{CII}}(k)$ that is more dominated by the small-scale
Poisson noise and a stronger shape evolution. Considering a metric of the
power spectrum shape contrast,
$\mathcal{X}(k_{1},z_{1},k_{2},z_{2})=\Delta^{2}(k_{1},z_{2})/\Delta^{2}(k_{2},z_{2})-\Delta^{2}(k_{1},z_{1})/\Delta^{2}(k_{2},z_{1})$,
which characterizes the change in the dominance of small-scale Poisson noise
in the power spectrum between two redshifts $z_{1}$ and $z_{2}$, we find
$\mathcal{X}_{\mathrm{CII}}(3\,h/\mathrm{Mpc},6,0.1\,h/\mathrm{Mpc},9)=20$,
40, 65, and 140 for the “weak”, “momentum”, “energy”, and “strong” models,
respectively. Here $k=3$ and $0.1\,h/\mathrm{Mpc}$ roughly correspond to the
smallest and largest scales accessed by our simulation, and a larger, positive
$\mathcal{X}$ indicates that from $z=6$ to 9 the “strong” model implies a
larger increase in the dominance of the Poisson noise contribution. Such a
correlation between $\mathcal{X}$ and the feedback strength is insensitive to
factors that only affect the power spectrum normalization, and thus marks a
potentially useful application of auto-correlation analysis to the
understanding of galaxy formation physics.

####  4.3.2 Information From Cross-Power Spectra

In practice, measurements of the auto-power spectrum are often unfortunately
complicated by a variety of astrophysical and instrumental effects. One of the
main obstacles is foreground contamination, which can overwhelm the target LIM
signal by several orders of magnitude. Even though a multitude of cleaning
techniques have been devised to remove foreground contamination of various
origins, cross-correlating signals with uncorrelated foregrounds still has its
unique advantages. Therefore, it is interesting to understand how cross-
correlations between different lines, especially the 21 cm line and nebular
lines tracing ionizing sources, may be leveraged to characterize the effect of
feedback in high-$z$ galaxy formation.

In the rightmost column of Figure 10, we compare different feedback models by
showing how their predicted 21 cm–[C II] cross-power spectra evolve with
redshift in their amplitude and shape. From the redshift evolution, the
strength of feedback determines how rapidly the cross-power amplitude evolves.
Moreover, thanks to the counteractive evolution of 21 cm and [C II] amplitudes
with redshift, the peaks intrinsic to the 21 cm contribution become broadened
for $\Delta^{2}_{\rm H\,\textsc{i},C\,\textsc{ii}}(z)$ compared with
$\Delta^{2}_{\rm H\,\textsc{i}}(z)$ (especially for the peak due to X-ray
heating), and the extent of the broadening depends on how long the
counteractive effect persists. Different feedback modes are therefore easier
to be distinguished by $\Delta^{2}_{\rm H\,\textsc{i},C\,\textsc{ii}}(z)$,
wherein less rapid evolution with broad, flattened, and later peaks in
redshift corresponding to weaker feedback. From the shape of the cross-power
spectrum, on the other hand, the most pronounced feature is the dependence of
the scale at which the cross-power changes sign on the feedback mode in play.
As will be shown next, both feedback and the physics of line emission affect
the interpretation of such a characteristic scale, which has been perceived as
an indicator of the typical size of ionized bubbles during the EoR.

Figure 11: Top: cross-correlation coefficients between the 21 cm line and
H$\alpha$, Ly$\alpha$, [C II], and CO lines derived from the maps simulated by
LIMFAST at $\langle x_{\rm H\textsc{i}}\rangle\approx 0.2$, assuming energy-
driven (Model II, blue set of points) and strong (Model IV, red set of points)
feedback. Bottom: the evolution of the comoving transition scale, at which
$r_{\mathrm{line\times 21cm}}(k_{\mathrm{trans}})=0$, with the mean IGM
neutral fraction. Scales inaccessible by our simulation outputs of limited
resolution are greyed out, and $k^{-1}_{\mathrm{trans}}$ values below which
are marked and interpreted as upper limits.

We further inspect the effect of feedback on the cross-correlation signals in
Figure 11 by comparing the cross-correlation coefficient $r_{1\times
2}(k)=P_{\mathrm{1\times 2}}(k)/\sqrt{P_{1}(k)P_{2}(k)}$ of the 21 cm signal
with a variety of emission-line tracers of galaxies, including H$\alpha$,
Ly$\alpha$, [C II], and CO lines. In particular, we focus on how the scale
dependence of $r(k)$ differs for different cross-correlations assuming
different feedback assumptions, especially the transition scale
$k_{\mathrm{trans}}$ where $r(k_{\mathrm{trans}})=0$. In the top panel of
Figure 11, $r(k)$ of different feedback models and line tracers when $\langle
x_{\rm H\textsc{i}}\rangle\approx 0.2$ are shown in different hues and tints,
respectively. The fact that, at a fixed stage (i.e., $\langle x_{\rm
H\textsc{i}}\rangle$) of reionization, stronger feedback predicts faster de-
correlation between 21 cm and nebular lines as $k$ increases (e.g., from
$0.1\,h/\mathrm{Mpc}$ to $1\,h/\mathrm{Mpc}$) is consistent with the more
“top-heavy” BSD skewed towards larger bubble sizes expected in this case.
Nonetheless, in the bottom panel of Figure 11, we show that
$k_{\mathrm{trans}}$ is only modestly sensitive to feedback (and thus the BSD)
for a given line tracer, although, as expected, it indeed traces the
macroscopic progress of reionization described by $\langle x_{\rm
H\textsc{i}}\rangle$.

Figure 12: Power spectra of auto-correlations (top row) of Ly$\alpha$, [C II],
and CO lines, as well as their cross-correlations with the H I 21 cm signal
(bottom row) during the EoR predicted by the three LIMFAST models assuming the
same momentum-driven feedback but different star formation laws.

Another noteworthy feature in Figure 11 is the discrepancies in the individual
cross-correlations for a given feedback mode. Unlike naively expected, there
are non-trivial differences in both $r(k)$ and $k_{\mathrm{trans}}$ among
different line tracers of galaxies. Lines like [C II] and CO tracing the
neutral ISM, whose luminosities evolve steeply with mass due to e.g., their
strong dependence on the gas metallicity, exhibit a modestly lower level of
(negative) correlation with the 21 cm line, when compared with tracers of the
ionized ISM less sensitive to metallicity, such as H$\alpha$ and Ly$\alpha$.
This, in turn, makes $k_{\mathrm{trans}}$ vary among different nebular lines
with essentially different effective bias factors for a fixed
feedback/reionization scenario. For example, $k_{\mathrm{trans}}^{-1}$, as a
proxy for the bubble size, can differ by more than 50% at $\langle x_{\rm
H\textsc{i}}\rangle\sim 0.3$ depending on whether H$\alpha$ or CO is cross-
correlated with the 21 cm line. Similar effects have been noted previously by
several other authors (Dumitru et al., 2019; Kannan et al., 2022b; Cox et al.,
2022), despite using less explicit formulations of the connection between
nebular line emission and galaxy formation. Finally, for Ly$\alpha$, we note
that a qualitatively different trend appears for $k_{\mathrm{trans}}^{-1}$ as
a function of $\langle x_{\rm H\textsc{i}}\rangle$, which is caused by the
additional diffuse component from recombinations in the diffuse ionized IGM
(see Figure 9) that can strongly modulate $k_{\mathrm{trans}}^{-1}$ especially
towards the end of the EoR.

###  4.4 Characterizing the Star Formation Law With LIM

Besides stellar feedback, the other way that the astrophysics of galaxy
formation can affect the luminosity–halo mass relation of nebular lines is
through the star formation law. In particular, because the star formation law
only alters the relative gas content of galaxies instead of the amount of star
formation, as illustrated in Figure 1, lines originating from the neutral ISM
are most sensitive to changes in the star formation law. We note, though, that
measurements of large-scale structure using LIM signals of star-formation
lines from H II regions can still be useful probes of the star formation law
across cosmic time (Sun, 2022).

Figure 12 shows the auto-power spectra of Ly$\alpha$, [C II], and CO lines and
their cross-power spectra with the 21 cm signal during the EoR. Clearly, the
statistics of Ly$\alpha$, whose luminosity simply scales with the star
formation rate, are little affected by using different forms of the star
formation law (Models Ia, Ib, and Ic). Even though in principle there can be a
small indirect effect through the different metallicities implied by different
star formation laws, it is barely visible from the comparison of Ly$\alpha$
power spectra. On the contrary, given the substantial difference in the
luminosity–halo mass (or SFR) relation caused by the gas mass dependence (see
Figure 5), [C II] and CO lines have power spectra varying significantly with
the assumed star formation law in both the shape and amplitude. Models
implying more efficient star formation out of the gas reservoir and therefore
a shallower luminosity–SFR relation, e.g., Model Ic, tend to yield auto-power
spectra less dominated by the Poisson noise and evolving less rapidly with
redshift — an unsurprising result given that large-scale fluctuations are
mainly contributed by the more abundant fainter sources, whereas the small-
scale Poisson fluctuations dominating are mainly contributed by the very rare
and bright sources. For reference, we find
$\mathcal{X}_{\mathrm{CII}}(3\,h/\mathrm{Mpc},6,0.1\,h/\mathrm{Mpc},9)=41$,
25, 17, and
$\mathcal{X}_{\mathrm{CO}}(3\,h/\mathrm{Mpc},6,0.1\,h/\mathrm{Mpc},9)=55$, 36,
26 for Model Ia, Ib, and Ic, respectively, suggesting that indeed Model Ic
assuming the FQH13 star formation law predicts power spectra the least Poisson
noise-dominated.

Figure 13: Cross-correlation coefficients between the 21 cm line and
Ly$\alpha$, [C II], and CO lines and derived from the maps simulated by
LIMFAST assuming the default (lines, Model Ia), KS (filled markers, Model Ib),
and FQH13 (empty markers, Model Ic) star formation law.

The cross-power spectrum between the 21 cm and [C II] or CO lines also
exhibits a clear dependence on the star formation law assumed, even though the
reionization scenario is largely independent of it. With a steeper star
formation law (i.e., more efficient gas to stellar mass conversion), the
intensity field of [C II] or CO becomes less dominated by bright sources and
therefore de-correlates with the 21 cm field at smaller scales, causing a
noticeable shape difference potentially useful for testing star formation law
models. Moreover, the FQH13 model (Model Ic) predicts the strongest redshift
evolution of the cross-power amplitude over $6<z<9$, again due to the weaker
counteractive evolution of 21 cm and [C II] or CO lines in the case of a
steeper star formation law.

In Figure 13, we show the cross-correlation coefficients, $r(k)$, between
Ly$\alpha$, [C II], CO lines and the 21 cm line at $z=6$ and $z=9$ under
various assumptions of the star formation law. Since the reionization scenario
is nearly insensitive to changes in the star formation law, any difference in
$r(k)$ shown in this figure is due to the nebular line intensity signal rather
than the 21 cm signal. Several interesting features are noteworthy. First, as
expected, $r_{\mathrm{Ly\alpha\times 21cm}}(k)$ remains almost unchanged in
different star formation law models because Ly$\alpha$ only depends on the
SFR. Second, similar to what is shown in Figure 11, for either [C II] or CO
line the level of (negative) correlation at a given scale depends moderately
on the star formation law assumed, with steeper star formation law yielding a
less rapid de-correlation as $k$ increases. Lastly, a change in the relative
order of $r(k)$ for Ly$\alpha$, [C II], and CO lines is observed by
contrasting Model Ib with Model Ic, which may be utilized for star formation
law model selection. Similar to stellar feedback, the star formation law also
serves as a source of complications in the interpretation of typical ionized
bubble size from $r(k)$ or $k_{\mathrm{trans}}$.

##  5 Discussion

In what follows, we compare the results and their implications from this work
against some previous literature, and discuss potential caveats and
limitations of our methods. In addition, we also outline several promising
directions to extend the current framework of LIMFAST in the future.

###  5.1 Comparison to Previous Work

A number of studies have previously studied and demonstrated the huge
potential of LIM observations targeting at different tracers for understanding
the cosmic dawn and reionization eras. By developing LIMFAST, we provide an
efficient modeling framework to self-consistently simulate a large number of
LIM signals during the EoR that have been investigated individually (or in
small subsets) before by different authors, such as H$\alpha$ (e.g., Heneka et
al., 2017; Silva et al., 2018; Heneka & Cooray, 2021), Ly$\alpha$ (e.g., Silva
et al., 2013; Pullen et al., 2014; Heneka et al., 2017), [C II] (e.g., Gong et
al., 2012; Chung et al., 2020; Sun et al., 2021b), [O III] (e.g., Padmanabhan
et al., 2021), and CO (e.g., Lidz et al., 2011; Mashian et al., 2015; Breysse
et al., 2022). As has been demonstrated in Paper I, the heterogeneous
assumptions made by different studies about physics of the ISM, star
formation, feedback, and the metal and dust content often make direct
comparisons between distinct line signals or distinct models of the same
signal challenging and difficult to interpret. While qualitative comparisons
may still reveal interesting astrophysical information, without quantitative
assessments of the discrepancies observed it is unlikely to reliably test and
compare different models against data. This urges the need to be able to
describe and forecast various target LIM signals during the EoR — usually
differing in both the natal phase of gas and the connection to galaxy
properties — with a unified picture of high-$z$ galaxy formation. Nonetheless,
as we show in Paper I, non-trivial offsets exist between our results and other
individual, line-specific models involving vastly varying assumptions of the
galaxy population and spectral line production. Thus, coherently modeling the
otherwise disconnected physical conditions of multiple emission lines and
galaxy evolution, using tools like LIMFAST, is imperative to understand and
exploit the multi-tracer LIM technique for studying the EoR.

On the usage of the scale $k_{\mathrm{trans}}$ at which the cross-correlation
coefficient between the 21 cm signal and a given galaxy tracer changes sign,
our findings are qualitatively similar to previous analyses by Lidz et al.
(2011), Dumitru et al. (2019), and most recently Kannan et al. (2022b). Put
briefly, the general redshift evolution of $k_{\mathrm{trans}}$ does reflect
the overall progress of the reionization as measured by $\langle x_{\rm
H\textsc{i}}\rangle$, but such evolution is complicated by uncertainties of
the source population that affect signals of both tracers being cross-
correlated. Specifically, at any given $\langle x_{\rm H\textsc{i}}\rangle$,
variations of our galaxy model in either feedback or the star formation law
can modulate $k_{\mathrm{trans}}$ through of the BSD and/or the effective bias
of the galaxy tracer. While previous analyses often adopt a sharp dichotomy of
halo emissivities in terms of ionizing photon production to distinguish
between reionization scenarios dominated by faint vs. bright sources (e.g.,
Dumitru et al., 2019; Kannan et al., 2022b), our model allows galaxies of
different luminosities to more smoothly impact both the neutral gas and line
intensity distributions in a consistent manner. Such smooth transitions in the
contribution from different sources to signatures of reionization are not
merely more realistic, but also essential for shedding light on how high-$z$
galaxies driving the reionization might be shaped by the balance between star
formation and feedback.

###  5.2 Limitations of the Galaxy Formation Model

In LIMFAST, we have implemented and leveraged the simple, quasi-equilibrium
model of high-$z$ galaxy formation described in Furlanetto et al. (2017) and
Furlanetto (2021) to study the impact of the astrophysics of galaxies on
various target LIM signals. Although it already represents an improvement over
the source modeling in the latest release of 21cmFAST Murray et al. (2020) in
aspects such as the physical connection between star formation and feedback
regulation, some intrinsic limitations of the method need to be noted and are
likely worthy of further exploration in future work.

A key assumption made in our galaxy formation model is that in the high-$z$
universe a quasi-equilibrium state can already be established by proto-
galaxies in the form of a settled disc where stars steadily form. Making this
assumption provides a neat way to describe the formation of EoR galaxies by
analogy to their low-$z$ counterparts, though one may question how valid such
a scenario can be in the highly dynamic and uncertain stage of early galaxy
formation. Recently studies, including a follow-up study to the Furlanetto
(2021) model, have shown that star formation might be highly bursty during the
early phase of galaxy formation, before some critical mass is reached and
stars can steadily form. For example, Furlanetto & Mirocha (2022) generalize
the quasi-equilibrium disc model by introducing a non-trivial perturbation
arising from the time delay between star formation and stellar feedback at
high redshifts. Numerical simulations also find strong evidence for strongly
time-variable star formation in early, low-mass galaxies before a
rotationally-supported ISM emerges from a rapid process of disc settling
(Gurvich et al., 2022), which turns out to be supported by Galactic
archaeology of the in situ, metal-poor component of the Milky Way’s stellar
halo (Belokurov & Kravtsov, 2022), indicating a potential requirement for
full, non-equilibrium approaches. Given the intimate connection between star
formation and spectral line emission in galaxies, as demonstrated in this
work, it is crucial to quantify in future studies the effects of highly time-
variable star formation on multi-tracer LIM observations of the EoR.

Even if the quasi-equilibrium model indeed approximates the formation and
evolution of high-$z$ star-forming galaxies well, it is admittedly simplistic
in many ways, some of which are closely related to subgrid modeling that will
be discussed in the next sub-section. One important simplification is
associated with the diversity of galaxy formation histories. As demonstrated
by Mirocha et al. (2021), simple subgrid, HAM-based models tend to produce
biased signatures of the reionization process, when compared against fully
numerical methods accounting for both halo mergers and the stochasticity of
the halo mass accretion rate. A hybrid or numerically-calibrated approach will
therefore be useful for further improvements in the model accuracy (see also
Section 5.4). Relatedly, we have also neglected the scatter in astrophysical
parameters of our galaxy formation model, which can impact LIM signals of
interest in a non-trivial way (Shekhar Murmu et al., 2021; Reis et al., 2022)
and therefore should be taken into account in future development of LIMFAST by
e.g., cell-level stochastic sampling of astrophysical parameters. Pop III
stars are another missing piece of the current model that can have non-trivial
effects on the EoR, whose physical properties and formation histories may be
studied either through their influence on the 21 cm signal (Mirocha et al.,
2018; Mebane et al., 2020; Qin et al., 2021b; Muñoz et al., 2022) or by
mapping the emission of nebular lines characteristic of Pop III stars, such as
the He II 1640 Å line (Visbal et al., 2015; Parsons et al., 2021). While
extensions of 21cmFAST-like, semi-numerical simulations have attempted to
self-consistently model the formation of Pop III and Pop II stars altogether
(e.g., Tanaka & Hasegawa, 2021; Muñoz et al., 2022), observational
constraints, either direct or indirect, are pivotal to the down-selection of
the poorly constrained model space (see e.g., Mirocha et al., 2018; Sun et
al., 2021a).

###  5.3 Uncertainties With Subgrid Astrophysics

We note that a range of simplifications and model assumptions are made for the
subgrid astrophysics of galaxy formation and evolution, which are essential
for the application of LIMFAST to the EoR science, but in the meantime serve
as important sources of uncertainty. For instance, the galaxy properties
captured by our quasi-equilibrium model are highly simplified, which in turn
limits how closely galaxy evolution and the production of the various kinds of
spectral line emission can be modeled coherently. In particular, several
physical conditions of the stellar population and the ISM must be specified
manually, such as the star formation history of galaxies, the gas density, the
interstellar radiation field strength, and the dust content and properties,
all of which are likely influential for the modeling of both galaxy evolution
and the LIM signals of interest (e.g., Lagache et al., 2018; Mirocha, 2020;
Mirocha et al., 2021; Yang et al., 2021).

Given all the aforementioned uncertainties, as well as those mentioned in
Section 5.3 about the galaxy formation model, accurately computing the line
emission and eventually applying the model to reverse-engineer the properties
of galaxies from upcoming cosmological surveys of the EoR will be a non-
trivial task. Insights from observations and detailed numerical simulations of
the mechanisms behind and the connections among different ingredients of
subgrid astrophysics, such as the co-evolution of gas, metals, and dust across
cosmic time (Li et al., 2019), the connection between the ionization parameter
and metallicity (Ji & Yan, 2022), and the effects of non-equilibrium
photoionization and metal cooling (Katz, 2022), will be extremely valuable for
better understanding and further improvements of the source modeling in semi-
numerical tools like LIMFAST, especially for the applications of interpreting
future observations.

###  5.4 Extension of the Current Framework

In Paper I and this work, we present the current structure and functionalities
of LIMFAST focusing on its capability of forward modeling the multi-tracer LIM
observations of the EoR. It is useful to note that the current framework may
be readily extended in various promising ways and applied to a broader range
of EoR studies, thanks to the modular nature of LIMFAST.

First, additional probes of the EoR can be incorporated into the same modeling
framework in a consistent manner similar to the existing ones. For instance,
several authors have demonstrated that semi-numerical simulations are ideal
tools for studying the kinetic Sunyaev-Zel’dovich (kSZ) effect from patchy
reionization and its synergy with the 21 cm signal for constraining the
reionization history (Battaglia et al., 2013; La Plante et al., 2020; Gorce et
al., 2022). Cross-correlating the kSZ signal derived from the simulated
ionization and velocity fields with line tracers of galaxies provides the
redshift information missing in kSZ measurements. Similar ideas can be applied
to other types broad-band, two-dimensional datasets such as the CMB lensing
and the cosmic near-infrared background, through the large-scale fluctuations
of which rich information about the population of ionizing sources may be
extracted (Helgason et al., 2016; Maniyar et al., 2022; Sun et al., 2021a;
Mirocha et al., 2022). Furthermore, as demonstrated already in the low-$z$
universe, three-dimensional Ly$\alpha$ forest tomography serves as a promising
probe of the large-scale distribution of the neutral IGM (Lee et al., 2018;
Newman et al., 2020), which can be ideally suited for studying the late stages
of the reionization process by itself or in combination with LIM datasets (Qin
et al., 2021a). It is interesting to implement these additional observables
into LIMFAST to quantitatively assess their potential for probing the EoR,
especially when jointly analyzed with LIM observations, and examine methods
required for overcoming observational challenges like foreground contamination
(e.g., Zhu et al., 2018; Gagnon-Hartman et al., 2021).

Besides taking into account extra probes of the EoR, it is also of interest to
extend LIMFAST further into the post-reionization universe ($0<z<5$). Galaxies
during this age of active assembly and evolution are not only interesting by
themselves but also important witnesses of the impact of reionization on
galaxy formation, which will be studied by a number of forthcoming LIM surveys
of galaxies at low-to-intermediate redshift, such as COMAP (Cleary et al.,
2021), EXCLAIM (Cataldo et al., 2020), SPHEREx (Doré et al., 2018), and TIM
(Vieira et al., 2020). That said, despite showing great promise, the low-$z$
extension of LIMFAST faces two main challenges. First, at lower redshift, the
halo occupation distribution (HOD) becomes more sophisticated due to the
increased population of satellite galaxies (Kravtsov et al., 2004; Bhowmick et
al., 2018; Behroozi et al., 2019), and quenching becomes a more and more
important process in galaxy formation and evolution (Tal et al., 2014; Brennan
et al., 2015; Donnari et al., 2021). Both factors call for more detailed
subgrid models for the luminosity–halo mass relation. Meanwhile, accurately
modeling the partitioning of mass into halos becomes more challenging at lower
redshift due to the increased importance of halo mergers. LIMFAST inherits the
formulation of large-scale structure and radiation field approximation from
21cmFAST, where the generation of halo source fields by a halo finding
algorithm is bypassed. To properly account for halo merger histories in the
low-$z$ extension, an explicit halo finding algorithm, with either an extended
dynamic range to resolve small halos at the cooling threshold, or an enhanced
subgrid modeling of halo source fields involving merger trees and a stochastic
population of simulation cells with unresolved halos, will be required at the
cost of extra RAM capacity and a slower speed (see discussions in e.g.,
Mesinger & Furlanetto, 2007; Mesinger et al., 2011, and references therein).

##  6 Conclusions

Using simulations generated by the LIMFAST code introduced in Paper I, we have
presented in this paper a unified picture of how the astrophysics of high-$z$
galaxy formation affect and therefore can be reveal by multi-tracer LIM
observations of the EoR. We investigate the impact of different stellar
feedback and star formation law prescriptions on a variety of signatures of
reionization, including the 21 cm signal and LIM signals of nebular emission
lines from the multi-phase ISM, such as H$\alpha$, Ly$\alpha$, [O III], [C
II], and CO. Our main findings can be summarized as follows:

  1. 1.

Because the cosmic star formation history is sensitive to feedback-regulated
star formation in individual galaxies, the efficiency of stellar feedback
directly impacts the history, geometry, and thereby the variety of
observational signatures of the reionization process. On the other hand,
besides a small indirect effect through the metallicity, the star formation
law only affects tracers of the neutral ISM of galaxies as indirect probes of
the reionization.

  2. 2.

The redshift evolution of multiple sky-averaged line signals already serves as
a useful probe of the astrophysics of high-$z$ galaxy formation. Timings of
the extrema in the 21 cm global signal are tightly connected to the feedback
efficiency through radiation fields scaling with the cosmic SFRD. Due to the
strong metallicity dependence of metal cooling lines like [C II], a comparison
between their sky-averaged signal evolution and that of hydrogen lines like
H$\alpha$ can inform the (cosmic mean) stellar feedback strength.

  3. 3.

Rich astrophysical information about the reionization and its driving sources
can be extracted from the auto- and cross-power spectra of multi-tracer LIM
data. Both feedback and the star formation law can modulate the shape and
amplitude of power spectra and their evolution across cosmic time. Power
spectral analyses combining multiple, complementary tracers therefore allows
cross-checks and the separation of effects due to the reionization itself and
those associated with galaxy formation and evolution.

  4. 4.

The cross-correlation between the 21 cm line and a spectral line tracer of
galaxies is particularly useful for tracing the overall progress of the EoR.
However, even though the transition scale $k_{\mathrm{trans}}$ roughly probes
the neutral fraction evolution, the exact interpretation and implications of
the cross-correlation are subject to complications due to astrophysics of
galaxy formation and the resulting properties of galaxies, and thus dependent
on the specific tracer considered. Multi-tracer LIM makes it possible to
better understand how LIM signals are influenced by astrophysical processes
such as feedback and the star formation law, on which the usage of
$k_{\mathrm{trans}}$ or the cross-correlation analysis in general is premised.

  5. 5.

By accessing a larger fraction of the faint galaxy population than individual
source detection, LIM surveys can use the inferred SFRD to offer more
sensitive tests for processes central to galaxy formation like the stellar
feedback. This makes LIM a highly complementary method for studying high-$z$
galaxy formation even in the era of new-generation telescopes.

In summary, there is great potential for multi-tracer LIM to transform our
understanding of cosmic reionization and the formation and evolution of
high-$z$ galaxies that drive the reionization process. In spite of the various
challenges that commonly exist in practice for different tracers, such as the
mismatch of scales and issues of foreground contamination (see the review
e.g., Liu & Shaw, 2020), careful coordination and optimization for future
multi-tracer synergies will eventually allow the invaluable astrophysical
information to be extracted and applied to tests of the galaxy formation
theory at high redshift. Reliable semi-numerical simulations like LIMFAST, in
its current and future forms, are essential tools for accurately modeling and
analyzing the vast amount of observational data to come.

Figure 14: The $L$–$M$ relations of H$\alpha$, [C II], [O III], and CO(1–0)
lines assuming different metallicity values. The solid curve assumes a varying
metallicity predicted by our galaxy model (Model Ia), whereas the dashed and
dotted curves assume fixed metallicity values of $Z=0.003$ and $Z=0.0001$,
respectively.

Figure 15: Same as Figure 10, but for CO(1–0), [O III], H$\alpha$, and
Ly$\alpha$ lines. Note the different format in which the auto-power spectra
are plotted as a function of $k$ to facilitate visual comparison.

We thank Adam Lidz, Fred Davies, Jordan Mirocha, Rahul Kannan, and Yuxiang Qin
for helpful discussions and comments that helped improve this paper, as well
as Mauro Stefanon for sharing the data of the stellar-to-halo mass ratios. We
acknowledge support from the JPL R&TD strategic initiative grant on line
intensity mapping. Part of this work was done at Jet Propulsion Laboratory,
California Institute of Technology, under a contract with the National
Aeronautics and Space Administration.

##  Appendix A Luminosity–Halo Mass Relation

The $L$–$M$ relation directly dictates the way the underlying matter density
field is traced by the line intensity map observed, as has been discussed in
previous studies (e.g., Kannan et al., 2022b). Varying behaviors of the
$L$–$M$ relation are therefore essential for understanding the properties and
statistics of the spectral line tracer in different variations of the galaxy
model.

Figure 14 shows a comparison of the $L$–$M$ relation of different lines,
assuming either a fixed metallicity or a varying metallicity as predicted by
our galaxy model. The similar shapes of curves with a fixed metallicity
suggests that the scaling with the SFR (e.g., H$\alpha$ and [O III]) or the
gas mass (e.g., [C II] and CO) barely affects the $L$–$M$ relation in the case
of the simple star formation law assumed for Model Ia (and Models II–IV). It
is really the metallicity dependence and evolution that result in different
$L$–$M$ relations of hydrogen and metal lines, which in turn lead to the
different effective bias factors of these spectral line tracers. It is also
interesting to note that the $L$–$M$ relation is coupled to the metallicity
not only through the metal content of the ISM, but also through the
metallicity dependence of the ISRF, whose synthetic spectrum from BPASS v2.1
(Eldridge et al., 2017) is supplied to the cloudy simulations (see Paper I).
For instance, a more ionizing ISRF for a lower metallicity produces modestly
brighter H$\alpha$ emission, and partially counteracts the effect of a more
metal-poor ISM towards lower halo masses for metal lines highly sensitive to
the ionizing radiation like [O III].

##  Appendix B Power Spectra of Nebular and 21 cm Lines With Varying Feedback

Supplementing Figure 10 which uses [C II] and 21 cm lines to exemplify the
ways feedback affects the redshift evolution and scale dependence of LIM power
spectra, we further in Figure 15 similar results for other nebular lines
considered in this work, including [O III], CO(1–0), H$\alpha$, and Ly$\alpha$
(sum of star formation and IGM contributions). Comparing the evolution with
redshift and scale of different tracers shows interesting (though in some
cases subtle) trends that inform about how these lines are sensitive to
different aspects of the galaxy evolution such as the gas content and
metallicity, which may be systematically probed by combining LIM surveys of
all these different line tracers.

## References

  * Battaglia et al. (2013) Battaglia, N., Natarajan, A., Trac, H., Cen, R., & Loeb, A. 2013, ApJ, 776, 83, doi: 10.1088/0004-637X/776/2/83
  * Becker et al. (2021) Becker, G. D., D’Aloisio, A., Christenson, H. M., et al. 2021, MNRAS, 508, 1853, doi: 10.1093/mnras/stab2696
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Belokurov & Kravtsov (2022) Belokurov, V., & Kravtsov, A. 2022, arXiv e-prints, arXiv:2203.04980.  https://arxiv.org/abs/2203.04980
  * Benson (2010) Benson, A. J. 2010, Phys. Rep., 495, 33, doi: 10.1016/j.physrep.2010.06.001
  * Bethermin et al. (2022) Bethermin, M., Gkogkou, A., Van Cuyck, M., et al. 2022, arXiv e-prints, arXiv:2204.12827.  https://arxiv.org/abs/2204.12827
  * Bhowmick et al. (2018) Bhowmick, A. K., Campbell, D., Di Matteo, T., & Feng, Y. 2018, MNRAS, 480, 3177, doi: 10.1093/mnras/sty2128
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440, doi: 10.1086/170520
  * 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
  * Brennan et al. (2015) Brennan, R., Pandya, V., Somerville, R. S., et al. 2015, MNRAS, 451, 2933, doi: 10.1093/mnras/stv1007
  * Breysse et al. (2022) Breysse, P. C., Chung, D. T., Cleary, K. A., et al. 2022, ApJ, 933, 188, doi: 10.3847/1538-4357/ac63c9
  * Cain et al. (2021) Cain, C., D’Aloisio, A., Gangolli, N., & Becker, G. D. 2021, ApJ, 917, L37, doi: 10.3847/2041-8213/ac1ace
  * Cataldo et al. (2020) Cataldo, G., Ade, P. A. R., Anderson, C. J., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11445, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1144524, doi: 10.1117/12.2561069
  * Chang et al. (2019) Chang, T.-C., Beane, A., Dore, O., et al. 2019, BAAS, 51, 282.  https://arxiv.org/abs/1903.11744
  * Chung et al. (2020) Chung, D. T., Viero, M. P., Church, S. E., & Wechsler, R. H. 2020, ApJ, 892, 51, doi: 10.3847/1538-4357/ab798f
  * Chung et al. (2022) Chung, D. T., Breysse, P. C., Cleary, K. A., et al. 2022, ApJ, 933, 186, doi: 10.3847/1538-4357/ac63c7
  * Cleary et al. (2021) Cleary, K. A., Borowska, J., Breysse, P. C., et al. 2021, arXiv e-prints, arXiv:2111.05927.  https://arxiv.org/abs/2111.05927
  * Cohen et al. (2017) Cohen, A., Fialkov, A., Barkana, R., & Lotem, M. 2017, MNRAS, 472, 1915, doi: 10.1093/mnras/stx2065
  * Cox et al. (2022) Cox, T. A., Jacobs, D. C., & Murray, S. G. 2022, MNRAS, 512, 792, doi: 10.1093/mnras/stac486
  * 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
  * Dayal & Ferrara (2018) Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1, doi: 10.1016/j.physrep.2018.10.002
  * De Looze et al. (2014) De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62, doi: 10.1051/0004-6361/201322489
  * Donnari et al. (2021) Donnari, M., Pillepich, A., Nelson, D., et al. 2021, MNRAS, 506, 4760, doi: 10.1093/mnras/stab1950
  * Doré et al. (2018) Doré, O., Werner, M. W., Ashby, M. L. N., et al. 2018, arXiv e-prints, arXiv:1805.05489.  https://arxiv.org/abs/1805.05489
  * Dumitru et al. (2019) Dumitru, S., Kulkarni, G., Lagache, G., & Haehnelt, M. G. 2019, MNRAS, 485, 3486, doi: 10.1093/mnras/stz617
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Fan et al. (2006) Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415, doi: 10.1146/annurev.astro.44.051905.092514
  * Faucher-Giguère et al. (2013) Faucher-Giguère, C.-A., Quataert, E., & Hopkins, P. F. 2013, MNRAS, 433, 1970, doi: 10.1093/mnras/stt866
  * Feng et al. (2017) Feng, C., Cooray, A., & Keating, B. 2017, ApJ, 846, 21, doi: 10.3847/1538-4357/aa7ff1
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385\.  https://arxiv.org/abs/1705.10877
  * Finkelstein et al. (2015) Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015, ApJ, 814, 95, doi: 10.1088/0004-637X/814/2/95
  * Furlanetto (2021) Furlanetto, S. R. 2021, MNRAS, 500, 3394, doi: 10.1093/mnras/staa3451
  * Furlanetto & Mirocha (2022) Furlanetto, S. R., & Mirocha, J. 2022, MNRAS, 511, 3895, doi: 10.1093/mnras/stac310
  * Furlanetto et al. (2017) Furlanetto, S. R., Mirocha, J., Mebane, R. H., & Sun, G. 2017, MNRAS, 472, 1576, doi: 10.1093/mnras/stx2132
  * Gagnon-Hartman et al. (2021) Gagnon-Hartman, S., Cui, Y., Liu, A., & Ravanbakhsh, S. 2021, MNRAS, 504, 4716, doi: 10.1093/mnras/stab1158
  * Gong et al. (2012) Gong, Y., Cooray, A., Silva, M., et al. 2012, ApJ, 745, 49, doi: 10.1088/0004-637X/745/1/49
  * Gorce et al. (2022) Gorce, A., Douspis, M., & Salvati, L. 2022, arXiv e-prints, arXiv:2202.08698.  https://arxiv.org/abs/2202.08698
  * Gullberg et al. (2015) Gullberg, B., De Breuck, C., Vieira, J. D., et al. 2015, MNRAS, 449, 2883, doi: 10.1093/mnras/stv372
  * Gurvich et al. (2022) Gurvich, A. B., Stern, J., Faucher-Giguère, C.-A., et al. 2022, arXiv e-prints, arXiv:2203.04321.  https://arxiv.org/abs/2203.04321
  * Harikane et al. (2018) Harikane, Y., Ouchi, M., Ono, Y., et al. 2018, PASJ, 70, S11, doi: 10.1093/pasj/psx097
  * Harikane et al. (2020) Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020, ApJ, 896, 93, doi: 10.3847/1538-4357/ab94bd
  * Harikane et al. (2022) Harikane, Y., Ono, Y., Ouchi, M., et al. 2022, ApJS, 259, 20, doi: 10.3847/1538-4365/ac3dfc
  * Hashimoto et al. (2019) Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71, doi: 10.1093/pasj/psz049
  * Helgason et al. (2016) Helgason, K., Ricotti, M., Kashlinsky, A., & Bromm, V. 2016, MNRAS, 455, 282, doi: 10.1093/mnras/stv2209
  * 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
  * Inoue et al. (2016) Inoue, A. K., Tamura, Y., Matsuo, H., et al. 2016, Science, 352, 1559, doi: 10.1126/science.aaf0714
  * Ji & Yan (2022) Ji, X., & Yan, R. 2022, A&A, 659, A112, doi: 10.1051/0004-6361/202142312
  * Kannan et al. (2022a) Kannan, R., Garaldi, E., Smith, A., et al. 2022a, MNRAS, 511, 4005, doi: 10.1093/mnras/stab3710
  * Kannan et al. (2022b) Kannan, R., Smith, A., Garaldi, E., et al. 2022b, MNRAS, doi: 10.1093/mnras/stac1557
  * Katz (2022) Katz, H. 2022, MNRAS, 512, 348, doi: 10.1093/mnras/stac423
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Kovetz et al. (2017) Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, arXiv e-prints, arXiv:1709.09066.  https://arxiv.org/abs/1709.09066
  * Kravtsov et al. (2004) Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004, ApJ, 609, 35, doi: 10.1086/420959
  * Krumholz et al. (2012) Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69, doi: 10.1088/0004-637X/745/1/69
  * La Plante et al. (2020) La Plante, P., Lidz, A., Aguirre, J., & Kohn, S. 2020, ApJ, 899, 40, doi: 10.3847/1538-4357/aba2ed
  * Lacey & Cole (1993) Lacey, C., & Cole, S. 1993, MNRAS, 262, 627, doi: 10.1093/mnras/262.3.627
  * Lagache et al. (2018) Lagache, G., Cousin, M., & Chatzikos, M. 2018, A&A, 609, A130, doi: 10.1051/0004-6361/201732019
  * Laporte et al. (2019) Laporte, N., Katz, H., Ellis, R. S., et al. 2019, MNRAS, 487, L81, doi: 10.1093/mnrasl/slz094
  * Lee et al. (2018) Lee, K.-G., Krolewski, A., White, M., et al. 2018, ApJS, 237, 31, doi: 10.3847/1538-4365/aace58
  * Li et al. (2019) Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425, doi: 10.1093/mnras/stz2684
  * 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
  * Maniyar et al. (2022) Maniyar, A. S., Schaan, E., & Pullen, A. R. 2022, Phys. Rev. D, 105, 083509, doi: 10.1103/PhysRevD.105.083509
  * Mashian et al. (2016) Mashian, N., Oesch, P. A., & Loeb, A. 2016, MNRAS, 455, 2101, doi: 10.1093/mnras/stv2469
  * Mashian et al. (2015) Mashian, N., Sternberg, A., & Loeb, A. 2015, J. Cosmology Astropart. Phys, 2015, 028, doi: 10.1088/1475-7516/2015/11/028
  * 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. (2018) Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, 856, 2, doi: 10.3847/1538-4357/aab0a7
  * Mason et al. (2019) Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947, doi: 10.1093/mnras/stz632
  * McGreer et al. (2015) McGreer, I. D., Mesinger, A., & D’Odorico, V. 2015, MNRAS, 447, 499, doi: 10.1093/mnras/stu2449
  * 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, MNRAS, 411, 955, doi: 10.1111/j.1365-2966.2010.17731.x
  * Mirocha (2020) Mirocha, J. 2020, MNRAS, 499, 4534, doi: 10.1093/mnras/staa3150
  * 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. (2015) Mirocha, J., Harker, G. J. A., & Burns, J. O. 2015, ApJ, 813, 11, doi: 10.1088/0004-637X/813/1/11
  * Mirocha et al. (2021) Mirocha, J., La Plante, P., & Liu, A. 2021, MNRAS, 507, 3872, doi: 10.1093/mnras/stab1871
  * 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 & Furlanetto (2013) Muñoz, J. A., & Furlanetto, S. R. 2013, MNRAS, 435, 2676, doi: 10.1093/mnras/stt1480
  * Muñoz et al. (2022) Muñoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657, doi: 10.1093/mnras/stac185
  * Murray et al. (2020) Murray, S., Greig, B., Mesinger, A., et al. 2020, The Journal of Open Source Software, 5, 2582, doi: 10.21105/joss.02582
  * Naidu et al. (2020) Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109, doi: 10.3847/1538-4357/ab7cc9
  * Newman et al. (2020) Newman, A. B., Rudie, G. C., Blanc, G. A., et al. 2020, ApJ, 891, 147, doi: 10.3847/1538-4357/ab75ee
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105, doi: 10.3847/1538-4357/aab03f
  * Oesch et al. (2016) Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129, doi: 10.3847/0004-637X/819/2/129
  * 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
  * Padoan & Nordlund (2002) Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870, doi: 10.1086/341790
  * Pagano et al. (2020) Pagano, L., Delouis, J. M., Mottet, S., Puget, J. L., & Vibert, L. 2020, A&A, 635, A99, doi: 10.1051/0004-6361/201936630
  * Park et al. (2019) Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933, doi: 10.1093/mnras/stz032
  * Parsons et al. (2021) Parsons, J., Mas-Ribas, L., Sun, G., et al. 2021, arXiv e-prints, arXiv:2112.06407.  https://arxiv.org/abs/2112.06407
  * Passot & Vázquez-Semadeni (1998) Passot, T., & Vázquez-Semadeni, E. 1998, Phys. Rev. E, 58, 4501, doi: 10.1103/PhysRevE.58.4501
  * 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
  * 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. (2021a) Qin, Y., Mesinger, A., Bosman, S. E. I., & Viel, M. 2021a, MNRAS, 506, 2390, doi: 10.1093/mnras/stab1833
  * Qin et al. (2021b) Qin, Y., Mesinger, A., Greig, B., & Park, J. 2021b, MNRAS, 501, 4748, doi: 10.1093/mnras/staa3408
  * Qin et al. (2020) Qin, Y., Poulin, V., Mesinger, A., et al. 2020, MNRAS, 499, 550, doi: 10.1093/mnras/staa2797
  * Qin et al. (2017) Qin, Y., Mutch, S. J., Poole, G. B., et al. 2017, MNRAS, 472, 2009, doi: 10.1093/mnras/stx1909
  * Reis et al. (2022) Reis, I., Barkana, R., & Fialkov, A. 2022, MNRAS, 511, 5265, doi: 10.1093/mnras/stac411
  * Robertson (2021) Robertson, B. E. 2021, arXiv e-prints, arXiv:2110.13160.  https://arxiv.org/abs/2110.13160
  * Robertson et al. (2015) Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19, doi: 10.1088/2041-8205/802/2/L19
  * Safarzadeh & Scannapieco (2016) Safarzadeh, M., & Scannapieco, E. 2016, ApJ, 832, L9, doi: 10.3847/2041-8205/832/1/L9
  * Schaan & White (2021) Schaan, E., & White, M. 2021, J. Cosmology Astropart. Phys, 2021, 068, doi: 10.1088/1475-7516/2021/05/068
  * Schaerer et al. (2020) Schaerer, D., Ginolfi, M., Béthermin, M., et al. 2020, A&A, 643, A3, doi: 10.1051/0004-6361/202037617
  * Shekhar Murmu et al. (2021) Shekhar Murmu, C., Olsen, K. P., Greve, T. R., et al. 2021, arXiv e-prints, arXiv:2110.10687.  https://arxiv.org/abs/2110.10687
  * Silva et al. (2018) Silva, B. M., Zaroubi, S., Kooistra, R., & Cooray, A. 2018, MNRAS, 475, 1587, doi: 10.1093/mnras/stx3265
  * 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
  * Stark (2016) Stark, D. P. 2016, ARA&A, 54, 761, doi: 10.1146/annurev-astro-081915-023417
  * Stefanon et al. (2021) Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29, doi: 10.3847/1538-4357/ac1bb6
  * Sun (2022) Sun, G. 2022, ApJ, 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. (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
  * Tacchella et al. (2018) Tacchella, S., Bose, S., Conroy, C., Eisenstein, D. J., & Johnson, B. D. 2018, ApJ, 868, 92, doi: 10.3847/1538-4357/aae8e0
  * Tal et al. (2014) Tal, T., Dekel, A., Oesch, P., et al. 2014, ApJ, 789, 164, doi: 10.1088/0004-637X/789/2/164
  * Tanaka & Hasegawa (2021) Tanaka, T., & Hasegawa, K. 2021, MNRAS, 502, 463, doi: 10.1093/mnras/stab072
  * Trac (2018) Trac, H. 2018, ApJ, 858, L11, doi: 10.3847/2041-8213/aabff0
  * Vallini et al. (2018) Vallini, L., Pallottini, A., Ferrara, A., et al. 2018, MNRAS, 473, 271, doi: 10.1093/mnras/stx2376
  * Vieira et al. (2020) Vieira, J., Aguirre, J., Bradford, C. M., et al. 2020, arXiv e-prints, arXiv:2009.14340.  https://arxiv.org/abs/2009.14340
  * Visbal et al. (2015) Visbal, E., Haiman, Z., & Bryan, G. L. 2015, MNRAS, 450, 2506, doi: 10.1093/mnras/stv785
  * Wardlow et al. (2017) Wardlow, J. L., Cooray, A., Osage, W., et al. 2017, ApJ, 837, 12, doi: 10.3847/1538-4357/837/1/12
  * Williams et al. (2018) Williams, C. C., Curtis-Lake, E., Hainline, K. N., et al. 2018, ApJS, 236, 33, doi: 10.3847/1538-4365/aabcbb
  * Yang et al. (2021) Yang, S., Somerville, R. S., Pullen, A. R., et al. 2021, ApJ, 911, 132, doi: 10.3847/1538-4357/abec75
  * Zhu et al. (2018) Zhu, H.-M., Pen, U.-L., Yu, Y., & Chen, X. 2018, Phys. Rev. D, 98, 043511, doi: 10.1103/PhysRevD.98.043511