#  NeutralUniverseMachine: An Empirical Model for the Evolution of H i and H2
Gas in the Universe

Hong Guo Shanghai Astronomical Observatory, Chinese Academy of Sciences,
Shanghai 200030, China<EMAIL_ADDRESS>Jing Wang Kavli Institute for
Astronomy and Astrophysics, Peking University, Beijing 100871, China  Michael
G. Jones Steward Observatory, University of Arizona, 933 N Cherry Ave.,
Tucson, AZ 85721, USA  Peter Behroozi Steward Observatory, University of
Arizona, 933 N Cherry Ave., Tucson, AZ 85721, USA  Division of Science,
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan

###### Abstract

Accurately modeling the cold gas content in the universe is challenging for
current theoretical models. We propose a new empirical model
https://halos.as.arizona.edu/UniverseMachine/DR1/Gas_Masses_NeutralUniverseMachine/
(catalog NeutralUniverseMachine) for the evolution of H i and H2 gas along
with dark matter halos based on the
http://halos.as.arizona.edu/UniverseMachine/DR1/ (catalog UniverseMachine)
catalog. It is able to accurately describe the observed H i and H2 mass
functions, molecular-to-atomic ratio, H i-halo mass relation, H i/H2-stellar
mass relations at $z\sim 0$, as well as the evolution of cosmic gas densities
$\rho_{\textrm{H {i}}}$ and $\rho_{\rm H_{2}}$ in $0<z<6$. The predictions
from our model include: (i) There is weak evolution of H i mass function in
$0<z<3$, but the evolution of H2 mass function is much stronger at the massive
end. (ii) The average H i and H2 masses at a given stellar mass decrease by
around 1 dex since $z=3$ for the star-forming galaxies, but the evolution for
the quenched galaxies is much weaker. (iii) Star-forming galaxies have varying
H i depletion time $\tau_{\textrm{H {i}}}$ from 0.1 Gyr to 10 Gyr, and the
dependence of $\tau_{\textrm{H {i}}}$ on stellar mass and redshift is much
stronger than those of H2 depletion time. The quenched galaxies have much
longer gas depletion time and weaker redshift evolution. (iv) The cosmic
baryon density associated with galaxies is dominated by stars for $z<1.2$ and
mainly contributed by H i gas at higher redshifts. (v) The H i bias gradually
increases with the redshift from 0.69 to 2.33 in $0<z<3$ and is consistent
with recent H i intensity mapping experiments.

##  1 Introduction

Cold neutral gas, both in the atomic and molecular phases, is a key component
of the cosmic baryon cycle (see e.g. Tumlinson et al., 2017; Péroux & Howk,
2020, for reviews). As fuel for star formation, cold gas is fully involved in
the complex baryon physics of gas accretion, heating, and cooling, and stellar
and active galactic nucleus (AGN) feedback in the galactic ecosystem (e.g.,
Tacconi et al., 2020; Saintonge & Catinella, 2022). The distribution and
evolution of cold gas thus provide essential information about the formation
histories of galaxies, as well as their coevolution with the host dark matter
halos.

The distribution of atomic gas in the local universe has been well mapped with
the H i 21 cm surveys, such as the H i Parkes All-Sky Survey (HIPASS; Barnes
et al., 2001; Meyer et al., 2004), the Arecibo Fast Legacy ALFA Survey
(ALFALFA; Giovanelli et al., 2005; Haynes et al., 2011), the HI Nearby Galaxy
Survey (THINGS; Walter et al., 2008), the GALEX Arecibo SDSS Survey (GASS;
Catinella et al., 2010) and the Apertif survey (Adams et al., 2022). At higher
redshifts, the detection of individual galaxies becomes increasingly difficult
because of faint signals and radio frequency interference. The H i spectra
stacking and intensity mapping techniques have been applied to infer the
cosmic H i density at $z<1$ (e.g., Lah et al., 2007; Delhaize et al., 2013;
Masui et al., 2013; Rhee et al., 2013, 2018; Kanekar et al., 2016; Bera et
al., 2019). At even higher redshifts of $z>1.5$, the cosmic H i density can be
observed through the damped Ly$\alpha$ (DLA) systems (Wolfe et al., 2005)
arising from background quasar absorption lines (e.g., Péroux et al., 2003;
Prochaska et al., 2005; Prochaska & Wolfe, 2009; Noterdaeme et al., 2009,
2012; Zafar et al., 2013; Crighton et al., 2015; Neeleman et al., 2016; Bird
et al., 2017).

The measurements of molecular hydrogen (H2) are commonly probed with surveys
of CO rotational lines over a large redshift range (see e.g., Tacconi et al.,
2020, for a review), such as FCRAO Extragalactic CO Survey (Young et al.,
1995), HERACLES ($z\sim 0$, Leroy et al., 2009), xCOLD GASS ($z\sim 0$,
Saintonge et al., 2017), ASPECS ($1<z<3.5$, Decarli et al., 2016, 2019), COLDz
(Riechers et al., 2019) PHIBSS2 ($0.48<z<5.25$, Lenkić et al., 2020), and
NOEMA HDF-N Survey (Boogaard et al., 2023). The molecular gas mass can also be
estimated with the far-IR dust continuum using reasonable molecular gas-to-
dust mass ratios (see e.g., Carilli & Walter, 2013; Berta et al., 2016;
Scoville et al., 2017; Magnelli et al., 2020).

Together with galaxy stellar properties, e.g., stellar mass ($M_{\ast}$) and
star formation rate (SFR), obtained from optical surveys such as SDSS (York et
al., 2000) and CANDELS (Grogin et al., 2011; Koekemoer et al., 2011), well-
defined scaling relations between gas mass and other galaxy properties have
been established with statistical and representative samples (see summaries in
e.g., Tacconi et al., 2018; Saintonge & Catinella, 2022). Although the H iand
H2 masses typically increase with $M_{\ast}$, the ratios of $f_{\textrm{H
{i}}}\equiv M_{\textrm{H {i}}}/M_{\ast}$ and $f_{\rm H_{2}}\equiv M_{\rm
H_{2}}/M_{\ast}$ are anticorrelated with $M_{\ast}$, i.e., galaxies become
progressively more gas-poor as their stellar mass grows. There are also
positive tight correlations between gas fractions and the specific SFR (${\rm
sSFR}\equiv{\rm SFR/M_{\ast}}$), with the $f_{\rm H_{2}}$–sSFR relation at
$z\sim 0$ showing a small intrinsic scatter of 0.12 dex (Saintonge et al.,
2017). In terms of molecular gas depletion time, $\tau_{\rm H_{2}}\equiv
M_{\rm H_{2}}/{\rm SFR}$, it is typically around $\sim 1$ Gyr in the local
universe and decreases slightly towards higher redshifts (Tacconi et al.,
2020).

Combining atomic and molecular gas with stars, we can quantify the total
baryon content associated with galaxies (e.g., Walter et al., 2020). The exact
contributions from different components provide valuable constraints to galaxy
formation and evolution models. However, it is still a long-standing question
how the different baryonic components of galaxies evolve with the mass
assembly histories of dark matter halos. Various theoretical models have been
proposed to understand the cosmic evolution of the cold gas reservoir, using
hydrodynamical simulations (e.g., Lagos et al., 2015; Crain et al., 2017;
Diemer et al., 2019; Davé et al., 2020) and semi-analytical models (e.g.,
Lagos et al., 2011; Fu et al., 2013; Popping et al., 2014; Kim et al., 2017;
Xie et al., 2017; Baugh et al., 2019; Chauhan et al., 2020; Spinelli et al.,
2020). By modeling the star formation and feedback mechanisms in simulations,
these models show reasonable agreement in reproducing some of the observed
cold gas properties, such as the cold gas mass functions. However, different
galaxy formation models have varying predictions of gas accretion and
consumption, causing large discrepancies when comparing other gas scaling
relations, e.g., H i-halo mass relation (Guo et al., 2020; Li et al., 2022a,
b; Dev et al., 2023; Rhee et al., 2023) and H i-stellar mass relation
(Catinella et al., 2018; Janowiecki et al., 2020; Guo et al., 2021).

Popping et al. (2015) proposed a semi-empirical approach based on the subhalo
abundance matching model of Behroozi et al. (2013a) that provides realistic
star formation histories for each subhalo in the simulations by fitting to the
observed stellar mass functions, cosmic SFRs and sSFRs. They estimated the
cold gas masses with the gas density–SFR relation of Bigiel et al. (2008) and
pressure-regulated H i-to-H2 transition model of Blitz & Rosolowsky (2006) by
assuming exponential distributions of gas disks. After calibrating with the
observed H i and H2 masses, it shows reasonable agreement with the local H i
and H2 mass functions for gas-rich galaxies. Although it still has difficulty
reproducing the distributions of gas-poor galaxies and cosmic evolution of
cold gas densities, such an approach shows the advantage of empirical models
in characterizing the realistic evolution of gas content in halos, without
introducing complex galaxy formation recipes.

Other simpler empirical models have been suggested to match some of the
observed gas properties. For example, conditional H i mass distributions based
on galaxy luminosity and color (Paul et al., 2018; Dutta & Khandai, 2021),
stellar mass and morphology (Calette et al., 2021a, b), or even a more complex
multi-parameter H imass estimator (Li et al., 2022a), show better agreement
with fitted observations and provide reasonable predictions for other gas
properties. Currently, all of these models are limited to $z\sim 0$, where a
wealth of direct H i observations are available to constrain the models. In
principle, it is possible to construct empirical models to coherently trace
the evolution of cold gas within the dark matter halos.

In this paper, we propose a new empirical model framework,
NeutralUniverseMachine, that can self-consistently predict the evolution of
both H i and H2 in a broad redshift range of $0<z<6$. Our work is based on the
empirical model of UniverseMachine (Behroozi et al., 2019) that successfully
matches many of the observed galaxy stellar properties, such as stellar mass
functions, SFRs, quenched fractions, and correlation functions over $0<z<10$.
The key parameters that we adopt from the UniverseMachine model are $M_{\ast}$
and SFR, which are essential in modeling the cold gas content. Our H i and H2
models based on UniverseMachine would then provide a comprehensive description
of the baryon content of galaxies.

The organization of this paper is as follows. In Section 2, we introduce our
modeling method, as well as the observational data used in the constraints. We
show the results in Section 3. The discussion and conclusions are presented in
Section 4 and Section 5, respectively. Throughout the paper, all masses are
expressed in units of $M_{\sun}$. We assume a flat $\Lambda$CDM cosmology of
$\Omega_{\rm m}=0.307$, $h=0.678$, $\Omega_{\rm b}=0.048$ and
$\sigma_{8}=0.823$, consistent with the Planck15 results (Planck Collaboration
et al., 2016).

##  2 Data and Method

###  2.1 UNIVERSEMACHINE Catalog

The UniverseMachine model parametrizes the probability distribution of galaxy
SFRs in halos as a function of $v_{\rm M_{peak}}$, $z$ and $\Delta v_{\rm
max}$, where $v_{\rm M_{peak}}$ is the maximum halo circular velocity ($v_{\rm
max}$) at the redshift of the peak halo mass ($M_{\rm peak}$) and $\Delta
v_{\rm max}$ is the relative change in $v_{\rm max}$ over the past halo
dynamical time. $v_{\rm M_{peak}}$ and $\Delta v_{\rm max}$ are used as
proxies for the halo mass and mass accretion rate, respectively. An SFR will
be assigned to each halo (or subhalo) according to the probability
distribution. Meanwhile, halos of larger $\Delta v_{\rm max}$ are assigned
higher SFRs, allowing for random scatters. The galaxy stellar masses can be
obtained from integrating SFRs along the halo merger trees. In this way, the
galaxy mass growth history is self-consistently calculated from the halo
assembly and star formation histories. The 44 model-free parameters are
determined from fitting to the observations of stellar mass functions, cosmic
SFRs and sSFRs, quenched fractions, UV luminosity functions, UV–stellar mass
relations and correlation functions. We refer the readers to Behroozi et al.
(2019) for more details.

In this paper, we use the public catalogs of
http://halos.as.arizona.edu/UniverseMachine/DR1/ (catalog UniverseMachine
DR1), which are run on the Bolshoi-Planck $N$-body simulation (Klypin et al.,
2016) with a box size of $250\,h^{-1}{\rm{Mpc}}$ on a side and a dark matter
particle mass resolution of $2.3\times 10^{8}M_{\sun}$. The cosmological
parameters are the same as our adopted cosmology. Dark-matter halos (and
subhalos) are identified using the Rockstar halo finder (Behroozi et al.,
2013b), while halo merger trees are constructed through the Consistent Trees
algorithm (Behroozi et al., 2013c).

We will use the following quantities in the UniverseMachine catalogs, $M_{\rm
vir}$ (halo or subhalo virial mass, defined as in Bryan & Norman 1998),
$M_{\rm peak}$ (the peak halo or subhalo virial mass over its entire merger
history), $v_{\rm M_{peak}}$ (as defined above), $z_{\rm form}$ (halo
formation time, defined as the redshift when the most massive progenitor
reaches $0.5M_{\rm peak}$), $R_{\rm vir}$ (halo virial radius, also defined in
Bryan & Norman 1998), $M_{\ast}$, and SFR. To match the observational
measurements, we adopted the galaxy properties of $M_{\ast}$ and SFR with
random and systematic uncertainties included (see the discussions in Section
3.5 of Behroozi et al. 2019).

It is suggested in Behroozi et al. (2019) that the inclusion of orphan
galaxies (i.e., tidally stripped subhalos with masses below the simulation
resolution) will improve the matching with observed clustering measurements
and produce the correct satellite evolution. We also include the orphan
galaxies from the UniverseMachine catalogs in our cold-gas model and trace
their merger trees to derive the corresponding measurements of $z_{\rm form}$.
Our modeling results would also be less affected by the simulation resolution
of Bolshoi-Planck after the inclusion of orphans.

###  2.2 NeutralUniverseMachine Model

As noted in Saintonge & Catinella (2022), the common practice in the H i
community is to calculate the H i mass using observed 21 cm fluxes, without
accounting for the contribution of helium and heavier elements. However, the
H2 mass measurements typically included an upward correction of 1.36 to obtain
all molecular gas masses. To avoid confusion and be consistent with the
literature, we still use $M_{\textrm{H {i}}}$ to represent H i mass only and
use $M_{\rm H_{2}}$ to indicate the total mass of the molecular gas in the
following sections.

####  2.2.1 H I Model

The H i mass of each halo or subhalo is determined as a function of its own
$M_{\rm vir}$, $z_{\rm form}$, SFR, and $z$, as follows.

$\displaystyle M_{\textrm{H {i}}}$ $\displaystyle=$ $\displaystyle\frac{\kappa
M_{\rm vir}}{\mu^{-\alpha}+\mu^{\beta}}\left(\frac{1+z}{1+z_{\rm
form}}\right)^{\gamma}\left(\frac{\rm SFR}{\rm SFR_{MS,obs}}\right)^{\lambda}$
(1) $\displaystyle\mu$ $\displaystyle=$ $\displaystyle M_{\rm vir}/M_{\rm
crit}$ (2) $\displaystyle\log\kappa$ $\displaystyle=$
$\displaystyle\kappa_{0}+\kappa_{1}z+\kappa_{2}z^{2}$ (3) $\displaystyle\log
M_{\rm crit}$ $\displaystyle=$ $\displaystyle M_{0}+M_{1}z+M_{2}z^{2},$ (4)

where SFR is in units of $M_{\sun}/{\rm yr}$ and ${\rm SFR_{MS,obs}}$ is the
best-fitting median (observed) SFR of star-forming galaxies defined in
Behroozi et al. (2019),

$\displaystyle{\rm SFR_{MS,obs}}$ $\displaystyle=$ $\displaystyle{\rm
SFR_{MS}}+0.041-0.044z/(1+z)+$ (5) $\displaystyle
0.314\exp\left(-\frac{(z-2)^{2}}{2}\right)$ $\displaystyle{\rm SFR_{MS}}$
$\displaystyle=$
$\displaystyle\epsilon\left[(v^{\alpha^{\prime}}+v^{\beta^{\prime}})^{-1}+\gamma^{\prime}\exp\left(-\frac{(\log
v)^{2}}{2\delta^{2}}\right)\right]$ (6) $\displaystyle v$ $\displaystyle=$
$\displaystyle v_{\rm M_{peak}}/V$ (7) $\displaystyle\log V$ $\displaystyle=$
$\displaystyle 2.151-1.658z/(1+z)$ (8) $\displaystyle+1.680\ln(1+z)-0.233z$
$\displaystyle\log\epsilon$ $\displaystyle=$ $\displaystyle
0.109-3.441z/(1+z)+$ (9) $\displaystyle 5.079\ln(1+z)-0.781z$
$\displaystyle\alpha^{\prime}$ $\displaystyle=$
$\displaystyle-5.598-20.731z/(1+z)+$ (10) $\displaystyle
13.455\ln(1+z)-1.321z$ $\displaystyle\beta^{\prime}$ $\displaystyle=$
$\displaystyle-1.911+0.395z/(1+z)+0.747z$ (11)
$\displaystyle\log\gamma^{\prime}$ $\displaystyle=$
$\displaystyle-1.699+4.206z/(1+z)-0.809z$ (12) $\displaystyle\delta$
$\displaystyle=$ $\displaystyle 0.055$ (13)

Equation (5) corrects for the systematic offset between observed and intrinsic
SFRs as in Behroozi et al. (2019). The functional formal of Equation (1) has
three parts, the double power law relation of $M_{\rm vir}$, the scaling with
the halo formation time $z_{\rm form}$ and the dependence on SFR. In summary,
we have 10 free model parameters for the H i gas in halos, $\kappa_{0}$,
$\kappa_{1}$, $\kappa_{2}$, $M_{0}$, $M_{1}$, $M_{2}$, $\alpha$, $\beta$,
$\gamma$, and $\lambda$.

Figure 1: Distributions of galaxies as a function of $M_{\ast}$ and SFR from
$z=0$ to $z=5$ in the UniverseMachine catalog. The color is coded by the
logarithmic number densities of galaxies, with redder colors for higher
densities. The solid black line in each panel indicates the positions of ${\rm
SFR_{MS,obs}}$.

In Figure 1, we show the distributions of galaxies as a function of $M_{\ast}$
and SFR from $z=0$ to $z=5$ in the UniverseMachine catalog. The black solid
line in each panel indicates the positions of ${\rm SFR_{MS,obs}}$. The
feature of bimodal distributions in SFR becomes more prominent at lower
redshifts. The SFR distribution in UniverseMachine is simply parameterized as
the sum of two lognormal distributions that mimic the bimodal distribution.

The double power law of $M_{\rm vir}$ is motivated by the observed H i-halo
mass relation of Guo et al. (2020), where H i mass increases in general with
the halo mass and the slope becomes steeper at the lower mass end, which is
also confirmed by a recent study of Dev et al. (2023). Physically, a galaxy’s
H i disk is typically much more extended than its stellar disk (Wang et al.,
2016) and thus more affected by the halo environment for low-mass galaxies
with low binding energy, as found in both central (Guo et al., 2021) and
satellite galaxies (Brown et al., 2017; Stevens et al., 2019). The redshift
evolution of normalization parameter $\kappa$ is to reproduce the change of
cosmic H i density at higher redshifts (Rhee et al., 2018; Walter et al.,
2020). We assume a similar functional form for the evolution of the
characteristic mass $M_{\rm crit}$. The parameters $\alpha$ and $\beta$ are
the slopes of the low and high mass halos, respectively. We assume no
evolution in these slopes, consistent with the trend found in hydrodynamical
simulations and semi-analytical models (e.g., Villaescusa-Navarro et al.,
2018; Spinelli et al., 2020; Chauhan et al., 2020).

The dependence of $M_{\textrm{H {i}}}$ on both the halo mass and the halo
formation time is found by Guo et al. (2017) using the spatial two-point
correlation function measurements of the H i-selected galaxies in ALFALFA.
They found that the clustering of H i-rich galaxies is significantly lower
than the predictions from simple abundance matching models, strongly favoring
the additional dependence on the halo formation time (see also Li et al.,
2022b). It is also confirmed by Stiskalek et al. (2021) using a different
abundance matching method. Late-forming halos tend to host more H i gas, which
is probably caused by the more recent accretion of cold gas or the new gas
reservoir brought in by recent wet mergers. We include the dependence of
formation time as $M_{\textrm{H {i}}}\propto[(1+z)/(1+z_{\rm
form})]^{\gamma}$, based on the results of Guo et al. (2017) and Li et al.
(2022b). In this way, the dependence becomes weaker at higher redshifts when
the halos are in the rapid growth phase with smaller differences between $z$
and $z_{\rm form}$ (Zhao et al., 2009).

As shown in Guo et al. (2021), the logarithmic growth/reduction of a galaxy’s
H i mass is proportional to its logarithmic offset from the star formation
main sequence (SFMS), $\log({\rm SFR/SFR_{MS}})$. For galaxies in the SFMS, we
can define the H i main sequence as,

$M_{\textrm{H {i}},{\rm MS}}=\frac{\kappa M_{\rm
vir}}{\mu^{-\alpha}+\mu^{\beta}}\left(\frac{1+z}{1+z_{\rm
form}}\right)^{\gamma}.$ (14)

Then the scaling of H i mass with SFR is simply $\log(M_{\textrm{H
{i}}}/M_{\textrm{H {i}},{\rm MS}})=\lambda\log({\rm SFR/SFR_{MS,obs}})$, where
$\lambda$ is found to be $\sim 0.4$ at $z\sim 0$ (Guo et al., 2021).

####  2.2.2 H2 Model

Tacconi et al. (2020) proposed a fitting formula for $M_{\rm H_{2}}$, as a
function of $M_{\ast}$, ${\rm SFR/SFR_{MS}}$ and $\log(1+z)$, which shows good
agreement with a compilation of $M_{\rm H_{2}}$ measurements in the
literature. We follow their functional form for the molecular gas as,

$\displaystyle M_{\rm H_{2}}$ $\displaystyle=$ $\displaystyle\zeta
M_{\ast}^{\nu}\left(\frac{\rm SFR}{\rm SFR_{MS,obs}}\right)^{\eta}$ (15)
$\displaystyle\log\zeta$ $\displaystyle=$
$\displaystyle\zeta_{0}+\zeta_{1}\ln(1+z)+\zeta_{2}\left[\ln(1+z)\right]^{2},$
(16)

The five free parameters are $\zeta_{0}$, $\zeta_{1}$, $\zeta_{2}$, $\nu$, and
$\eta$. As in the case of H i, for galaxies on the SFMS, their H2 main
sequence is described by $M_{\rm H_{2},MS}=\zeta M_{\ast}^{\nu}$. The star
formation law is then $\log(M_{\rm H_{2}}/M_{\rm H_{2},MS})=\eta\log({\rm
SFR/SFR_{MS,obs}})$. Although the effect of the environment on the molecular
gas content is still debated (e.g., Noble et al., 2017; Darvish et al., 2018;
Tadaki et al., 2019), it is generally found to be small. Interestingly, since
the molecular gas is distributed mostly in the stellar disks, it would be less
affected by the halo environment than in the case of H i. Thus, our H2 main
sequence is defined on the basis of stellar mass rather than halo mass.

###  2.3 Observational Measurements

In order to fully constrain the parameters of the cold gas model, we have
included recent observations of the H i mass function, H i-halo mass relation,
H i-stellar mass relations, H2 mass function, H2-stellar mass relation, and
molecular-to-atomic mass ratio for $z\sim 0$, as summarized in Table 1. For
higher redshift observations, we used the collections of cosmic H i density
and molecular gas density measurements in Walter et al. (2020), as well as the
H i-stellar mass relation of Chowdhury et al. (2022a) at $z\sim 1.1$. The
details are described as follows.

Table 1: Observational Constraints Measurements | References | Redshifts
---|---|---
H i Mass Function | This work (Appendix A) | $z\sim 0$
H i-Halo Mass Relation | Guo et al. (2020) | $z\sim 0$
H i-Stellar Mass Relation | Guo et al. (2021) | $z\sim 0$
H2 Mass Function | Fletcher et al. (2021) | $z\sim 0$
H2-Stellar Mass Relation | Saintonge et al. (2017) | $z\sim 0$
H2-to-H i mass ratio | Catinella et al. (2018) | $z\sim 0$
H i-Stellar Mass Relation | Chowdhury et al. (2022a) | $z\sim 1.1$
cosmic H i density | Walter et al. (2020) | $0<z<5$
cosmic H2 density | Walter et al. (2020) | $0<z<6$

####  2.3.1 Local Universe

The H i mass function (HIMF), $\phi(M_{\textrm{H {i}}})$, in the local
universe has been measured from the final catalog of ALFALFA (Haynes et al.,
2018) covering $\sim 6900\leavevmode\nobreak\ {\rm deg}^{2}$ over the redshift
range of $0<z<0.05$ in Jones et al. (2018). It describes the average number
densities of galaxies in given H i mass bins. By including an effective volume
weight of each galaxy in the catalog, the sample selections of ALFALFA and the
influence of local large-scale structures can be robustly taken into account
in the final measurements of $\phi(M_{\textrm{H {i}}})$ (Martin et al., 2010).
The H i observations of ALFALFA depend on both the H i flux and the line
profile width. After selecting H i targets above the 50% completeness limit
Jones et al. (2018) measured the HIMF using the two-dimensional stepwise
maximum likelihood (2DSWML) method (see e.g., Zwaan et al., 2005; Martin et
al., 2010). They have approximated the completeness limit as a sharp cut-off
at the location of 50% completeness (Eq. 5 of Haynes et al., 2011). However,
given the large number of sources now available in the ALFALFA catalog, the
missing sources between 50% and 100% completeness could potentially lead to an
overestimation of the effective volumes for low-$M_{\textrm{H {i}}}$ sources
and thus underestimate low-mass bins of the HIMF. We have therefore decided to
remeasure the HIMF following the same steps of Jones et al. (2018), but only
consider sources above the 90% completeness limit (Eq. 4 of Haynes et al.,
2011). It avoids the need to heavily weight the H i sources with low
completeness. We verify that changing from 90% completeness to even higher
thresholds does not have any visible effect on HIMF. We estimate the errors of
$\phi(M_{\textrm{H {i}}})$ using the fractional errors from 64 mock
observations of the ALFALFA final sample, as will be detailed in Section 3.2.
We list the corrected HIMF measurements in Appendix A.

The H i-halo mass relation, $\langle M_{\textrm{H {i}},{\rm tot}}|M_{\rm
vir}\rangle$, describes the total H i mass (i.e. including the H i gas of all
member galaxies) in halos of different masses. By selecting halos in the
overlapping regions between ALFALFA and the galaxy groups constructed from
SDSS DR7 (Lim et al., 2017), Guo et al. (2020) stacked the total H i spectra
for halos in different $M_{\rm vir}$ bins, using the aperture size determined
from the halo virial radius with a minimum diameter of $8\arcmin$, which is
about twice the ALFALFA beam width. The average H i mass is then obtained by
dividing the total H i mass in each $M_{\rm vir}$ bin by the corresponding
number of galaxy groups. The uncertainties in $\langle M_{\textrm{H {i}},{\rm
tot}}\rangle$ were determined from bootstrapping the stacking samples.

The H i-stellar mass relation, $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$,
provides additional constraints to the correlation between the stellar and gas
components of galaxies. By applying the H i spectra stacking technique, Guo et
al. (2021) accurately measured $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$
for star-forming and quenched _central_ galaxies (i.e., divided by cuts in
SFR) using the same galaxy sample as in Guo et al. (2020). This is especially
important for quenched galaxies, most of which are usually without individual
21 cm detections (Catinella et al., 2018). The use of $\langle M_{\textrm{H
{i}}}(M_{\ast})\rangle$ for star-forming and quenched central galaxies will
also constrain the power-law index $\eta$ of SFR in Equation (1).

The H2 mass function (H2MF), $\phi(M_{\rm H_{2}})$, at $z\sim 0$ is robustly
measured by Fletcher et al. (2021) using a large set of H2 measurements (532
galaxies with a detection rate of 63%) in xCOLD GASS (Saintonge et al., 2017).
We use their $\phi(M_{\rm H_{2}})$ values for $8<\log(M_{\rm
H_{2}}/M_{\sun})<10.16$, where the measurements are more complete. When
integrating over the best-fitting Schechter function, they found a cosmic H2
density of $\rho_{\rm H_{2}}=1.04\times 10^{7}M_{\sun}/{\rm Mpc^{3}}$.

We also utilize the H2-stellar mass relation, in terms of the molecular gas
fraction $f_{\rm H_{2}}(M_{\ast})$ measured by Saintonge et al. (2017) using
xCOLD GASS. We adopt their stacking measurements of $\log\langle f_{\rm
H_{2}}(M_{\ast})\rangle$, which properly include the contribution from
galaxies without detection. We adopt their measurements for the main-sequence
galaxies, defined as $|\log({\rm SFR}/{\rm SFR_{MS}})|<0.4$. For fair
comparisons, we would apply the same SFMS definition of Saintonge et al.
(2016) (their Eq. 5) in our model when fitting to $f_{\rm H_{2}}(M_{\ast})$ of
the main sequence galaxies.

To connect the H2 model with that of H i, it is important to have constraints
from the molecular-to-atomic mass ratio (i.e. $R_{\rm mol}\equiv M_{\rm
H_{2}}/M_{\textrm{H {i}}}$). We adopt the weighted average relation between
$\langle\log R_{\rm mol}\rangle$ and $M_{\ast}$ from Catinella et al. (2018)
(their Table 3). We note that the H2 measurements in Catinella et al. (2018)
were obtained from xCOLD GASS, but without the correction for helium and
heavier elements. Therefore, we weighted their $R_{\rm mol}$ measurements by a
factor of 1.36 to be consistent with our definition.

The above observational measurements at $z\sim 0$ come from the ALFALFA,
xGASS, and xCOLD GASS surveys. Both ALFALFA and xGASS data were observed with
the Arecibo telescope (Catinella et al., 2010). The use of these self-
consistent measurements will then provide tight constraints on our model. We
also note that the H i measurements adopted here do not include the self-
absorption correction, which is quite uncertain and still under debate (Jones
et al., 2018).

####  2.3.2 High Redshifts

There are fewer robust gas scaling relations available at higher redshifts,
due to the lack of large statistical samples, especially for the H i
observations. To anchor the correct H i mass distribution, we adopt the latest
H i stacking measurements of Chowdhury et al. (2022a) using the Giant
Metrewave Radio Telescope (GMRT) for star-forming galaxies in two different
stellar mass bins at $0.74<z<1.45$ ($\langle z\rangle\sim 1.1$).

It is much easier to estimate the cosmic H i and H2 densities than to obtain
robust scaling relations. Since our model only includes 6 parameters
($\kappa_{1}$, $\kappa_{2}$, $M_{1}$, $M_{2}$, $\zeta_{1}$, $\zeta_{2}$) to
describe the redshift dependence, they can be well constrained with the cosmic
gas density measurements at various redshifts. We adopt the collections of
various cosmic H i density ($\rho_{\textrm{H {i}}}$) and H2 density
($\rho_{\rm H_{2}}$) measurements in the literature from Walter et al. (2020)
in the redshift range of $0<z<6$. However, we note that their H i densities
(their Table 2) have included the helium contribution. To match our
definition, we have corrected for the different cosmologies and removed the
helium contribution using the original measurements in the literature.

As summarized in Walter et al. (2020), the cosmic H i densities can be derived
by measuring the 21 cm emission in the local universe, stacking H i spectra
and 21 cm intensity mapping at intermediate redshifts, and using DLA systems
at $z>1.5$. The uncertainties of these methods become increasingly larger at
higher redshifts, especially for DLAs. DLA measurements are generally limited
to systems above an H i column density threshold of $\log N({\textrm{H
{i}}})>20.3\,{\rm cm}^{-2}$ (e.g., Péroux et al., 2003; O’Meara et al., 2007).
The ignorance of lower $N({\textrm{H {i}}})$ systems will cause an
underestimate of H i density by 10–20% (Zafar et al., 2013; Berg et al., 2019;
Péroux & Howk, 2020).

On the other hand, the DLAs probe H i gas in and around galaxies. In fact,
there is an increased contribution of H i gas outside the galaxies (i.e., the
intergalactic medium) toward high redshifts (Péroux & Howk, 2020). Using the
hydrodynamical simulation of IllustrisTNG (Nelson et al., 2019), Villaescusa-
Navarro et al. (2018) found that about 80% and 90% of the contributions to
cosmic H i densities at $z=5$ are still from galaxies and within halos,
respectively. It means that our model will also likely miss 10–20% of H i
signals at high redshifts, counterbalancing the underestimation of H i
densities from DLAs. Thus, DLA measurements at high redshifts will still
provide reasonable constraints to the H i gas within halos.

The uncertainties associated with the cosmic molecular gas density inferred
from CO line luminosity lie mainly in the systematics of the CO to H2
conversion factor (commonly denoted as $\alpha_{\rm CO}$) (see discussions in
Tacconi et al., 2020). $M_{\rm H_{2}}$ measurements in xCOLD GASS were derived
with a conversion factor that depends on both the metallicity and the offset
from the SFMS (Accurso et al., 2017). The median value of $\alpha_{\rm CO}$ is
around $3.3M_{\sun}{\rm(K\,km\,s^{-1}\,pc^{2})}^{-1}$ in the xCOLD GASS
sample. But other CO measurements, as listed in Table 3 of Walter et al.
(2020) adopted a constant $\alpha_{\rm
CO}=3.6M_{\sun}{\rm(K\,km\,s^{-1}\,pc^{2})}^{-1}$. It will cause a small level
of inconsistency among the CO measurements, but the large uncertainties for
$\rho_{\rm H_{2}}$ at high redshifts make it less an issue. Using the dust
continuum to infer $\rho_{\rm H_{2}}$ also requires the assumptions of a
metallicity-dependent gas-to-dust mass ratio, dust temperature, and
emissivity, causing the systematic uncertainties of $\rho_{\rm H_{2}}$ in the
literature.

There are still some H i and H2 gas scaling relation measurements available at
high redshift in the literature. However, their sample selections would vary
from one to another. We will use these measurements as consistency check of
our model predictions in the following sections, rather than apply the
different selection cuts to constrain our model.

###  2.4 Model Fitting

To fit the observational H i and H2 measurements using our 15-parameter
NeutralUniverseMachine model, we apply the Bayesian inference tool of
MultiNest (Feroz et al., 2009) to explore the parameter space. For each set of
model parameters in the Monte Carlo Markov Chain (MCMC), we can generate the
corresponding H i and H2 masses for each halo (central galaxy) and subhalo
(satellite galaxy) following Equations (1) and (15). To fit the observational
data of $\phi(M_{\textrm{H {i}}})$, $\langle M_{\textrm{H {i}},{\rm
tot}}|M_{\rm vir}\rangle$, $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$,
$\phi(M_{\rm H_{2}})$, $f_{\rm H_{2}}(M_{\ast})$ and $R_{\rm mol}$, we adopt
the same binning schemes as in the references (Table 1) to avoid the binning
effect. To be consistent with the observations using H i spectra stacking, we
calculate the average H i mass in $\langle M_{\textrm{H {i}},{\rm tot}}|M_{\rm
vir}\rangle$ and $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$ as
$\Sigma_{i}M_{{\textrm{H {i}}},i}/N$, where $N$ is the total number of halos
(or galaxies) in each bin. In order to fit the observed H2MF, we have added a
Gaussian random scatter of $0.2$ dex to $M_{\rm H_{2}}$ for each galaxy as
suggested in Saintonge & Catinella (2022), mimicking measurement errors.

To match the observations of $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$ for
star-forming and quenched galaxies in Guo et al. (2021), we adopt the same SFR
cut of $\log({\rm SFR_{cut}}/{\rm
yr}^{-1}M_{\sun})=0.65\log(M_{\ast}/M_{\sun})-7.25$ to separate the two
populations. To fit the $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$ relation
at $z\sim 1.1$ for star-forming galaxies from Chowdhury et al. (2022a), we
define our star-forming population using $\log({\rm SFR/SFR_{MS,obs}})>-1$. As
shown in Figure 1, it provides a reasonable cut to separate star-forming and
quenched galaxies.

Other observables are straightforward to measure in the UniverseMachine
catalogs. There are 149 redshift outputs between $z=0$ and $z=6$ in the
http://halos.as.arizona.edu/UniverseMachine/DR1/ (catalog UniverseMachine DR1)
of the Bolshoi-Planck simulation. We calculate $\rho_{\textrm{H {i}}}$ and
$\rho_{\rm H_{2}}$ for all redshift outputs in each run of MCMC, and
interpolate within the outputs to derive gas densities at the observed
redshifts.

The $z\sim 0$ measurements are self-consistent with each other and the
uncertainties are also small. But the high-redshift gas density measurements
suffer from large errors and systematic uncertainties, as well as
inconsistency between different sets of observations. It is more practical and
reliable to first constrain the nine redshift-independent model parameters
using $z\sim 0$ measurements. The remaining 6 parameters ($\kappa_{1}$,
$\kappa_{2}$, $M_{1}$, $M_{2}$, $\zeta_{1}$, $\zeta_{2}$) can then be better
constrained with the cosmic gas densities at higher redshifts. Taking into
account the large uncertainties in the H i and H2 densities, we emphasize that
our best-fitting models at these high redshifts are based on the available
data sets and may be improved with future observations.

The likelihood of each MCMC run is proportional to $\exp(-\chi^{2}/2)$, where
the total $\chi^{2}$ is determined as the sum of $\chi^{2}$ for each set of
measurements in Table 1. For different sets of measurements, the covariance
between the data points is hard to quantify, and we simply use the error of
each data point to calculate $\chi^{2}$. Minimizing the total $\chi^{2}$ seems
to bias the fit toward data sets with more data points (e.g., the HIMF).
Intuitively, this corresponds to weighting the data sets by information
content, which is the optimal way to use the available data. In cases where
one data set is not fit well, this could be due to some combination of: (1) a
fitting function that is less flexible than needed, (2) unmodeled systematics,
instrumental effects or analysis assumptions that cause tension between data
sets, and (3) statistical fluctuations. In the first two cases, the
statistically valid solution is to add parameters (to the fitting function and
to the nuisance parameter set, respectively) to better account for the data,
as in Behroozi et al. (2019). In the latter case, the statistically valid
solution is to weight the points by the total $\chi^{2}$.

##  3 Results

###  3.1 Best-fitting Model

Our best-fitting model is determined from the MCMC run with the maximum
likelihood, which is capable of reproducing all the observations in Table 1.
The observational data, as well as the best-fitting models, are shown for
$\phi(M_{\textrm{H {i}}})$ (Fig. 2, left panel), $\langle M_{\textrm{H
{i}},{\rm tot}}|M_{\rm vir}\rangle$ (Fig. 2, middle panel), $\langle
M_{\textrm{H {i}}}(M_{\ast})\rangle$ (Fig. 2, right panel), $\phi(M_{\rm
H_{2}})$ (Fig. 3, left panel), $f_{\rm H_{2}}(M_{\ast})$ (Fig. 3, middle
panel), $R_{\rm mol}$ (Fig. 3, right panel), $\rho_{\textrm{H {i}}}(z)$ (Fig.
4, left panel), $\rho_{\rm H_{2}}(z)$ (Fig. 4, right panel), and $\langle
M_{\textrm{H {i}}}(M_{\ast})\rangle$ at $z\sim 1.1$ (Fig. 9, right panel). The
observational measurements are shown as the points with errors, whereas our
model predictions are displayed as solid lines.

The best-fitting model parameters are,

$\displaystyle\log\kappa=-0.972^{+0.012}_{-0.062}-0.180^{+0.046}_{-0.126}z+0.053^{+0.067}_{-0.026}z^{2}$
(17) $\displaystyle\log M_{\rm
crit}=10.832^{+0.045}_{-0.016}+0.835^{+0.068}_{-0.253}z-0.246^{+0.134}_{-0.044}z^{2}$
(18) $\displaystyle\alpha=1.346^{+0.162}_{-0.071},\ \
\beta=0.604^{+0.011}_{-0.017}$ (19)
$\displaystyle\gamma=2.233^{+0.037}_{-0.137},\ \
\lambda=0.433^{+0.010}_{-0.007}$ (20)
$\displaystyle\nu=0.921^{+0.009}_{-0.016},\quad\eta=0.896^{+0.036}_{-0.021}$
(21)
$\displaystyle\log\zeta=-0.384^{+0.166}_{-0.094}+1.420^{+0.084}_{-0.071}\ln(1+z)$
$\displaystyle\qquad\quad\ -0.425^{+0.059}_{-0.069}\left[\ln(1+z)\right]^{2},$
(22)

where $M_{\rm crit}$ is in units of $M_{\sun}$. In Appendix B, we display the
density distributions of the model parameters in the MCMC chains.

Our best-fitting model shows good agreement with the HIMF measurements
$\phi(M_{\textrm{H {i}}})$ in the left panel of Figure 2. The HIMF is
typically fitted with a Schechter function in the literature (Zwaan et al.,
2005; Martin et al., 2010; Jones et al., 2018). But it is naturally explained
in our halo model by integrating the halo mass function weighted by the H i
mass. We note that our halo-based H i model would be affected by the dark
matter particle mass resolution of the Bolshoi-Planck simulation ($2.3\times
10^{8}M_{\sun}$). The halo mass function is only accurate for halos with
$M_{\rm vir}>10^{10}M_{\sun}$ (Klypin et al., 2016), which corresponds to an
average H i mass limit of $10^{6.9}M_{\sun}$ for our best-fitting model. We
indicate this mass limit as the vertical dotted line in the left panel of
Figure 2. Galaxies with $M_{\textrm{H {i}}}$ below this limit should be used
with caution.

Our predicted H i-halo mass relation $\langle M_{\textrm{H {i}},{\rm
tot}}|M_{\rm vir}\rangle$ is shown as the solid line in the middle panel. It
seems to slightly overestimate the total H i mass in halos of $M_{\rm vir}\sim
10^{11.25}M_{\sun}$. As will be shown in Section 3.3.1, it is likely caused by
both the sample selection of SDSS galaxies and the halo mass uncertainties in
the SDSS group catalog. As shown in Fig. 10 of Lim et al. (2017), the
estimates of the halo mass from the group finder are biased high for the halo
of $M_{\rm vir}<10^{11}M_{\sun}$. We note that the minimal halo mass in the
group catalog is $10^{11.14}M_{\sun}$, which is potentially overestimated for
those low-mass halos.

In the right panel of Figure 2, the H i-stellar mass relations are shown for
both star-forming (blue symbols) and quenched _central_ galaxies (red
symbols). As expected, the best-fitting power index $\lambda$ is consistent
with the observation ($\sim 0.4$). The fitting to quenched galaxies is
somewhat worse at the low-mass end, mainly due to the slight differences
between the SFR distributions of SDSS and UniverseMachine model. For
comparison, we also display the $\langle M_{\textrm{H {i}}}(M_{\ast})\rangle$
measurements for all central galaxies (gray symbols), which are not included
in observational constraints but are still well reproduced by our model
(dotted line), because it is simply the summation of $M_{\textrm{H {i}}}$ in
the star-forming and quenched populations weighted by the corresponding
numbers of galaxies.

Figure 2: Comparisons between the observed measurements and the best-fitting
model predictions, for HIMF (top left panel), H i-halo mass relation (top
right panel), H i-stellar mass relations of central galaxies (bottom left
panel) and projected 2PCF (bottom right panel) at $z\sim 0$. The observational
data are shown as points with errors, while the best-fitting models are
represented by solid lines. The estimated H i mass limit of our halo-based
model is shown as the vertical dotted line in the left panel (see text for
details). The measurements of H i-halo mass relation with confusion correction
in Guo et al. (2020) are shown as the open blue circles in the top right
panel. The H i-stellar mass relations of central galaxies are shown for both
star-forming (blue points and blue line) and quenched galaxies (red points and
red line). We also show for comparison the $\langle M_{\textrm{H
{i}}}(M_{\ast})\rangle$ measurements for all central galaxies (gray points),
which are not used in the model fittings. The corresponding prediction of the
best-fitting model is shown as the gray dotted line.  Figure 3: Comparisons
between the best-fitting models and the observed measurements, for H2MF
$\phi(M_{\rm H_{2}})$ (left panel), H2 mass fraction $f_{\rm H_{2}}(M_{\ast})$
(middle panel) and H2-to-H i mass ratio $R_{\rm mol}$ (right panel) at $z\sim
0$. The best-fitting model is shown as the solid line in each panel. In
addition to the observational constraints of H2MF in Fletcher et al. (2021),
we also display in the left panel the measurements of Keres et al. (2003)
(crosses) and Andreani et al. (2020) (diamonds). The $R_{\rm mol}$
measurements of galaxies in the xCOLD GASS sample are shown as the gray
symbols in the right panel. The nondetections in the sample are displayed as
the open circles, using the upper limits of the mass measurements. Figure 4:
Comparison between observed cosmic gas densities and the best-fitting models
for H i (left) and H2 (right). The measurements from H i 21 cm emission and
intensity mapping are shown as red symbols in the left panel, while those from
DLA are shown as blue symbols. Similarly, the H2 densities of the CO line
luminosity and dust continuum are shown as red and blue symbols in the right
panel, respectively. The best-fitting models are shown as solid lines, while
the dotted lines indicate the model predictions without redshift evolution
(i.e., $\kappa_{1}$, $\kappa_{2}$, $M_{1}$, $M_{2}$, $\zeta_{1}$, and
$\zeta_{2}$ set as 0).

In the left panel of Figure 3, we show the best-fitting model for the H2MF
$\phi(M_{\rm H_{2}})$. For comparisons, we also display the H2MF measurements
at $z\sim 0$ from Keres et al. (2003) (crosses) and Andreani et al. (2020)
(diamonds). While the H2MF measurements of Keres et al. (2003) were estimated
from the CO luminosity function based on the FCRAO Extragalactic CO Survey,
the measurements of Andreani et al. (2020) were derived from the bivariate
K-band luminosity-H2 mass function of the Herschel Reference Survey. The
different measurements are consistent with each other within errors. But the
H2MF is most robustly measured with the xCOLD GASS sample (Fletcher et al.,
2021), where sample selection bias is carefully taken into account (Saintonge
& Catinella, 2022). Our model agrees with the measurements of Fletcher et al.
(2021), but predicts a somewhat higher $\phi(M_{\rm H_{2}})$ for $M_{\rm
H_{2}}<10^{8}M_{\sun}$ compared to Andreani et al. (2020). Since the H2MF
measurements at these low-mass ranges are quite uncertain, the discrepancies
can be verified with future surveys.

Reasonable agreement between data and the model is seen for the H2-stellar
mass relation, as well as the average H2-to-H i mass ratio $R_{\rm mol}$ in
the middle and right panels of Figure 3. The H2 fraction $f_{\rm H_{2}}$ of
star-forming galaxies is decreasing with $M_{\ast}$ from $13\%$ to $2\%$ in
the mass range probed by xCOLD GASS. Since the UniverseMachine model
accurately fits the observed galaxy stellar mass function with
$M_{\ast}>10^{7}M_{\sun}$ and the average H2 fraction is approaching a
constant of $13\%$ for low-mass galaxies, the mass resolution effect of
Bolshoi-Planck simulation is then negligible for $M_{\rm
H_{2}}>10^{6.1}M_{\sun}$, which is well below the current detection limit of
H2 gas.

The average H2-to-H i mass ratio $R_{\rm mol}$ has only a weak dependence on
$M_{\ast}$. The $R_{\rm mol}$ measurements of galaxies in the xCOLD GASS
sample are shown as gray symbols. The non-detections in the xCOLD GASS survey
are shown as the open circles using the upper limits provided in the catalog.
The small level of inconsistency between $R_{\rm mol}$ of our best-fitting
model and those of Catinella et al. (2018) is likely due to the treatment of
non-detections in xCOLD GASS, since only the upper limits for non-detections
are used in the estimation of $R_{\rm mol}$. The scatter around the mean
$R_{\rm mol}$ is still quite large and varies by two orders of magnitude. Our
model also predicts a similarly large scatter of $\sim 0.5$ dex for $R_{\rm
mol}$, consistent with the results of Catinella et al. (2018).

The predictions of the best-fitting model for $\rho_{\textrm{H {i}}}(z)$ and
$\rho_{\rm H_{2}}(z)$ are shown in Figure 4 as solid lines. The measurements
from H i 21 cm emission and intensity mapping are shown as red symbols in the
left panel of Figure 4, while those from DLAs are shown as blue symbols.
Similarly, the H2 densities from the CO line luminosity and the dust continuum
are shown as red and blue symbols in the right panel, respectively. To see the
effects of the redshift evolution of the model parameters in Equations (1) and
(15), we also show the model predictions without redshift evolution as green
dotted lines (that is, $\kappa_{1}$, $\kappa_{2}$, $M_{1}$, $M_{2}$,
$\zeta_{1}$ and $\zeta_{2}$ set as 0). It shows the importance of including
the redshift evolution, otherwise $\rho_{\textrm{H {i}}}$ would quickly
increase to the peak at around $z\sim 1.8$ and significantly decrease towards
higher redshifts, following the evolution of cosmic SFR densities (Behroozi et
al., 2019).

Our best-fitting model shows that $\rho_{\textrm{H {i}}}$ has a fast increase
by $\sim 0.25$ dex from $z=0$ to $z=1$, but it remains flat in $1<z<3$. The
increase and decrease of $\rho_{\textrm{H {i}}}$ in $3<z<5$ is caused by our
functional form of $\kappa$ and $M_{\rm crit}$ to fit the available data. More
complex functional forms as in Behroozi et al. (2019) might be necessary to
explore the behavior at these high redshifts, but they cannot be well
constrained by the current data. For even higher redshifts of $z>5$, our
current model is not accurate enough, because there is an increasing fraction
of H i gas outside the halo boundary (Villaescusa-Navarro et al., 2018). Our
halo-based model may not encompass the correct amount of H i gas at these
redshifts. As will be shown in the following section, the simulation
resolution of Bolshoi-Planck is also not able to resolve the ambient H i gas
in low-mass halos below $10^{10}M_{\sun}$ at these redshifts. Without the
redshift evolution of $\zeta$ in Eq. (15), $\rho_{\rm H_{2}}$ is decreasing
with redshift (green dotted line), following the trend of stellar mass
accretion. The molecular fraction of the cosmic neutral gas density is highest
around $z\sim 1.5$ and quickly decreases for both lower and higher redshifts.

We emphasize that the current measurements of $\rho_{\textrm{H {i}}}$ and
$\rho_{\rm H_{2}}$ from DLAs at high redshifts (especially $z>3$) have large
uncertainties. Our best-fitting model is based on all these available
measurements and will be improved with future cold gas surveys. Moreover,
other tracers of H i densities (e.g., the [C II] emission) have been suggested
to provide additional constraints to the H i gas within galaxies (Heintz et
al., 2021, 2022), which will be useful to further understand the neutral gas
distribution in and around galaxies.

###  3.2 Mock Observation

Since the ALFALFA sample only covers a small volume in the local universe, it
is important to investigate the effect of cosmic variance on the available
measurements (see e.g., Chauhan et al., 2019; Chen et al., 2019). Both HIMF
and H i clustering would be affected by the cosmic variance effect (Li et al.,
2022a). The H i clustering measurements would be underestimated on large
scales due to the integral constraint effect (see discussions in Section 3.3
of Guo et al., 2017). This is more severe for galaxies with lower
$M_{\textrm{H {i}}}$, as they are observed in smaller volumes due to lower H i
fluxes.

To fully account for this effect, we construct realistic mock catalogs as in
observations to measure the HIMF and H i clustering by adopting the best-
fitting model parameters. We first transform the Cartesian coordinates in the
UniverseMachine catalogs into celestial coordinates and then apply the same
geometry as in ALFALFA. The galaxy redshift is determined by including the
distortion from the LOS velocity of each subhalo. The mock line width $W_{50}$
of the H i profile for each galaxy is estimated to be $W_{50}=2V_{\rm
max}\sin(i)$, where $V_{\rm max}$ is the maximum circular velocity of each
halo/subhalo and $i$ is a random inclination angle. It is found in Chauhan et
al. (2019) (their Fig. 5) that $W_{50}$ has a strong correlation with $V_{\rm
max}$. But our results are not very sensitive to the details of the modeling
$W_{50}$. We further apply the 90% completeness limit of ALFALFA for the H i
flux at a given $W_{50}$ as in Eq. (5) of Haynes et al. (2011). We divide the
Bolshoi-Planck simulation into 64 sub-boxes. By applying the periodic boundary
condition, we then construct 64 mock catalogs by placing an observer at the
center of each sub-box.

To measure the HIMF in each mock, we apply the 2DSWML algorithm as in
observation to derive the effective volumes for mock galaxies (Jones et al.,
2018). The HIMF $\phi(M_{\textrm{H {i}}})$ is then simply calculated as the
sum of inverse volumes of galaxies in each $M_{\textrm{H {i}}}$ bin. The
resulting $\phi(M_{\textrm{H {i}}})$ for the 64 mocks are shown as the blue
lines in the left panel of Figure 5. The median HIMF is shown as the black
thick line, which agrees with our corrected ALFALFA HIMF (shown as red open
circles) as expected. The scatters among the different mock HIMFs are
relatively small, demonstrating that the cosmic variance effect is
significantly suppressed using the 2DSWML method. For comparison, we also
display the original ALFALFA HIMF of Jones et al. (2018) using the 50%
completeness cut as the yellow pluses. Although the two ALFALFA measurements
agree with each other for $M_{\textrm{H {i}}}>10^{9.5}M_{\sun}$, the low-mass
end HIMF is enhanced by 50% using the 90% completeness cut. But the Jones et
al. (2018) measurements are still within the range of mock variations. Using
the comprehensive H i mass estimator of Li et al. (2012), Li et al. (2022a)
also estimated the HIMF using the optical SDSS galaxies and correct for the
cosmic variance effect using the method of Chen et al. (2019). Their HIMF
(green crosses) agrees with our measurements for $M_{\textrm{H
{i}}}>10^{8.4}M_{\sun}$, but slightly higher for smaller $M_{\textrm{H {i}}}$.
This difference is likely due to the fact that the underdense distribution of
galaxies in the local universe would cause the galaxy stellar and H i mass
functions to be underestimated at the low mass end, as noted in Chen et al.
(2019) and Li et al. (2022a).

Figure 5: Left: Measurements of HIMF from 64 mock observations (blue lines).
The median of the mocks is shown as the black line, whereas our corrected
ALFALFA measurements are shown as the red open circles. For comparison, we
also show the original ALFALFA measurements of Jones et al. (2018) (yellow
pluses) using the 50% completeness cut and the HIMF estimates (green crosses)
of Li et al. (2022a) applying the H i mass estimator on the SDSS optical
galaxies. Right: Projected H i clustering measurements $w_{\rm p}(r_{\rm p})$
for the mock observations (blue lines) and their median values (black line).
The corresponding measurements for the ALFALFA sample with galaxies above the
90% completeness limit are shown as the red open circles.

In the right panel of Figure 5, we show the H i clustering measurements of the
64 mocks (blue lines), as well as the median measurement (black line). The
observed H i clustering measurements of ALFALFA are shown as red open circles.
Based on the method of Guo et al. (2017), Li et al. (2022b) measured the
projected two-point correlation function (2PCF), $w_{\rm p}(r_{\rm p})$, for
the H i-selected galaxies in ALFALFA, where $r_{\rm p}$ is the projected
separation of galaxy pairs and is measured in logarithmic bins of $\Delta\log
r_{\rm p}=0.2$ covering $0.1\,h^{-1}{\rm{Mpc}}$–$25.1\,h^{-1}{\rm{Mpc}}$. As
detailed in Li et al. (2022b), $w_{\rm p}(r_{\rm p})$ was obtained by
integrating the 3D 2PCF $\xi(r_{\rm p},r_{\rm\pi})$ along the line of sight
(LOS) to a distance of $r_{\rm\pi,max}=20\,h^{-1}{\rm{Mpc}}$ with $r_{\rm\pi}$
being the pair separation along LOS, i.e.,

$w_{\rm p}(r_{\rm p})=2\int_{0}^{r_{\rm\pi,max}}\xi(r_{\rm
p},r_{\rm\pi})dr_{\rm\pi},$ (23)

where $\xi(r_{\rm p},r_{\rm\pi})$ is measured using the Landy-Szalay estimator
(Landy & Szalay, 1993). To measure the clustering of H i gas and correct for
the selection bias of the ALFALFA galaxies, each galaxy pair (between galaxies
$i$ and $j$) should be weighted by $M_{{\textrm{H {i}}},i}M_{{\textrm{H
{i}}},j}/V_{ij}$, where $V_{ij}={\rm min}(V_{{\rm eff},i},V_{{\rm eff},j})$,
with $V_{{\rm eff},i}$ and $V_{{\rm eff},j}$ being the effective volumes
accessible to the galaxy pair. In this paper, we update their measurements by
using the ALFALFA galaxies above the 90% completeness limit. The effect of
completeness cuts on the H i clustering measurements is relatively minor.

The median $w_{\rm p}(r_{\rm p})$ measurements of the mock observations agree
well with the observed ALFALFA H i clustering. It demonstrates that the H i
clustering can be recovered with the information of the HIMF and H i-halo mass
relation. However, the cosmic variance effect is still very significant for
clustering measurements even with effective volume weights, necessitating
future H i surveys of much larger volumes. The best-fitting value of $\gamma$
($2.233^{+0.037}_{-0.137}$) strongly favors the dependence of H i clustering
on the halo formation time, confirms the findings of Guo et al. (2017) and Li
et al. (2022b). The $M_{\textrm{H {i}}}$–$z_{\rm form}$ relation in our model
is in agreement with the result of Guo et al. (2017) (their Fig. 12), where
the $z_{\rm form}$ dependence is constrained with the clustering measurements
in their extended subhalo abundance matching model.

###  3.3 Comparison with Literature

By fitting the observational constraints listed in Table 1, our model is able
to explain a set of important gas scaling relations at various redshifts.
However, comparison to other measurements not used in the model fitting will
provide an independent opportunity to further verify our model.

####  3.3.1 H I Measurements for Satellites

Figure 6: Left: Comparison between the observed $M_{\textrm{H {i}},{\rm
tot}}$–$M_{\rm vir}$ relation and the model predictions for halos of different
richness values ($N_{\rm g}$). The predictions of the model are shown as solid
lines, while the observational measurements of Guo et al. (2020) are displayed
as symbols of different colors. Right: Corresponding measurements for the
total H i masses of all satellite galaxies ($M_{\textrm{H {i}},{\rm sat}}$) in
each halo.  Figure 7: Differences in H i fraction with respective to the
values at virial radius, $\Delta\log(M_{\textrm{H {i}}}/M_{\ast})$, as a
function of halo-centric distance $r_{\rm cen}$ scaled by the halo virial
radius ($R_{\rm vir}$). We investigate the H i fraction in the same four
stellar mass bins as in Zhang et al. (2013), represented by lines of different
colors.

Guo et al. (2020) found that the total $M_{\textrm{H {i}},{\rm tot}}$ in halos
depends on both the halo mass and halo richness, with halos of higher richness
possessing larger amount of H i gas. They explain this behavior as the
dependence of H i mass on the halo assembly history. In the left panel of
Figure 6, we compare our model predictions (solid lines) for halos of
different richness values ($N_{\rm g}$) with measurements of Guo et al. (2020)
(symbols of different colors). We follow Guo et al. (2020) by applying a
stellar mass threshold of $M_{\ast}>10^{9}M_{\sun}$ to the UniverseMachine
catalog, to mimic the optical flux limit of SDSS observations. The halo
richness is then determined by calculating the number of galaxies (including
the orphans) with $M_{\ast}>10^{9}M_{\sun}$ for each halo. We compare our
model predictions to the observed measurements of Guo et al. (2020) after
applying the confusion correction. Our model is in remarkably good agreement
with observation. The slight overestimation of $M_{\textrm{H {i}},{\rm tot}}$
at $\log(M_{\rm vir}/M_{\sun})\sim 11.25$ in Figure 2 is improved, because
those halos with $N_{\rm g}=0$ (i.e., all member galaxies have masses smaller
than $10^{9}M_{\sun}$) also form later than halos hosting more massive
galaxies and thus have higher H i masses. The trend of $M_{\textrm{H {i}},{\rm
tot}}$ with halo richness is well reproduced. It was shown in Wechsler et al.
(2006) that the halos with later formation time have higher richness. Our
result further confirms that the dependence of H i mass on halo richness
originates from the dependence on halo formation time.

In our model constraints of H i content at $z\sim 0$, we include the H i-halo
mass relation for the total H i gas and the H i-stellar mass relation for
central galaxies. We do not explicitly constrain the H i content for satellite
galaxies. It is helpful to investigate how well the model works for the
satellite galaxies. In the right panel of Figure 6, we show the measurements
of the total H i mass contributed by the satellite galaxies (symbols of
different colors for halos of different $N_{\rm g}$), measured as the
subtraction between $M_{\textrm{H {i}},{\rm tot}}$ and $M_{\textrm{H {i}}}$ of
central galaxies in Guo et al. (2020). Our model predictions are shown as the
corresponding solid lines. Measurements and model predictions agree with each
other for $N_{\rm g}\geq 1$. The satellite H i masses for halos of higher
$N_{\rm g}$ are slightly underestimated in our model, but the agreement is
still reasonable considering the large errors.

Zhang et al. (2013) quantified the effect of H i depletion for cluster
galaxies using groups and clusters identified in SDSS. They found a smaller H
i fraction $M_{\textrm{H {i}}}/M_{\ast}$ for satellite galaxies with smaller
cluster-centric distances (their Fig. 5) (see also, Wang et al., 2020, 2021a).
This effect is stronger for lower-mass galaxies. The changes of $M_{\textrm{H
{i}}}/M_{\ast}$ with cluster-centric distance are around 0.37 dex and 0.68 dex
for $M_{\ast}\sim 10^{10.6}M_{\sun}$ and $M_{\ast}\sim 10^{9.8}M_{\sun}$,
respectively. In our model, the gas depletion effect is naturally included
along with the decrease of subhalo mass after infall. While the cluster-
centric distances in Zhang et al. (2013) suffer from the redshift-space
distortion effect, we can compare the relative change of $M_{\textrm{H
{i}}}/M_{\ast}$ for satellite galaxies from outer parts to inner parts.

In Figure 7, we show the differences in H i fraction with respective to the
values at the virial radius as a function of halo-centric distance ($r_{\rm
cen}$) scaled by the virial radius ($R_{\rm vir}$). We adopt the same four
stellar mass bins as in Zhang et al. (2013). We find similar levels of
decrease in $M_{\textrm{H {i}}}/M_{\ast}$ between our model and Zhang et al.
(2013) for all stellar mass bins, further confirming that the H i depletion of
satellite galaxies is reasonably reproduced in our empirical model. The
distributions of H i gas for satellite galaxies in halos of different masses
can also be investigated using the conditional HIMF as in Li et al. (2022a).
We will compare the predictions of the conditional HIMF in our future work.

####  3.3.2 HIMF in Groups

Figure 8: Comparison between the observed HIMFs for galaxies in groups of
different halo masses from Jones et al. (2020) (symbols) (updated using
ALFALFA galaxies above the 90% completeness limit) and the model predictions
(solid lines).

Jones et al. (2020) measured the HIMF for galaxies residing in groups by
making use of the overlapping regions between the ALFALFA survey and SDSS
group catalogs. In this paper, we have updated their measurements by using the
ALFALFA galaxies above 90% completeness limit. However, the calculation of
HIMF in Jones et al. (2020) is relative to the volume covered by the groups.
For fair comparisons with the HIMF shown in Figure 2, we normalized their HIMF
measurements using the volume of the entire galaxy sample. The galaxies in
Jones et al. (2020) were separated into field galaxies and group galaxies of
different halo masses. But the definition of field galaxies is somewhat
ambiguous and mainly refers to isolated galaxies in low-mass halos, which
might be affected by the flux limit of the SDSS sample. Therefore, we combine
their HIMF measurements of field galaxies with those in halos of $M_{\rm
vir}<10^{13}h^{-1}M_{\sun}$, and compare them to the HIMF measurements for
galaxies in more massive halos, shown as the colored symbols in Figure 8. Our
model predictions are shown as solid lines with the corresponding colors.

We find reasonable agreement between observation and model for different halo
mass bins. But our model predictions of HIMF in massive halos of
$13<\log(M_{\rm vir}/h^{-1}M_{\sun})<14$ are slightly underestimated. Since
more than $95\%$ of the H i-selected galaxies in these halos are satellites.
The misidentification of central and satellite galaxies in the SDSS group
catalog, as well as the errors in the halo mass estimates, could easily cause
the level of discrepancy seen in the figure (Campbell et al., 2015). We will
explore more accurate models for satellite galaxies in our future work. But
the overall trend of HIMF with the halo mass is reasonably reproduced in our
model. Jones et al. (2020) pointed out that the quick drop-off in the measured
HIMFs at the low mass end in massive halos is caused by the small volume of
the ALFALFA survey. The effect is more severe for larger halos. As confirmed
by our model prediction, the low-mass slopes of the HIMF measurements in
different halos are quite similar. But HIMFs in massive halos have much
shallower slopes at the massive end.

####  3.3.3 H I and H2 Measurements at High Redshifts

Figure 9: Comparison between observed data and best-fitting models for the
SFR–$M_{\ast}$ (left) and H i-$M_{\ast}$ relations (right) at $z\sim 0.37$ and
$z\sim 1.1$. The fitted H i-stack measurements of Chowdhury et al. (2022a) are
shown as filled cycles, while the other measurements of Chowdhury et al.
(2020), Rhee et al. (2018), Bera et al. (2019) and Sinigaglia et al. (2022)
are shown as open symbols. We display the best-fitting model at $z=0$ as the
black solid lines for comparison.

In the right panel of Figure 9, we show our best-fitting model of
$M_{\textrm{H {i}}}$-$M_{\ast}$ relation at $z\sim 0.37$ and $z\sim 1.1$ for
star-forming galaxies as the blue and red solid lines, respectively. The
fitted H i-stacking measurements of Chowdhury et al. (2022a) are shown as
filled circles. For comparison, we also display the H i-stacking measurements
of Chowdhury et al. (2020) at $z\sim 1.1$ and those of Rhee et al. (2018),
Bera et al. (2019), and Sinigaglia et al. (2022) at $z\sim 0.37$ as open
symbols. Except for the measurements of Sinigaglia et al. (2022) using the
MeerKAT radio telescope, all other measurements were observed with GMRT. For
comparison, the model predictions at $z=0$ are shown as the black solid line.
Although our model predicts a consistent slope of the H i-stellar mass
relation at $0<z<1$, the observed $M_{\textrm{H {i}}}$ measurements at $z\sim
0.37$ seem to have a weak dependence on $M_{\ast}$, with $M_{\textrm{H
{i}}}\sim 10^{9.7}M_{\sun}$.

We note that all these measurements were made for the star-forming galaxies,
as only these galaxies can have reliable measurements of H i gas even with the
stacking method. We show in the left panel of Figure 9 the available SFR
measurements from the corresponding references using the same symbols, and the
positions of the SFMS in the best-fitting models are shown as solid lines. In
Sinigaglia et al. (2022) it is claimed that the slope of the H i-stellar mass
relation is becoming flat at the massive end at $z\sim 0.37$, compared to the
measurements of $z=0$. However, Chowdhury et al. (2022b) found no evolution in
slope for the H i-stellar mass relation by comparing measurements between
$z=0$ and $z\sim 1.1$ (see also Bera et al., 2023). We emphasize that
$M_{\textrm{H {i}}}$ depends on both $M_{\ast}$ and SFR. It is still important
to check their SFR distributions when comparing the measurements with the
best-fitting models. Part of the decrease in $M_{\textrm{H {i}}}$ at
$M_{\ast}\sim 10^{11}M_{\sun}$ seen for Sinigaglia et al. (2022) can be
attributed to the lower SFR values relative to SFMS. High redshift
measurements of H i-stellar mass relation and HIMF (e.g., Bera et al., 2022)
are still scarce at the moment, but we expect the upcoming H i surveys to
provide more insight.

Figure 10: Comparison between observed H2 gas masses of Wang et al. (2022)
using the A3COSMOS catalog (symbols) and our model predictions (solid lines)
for main sequence galaxies. The measurements in different stellar mass bins
are shown as symbols and lines of different colors. For comparison, the
average H2 masses derived from the compilation of measurements in the
literature (Tacconi et al., 2020) in the mass range of
$11<\log(M_{\ast}/M_{\sun})<12$ are also shown as the gray circles.

Using a mass-complete sample of main sequence galaxies in the archival Atacama
Large Millimeter/submillimeter Array (ALMA) observations of COSMOS field
(A3COSMOS; Liu et al., 2019), Wang et al. (2022) measured the redshift
evolution of stacked molecular gas mass in different stellar mass bins. In
Figure 10, we show the comparison between their measurements (symbols) and our
model predictions (solid lines) for main sequence galaxies ($|\log({\rm
SFR}/{\rm SFR_{MS,obs}})|<0.5$). Measurements in different stellar mass bins
are shown as different colors.

We find good agreement between their measurements and our model predictions in
$0.5<z<3.5$. For comparison, we also derive the average H2 masses for the most
massive main-sequence galaxies ($11<\log(M_{\ast}/M_{\sun})<12$) using the
compilation of H2 measurements from Tacconi et al. (2020), as measurements for
these massive galaxies are more complete and representative. They are shown as
gray circles in Figure 10, also in agreement with our model predictions.

It is interesting that at a given stellar mass bin, $M_{\rm H_{2}}$ is
generally increasing with redshift, but the trend is becoming flatter at
$z>2$. Combining with the fact that there are significantly less massive
galaxies at higher redshifts (e.g., Figure 2 of Behroozi et al., 2019), the
cosmic molecular gas density $\rho_{\rm H_{2}}$ would then reach the peak at
round $z\sim 2$. It would be more intriguing to trace the evolution of gas
densities along with the galaxy stellar mass growth, as will be shown in the
following section.

###  3.4 Model Predictions

With the best-fitting model that is able to explain most H i and H2
observations, we can make valuable predictions about the evolution of the
properties of cold gas for the upcoming surveys.

####  3.4.1 Cosmic Baryon Budget

Figure 11: Fractional contributions of different baryon components to the
total cosmic baryon density, $\rho_{\rm dm}\Omega_{\rm b}/\Omega_{\rm m}$.
Contributions from stars, H i and H2 gas are shown as the orange, blue, and
pink lines, respectively. The overall contribution of $\rho_{\rm
baryon}=\rho_{\rm star}+1.36\rho_{\textrm{H {i}}}+\rho_{\rm H_{2}}$ is
displayed as the shaded area. Figure 12: Similar to Figure 11, but for the
ratios between the masses of different baryon components and the halo mass
$M_{\rm vir}$ at six typical redshifts.

Since our model includes the stellar, H i and H2 gas components, we are able
to quantify the evolution of cosmic baryon density associated with galaxies,
i.e., $\rho_{\rm baryon}=\rho_{\rm star}+1.36\rho_{\textrm{H {i}}}+\rho_{\rm
H_{2}}$. We have multiplied the H i density by a factor of 1.36 to account for
the contribution of helium and heavier elements. We can normalize the cosmic
densities with $\rho_{\rm dm}\Omega_{\rm b}/\Omega_{\rm m}$ ($\rho_{\rm dm}$
is the dark matter density) to quantify the fractional contributions of
different components to the total baryon density, as shown in Figure 11.
Different components are represented as lines of different colors, while the
overall contribution of $\rho_{\rm baryon}$ is displayed as the shaded area.

It is clear that the baryon densities associated with galaxies only account
for about 5% of the total baryon budget at $z=0$ and this fraction is slowly
decreasing at higher redshifts. Since our halo-based H i model has already
included the contribution of H i gas in the circumgalactic medium (CGM), the
majority baryons are in the form of hot gas distributed in the intergalactic
medium (IGM) and CGM (see more dicussions in Tumlinson et al., 2017).

We note that $\rho_{\rm baryon}$ is dominated by H i gas at $z>1.2$ and by
stars thereafter. The contribution of $\rho_{\rm H_{2}}$ to $\rho_{\rm
baryon}$ is generally less than $23\%$ and it is even dominating over
$\rho_{\rm star}$ at $z>2.6$. At low redshifts of $z<1.2$, the decrease of H i
and H2 gas is roughly comparable to stellar mass growth, which indicates the
persistent conversion of cold gas into stars in a quasi-steady state (Lilly et
al., 2013).

To better understand the roles of different components, we show in Figure 12
the ratios between the masses of different components and $M_{\rm vir}$, as
well as the combined contribution of stars and cold gas, $M_{\rm
baryon}=M_{\ast}+1.36M_{\textrm{H {i}}}+M_{\rm H_{2}}$, similarly as in Figure
11. In the calculation of mass for each component, we have included the
contribution from all galaxies in the halos. The ratios of $M_{\ast}/M_{\rm
vir}$ and $M_{\rm H_{2}}/M_{\rm vir}$ peak around $M_{\rm vir}\sim
10^{12}M_{\sun}$ and slightly shift toward more massive halos at higher
redshifts. The two ratios also have similar shapes at all redshifts, arising
from the tight correlation between $M_{\rm H_{2}}$ and $M_{\ast}$ in Eq. 15.

However, the peak of $M_{\textrm{H {i}}}/M_{\rm vir}$ occurs around $M_{\rm
vir}\sim 10^{11.4}M_{\sun}$ for $z<3$ and quickly moves to halos of lower mass
at higher redshifts. It highlights the importance of using high-resolution
simulations when modeling the H i gas at high redshifts. With the particle
mass resolution of $2.3\times 10^{8}M_{\sun}$, the Bolshoi-Planck simulation
is unable to correctly sample halos of $M_{\rm vir}<10^{10}M_{\sun}$ (i.e.,
with less than 50 particles). Our current model using Bolshoi-Planck
simulation is less accurate for $z\sim 5$, where the contribution of H i gas
to halos of $M_{\rm vir}<10^{10}M_{\sun}$ is significant.

At higher redshifts, particularly above $z=2$, we only have constraints on the
cosmic H i and H2 densities, as opposed to H i or H2 as a function of galaxy
mass. Therefore, it is useful to examine whether the relationship between the
gas mass and the halo mass in Fig. 12 continues to be reasonable at $z>2$. At
$z=3$–$5$, the H2 mass scales roughly proportionally to the stellar mass,
which is as expected for galaxies on a star-forming main sequence with sSFR
roughly independent of stellar mass. The behavior of H i is more interesting,
with the peak of $M_{\textrm{H {i}}}/M_{\rm vir}$ dramatically shifting to
lower halo masses, e.g., at $z=4$–$5$. We tested many alternate forms for the
redshift scaling of H i mass with the halo mass, but all forms that matched
the observational data required a similarly rapid shifting of the peak of
$M_{\textrm{H {i}}}/M_{\rm vir}$ to lower masses. Fig. 12 provides a simple
explanation for this behavior: even with the rapid change in the shape of the
$M_{\textrm{H {i}}}$–$M_{\rm vir}$ relation, the total baryon fraction sits
relatively close to the cosmic baryon fraction of 0.16. Therefore, to match
the high observed cosmic density of H i in DLAs, any less rapid shift in the
peak of $M_{\textrm{H {i}}}/M_{\rm vir}$ would result in halos where the total
mass in H i, H2, and stars exceeded the cosmic baryon fraction. Put another
way, the shape of the $M_{\textrm{H {i}}}$–$M_{\rm vir}$ relation at $z>2$ is
constrained above by the total cosmic baryon fraction and below by the high
observed number density of DLAs at high redshifts. Physically, we can
interpret this as cooling becoming more and more efficient in high-redshift
halos of all masses, such that a larger fraction of the available baryons in
the halos are cooling to H i and H2. High star formation rates in massive
halos (e.g., $M_{\rm vir}\sim 10^{12}M_{\sun}$) mean that most of the gas in
those halos is in the form of H2 rather than H i, and lower star formation
rates in less massive halos mean that most of the gas in the lower mass halos
will be in H i instead of H2.

####  3.4.2 Evolution of Cold Gas Mass Functions and Scaling Relations

In Figure 13, we show the predictions of the best fit model for the evolution
of H i (left panel) and H2 mass functions (right panel) from $z=3$ to $z=0$,
shown as lines of different colors. We have included the predictions of the
model at $z=0.5$ that could be compared with observations in the upcoming H i
surveys probing similar redshifts. The evolution of HIMF is relatively weak,
and the trends at low- and high-mass ends vary with redshift. But generally we
have a higher chance of detecting H i-rich galaxies in the future deep H i
surveys than those in the local universe (see e.g., Xi et al., 2021;
Ponomareva et al., 2023).

There is also weak evolution of H2MF for $M_{\rm H_{2}}<10^{9}M_{\sun}$, but
the evolution is much stronger at the massive end. There would be many more
galaxies with a large H2 reservoir at high redshifts. For comparison, we show
the collected measurements of H2MF at $2<z<2.5$ from Tacconi et al. (2020),
which agrees with our model predictions at $z\sim 2.5$ (shown as the dotted
line).

Figure 13: Best-fitting model predictions for the H i (left panel) and H2 mass
functions (right panel) at $0<z<3$, shown as lines of different colors. For
comparison, The collected measurements of H2MF at $2<z<2.5$ from Tacconi et
al. (2020) are shown as open circles in the right panel. Our model prediction
of H2MF at $z\sim 2.5$ is shown as the dotted line in the right panel. Figure
14: Best-fitting model predictions for the $M_{\textrm{H {i}}}$–$M_{\ast}$
(left panel) and $M_{\rm H_{2}}$–$M_{\ast}$ relations (right panel) at
$0<z<3$, shown as lines of different colors. The measurements of star-forming
and quenched galaxies are displayed as solid and dashed lines, respectively.
Figure 15: Average H i (left panel) and H2 masses (right panel) ($\langle\log
M_{\textrm{H {i}}}\rangle$ and $\langle\log M_{\rm H_{2}}\rangle$) as a
function of $M_{\ast}$ and $z$. The values of $\langle\log M_{\textrm{H
{i}}}\rangle$ and $\langle\log M_{\rm H_{2}}\rangle$ are color coded as
labeled.

In Figure 14, we show the scaling relations of $M_{\textrm{H {i}}}$–$M_{\ast}$
(left panel) and $M_{\rm H_{2}}$–$M_{\ast}$ (right panel) for the star-forming
(solid lines) and quenched galaxies (dashed lines) at various redshifts as in
Figure 13. We show in Figure 14 the model predictions of $\log\langle
M_{\textrm{H {i}}}\rangle$ and $\log\langle M_{\rm H_{2}}\rangle$. Star-
forming and quenched galaxies are separated by the demarcation line of
$\log({\rm SFR}/{\rm SFR_{MS,obs}})=-1$. The measurements at $z=0$ are
slightly different from the results in the bottom left panel of Figure 2, due
to the different selection cuts of star-forming and quenched galaxies.

The H i mass is generally increasing with $M_{\ast}$ and redshift, but the
slope of the $M_{\textrm{H {i}}}$–$M_{\ast}$ relation for massive galaxies
becomes shallower at higher redshifts. Since in our model the slopes $\alpha$
and $\beta$ do not vary with redshift, the shallower slope is caused by the
slope changes of the stellar-halo mass relation (Behroozi et al., 2019). The
offsets of $M_{\textrm{H {i}}}$ between star-forming and quenched galaxies
become larger at higher redshifts. It means that the massive end of HIMF is
mainly contributed by star-forming galaxies at high redshifts. It is then
reasonable to target star-forming galaxies for H i surveys at high redshifts.

The evolution of the $M_{\rm H_{2}}$–$M_{\ast}$ relation is similar to the
case of H i. But the dependence of $M_{\rm H_{2}}$ on $M_{\ast}$ for quenched
galaxies has almost no evolution, with a slope of $M_{\rm H_{2}}\propto
M_{\ast}$. The variation of $M_{\rm H_{2}}$ with redshift is larger for star-
forming galaxies, but there is a very weak evolution of $M_{\rm
H_{2}}$–$M_{\ast}$ relation for $z>2$. The change of $\rho_{\rm H_{2}}$ at
these high redshifts is due to the evolution of the stellar mass function.

For a complete description of the redshift evolution of the $M_{\textrm{H
{i}}}$–$M_{\ast}$ and $M_{\rm H_{2}}$–$M_{\ast}$ relations, we show in Figure
15 the average H i (left panel) and H2 masses (right panel), $\langle\log
M_{\textrm{H {i}}}\rangle$ and $\langle\log M_{\rm H_{2}}\rangle$, as a
function of $M_{\ast}$ and $z$. As shown in Figure 14, the high mass end slope
of the $M_{\textrm{H {i}}}$–$M_{\ast}$ relation is significantly flatter at
$z\sim 3$, slightly before the peak of cosmic SFR density, likely indicating
that the consumption of H i gas is much higher than the cooling and cold gas
accretion for these massive galaxies. The evolutionary trend for the $M_{\rm
H_{2}}$–$M_{\ast}$ relation has a roughly constant slope, originating from the
H2 model in Eq. (15).

Figure 16: Similar to Figure 14, but for the depletion time scales of H i
(left panel) and H2 (right panel). For comparison, we show the measurements of
depletion time scales for star-forming galaxies from Guo et al. (2021) at
$z\sim 0$ as the open circles.

Another important parameter is the cold gas depletion time scale $\tau$.
Similarly to Figure 14, we show in Figure 16 the depletion time scales for H i
($\log\tau_{\textrm{H {i}}}\equiv\log\langle M_{\textrm{H
{i}}}\rangle-\langle\log{\rm SFR}\rangle$, left panel) and H2 ($\log\tau_{\rm
H_{2}}\equiv\log\langle M_{\rm H_{2}}\rangle-\langle\log{\rm SFR}\rangle$,
right panel). We note that the H i and H2 masses are measured as stacked
averages in each mass bin for direct comparisons with the H i observations of
Guo et al. (2021) at $z\sim 0$ (shown as open circles). Massive galaxies
typically have lower $\tau_{\textrm{H {i}}}$ and the dependence on $M_{\ast}$
is stronger at higher redshifts. The H i depletion time scale for star-forming
galaxies varies from 0.1 Gyr to 10 Gyr, but that of the quenched galaxies is
significantly longer. These quenched galaxies can hardly deplete their H i
reservoir (albeit small) with low SFRs.

The dependences of the H2 depletion time scale $\tau_{\rm H_{2}}$ on
$M_{\ast}$ and redshift are much weaker. The star-forming and quenched
galaxies have similar depletion time scales around 0.3–3 Gyr and it is
slightly decreasing with redshift, as also shown in Tacconi et al. (2020)
(their Figure 3). Similar levels of $\tau_{\rm H_{2}}$ were also found in the
empirical model of Padmanabhan & Loeb (2020). The smallest $\tau_{\rm H_{2}}$
is found in galaxies of $\sim 10^{10.5}M_{\sun}$ at $z=2$, indicating the fast
consumption of molecular gas for the formation of stars. We note that the
lines in Figure 16 are the average values of the depletion time scales. The
scatters for $\log\tau_{\textrm{H {i}}}$ and $\log\tau_{\rm H_{2}}$ for
individual galaxies are around 0.8 dex and 0.4 dex, respectively.

####  3.4.3 Cold Gas Accretion Histories

Figure 17: Best-fitting model predictions for the growth histories of central
galaxies for the stellar mass (top left panel), SFR (top right panel, in units
of $M_{\sun}/{\rm Gyr}$), H i mass (bottom left panel) and H2 mass (bottom
right panel). The blue, orange and pink lines are for halos with
$10.8<\log(M_{\rm vir}/M_{\sun})<11.2$, $11.8<\log(M_{\rm vir}/M_{\sun})<12.2$
and $12.8<\log(M_{\rm vir}/M_{\sun})<13.2$ at $z=0$, respectively. The
measurements for star-forming galaxies are shown as solid lines, while those
of quenched galaxies are displayed as dotted lines.

Using the halo merger trees of the Bolshoi-Planck simulation, we are able to
trace the evolution history of the cold gas content for individual halos. In
Figure 17, we show the growth histories of central galaxies for the stellar
mass (top left panel), SFR (top right panel), H i mass (bottom left panel),
and H2 mass (bottom right panel). The blue, orange and pink lines are for
halos with $10.8<\log(M_{\rm vir}/M_{\sun})<11.2$, $11.8<\log(M_{\rm
vir}/M_{\sun})<12.2$ and $12.8<\log(M_{\rm vir}/M_{\sun})<13.2$ at $z=0$,
respectively. The measurements for star-forming galaxies are shown as solid
lines, while those of quenched galaxies are displayed as dotted lines.

The stellar mass growth histories for halos of different masses are quite
similar. In low-mass halos of $M_{\rm vir}\sim 10^{11}M_{\sun}$, the quenched
galaxies have $\sim 0.3$ dex higher $M_{\ast}$ than their star-forming
counterparts. For more massive halos, there are no strong differences in the
mass growth histories of the two populations. But the evolution of SFR varies
significantly in halos of different masses. The SFR history peaks are shifting
from low to high redshifts, with increasing $M_{\rm vir}$, as also seen in
Figure 16 of Behroozi et al. (2019). Significant differences between star-
forming and quenched galaxies only happen since $z=0.5$ (i.e., within $\sim$5
Gyr of $z=0$). Along with the general decrease of SFR toward $z=0$ after
reaching the peaks, the stellar mass growth becomes slower, shown as the
shallower curves in the top left panel. Thus, even though the quenched
galaxies have much smaller SFRs at $z=0$ than the star-forming ones, their
final stellar masses are still quite similar. The higher $M_{\ast}$ for
quenched galaxies at $z=0$ in low-mass halos is caused by their higher SFRs at
$z>0.5$.

When we compare the cold gas accretion histories shown in the bottom panels of
Figure 17, it is clear that the evolution of molecular gas closely follows the
star formation histories, as expected from the tight correlation in Equation
(15). The evolution of the H i gas is also similar to that of molecular gas,
but its dependence on star formation histories is much weaker, especially for
massive halos. In halos of $M_{\rm vir}>10^{12}M_{\sun}$, most of the H i gas
has already been accreted much earlier than the redshifts of peak star
formation histories, with a fast growing phase at $z>4$. The H i gas will then
be converted to H2 and form stars. The star-forming galaxies in low-mass halos
of $M_{\rm vir}\sim 10^{11}M_{\sun}$ are still experiencing rapid accretion of
H igas at low redshifts, leading to increasing SFRs toward $z=0$.

The time scale required to transform H i to H2 will increase the time lag
between the peaks of H i and H2 accretion histories, which is about 0.8 Gyr
for star-forming galaxies and 1 Gyr for quenched galaxies. The H i gas starts
to be depleted well before the star formation quenching occurs. The sharp
decrease of $M_{\textrm{H {i}}}$ for quenched galaxies at low redshifts could
be caused by additional physical mechanisms, such as stellar and AGN feedback
(Guo et al., 2022; Ma et al., 2022).

It is interesting that the SFRs of star-forming galaxies in massive halos of
$M_{\rm vir}\sim 10^{13}M_{\sun}$ will decrease from $z=2$ to $z=0.5$, but
increase again since $z=0.5$. This is mainly due to the selection effect of
star-forming galaxies. As most of the galaxies in these massive halos have low
SFRs, selecting star-forming galaxies will tend to pick out those that have
undergone recent rejuvenation effects, leading to the upturn seen in the star
formation history (see also, Fig. 16 of Behroozi et al., 2019). Yet, the SFRs
for these galaxies are still below the SFMS and lying in the ‘green valley’
region, consistent with the recent observations of rejuvenation events (Chauke
et al., 2019). As seen in the bottom panels of Figure 17, the rejuvenation is
accompanied by efficient accretion of cold gases in massive halos in the
recent 5 Gyr.

####  3.4.4 H I Bias

Figure 18: Best-fitting model predictions for the H i bias (left panel) and
the combination of $\Omega_{\textrm{H {i}}}b_{\textrm{H {i}}}$ (right panel)
at $0<z<3$. The solid and dotted lines in the right panel are for the model
predictions in the real and redshift spaces, respectively. We also show for
comparison the measurements of Wolz et al. (2022) and Amiri et al. (2023) as
different symbols.

Our halo-based H i model can be used to predict the evolution of H i bias, as
well as making predictions for the future H i intensity mapping surveys. The
power spectrum measurements of H i intensity mapping experiments constrain the
product of H i bias ($b_{\textrm{H {i}}}$) and cosmic H i abundance
($\Omega_{\textrm{H {i}}}\equiv\rho_{\textrm{H {i}}}/\rho_{\rm c}$), where
$\rho_{\rm c}$ is the critical density at $z=0$. Cross-correlations between H
i signals and optical galaxy samples have frequently been applied in intensity
mapping measurements (Chang et al., 2010; Masui et al., 2013; Switzer et al.,
2013; Amiri et al., 2023; Wolz et al., 2022; Cunnington et al., 2023), because
autocorrelation measurements would suffer from more severe systematic
uncertainties.

In our model, we can calculate the H i bias as,

$b_{\textrm{H {i}}}(r)=\sqrt{\frac{\xi_{\textrm{H {i}}}(r)}{\xi_{\rm
dm}(r)}},$ (24)

where $r$ is the separation between galaxy pairs. $\xi_{\textrm{H {i}}}(r)$ is
the H i mass weighted real-space galaxy two-point correlation function and
$\xi_{\rm dm}(r)$ is the corresponding measurement for dark matter. To
calculate the large-scale linear H i bias, we use the correlation functions in
$3.25\,h^{-1}{\rm{Mpc}}<r<12.92\,h^{-1}{\rm{Mpc}}$ to derive the average
$b_{\textrm{H {i}}}$.

In the left panel of Figure 18, we show the model prediction of $b_{\textrm{H
{i}}}$ as the solid line. The H i bias gradually increases from $0.69$ at
$z=0$ to 2.33 at $z=3$, which is in agreement with the semi-analytical model
prediction of Wang et al. (2021b). The critical value of $b_{\textrm{H
{i}}}=1$ is reached at $z\sim 0.85$. As shown in Guo et al. (2017), the
clustering amplitudes of H i-selected galaxies will increase with
$M_{\textrm{H {i}}}$. The overall H i bias then depends on the galaxy
population hosting the H i gas. Therefore, correctly modeling the H i-halo
mass relation is essential to obtain the accurate H i bias.

Our model prediction for the combined value of $\Omega_{\textrm{H
{i}}}b_{\textrm{H {i}}}$ probed by the H i intensity mapping experiments is
shown in the right panel of Figure 18 as the solid line. We also show for
comparison the recent measurements of Wolz et al. (2022) and Amiri et al.
(2023) as symbols of different colors.

Figure 19: Host halo mass probability distributions of a random H i atom (left
panel) and a random H2 molecule (right panel) in the universe. We show the
logarithmic probabilities $\log P_{\textrm{H {i}}}(M_{\rm vir}|z)$ and $\log
P_{\rm H_{2}}(M_{\rm vir}|z)$ as contours of different colors. See the text
for details.

Wolz et al. (2022) cross-correlated the H i intensity mapping data from the
Green Bank Telescope (GBT) with the optical galaxies from WiggleZ Dark Energy
Survey (Drinkwater et al., 2010) in the redshift range of $0.6<z<1$, based on
the work of Masui et al. (2013). They found that $\Omega_{\textrm{H
{i}}}b_{\textrm{H {i}}}r_{\textrm{H {i}}}=(0.58\pm 0.14)\times 10^{-3}$ on the
effective scale $k_{\rm eff}=0.31h\,{\rm Mpc}^{-1}$, where $r_{\textrm{H
{i}}}$ is the galaxy-H i cross-correlation coefficient. The main uncertainty
lies in the estimate of $r_{\textrm{H {i}}}$, which depends on both the H i
and the galaxy samples. We adopt their estimate of $r_{\textrm{H {i}}}=0.9$
and derive the combination of $\Omega_{\textrm{H {i}}}b_{\textrm{H
{i}}}=(0.64\pm 0.16)\times 10^{-3}$ (consistent with the measurements of
Switzer et al. 2013), which shows good agreement with the solid line.

Amiri et al. (2023) made the first H i intensity mapping measurements with the
interferometer from the Canadian Hydrogen Intensity Mapping Experiment (CHIME;
CHIME Collaboration et al., 2022). By cross-correlating the H i data with
luminous red galaxies (LRGs), emission line galaxies (ELGs) and quasars (QSOs)
from the eBOSS sample (Dawson et al., 2016), they can constrain the effective
H i clustering amplitude $\mathcal{A}_{\textrm{H {i}}}\equiv
10^{3}\Omega_{\textrm{H {i}}}(b_{\textrm{H {i}}}+\langle f\mu^{2}\rangle)$,
where the term $\langle f\mu^{2}\rangle$ is to account for the redshift-space
distortion effect in observation, with $f$ the linear growth rate and $\mu$
the cosine of angle between the wavevector $\mathbf{k}$ and line of sight.
They measured $\mathcal{A}_{\textrm{H {i}}}$ for LRGs ($\langle
z\rangle=0.84$), ELGs ($\langle z\rangle=0.96$) and QSOs ($\langle
z\rangle=1.20$) that provide constraints at higher redshifts than previous
measurements. For fair comparisons, we present our model prediction of
$\Omega_{\textrm{H {i}}}b_{\textrm{H {i}}}$ for H i bias measured in the
redshift space as dotted line in the right panel of Figure 18. The redshift-
space distortion effect would increase the measured $b_{\textrm{H {i}}}$ by
about 20%. The measurements of Amiri et al. (2023) are slightly higher than
the predictions of our model, but generally in agreement within the errors. As
discussed in Amiri et al. (2023), the higher amplitude of
$\mathcal{A}_{\textrm{H {i}}}$ is likely caused by the degeneracy between
$\mathcal{A}_{\textrm{H {i}}}$ and their model parameters that describe the
Fingers-of-God effect on non-linear scales.

To better understand the distributions of H i and H2 in the universe, we show
in Figure 19 the probability distributions of H i and H2 gas in halos of
different masses, i.e., $P_{\textrm{H {i}}}(M_{\rm vir}|z)$ (left panel) and
$P_{\rm H_{2}}(M_{\rm vir}|z)$ (right panel). At each redshift, we calculate
the probabilities as follows,

$\displaystyle P_{\textrm{H {i}}}(M_{\rm vir}|z)$ $\displaystyle=$
$\displaystyle\frac{M_{\textrm{H {i}}}dn}{d\log M_{\rm vir}}$ (25)
$\displaystyle P_{\rm H_{2}}(M_{\rm vir}|z)$ $\displaystyle=$
$\displaystyle\frac{M_{\rm H_{2}}dn}{d\log M_{\rm vir}},$ (26)

where $dn/d\log M_{\rm vir}$ is the normalized halo mass function. Therefore,
$P_{\textrm{H {i}}}(M_{\rm vir}|z)$ and $P_{\rm H_{2}}(M_{\rm vir}|z)$
indicate the host halo mass distributions of any given H i atom and H2
molecule, respectively. As shown in Figure 19, the majority of H i gas is
distributed in halos of $10^{11}$–$10^{12}M_{\sun}$ in $0<z<3$. The
probability $P_{\textrm{H {i}}}(M_{\rm vir}|z=0)$ peaks at $M_{\rm vir}\sim
10^{11.3}M_{\sun}$, corresponding to a halo bias of $b\sim 0.7$ (Guo et al.,
2017), which is in agreement with the H i bias in Figure 18. At higher
redshifts, the peak halo mass of the H i gas gradually moves to lower values.
Our model based on Bolshoi-Planck simulation is not accurate enough at $z>5$,
as the peak halo mass drops below $10^{10}M_{\sun}$, as also indicated in
Figure 12.

The peak halo mass for the H2 molecule is higher than that of H i. It is
around $M_{\rm vir}\sim 10^{12}M_{\sun}$ at $0<z<4$, which is consistent with
the host-halo distributions of the cosmic SFR (Behroozi et al., 2019, their
Figure 13). We would expect that the H2 bias is higher than the H ibias, which
can be verified in future H2 surveys. The resolution of Bolshoi-Planck
simulation is also good enough to model the H2 gas.

##  4 Discussion

Besides the halo properties of $M_{\rm vir}$ and $z_{\rm form}$ that can be
easily measured in the $N$-body simulations, our cold gas model requires
additional information of $M_{\ast}$ and ${\rm SFR}$. All these galaxy and
halo properties can be found in the UniverseMachine catalog, making it a
perfect starting point to apply the gas model. However, the
NeutralUniverseMachine model can also be applied to other theoretical models
that include relevant information, e.g., semi-analytical models and
hydrodynamical simulations. The parameter ${\rm SFR_{MS,obs}}$ in our model
can be replaced with the corresponding SFMS in different models. The exact
functional form of ${\rm SFR_{MS,obs}}$ as in Eq. 5 is not essential for cold
gas, since only the offset from SFMS is used to scale the H i and H2 masses.

Although our empirical model lacks detailed baryon physics, it has the
advantage of correctly describing the observed cold gas properties, and the
model constraints can be further improved with future observations. Since cold
gas serves as the fuel for galactic star formation, it is very sensitive to
the complicated physical processes involved in galaxy formation and evolution.
Matching the current observational measurements of H i and H2 gas is still
challenging for theoretical models.

The first step is to match the observed HIMF and H2MF. As shown in Davé et al.
(2020), current hydrodynamical simulations can reasonably fit the $z=0$ HIMF
for $M_{\textrm{H {i}}}>10^{10}M_{\sun}$, but the discrepancies are greater
for smaller H i masses (e.g., Crain et al., 2017; Diemer et al., 2019). The
situation is similar when comparing the semi-analytical models that are based
on the $N$-body simulations and incorporate simple recipes to describe the
baryon physics in the galaxy formation and evolution (e.g., Fu et al., 2013;
Popping et al., 2014; Kim et al., 2017; Baugh et al., 2019; Chauhan et al.,
2020; Spinelli et al., 2020). The comparison with observed H2MF is relatively
better (see e.g., Lagos et al., 2015; Diemer et al., 2019; Davé et al., 2020),
since the molecular gas is determined mainly by stellar mass and SFR.

The H i-halo mass and H i-stellar mass relations provide more stringent tests
of theoretical models. The HIMF and H2MF set the probability distributions of
cold gas in the universe, while the scaling relations with stellar and halo
masses are more sensitive to the evolution of the baryon cycle. The H i
content in massive halos of $M_{\rm vir}>10^{12}M_{\sun}$ varies significantly
in different hydrodynamical simulations and semi-analytical models (e.g.,
Villaescusa-Navarro et al., 2018; Baugh et al., 2019; Obuljen et al., 2019;
Chauhan et al., 2020; Spinelli et al., 2020; Li et al., 2022b), due to the
various levels of gas accretion, heating and cooling. Comparing the gas
scaling relations for star-forming and quenched galaxies also helps
distinguish the dependence of cold gas on the global SFR (Ma et al., 2022).

The inclusion of halo formation time $z_{\rm form}$ is an essential component
of our empirical model, which was not considered in previous theoretical
models. The H i clustering measurements $w_{\rm p}(r_{\rm p})$ provide further
verification, as they are sensitive to the parameter $\gamma$ related to the
halo formation time. Accurate measurements of H i-halo and H i-stellar mass
relations of stacked H i signals have already placed tight constraints on
$\gamma$. If we fix $\gamma=0$, i.e., ignoring the dependence of $M_{\textrm{H
{i}}}$ on $z_{\rm form}$, the H i content in massive halos would be
overestimated and the shape of the H i-stellar mass relation would also be
inconsistent with observation. Furthermore, without dependence on $z_{\rm
form}$, the H i content in halos of different richness would be quite similar
to each other, different from the results shown in Figure 6.

Overall, the empirical model is very powerful in accurately capturing the
dependence of cold gas content on the various galaxy and halo properties,
which is not a simple task for the methods of hydrodynamical simulations and
semi-analytical models. But it is still important to understand the physical
processes that determine the cold gas content. The semi-empirical approach
proposed by Popping et al. (2015) shows some success in the attempt to combine
the empirical model of galaxy star formation histories with a physically
motivated cold gas model. Despite the failure to describe the observed HIMF at
$z\sim 0$ for $M_{\textrm{H {i}}}<10^{10}M_{\sun}$ (caused by the
overcorrection for the number of low-mass galaxies in Behroozi et al. 2013a)
and the cosmic evolution of $\rho_{\textrm{H {i}}}$, it is still worth the
effort to explore such an approach by employing more accurate physical models
in the future.

##  5 Conclusions

In this paper, we propose a new empirical model that is capable of accurately
describing the various statistics for the H i and H2 gas content in the
redshift range of $0<z<6$. The functional form of the empirical model is
motivated by the various observations describing the scaling relations between
cold gas mass and properties of their host galaxies and halos. Our results are
summarized as follows.

(i) Our empirical model can accurately describe the H i and H2 mass functions,
molecular-to-atomic mass ratio, H i-halo mass relation, H i\- and H2-stellar
mass relations, and H i clustering measurements at $z\sim 0$. Higher redshift
measurements of the H i-stellar mass relations, as well as the cosmic gas
densities of $\rho_{\textrm{H {i}}}$ and $\rho_{\rm H_{2}}$, are also well
reproduced. Our best-fitting model is further verified with a few sets of H i
and H2 measurements not used in the modeling constraints and shows good
agreement with all these different observations.

(ii) There is only weak evolution in the HIMF from $z=0$ to $z=3$. However,
the evolution of H2MF is significantly larger at the massive end and smaller
for $M_{\rm H_{2}}<10^{9}M_{\sun}$. The average $M_{\textrm{H {i}}}$ and
$M_{\rm H_{2}}$ increase by around 1 dex from $z=0$ to $z=3$ for star-forming
galaxies, but there is much weaker evolution for the quenched population.

(iii) The H i gas depletion time $\tau_{\textrm{H {i}}}$ generally decreases
with increasing stellar mass, and varies from 0.1 Gyr to 10 Gyr for the star-
forming galaxies. The H2 gas depletion time $\tau_{\rm H_{2}}$ has a weaker
dependence on the redshift and stellar mass. The quenched galaxies have much
longer H i gas depletion time, varying from 10 Gyr to 200 Gyr, i.e. they are
not likely to fully deplete their H i reservoir with the low SFRs.

(iv) From the growth histories of the galaxy stellar mass, SFR, H i and H2
masses in halos of different masses, we find that $M_{\rm H_{2}}$ closely
trace the evolution of the SFR, but the correlation between $M_{\textrm{H
{i}}}$ and SFR is weaker. There is also an apparent time lag between the
evolution trends of $M_{\textrm{H {i}}}$ and SFR, with $M_{\textrm{H {i}}}$
reaching the peaks earlier than SFR.

(v) The cosmic baryon density associated with galaxies is dominated by stars
for $z<1.2$, and mainly contributed by H i gas at higher redshifts. But they
only account for less than 5% of the total baryon budget. The ratios of
$M_{\ast}/M_{\rm vir}$ and $M_{\rm H_{2}}/M_{\rm vir}$ closely follow each
other and reach the peaks around $M_{\rm vir}\sim 10^{12}M_{\sun}$, while
$M_{\textrm{H {i}}}/M_{\rm vir}$ peaks around $M_{\rm vir}\sim
10^{11.4}M_{\sun}$ for $z<3$ and shifts to lower-mass halos at higher
redshifts.

(vi) Our model can predict the evolution of H i clustering in the universe.
The H i bias $b_{\textrm{H {i}}}$ is gradually evolving from $0.69$ at $z=0$
to $2.33$ at $z=3$. The combined value $\Omega_{\textrm{H {i}}}b_{\textrm{H
{i}}}$ increases from $0.24\times 10^{-3}$ to $1.81\times 10^{-3}$ in the same
redshift range and shows good agreement with recent H i intensity mapping
measurements.

We thank the anonymous reviewer for helpful suggestions that significantly
improve the presentation of this paper. This work is supported by the National
SKA Program of China (grant No. 2020SKA0110100), National Science Foundation
of China (Nos. 11922305, 11833005, 12073002, 11721303, 12011530159) and the
science research grants from the China Manned Space Project with NOs. CMS-
CSST-2021-A02 and CMS-CSST-2021-B01. We thank Toby Brown, Nissim Kanekar, and
Kasper E. Heintz for helpful discussions. We acknowledge the use of the High
Performance Computing Resource in the Core Facility for Advanced Research
Computing at the Shanghai Astronomical Observatory.

##  Appendix A ALFALFA H I Mass Function with 90% Completeness Cut

Table 2: ALFALFA H i Mass Function with 90% Completeness $\log M_{\textrm{H {i}}}$ | $\phi(M_{\textrm{H {i}}})$ | Error on $\phi(M_{\textrm{H {i}}})$ | $N_{\rm gal}$
---|---|---|---
$M_{\sun}$ | $\rm{Mpc^{-3}dex^{-1}}$ | $\rm{Mpc^{-3}dex^{-1}}$ |
7.1 | $8.509\times 10^{-2}$ | $2.038\times 10^{-2}$ | 41
7.3 | $7.433\times 10^{-2}$ | $2.450\times 10^{-2}$ | 52
7.5 | $8.045\times 10^{-2}$ | $2.298\times 10^{-2}$ | 106
7.7 | $5.629\times 10^{-2}$ | $1.252\times 10^{-2}$ | 134
7.9 | $4.705\times 10^{-2}$ | $0.943\times 10^{-2}$ | 172
8.1 | $4.581\times 10^{-2}$ | $0.586\times 10^{-2}$ | 245
8.3 | $3.851\times 10^{-2}$ | $0.395\times 10^{-2}$ | 308
8.5 | $3.214\times 10^{-2}$ | $0.279\times 10^{-2}$ | 362
8.7 | $3.049\times 10^{-2}$ | $0.252\times 10^{-2}$ | 568
8.9 | $2.412\times 10^{-2}$ | $0.193\times 10^{-2}$ | 827
9.1 | $2.346\times 10^{-2}$ | $0.153\times 10^{-2}$ | 1545
9.3 | $1.707\times 10^{-2}$ | $0.109\times 10^{-2}$ | 2174
9.5 | $1.132\times 10^{-2}$ | $0.072\times 10^{-2}$ | 2679
9.7 | $8.335\times 10^{-3}$ | $0.464\times 10^{-3}$ | 3479
9.9 | $5.079\times 10^{-3}$ | $0.297\times 10^{-3}$ | 3595
10.1 | $2.687\times 10^{-3}$ | $0.162\times 10^{-3}$ | 2795
10.3 | $8.495\times 10^{-4}$ | $0.816\times 10^{-4}$ | 1058
10.5 | $1.687\times 10^{-4}$ | $0.202\times 10^{-4}$ | 219
10.7 | $2.386\times 10^{-5}$ | $0.357\times 10^{-5}$ | 31

We show in Table 2 the corrected ALFALFA HIMF using galaxies above the 90%
completeness cut. The number of galaxies in each H i mass bin is also listed
in the last column. Even by applying the stricter cut, there are still fair
amount of galaxies in the ALFALFA final sample to achieve an accurate
estimation of the HIMF. To better compare with the previous ALFALFA HIMF of
Jones et al. (2018), we also fit a Schechter function to our measurements,

$\phi(M_{\textrm{H {i}}})=\ln(10)\phi_{\rm s}\left(\frac{M_{\textrm{H
{i}}}}{M_{\rm s}}\right)^{\alpha_{\rm s}+1}\exp(-\frac{M_{\textrm{H
{i}}}}{M_{\rm s}}).$ (A1)

Our best-fitting parameters are $\alpha_{\rm s}=-1.30\pm 0.02$, $\log(M_{\rm
s}/M_{\sun})=9.91\pm 0.01$, $\phi_{\rm s}=5.93\pm 0.29\times 10^{-3}{\rm
Mpc}^{-3}{\rm dex}^{-1}$. The corresponding parameters in Jones et al. (2018)
are $\alpha_{\rm s}=-1.25\pm 0.02$, $\log(M_{\rm s}/M_{\sun})=9.94\pm 0.01$,
$\phi_{\rm s}=4.5\pm 0.2\times 10^{-3}{\rm Mpc}^{-3}{\rm dex}^{-1}$. In the
corrected HIMF, the low-mass end slope $\alpha_{\rm s}$ is becoming slightly
steeper and the ‘knee’ mass is shifting to lower values.

Figure 20: Comparisons of Schechter function fittings of Martin et al. (2010)
(black line), Jones et al. (2018) (blue line) and our results (red line). The
corrected HIMF for the ALFALFA final sample is also shown as yellow symbols.

In Figure 20, we show comparisons of the Schechter function fittings between
our results (red line) and that of Jones et al. (2018) (blue line), as well as
the measurements of Martin et al. (2010) (black line) using the previous
ALFALFA release covering 40% of the final sample. Martin et al. (2010) also
applied the 2DSWML method to derive the HIMF, but imposed a strict 100%
completeness threshold, similar to our 90% completeness cut.

Our model fitting is in good agreement with that of Martin et al. (2010) (see
also, Oman, 2022), while the fitting from Jones et al. (2018) is
systematically lower by around 50% for $M_{\textrm{H {i}}}<10^{9}M_{\sun}$. It
emphasizes the importance of selecting complete samples or correcting for
incompleteness when applying the 2DSWML method. However, we also note that
there is still a large difference ($\sim 30\%$) between the HIMF measurements
for the ALFALFA spring and fall regions at the low-mass end, even with the
2DSWML method (see e.g., Fig. 3 of Jones et al., 2018). Larger samples in
future H i surveys are essential to fully reduce the impact of cosmic
variance.

The cosmic H i abundance $\Omega_{\textrm{H {i}}}$ can be obtained from
integrating the Schechter function (Martin et al., 2010; Jones et al., 2018),

$\Omega_{\textrm{H {i}}}=\frac{1}{\rho_{\rm c}}\int M_{\textrm{H
{i}}}\phi(M_{\textrm{H {i}}})dM_{\textrm{H {i}}}=\frac{\phi_{\rm s}}{\rho_{\rm
c}}M_{\rm s}\Gamma(\alpha_{\rm s}+2),$ (A2)

which gives $\Omega_{\textrm{H {i}}}=(4.55\pm 0.29)\times 10^{-4}$ for our
best-fitting parameters assuming $h=0.7$, which is about 29% higher than the
corresponding value ($3.5\times 10^{-4}$) of Jones et al. (2018) before the
correction of H i self-absorption.

##  Appendix B Probability Distributions of Cold Gas Model Parameters

Figure 21: Pairwise Density distributions of the NeutralUniverseMachine model
parameters of the MCMC runs. We separate the model parameters for $z=0$
(bottom left) and higher redshifts (top right), as they are not jointly
fitted. The best-fitting model parameters are shown as the black plus in each
panel. The distributions are color-coded by the number of galaxies in each
grid, with redder colors for higher densities.

In Figure 21, we show the pairwise density distributions of the
NeutralUniverseMachine model parameters from the MCMC runs. The model
parameters for $z=0$ and higher redshifts are shown in the bottom left and top
right, respectively. Correlations are generally weak between different
parameters. But there are relatively strong correlations for the
$\kappa_{0}$–$\gamma$, $\zeta_{0}$–$\nu$, $\kappa_{1}$–$\kappa_{2}$,
$M_{1}$–$M_{2}$, and $\zeta_{1}$–$\zeta_{2}$ pairs. The correlation between
$\kappa_{0}$ and $\gamma$ is caused by the small errors in the HIMF and H
i-halo mass relation. The distributions between $\zeta_{0}$ and $\nu$ can be
well described by a tight linear relation of $\zeta_{0}=9.40-10.63\nu$, which
is mainly constrained by small errors in the stacked H2-stellar mass relations
of Saintonge et al. (2017). Similarly, the strong correlations in
$\kappa_{1}$–$\kappa_{2}$, $M_{1}$–$M_{2}$, and $\zeta_{1}$–$\zeta_{2}$ pairs
are caused by the stacked H i measurements of Chowdhury et al. (2022a) at
$z\sim 1.1$. Future surveys of more accurate H i and H2 density measurements
will provide tighter constraints on the model parameters, as well as on the
functional form of the redshift evolution.

## References

  * Accurso et al. (2017) Accurso, G., Saintonge, A., Catinella, B., et al. 2017, MNRAS, 470, 4750, doi: 10.1093/mnras/stx1556
  * Adams et al. (2022) Adams, E. A. K., Adebahr, B., de Blok, W. J. G., et al. 2022, A&A, 667, A38, doi: 10.1051/0004-6361/202244007
  * Amiri et al. (2023) Amiri, M., Bandura, K., Chen, T., et al. 2023, ApJ, 947, 16, doi: 10.3847/1538-4357/acb13f
  * Andreani et al. (2020) Andreani, P., Miyamoto, Y., Kaneko, H., et al. 2020, A&A, 643, L11, doi: 10.1051/0004-6361/202038675
  * Barnes et al. (2001) Barnes, D. G., Staveley-Smith, L., de Blok, W. J. G., et al. 2001, MNRAS, 322, 486, doi: 10.1046/j.1365-8711.2001.04102.x
  * Baugh et al. (2019) Baugh, C. M., Gonzalez-Perez, V., Lagos, C. d. P., et al. 2019, MNRAS, 483, 4922, doi: 10.1093/mnras/sty3427
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013a, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013b, ApJ, 762, 109, doi: 10.1088/0004-637X/762/2/109
  * Behroozi et al. (2013c) Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013c, ApJ, 763, 18, doi: 10.1088/0004-637X/763/1/18
  * Bera et al. (2019) Bera, A., Kanekar, N., Chengalur, J. N., & Bagla, J. S. 2019, ApJ, 882, L7, doi: 10.3847/2041-8213/ab3656
  * Bera et al. (2022) —. 2022, ApJ, 940, L10, doi: 10.3847/2041-8213/ac9d32
  * Bera et al. (2023) —. 2023, ApJ, 950, L18, doi: 10.3847/2041-8213/acd0b3
  * Berg et al. (2019) Berg, T. A. M., Ellison, S. L., Sánchez-Ramírez, R., et al. 2019, MNRAS, 488, 4356, doi: 10.1093/mnras/stz2012
  * Berta et al. (2016) Berta, S., Lutz, D., Genzel, R., Förster-Schreiber, N. M., & Tacconi, L. J. 2016, A&A, 587, A73, doi: 10.1051/0004-6361/201527746
  * Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846, doi: 10.1088/0004-6256/136/6/2846
  * Bird et al. (2017) Bird, S., Garnett, R., & Ho, S. 2017, MNRAS, 466, 2111, doi: 10.1093/mnras/stw3246
  * Blitz & Rosolowsky (2006) Blitz, L., & Rosolowsky, E. 2006, ApJ, 650, 933, doi: 10.1086/505417
  * Boogaard et al. (2023) Boogaard, L. A., Decarli, R., Walter, F., et al. 2023, ApJ, 945, 111, doi: 10.3847/1538-4357/acb4f0
  * Brown et al. (2017) Brown, T., Catinella, B., Cortese, L., et al. 2017, MNRAS, 466, 1275, doi: 10.1093/mnras/stw2991
  * Bryan & Norman (1998) Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80, doi: 10.1086/305262
  * Calette et al. (2021a) Calette, A. R., Avila-Reese, V., Rodríguez-Puebla, A., Lagos, C. d. P., & Catinella, B. 2021a, MNRAS, 505, 304, doi: 10.1093/mnras/stab1282
  * Calette et al. (2021b) Calette, A. R., Rodríguez-Puebla, A., Avila-Reese, V., & Lagos, C. d. P. 2021b, MNRAS, 506, 1507, doi: 10.1093/mnras/stab1788
  * Campbell et al. (2015) Campbell, D., van den Bosch, F. C., Hearin, A., et al. 2015, MNRAS, 452, 444, doi: 10.1093/mnras/stv1091
  * Carilli & Walter (2013) Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105, doi: 10.1146/annurev-astro-082812-140953
  * Catinella et al. (2010) Catinella, B., Schiminovich, D., Kauffmann, G., et al. 2010, MNRAS, 403, 683, doi: 10.1111/j.1365-2966.2009.16180.x
  * Catinella et al. (2018) Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875, doi: 10.1093/mnras/sty089
  * Chang et al. (2010) Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Nature, 466, 463, doi: 10.1038/nature09187
  * Chauhan et al. (2019) Chauhan, G., Lagos, C. d. P., Obreschkow, D., et al. 2019, MNRAS, 488, 5898, doi: 10.1093/mnras/stz2069
  * Chauhan et al. (2020) Chauhan, G., Lagos, C. d. P., Stevens, A. R. H., et al. 2020, MNRAS, 498, 44, doi: 10.1093/mnras/staa2251
  * Chauke et al. (2019) Chauke, P., van der Wel, A., Pacifici, C., et al. 2019, ApJ, 877, 48, doi: 10.3847/1538-4357/ab164d
  * Chen et al. (2019) Chen, Y., Mo, H. J., Li, C., et al. 2019, ApJ, 872, 180, doi: 10.3847/1538-4357/ab0208
  * CHIME Collaboration et al. (2022) CHIME Collaboration, Amiri, M., Bandura, K., et al. 2022, ApJS, 261, 29, doi: 10.3847/1538-4365/ac6fd9
  * Chowdhury et al. (2022a) Chowdhury, A., Kanekar, N., & Chengalur, J. N. 2022a, ApJ, 931, L34, doi: 10.3847/2041-8213/ac6de7
  * Chowdhury et al. (2022b) —. 2022b, ApJ, 941, L6, doi: 10.3847/2041-8213/ac9d8a
  * Chowdhury et al. (2020) Chowdhury, A., Kanekar, N., Chengalur, J. N., Sethi, S., & Dwarakanath, K. S. 2020, Nature, 586, 369, doi: 10.1038/s41586-020-2794-7
  * Crain et al. (2017) Crain, R. A., Bahé, Y. M., Lagos, C. d. P., et al. 2017, MNRAS, 464, 4204, doi: 10.1093/mnras/stw2586
  * Crighton et al. (2015) Crighton, N. H. M., Murphy, M. T., Prochaska, J. X., et al. 2015, MNRAS, 452, 217, doi: 10.1093/mnras/stv1182
  * Cunnington et al. (2023) Cunnington, S., Li, Y., Santos, M. G., et al. 2023, MNRAS, 518, 6262, doi: 10.1093/mnras/stac3060
  * Darvish et al. (2018) Darvish, B., Scoville, N. Z., Martin, C., et al. 2018, ApJ, 860, 111, doi: 10.3847/1538-4357/aac836
  * Davé et al. (2020) Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146, doi: 10.1093/mnras/staa1894
  * Dawson et al. (2016) Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44, doi: 10.3847/0004-6256/151/2/44
  * Decarli et al. (2016) Decarli, R., Walter, F., Aravena, M., et al. 2016, ApJ, 833, 69, doi: 10.3847/1538-4357/833/1/69
  * Decarli et al. (2019) Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, ApJ, 882, 138, doi: 10.3847/1538-4357/ab30fe
  * Delhaize et al. (2013) Delhaize, J., Meyer, M. J., Staveley-Smith, L., & Boyle, B. J. 2013, MNRAS, 433, 1398, doi: 10.1093/mnras/stt810
  * Dev et al. (2023) Dev, A., Driver, S. P., Meyer, M., et al. 2023, MNRAS, 523, 2693, doi: 10.1093/mnras/stad1575
  * Diemer et al. (2019) Diemer, B., Stevens, A. R. H., Lagos, C. d. P., et al. 2019, MNRAS, 487, 1529, doi: 10.1093/mnras/stz1323
  * Drinkwater et al. (2010) Drinkwater, M. J., Jurek, R. J., Blake, C., et al. 2010, MNRAS, 401, 1429, doi: 10.1111/j.1365-2966.2009.15754.x
  * Dutta & Khandai (2021) Dutta, S., & Khandai, N. 2021, MNRAS, 500, L37, doi: 10.1093/mnrasl/slaa178
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601, doi: 10.1111/j.1365-2966.2009.14548.x
  * Fletcher et al. (2021) Fletcher, T. J., Saintonge, A., Soares, P. S., & Pontzen, A. 2021, MNRAS, 501, 411, doi: 10.1093/mnras/staa3025
  * Fu et al. (2013) Fu, J., Kauffmann, G., Huang, M.-l., et al. 2013, MNRAS, 434, 1531, doi: 10.1093/mnras/stt1117
  * Giovanelli et al. (2005) Giovanelli, R., Haynes, M. P., Kent, B. R., et al. 2005, AJ, 130, 2598, doi: 10.1086/497431
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35, doi: 10.1088/0067-0049/197/2/35
  * Guo et al. (2020) Guo, H., Jones, M. G., Haynes, M. P., & Fu, J. 2020, ApJ, 894, 92, doi: 10.3847/1538-4357/ab886f
  * Guo et al. (2022) Guo, H., Jones, M. G., & Wang, J. 2022, ApJ, 933, L12, doi: 10.3847/2041-8213/ac794f
  * Guo et al. (2021) Guo, H., Jones, M. G., Wang, J., & Lin, L. 2021, ApJ, 918, 53, doi: 10.3847/1538-4357/ac062e
  * Guo et al. (2017) Guo, H., Li, C., Zheng, Z., et al. 2017, ApJ, 846, 61, doi: 10.3847/1538-4357/aa85e7
  * Haynes et al. (2011) Haynes, M. P., Giovanelli, R., Martin, A. M., et al. 2011, AJ, 142, 170, doi: 10.1088/0004-6256/142/5/170
  * Haynes et al. (2018) Haynes, M. P., Giovanelli, R., Kent, B. R., et al. 2018, ApJ, 861, 49, doi: 10.3847/1538-4357/aac956
  * Heintz et al. (2021) Heintz, K. E., Watson, D., Oesch, P. A., Narayanan, D., & Madden, S. C. 2021, ApJ, 922, 147, doi: 10.3847/1538-4357/ac2231
  * Heintz et al. (2022) Heintz, K. E., Oesch, P. A., Aravena, M., et al. 2022, ApJ, 934, L27, doi: 10.3847/2041-8213/ac8057
  * Janowiecki et al. (2020) Janowiecki, S., Catinella, B., Cortese, L., Saintonge, A., & Wang, J. 2020, MNRAS, 493, 1982, doi: 10.1093/mnras/staa178
  * Jones et al. (2018) Jones, M. G., Haynes, M. P., Giovanelli, R., & Moorman, C. 2018, MNRAS, 477, 2, doi: 10.1093/mnras/sty521
  * Jones et al. (2020) Jones, M. G., Hess, K. M., Adams, E. A. K., & Verdes-Montenegro, L. 2020, MNRAS, 494, 2090, doi: 10.1093/mnras/staa810
  * Kanekar et al. (2016) Kanekar, N., Sethi, S., & Dwarakanath, K. S. 2016, ApJ, 818, L28, doi: 10.3847/2041-8205/818/2/L28
  * Keres et al. (2003) Keres, D., Yun, M. S., & Young, J. S. 2003, ApJ, 582, 659, doi: 10.1086/344820
  * Kim et al. (2017) Kim, H.-S., Wyithe, J. S. B., Baugh, C. M., et al. 2017, MNRAS, 465, 111, doi: 10.1093/mnras/stw2779
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340, doi: 10.1093/mnras/stw248
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36, doi: 10.1088/0067-0049/197/2/36
  * Lagos et al. (2011) Lagos, C. D. P., Baugh, C. M., Lacey, C. G., et al. 2011, MNRAS, 418, 1649, doi: 10.1111/j.1365-2966.2011.19583.x
  * Lagos et al. (2015) Lagos, C. d. P., Crain, R. A., Schaye, J., et al. 2015, MNRAS, 452, 3815, doi: 10.1093/mnras/stv1488
  * Lah et al. (2007) Lah, P., Chengalur, J. N., Briggs, F. H., et al. 2007, MNRAS, 376, 1357, doi: 10.1111/j.1365-2966.2007.11540.x
  * Landy & Szalay (1993) Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64, doi: 10.1086/172900
  * Lenkić et al. (2020) Lenkić, L., Bolatto, A. D., Förster Schreiber, N. M., et al. 2020, AJ, 159, 190, doi: 10.3847/1538-3881/ab7458
  * Leroy et al. (2009) Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137, 4670, doi: 10.1088/0004-6256/137/6/4670
  * Li et al. (2012) Li, C., Kauffmann, G., Fu, J., et al. 2012, MNRAS, 424, 1471, doi: 10.1111/j.1365-2966.2012.21337.x
  * Li et al. (2022a) Li, X., Li, C., Mo, H. J., Xiao, T., & Wang, J. 2022a, ApJ, 941, 48, doi: 10.3847/1538-4357/ac9ccb
  * Li et al. (2022b) Li, Z., Guo, H., & Mao, Y. 2022b, arXiv e-prints, arXiv:2207.10414.  https://arxiv.org/abs/2207.10414
  * Lilly et al. (2013) Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119, doi: 10.1088/0004-637X/772/2/119
  * Lim et al. (2017) Lim, S. H., Mo, H. J., Lu, Y., Wang, H., & Yang, X. 2017, MNRAS, 470, 2982, doi: 10.1093/mnras/stx1462
  * Liu et al. (2019) Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235, doi: 10.3847/1538-4357/ab578d
  * Ma et al. (2022) Ma, W., Liu, K., Guo, H., et al. 2022, ApJ, 941, 205, doi: 10.3847/1538-4357/aca326
  * Magnelli et al. (2020) Magnelli, B., Boogaard, L., Decarli, R., et al. 2020, ApJ, 892, 66, doi: 10.3847/1538-4357/ab7897
  * Martin et al. (2010) Martin, A. M., Papastergis, E., Giovanelli, R., et al. 2010, ApJ, 723, 1359, doi: 10.1088/0004-637X/723/2/1359
  * Masui et al. (2013) Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJ, 763, L20, doi: 10.1088/2041-8205/763/1/L20
  * Meyer et al. (2004) Meyer, M. J., Zwaan, M. A., Webster, R. L., et al. 2004, MNRAS, 350, 1195, doi: 10.1111/j.1365-2966.2004.07710.x
  * Neeleman et al. (2016) Neeleman, M., Prochaska, J. X., Ribaudo, J., et al. 2016, ApJ, 818, 113, doi: 10.3847/0004-637X/818/2/113
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Noble et al. (2017) Noble, A. G., McDonald, M., Muzzin, A., et al. 2017, ApJ, 842, L21, doi: 10.3847/2041-8213/aa77f3
  * Noterdaeme et al. (2009) Noterdaeme, P., Petitjean, P., Ledoux, C., & Srianand, R. 2009, A&A, 505, 1087, doi: 10.1051/0004-6361/200912768
  * Noterdaeme et al. (2012) Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1, doi: 10.1051/0004-6361/201220259
  * Obuljen et al. (2019) Obuljen, A., Alonso, D., Villaescusa-Navarro, F., Yoon, I., & Jones, M. 2019, MNRAS, 486, 5124, doi: 10.1093/mnras/stz1118
  * Oman (2022) Oman, K. A. 2022, MNRAS, 509, 3268, doi: 10.1093/mnras/stab3164
  * O’Meara et al. (2007) O’Meara, J. M., Prochaska, J. X., Burles, S., et al. 2007, ApJ, 656, 666, doi: 10.1086/510711
  * Padmanabhan & Loeb (2020) Padmanabhan, H., & Loeb, A. 2020, MNRAS, 496, 1124, doi: 10.1093/mnras/staa1565
  * Paul et al. (2018) Paul, N., Choudhury, T. R., & Paranjape, A. 2018, MNRAS, 479, 1627, doi: 10.1093/mnras/sty1539
  * Péroux & Howk (2020) Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363, doi: 10.1146/annurev-astro-021820-120014
  * Péroux et al. (2003) Péroux, C., McMahon, R. G., Storrie-Lombardi, L. J., & Irwin, M. J. 2003, MNRAS, 346, 1103, doi: 10.1111/j.1365-2966.2003.07129.x
  * 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
  * Ponomareva et al. (2023) Ponomareva, A. A., Jarvis, M. J., Pan, H., et al. 2023, MNRAS, 522, 5308, doi: 10.1093/mnras/stad1249
  * Popping et al. (2014) Popping, G., Somerville, R. S., & Trager, S. C. 2014, MNRAS, 442, 2398, doi: 10.1093/mnras/stu991
  * Popping et al. (2015) Popping, G., Caputi, K. I., Trager, S. C., et al. 2015, MNRAS, 454, 2258, doi: 10.1093/mnras/stv2136
  * Prochaska et al. (2005) Prochaska, J. X., Herbert-Fort, S., & Wolfe, A. M. 2005, ApJ, 635, 123, doi: 10.1086/497287
  * Prochaska & Wolfe (2009) Prochaska, J. X., & Wolfe, A. M. 2009, ApJ, 696, 1543, doi: 10.1088/0004-637X/696/2/1543
  * Rhee et al. (2018) Rhee, J., Lah, P., Briggs, F. H., et al. 2018, MNRAS, 473, 1879, doi: 10.1093/mnras/stx2461
  * Rhee et al. (2013) Rhee, J., Zwaan, M. A., Briggs, F. H., et al. 2013, MNRAS, 435, 2693, doi: 10.1093/mnras/stt1481
  * Rhee et al. (2023) Rhee, J., Meyer, M., Popping, A., et al. 2023, MNRAS, 518, 4646, doi: 10.1093/mnras/stac3065
  * Riechers et al. (2019) Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7, doi: 10.3847/1538-4357/aafc27
  * Saintonge & Catinella (2022) Saintonge, A., & Catinella, B. 2022, ARA&A, 60, 319, doi: 10.1146/annurev-astro-021022-043545
  * Saintonge et al. (2016) Saintonge, A., Catinella, B., Cortese, L., et al. 2016, MNRAS, 462, 1749, doi: 10.1093/mnras/stw1715
  * Saintonge et al. (2017) Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22, doi: 10.3847/1538-4365/aa97e0
  * Scoville et al. (2017) Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150, doi: 10.3847/1538-4357/aa61a0
  * Sinigaglia et al. (2022) Sinigaglia, F., Rodighiero, G., Elson, E., et al. 2022, ApJ, 935, L13, doi: 10.3847/2041-8213/ac85ae
  * Spinelli et al. (2020) Spinelli, M., Zoldan, A., De Lucia, G., Xie, L., & Viel, M. 2020, MNRAS, 493, 5434, doi: 10.1093/mnras/staa604
  * Stevens et al. (2019) Stevens, A. R. H., Diemer, B., Lagos, C. d. P., et al. 2019, MNRAS, 483, 5334, doi: 10.1093/mnras/sty3451
  * Stiskalek et al. (2021) Stiskalek, R., Desmond, H., Holvey, T., & Jones, M. G. 2021, MNRAS, 506, 3205, doi: 10.1093/mnras/stab1845
  * Switzer et al. (2013) Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, MNRAS, 434, L46, doi: 10.1093/mnrasl/slt074
  * Tacconi et al. (2020) Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157, doi: 10.1146/annurev-astro-082812-141034
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179, doi: 10.3847/1538-4357/aaa4b4
  * Tadaki et al. (2019) Tadaki, K.-i., Kodama, T., Hayashi, M., et al. 2019, PASJ, 71, 40, doi: 10.1093/pasj/psz005
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389, doi: 10.1146/annurev-astro-091916-055240
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Genel, S., Castorina, E., et al. 2018, ApJ, 866, 135, doi: 10.3847/1538-4357/aadba0
  * Walter et al. (2008) Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563, doi: 10.1088/0004-6256/136/6/2563
  * Walter et al. (2020) Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111, doi: 10.3847/1538-4357/abb82e
  * Wang et al. (2020) Wang, J., Catinella, B., Saintonge, A., et al. 2020, ApJ, 890, 63, doi: 10.3847/1538-4357/ab68dd
  * Wang et al. (2016) Wang, J., Koribalski, B. S., Serra, P., et al. 2016, MNRAS, 460, 2143, doi: 10.1093/mnras/stw1099
  * Wang et al. (2021a) Wang, J., Staveley-Smith, L., Westmeier, T., et al. 2021a, ApJ, 915, 70, doi: 10.3847/1538-4357/abfc52
  * Wang et al. (2022) Wang, T.-M., Magnelli, B., Schinnerer, E., et al. 2022, A&A, 660, A142, doi: 10.1051/0004-6361/202142299
  * Wang et al. (2021b) Wang, Z., Chen, Y., Mao, Y., et al. 2021b, ApJ, 907, 4, doi: 10.3847/1538-4357/abcb8a
  * Wechsler et al. (2006) Wechsler, R. H., Zentner, A. R., Bullock, J. S., Kravtsov, A. V., & Allgood, B. 2006, ApJ, 652, 71, doi: 10.1086/507120
  * Wolfe et al. (2005) Wolfe, A. M., Gawiser, E., & Prochaska, J. X. 2005, ARA&A, 43, 861, doi: 10.1146/annurev.astro.42.053102.133950
  * Wolz et al. (2022) Wolz, L., Pourtsidou, A., Masui, K. W., et al. 2022, MNRAS, 510, 3495, doi: 10.1093/mnras/stab3621
  * Xi et al. (2021) Xi, H., Staveley-Smith, L., For, B.-Q., et al. 2021, MNRAS, 501, 4550, doi: 10.1093/mnras/staa3931
  * Xie et al. (2017) Xie, L., De Lucia, G., Hirschmann, M., Fontanot, F., & Zoldan, A. 2017, MNRAS, 469, 968, doi: 10.1093/mnras/stx889
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Young et al. (1995) Young, J. S., Xie, S., Tacconi, L., et al. 1995, ApJS, 98, 219, doi: 10.1086/192159
  * Zafar et al. (2013) Zafar, T., Péroux, C., Popping, A., et al. 2013, A&A, 556, A141, doi: 10.1051/0004-6361/201321154
  * Zhang et al. (2013) Zhang, W., Li, C., Kauffmann, G., & Xiao, T. 2013, MNRAS, 429, 2191, doi: 10.1093/mnras/sts490
  * Zhao et al. (2009) Zhao, D. H., Jing, Y. P., Mo, H. J., & Börner, G. 2009, ApJ, 707, 354, doi: 10.1088/0004-637X/707/1/354
  * Zwaan et al. (2005) Zwaan, M. A., Meyer, M. J., Staveley-Smith, L., & Webster, R. L. 2005, MNRAS, 359, L30, doi: 10.1111/j.1745-3933.2005.00029.x