# Conditional Hi mass functions and the Hi-to-halo mass relation in the local
Universe

Xiao Li Department of Astronomy, Tsinghua University, Beijing 100084, China
<EMAIL_ADDRESS>Cheng Li Department of Astronomy, Tsinghua
University, Beijing 100084, China<EMAIL_ADDRESS>H.J. Mo Department
of Astronomy, University of Massachusetts Amherst, MA 01003, USA  Ting Xiao
Department of Physics, Zhejiang University, Hangzhou, Zhejiang 310027, China
Jing Wang Kavli Institude for Astronomy and Astrophysics, Peking University,
Beijing 100871, China

###### Abstract

We present a new Hi mass estimator which relates the logarithm of the Hi-to-
stellar mass ratio to a linear combination of four galaxy properties: stellar
surface mass density, color index $u-r$, stellar mass and concentration index,
with the scatter of individual galaxies around the mean Hi mass modeled with a
Gaussian distribution function. We calibrate the estimator using the xGASS
sample, including both Hi detection and non-detection, and constrain the model
parameters through Bayesian inferences. Tests with mock catalogs demonstrate
that our estimator provides unbiased Hi masses for optical samples like the
SDSS, thus suitable for statistical studies of Hi gas contents in galaxies and
dark matter halos. We apply our estimator to the SDSS spectroscopic sample to
estimate the Hi mass function (HIMF) of local galaxies, as well as the
conditional Hi mass function (CHIMF) in galaxy groups and the Hi-halo mass
(HIHM) relation. Our HIMF agrees with the ALFALFA measurements at
$\textrm{M}_{\textsc{Hi}}$$\gtrsim 5\times 10^{9}$$\textrm{M}_{\odot}$, but
with higher amplitude and a steeper slope at lower masses. We show that this
discrepancy is caused primarily by the cosmic variance which is corrected for
the SDSS sample but not for the ALFALFA. The CHIMFs for all halo masses can be
described by a single Schechter function, and this is true for red, blue and
satellite galaxies. For central galaxies the CHIMFs show a double-Gaussian
profile, with the two components contributed by the red and blue galaxies,
respectively. The total Hi mass in a group increases monotonically with halo
mass. The Hi mass of central galaxies in galaxy groups increases rapidly with
halo mass only at $\textrm{M}_{\textrm{h}}$$\lesssim
10^{12}$$\textrm{M}_{\odot}$, while the mass dependence becomes much weaker at
higher halo masses. The observed Hi-halo mass relation is not reproduced by
current hydrodynamic simulations and semi-analytic models of galaxy formation.

Neutral hydrogen clouds (690), Galaxy dark matter halos (1880)

##  1 Introduction

Galaxies are believed to form at the center of dark matter halos through gas
cooling and condensation (White & Rees, 1978; Mo et al., 2010). Theoretical
studies of the Hi-to-halo mass (HIHM) relation based on hydrodynamic
simulations or halo-based models in recent years (Barnes & Haehnelt, 2014; Guo
et al., 2017; Villaescusa-Navarro et al., 2018; Paul et al., 2018; Obuljen et
al., 2019; Davé et al., 2020; Chauhan et al., 2020, 2021; Li et al., 2022)
have suggested that the majority of the cold gas in the Universe, mostly in
neutral (atomic) hydrogen (Hi) and molecular hydrogen (${\rm H}_{2}$), is
expected to be in individual galaxies, with a small amount in the circum-
galactic medium (CGM) and intergalactic medium (IGM). On the other hand, both
theoretical studies and observations have shown that some of the cold gas may
be thrown out of galaxies due to galactic winds driven by stellar and AGN
feedback or environmental effects such as tidal stripping and ram-pressure
stripping in massive halos. Clearly, measurements of Hi gas mass for large
samples of galaxies down to small gas fractions and over large sky coverage
and redshift ranges are needed to provide stringent constraints on models of
galaxy formation.

Large surveys of Hi 21cm emission of galaxies have become available only in
the past two decades, such as the Hi Parkes All-Sky Survey (HIPASS; Meyer et
al., 2004; Wong et al., 2006) and the Arecibo Legacy Fast ALFA (ALFALFA;
Giovanelli et al., 2005). These surveys have detected Hi emission from tens of
thousands of galaxies up to $z\sim 0.06$, covering large areas in the sky. As
expected, most (if not all) cases of the Hi detection are associated with
optically-identified galaxies in the same sky position and redshift (Haynes et
al., 2018). Both HIPASS and ALFALFA have yielded consistent estimates of the
Hi mass function (HIMF) of the local galaxy population, across a mass range
over nearly five orders of magnitudes from $\textrm{M}_{\textsc{Hi}}$$\sim
10^{6}M_{\odot}$ to $\textrm{M}_{\textsc{Hi}}$$\sim 10^{11}M_{\odot}$ (Zwaan
et al., 2005; Jones et al., 2018). Compared to optical surveys, however,
current Hi surveys are still shallow and biased to gas-rich galaxies. There
have been recent efforts to measure the total Hi mass in dark matter halos,
e.g. by stacking ALFALFA data-cubes of galaxy groups of different halo mass
(Guo et al., 2020) or through deep observations of Hi emission of individual
galaxies (Rhee et al. in preparation). Guo et al. (2020) found a positive
correlation between the total Hi mass and the dark matter halo mass, with a
significant excess of Hi mass at $M_{h}\sim 10^{12}\textrm{M}_{\odot}$, while
Rhee et al. (in preparation) observed a nearly flat relation. Chauhan et al.
(2021) analyzed possible reasons for the discrepancy between these two
studies, and showed that these measurements may not be able to obtain the
intrinsic Hi-halo mass relation correctly, owing to uncertainties in the halo
mass estimate and limitations of the stacking method.

Hi samples with depths and sizes comparable to optical surveys, such as the
Sloan Digital Sky Survey (SDSS; York et al., 2000), are difficult to obtain in
general. As an alternative way to proceed, attempts have been made to estimate
Hi masses for optically-selected galaxies using their spectral and/or
photometric properties. For instance, Tremonti et al. (2004) estimated surface
gas mass densities for SDSS galaxies from surface densities of star formation
rate (given by H$\alpha$ luminosities) by inverting the Kennicutt-Schmidt law
(Schmidt, 1963; Kennicutt, 1998). Meanwhile, Kannappan (2004) found that the
Hi-to-stellar mass ratio ($\log(\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast})$)
of local galaxies is linearly correlated with optical (e.g. $u-r$) and
optical-to-near IR (e.g. $u-K$) colors, with a scatter of $\sim 0.4$ dex.
Motivated by these studies, Zhang et al. (2009) proposed an estimator for
$\log(\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast})$ based on a linear
combination of the $g-r$ color and $i$-band surface brightness $\mu_{i}$, with
a scatter of $\sim 0.3$ dex. They applied the estimator to the SDSS galaxy
sample to investigate the variation of Hi gas mass fraction on the stellar
mass-metallicity relation. A follow-up study by the team of the GALEX Arecibo
SDSS Survey (GASS; Catinella et al., 2010) calibrated a similar estimator to
link $\log(\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast})$ with $NUV-r$ and
surface stellar mass density $\mu_{\ast}$, which also has a scatter of $\sim
0.3$ dex. Li et al. (2012, hereafter L12) found that the two-parameter
estimators from Zhang et al. (2009) and Catinella et al. (2010) significantly
underestimate the Hi mass for the ALFALFA sample, which is expected to be
biased to Hi-rich galaxies. They proposed a new estimator with four galaxy
properties: surface stellar mass density $\mu_{\ast}$, $NUV-r$, stellar mass
$M_{\ast}$ and $g-i$ color gradient $\Delta_{g-i}$. With a smaller scatter of
$\sim 0.2$ dex, this estimator was applied to the SDSS sample to study both
the Hi mass-dependence of the galaxy clustering (L12) and environmental
effects on the gas depletion in clusters of galaxies (Zhang et al., 2013).

Using the REsolved Spectroscopy of a Local VolumE survey (RESOLVE; Kannappan
et al. in prep.), Eckert et al. (2015) calibrated a relation between the Hi
gas fraction, color index and axis ratio of galaxies. More recently, Zu (2020,
hereafter Z20) developed a two-parameter Hi mass estimator by linearly
combining stellar mass $M_{\ast}$ and $(g-r)$, and used it to study the mass-
metallicity relation of SDSS galaxies. To avoid the Malmquist bias produced by
limited detection depths of current Hi surveys, Z20 designed a likelihood
model to account for the detection probability of ALFALFA and constrained the
model parameters using Bayesian inferences. Lu et al. (2020) also used a
linear combination of $M_{\ast}$ and $(g-r)$ to estimate the Hi mass of SDSS
galaxies, in an attempt to estimate the Hi gas content in dark matter halos.
In addition to these linear estimators, a few nonlinear models have also been
obtained to estimate Hi masses of local galaxies, taking advantage of the
technique of machine learning (e.g. Teimoorinia et al., 2017; Rafieferantsoa
et al., 2018).

In this paper we extend previous studies by developing a new estimator to
predict the Hi gas content of galaxies in a large optical sample. Our
estimator makes improvements in the following aspects. First, we calibrate our
Hi mass estimator using the xGASS sample (Catinella et al., 2018), and we
include galaxies both with and without Hi detection in the calibration to
reduce bias. Recent studies have mostly used the ALFALFA sample, which is
shallow and biased to Hi-rich galaxies. Galaxies without HI detection contain
useful information about gas-poor galaxies, and so should not be ignored in
the calibration in order to obtain an unbiased estimator. Second, we follow
L12 to use a linear model using four galaxy properties for the estimator, but
we also take into account scatter in the model. This scatter includes not only
the intrinsic variance of individual galaxies but also uncertainties in the
measurements of the Hi mass and optical parameters. Finally, we follow Z20 to
constrain model parameters of our estimator using the Bayesian framework,
which allows us to explore the model parameter space efficiently, and to
better understand the correlations among model parameters. We use mock
catalogs to mimic the selection effects of the ALFALFA and xGASS samples and
to demonstrate that our estimator is unbiased for SDSS-like samples. We thus
can use our estimator to estimate the HIMF of local galaxies, the conditional
HIMF and the total Hi mass as a function of dark matter halo mass.

We organize our paper as follows. In section 2 we describe the Hi and optical
samples used in our analysis. In section 3 we present our Hi estimator and
describe the calibration procedure and test results. In section 4 we apply our
Hi estimator to the SDSS sample to estimate Hi contents of galaxies and their
host halos. We discuss and summarize our results in section 5 and section 6.
Throughout this paper we assume a flat $\Lambda$CDM cosmology with
$H_{0}=100h{\rm\ km\ s^{-1}Mpc^{-1}}$, $h=0.7$, $\Omega_{m}=0.3$, and
$\Omega_{\Lambda}=0.7$.

##  2 Data

###  2.1 The xGASS sample

We use the xGASS sample to calibrate our Hi mass estimator. The xGASS is an
extension of the GASS survey, an Hi 21cm survey observed with the Arecibo
telescope for a sample of galaxies with redshift $0.025<z<0.05$ and a flat
stellar mass distribution in the range
$10^{10}\textrm{M}_{\odot}<\textrm{M}_{\ast}<10^{11.5}\textrm{M}_{\odot}$. The
sample galaxies are randomly selected from a parent sample of $\sim$1,2000
galaxies located in the overlapping region among SDSS data release six (DR6;
Adelman-McCarthy et al., 2008), Galaxy Evolution Explorer (GALEX; Martin et
al., 2005) Medium Imaging Survey, and the ALFALFA survey footprint. Each
galaxy is observed with Arecibo until a significant Hi emission line is
detected, or the Hi to stellar mass ratio reaches an upper limit of
$\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast}\sim 1.5\%$ (see C10 for details).
The xGASS extends the survey down to a lower mass limit of
$10^{9}\textrm{M}_{\odot}$, by further observing a sample of galaxies with
$9.0<\log($$\textrm{M}_{\textsc{Hi}}$$/$$\textrm{M}_{\odot}$$)<10.2$ and111For
simplicity, we’ll denote $\log_{10}$ as $\log$ throughout this paper
$0.01<z<0.02$. Here we use the combined xGASS sample constructed by Catinella
et al. (2018), which includes 1179 galaxies obtained by GASS and xGASS and
supplemented with Hi-rich galaxies selected from the ALFALFA $\alpha.70$
sample that are not included in GASS/xGASS.

From the combined xGASS sample we select a subset of galaxies with good
measurements. A galaxy is included in our sample if all the following
requirements are met: (i) $\tt{HIconf\\_flag}\leq 0$, (ii) $\tt{HI\\_flag}\leq
1$, and (iii) $\sigma_{u-r}<0.3$. Here $\tt{HIconf\\_flag}$ is the confusion
flag provided in the combined xGASS catalog, with $\tt{HIconf\\_flag}=1$ for
certain confusion, $>0$ for a small companion, $0$ for no confusion and $-99$
for non-detection. Thus, the first criterion selects both detections with no
confusion and non-detections from the xGASS sample. The $\tt{HI\\_flag}$ is
the quality flag for xGASS detections, with $\tt{HI\\_flag}\leq 1$ for good
quality, $2$ for marginal, $5$ for confused and $3$ for both marginal and
confused observations. The error of $u-r$ color, $\sigma_{u-r}$, is calculated
through error propagation based on the inverse variance of the $u-$ and
$r-$band absolute magnitudes provided in the NASA Sloan Atlas (NSA; Blanton et
al., 2011). These restrictions yield a final sample, which contains 625
galaxies each with an Hi mass measurement, and 358 galaxies each with only an
Hi mass upper limit. For those with an Hi gas mass measurement,
$\textrm{M}_{\textsc{Hi}}$ is self-absorption corrected using $\Delta\log
M_{\rm{\textsc{Hi}}}=(0.13\pm 0.03)\log(a/b)$ (Jones et al., 2018), where
$a/b$ is the major-to-minor axis ratio. We refer to this sample as SampleX in
the rest of this paper.

###  2.2 The ALFALFA galaxy sample

The Arecibo Legacy Fast Arecibo L-band Feed Array Survey (ALFALFA; Giovanelli
et al., 2005) is a blind survey of Hi 21 cm emission over $\sim 7000\ {\rm
deg^{2}}$ of the sky and up to redshift $z\sim 0.06$. The final data release,
the ALFALFA $\alpha.100$ catalog (Haynes et al., 2018), consists of $\sim
31,500$ extragalactic Hi line sources. Two sky areas are covered by ALFALFA:
one in the northern Galactic hemisphere and one in the southern Galactic
hemisphere. Since our optical sample (see below) is mainly in the northern
Galactic hemisphere, we select from ALFALFA $\alpha.100$ a subsample that
covers the sky of $138^{\circ}<\alpha<232^{\circ},\
0^{\circ}<\delta<36.5^{\circ}$ and the redshift of $0.00<z<0.05$. The upper
limit in $z$ is set to avoid the contamination of RFI (Haynes et al., 2018).
We use ALFALFA sources flagged with ${\tt Code\ 1}$, i.e those with ${\rm
SNR\gtrsim 6.5}$, together with those flagged with ${\tt Code\ 2}$, i.e. those
with ${\rm SNR\ (\lesssim 6.5)}$ but matched with optical counterparts that
have redshift consistent with the observed 21cm line. Our ALFALFA galaxy
sample consists of 16,400 galaxies. In what follows this sample is referred to
as SampleA.

###  2.3 The observed HI mass function

Figure 1: The HIMF of ALFALFA and xGASS survey. The magenta errorbars are the
ALFALFA 100 percent HIMF datapoints (Jones et al., 2018). The magenta dashed
line is the best fit to the data. The black open circle with errorbar is the
xGASS HIMF derived in this work (see subsection 2.3).

A blind Hi survey, such as the ALFALFA, aims to detect all sources above some
detection limit, e.g. in 21cm flux. Such a survey is expected to be biased
against galaxies with low Hi mass. Thus, unlike an Hi survey which starts from
a given optical sample, such as xGASS described above, a blind survey does not
provide a fair sample to study the relationship between Hi and optical
properties of galaxies. However, being selected on the basis of 21cm flux,
blind surveys can provide complete samples to obtain summary statistics, such
as the Hi mass function (HIMF), above certain Hi mass. Such statistics are not
direct measurements of the relationship between Hi and optical properties,
they nevertheless provide constraints on the relationship. In our analysis, we
will use the HIMF obtained from ALFALFA as an additional constraint on our
model.

As a check of consistency between xGASS and ALFALFA, we show in Figure 1 the
HIMFs obtained from these two samples. The ALFALFA measurements are adopted
from Jones et al. (2018) and shown as the magenta error bars, together with
the fitting to a Schechter function shown by the magenta dashed line. The HIMF
for xGASS is estimated by us using the xGASS representative sample (SampleX).
To that end, we assign a weight $w_{i1}w_{i2}$ to the $i$-th galaxy in
SampleX, where $w_{i1}$ accounts for the mass-dependent selection of SampleX
with respect to the SDSS galaxy sample (SampleS), and $w_{i2}$ corrects for
the selection effect of SampleS using the $1/V_{\rm max}$ weighting scheme
(see subsection 4.1 for more details). We only obtain the xGASS HIMF at $\log
M_{\rm{\textsc{Hi}}}/M_{\odot}>9.5$ to avoid influences of non-detection and
incompleteness at the low-mass end. The errors of the HIMF are estimated by
bootstrap resampling of the galaxies in SampleX. As one can see from the
figure, at $\log M_{\rm{\textsc{Hi}}}/M_{\odot}>9.5$, both the xGASS and
ALFALFA measurements are consistent with each other. To ensure consistency, we
only use the ALFALFA measurements of the HIMF at $\log
M_{\rm{\textsc{Hi}}}/M_{\odot}>9.5$ as additional constraints on our Hi
estimator.

###  2.4 The SDSS galaxy sample

We select our optical galaxy sample from the New York University Value Added
Galaxy Catalog (NYU-VAGC)222http://sdss.physics.nyu.edu/vagc/ constructed by
Blanton et al. (2005) from the SDSS spectroscopic galaxy sample. We start with
the VAGC post-redshift sample bbright0, which contains 535,192 galaxies with
spectroscopically measured redshifts and brighter than $r=17.6$ mag, where $r$
is the $r$-band apparent Petrosian magnitude corrected for Galactic
extinction. We restrict ourselves to the redshift range of $0.003<z<0.05$,
where the upper redshift limit is set to match the Hi galaxy sample and the
lower limit is set to avoid large uncertainties in distances produced by
peculiar motions of nearby galaxies. This gives us a sample of 86,487
galaxies. We refer to this sample as SampleS in the rest of this paper. We
will apply our Hi mass estimator to this sample and derive the statistical
properties of Hi gas contents for both galaxies and their host dark matter
halos (see section 4).

###  2.5 The NSA and galaxy properties

To design our Hi mass estimator, we obtain, for each galaxy in the xGASS,
ALFALFA and SDSS samples described above, the following set of properties from
the NSA333http://www.nsatlas.org/:

  * •

$\textrm{M}_{\ast}$: stellar mass in units of solar mass, estimated with ${\tt
kcorrect\ v4\\_2}$,444http://kcorrect.org (Blanton & Roweis, 2007) which
performs SED fitting to the SDSS photometric data based on the stellar
spectral templates from Bruzual & Charlot (2003) using the stellar initial
mass function of Chabrier (2003) and the Padova 1994 isochrones, and the
ionized gas emission spectral templates of Kewley et al. (2001).

  * •

$\mu_{\ast}$: the surface stellar mass density, defined as $M_{*}/(2\pi
R_{50}^{2})$, where $R_{50}$, in units of kpc, is the radius enclosing a half
of the total light in $r$-band.

  * •

$u-r$: the $u-r$ color index, given by the K-corrected absolute magnitude
difference between $u$\- and $r$-band using elliptical Petrosian photometry of
SDSS images.

  * •

$NUV-r$: the $NUV-r$ color index, where $NUV$ is the K-corrected magnitude in
the $NUV$-band using elliptical Petrosian photometry of GALEX images.

  * •

$R_{90}/R_{50}$: the concentration index defined as the ratio between $R_{90}$
(the radius enclosing 90% of the total light in $r$-band) and $R_{50}$.

About 5% of the galaxies in SampleS do not have counterparts in the NSA. We
have examined the effects of these missing galaxies on our results and found
that they are negligible. Figure 2 displays the distribution of the SampleX,
SampleA and SampleS in the plane of $NUV-r$ color and stellar mass. Galaxies
with Hi detection and non-detection in SampleX are plotted as blue dots and
red crosses respectively, respectively. The cyan contours represent the
distribution of SampleA. The grey 2D histogram and dots indicate SampleS. We
can see that SampleA galaxies are mainly located in the star-forming sequence
with $NUV-r\lesssim 4$. Compared to SampleA, SampleX covers a larger region in
parameter space (although it is limited to
$\textrm{M}_{\ast}$$>10^{9}$$\textrm{M}_{\odot}$ and with a smaller sample
size), thus can better represent the general galaxy population. Therefore, we
use SampleX to calibrate our Hi mass estimator.

Figure 2: The distribution of our galaxy samples in NUV-r color and stellar
mass diagram. Blue points are Hi-detected galaxies in SampleX and red crosses
indicate Hi non-detections in SampleX. Cyan contours indicate the distribution
of SampleA galaxies. From inside to outside, the contours include 16%, 50% and
84% of the SampleA respectively. The grey 2D histogram and dots show the
distribution of SampleS galaxies.

##  3 Predicting the Hi mass of galaxies

###  3.1 The Hi mass estimator

Consider a sample of $N_{g}$ galaxies indexed by $i$, each of which has a set
of observational parameters, $\bm{x}_{i}$, in addition to the Hi observation,
$y_{i}$. We aim to calibrate a general relation between $y$ and $\bm{x}$ so
that it can be used to reliably predict $y$ from $\bm{x}$ for galaxies in
other samples. We model this relation as

$y(\bm{x})=y_{0}(\bm{x})+\delta y(\sigma),$ (1)

where $y_{0}$ gives the mean value of $y$ for galaxies of given $\bm{x}$, and
$\delta y$ quantifies the scatter of individual galaxies around the mean.
$\delta y$ is assumed to be a random variable that follows a Gaussian
distribution function with a zero mean and a dispersion of $\sigma$. We model
the mean of $y$ as a linear combination of $\bm{x}$:

$y_{0}(\bm{x})=\bm{\beta}^{T}\bm{x}+e,$ (2)

where $\bm{\beta}$ is a vector of model parameters that has the same number of
elements as $\bm{x}$, and $e$ is a constant. In this work we choose $y$ to be
the Hi-to-stellar mass ratio in logarithmic scale,
$y=\log(\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast})$, and We use SampleX to
calibrate model parameters. For $\bm{x}$, we adopt the following four
parameters: the surface stellar mass density $\mu_{\ast}$, the color index
$u-r$, the stellar mass $\textrm{M}_{\ast}$ and the concentration index
$R_{90}/R_{50}$ (see § 2.5). Thus

$y=\log\frac{\textrm{M}_{\textsc{Hi}}}{\textrm{M}_{\ast}},~{}~{}~{}\bm{x}=\left(\begin{aligned}
&\log\mu_{\ast}\\\ &u-r\\\ &\log\textrm{M}_{\ast}\\\
&\log\frac{R_{90}}{R_{50}}\end{aligned}\right),~{}~{}~{}\bm{\beta}=\left(\begin{aligned}
&a\\\ &b\\\ &c\\\ &d\end{aligned}\right).$ (3)

It is not straightforward to design an appropriate functional form for
$\sigma$. Our tests showed that the scatter actually depends on
$M_{\textsc{Hi},0}$, which is the mean value of Hi mass given by the Hi
estimator ($M_{\textsc{Hi},0}=y_{0}+\log\textrm{M}_{\ast}$), and the
dependence can be described by a piecewise linear function:

$\sigma(M_{\textsc{Hi},0})=\left\\{\begin{aligned} &|c_{a}\log
M_{\textsc{Hi},0}+c_{b}|,&\textrm{if}~{}M_{\textsc{Hi},0}\geq
M_{\textsc{Hi},0,t}\\\ &|c_{a}\log
M_{\textsc{Hi},0,t}+c_{b}|,&\textrm{if}~{}M_{\textsc{Hi},0}<M_{\textsc{Hi},0,t}\end{aligned}\right.$
(4)

where $c_{a}$ and $c_{b}$ are model parameters to be determined, and
$M_{\textsc{Hi},0,t}$ is an empirically determined cut in $M_{\textsc{Hi},0}$
below which the scatter is a constant. Our model thus has a total of 7 free
parameters: $a,\ b,\ c,\ d,\ e,\ c_{a},\ c_{b}$, to be determined. The four-
parameter model for the mean relation (Equation 2) is motivated by that in L12
where the Hi-to-stellar mass ratio of galaxies was considered as a function of
$\mu_{\ast}$, $NUV-r$, $\textrm{M}_{\ast}$ and the radial gradient of $g-i$.
Here we replace the $g-i$ color gradient by the concentration index which is
more widely used and easier to measure. We also replace $NUV-r$ by $u-r$ for
similar reasons. We note that our model and that of L12 lead to similar
results in terms of estimating the mean relation between Hi gas mass and
galaxy properties. We have attempted to include more parameters or to replace
some of the parameters, and found no significant improvements. As pointed out
in the introduction, our Hi mass estimator is an improvement on earlier
estimators because of the inclusion of both the scatter in the model and non-
detections in the calibration sample.

###  3.2 Bayesian inferences of model parameters

We make Bayesian inferences for the model parameters of our Hi mass estimator.
In what follows we denote the set of model parameters as
$\bm{\theta}\equiv\\{a,b,c,d,e,c_{a},c_{b}\\}$. We use a combination of two
sets of observational data. The first set includes measurements of the four
galaxy properties and the Hi-to-stellar mass ratio for all individual galaxies
in SampleX, and is denoted as $\bm{D}\equiv\\{\bm{x}_{i},y_{i}\\}$
($i=1,...,N_{g}$), where $N_{g}$ is the number of galaxies in the sample. The
second set is the HIMF estimated from ALFALFA, and is denoted as
$\Phi_{\textsc{Hi}}\equiv\Phi(\textrm{M}_{\textsc{Hi}})=\\{\Phi_{\textsc{Hi},j}\\}$
($j=1,...,N_{m}$), with $N_{m}$ the number of mass bins of the HIMF. Applying
Bayes Theorem, we can write the likelihood of $\bm{\theta}$ given both
$\bm{D}$ and $\Phi_{\textsc{Hi}}$ as

$\displaystyle P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})$
$\displaystyle=\frac{P(\Phi_{\textsc{Hi}},\bm{D}|\bm{\theta})\cdot
P(\bm{\theta})}{P(\Phi_{\textsc{Hi}},\bm{D})}$ (5)
$\displaystyle=\frac{P(\Phi_{\textsc{Hi}}|\bm{\theta})\cdot
P(\bm{D}|\bm{\theta})\cdot P(\bm{\theta})}{P(\Phi_{\textsc{Hi}})\cdot
P(\bm{D})}$ $\displaystyle=\frac{P(\Phi_{\textsc{Hi}}|\bm{\theta})\cdot
P(\bm{\theta}|\bm{D})}{P(\Phi_{\textsc{Hi}})},$

where in the second line we have assumed that the measurements of $D$ and
$\Phi_{\textsc{Hi}}$ are independent. The third line follows from the Bayes
relation,

$P(\bm{\theta}|\bm{D})=\frac{P(\bm{D}|\bm{\theta})\cdot
P(\bm{\theta})}{P(\bm{D})}$ (6)

where $P(\bm{\theta}|\bm{D})$ is the posterior distribution of $\bm{\theta}$
inferred from $\bm{D}$ alone. The above equations show that we can derive
$P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})$ in two steps. First, we obtain
$P(\bm{\theta}|\bm{D})$ from $\bm{D}$ using an assumed prior distribution
$P(\bm{\theta})$. Second, we infer $P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})$
from $\Phi_{\textsc{Hi}}$ using $P(\bm{\theta}|\bm{D})$ as the prior.

We start by deriving $P(\bm{\theta}|\bm{D})$ using Equation 6. Assuming that
galaxies are sampled with an inhomogeneous Poisson process, we can write
Equation 6 in logarithmic form as

$\displaystyle\ln P(\bm{\theta}|\bm{D})=$ $\displaystyle\ln
P(\bm{D}|\bm{\theta})+\ln P(\theta)-\ln P(\bm{D})$ (7) $\displaystyle=$
$\displaystyle\sum_{i\in\textrm{S}_{\textrm{d}}}\ln
P_{\textrm{d},i}(y_{i}|\bm{\theta},\bm{x}_{i})+\sum_{i\in\textrm{S}_{\textrm{n}}}\ln
P_{\textrm{n},i}(y_{i}|\bm{\theta},\bm{x}_{i})$
$\displaystyle+\textrm{const.}$

Here we have assumed a flat prior distribution for all model parameters, so
that $\ln P(\bm{\theta})$ and the data distribution $P(\bm{D})$ both become
constant and can be ignored. In the second equation $\textrm{S}_{\textrm{d}}$
and $\textrm{S}_{\textrm{n}}$ denote the subsets of galaxies in SampleX with
and without Hi detection, respectively, while $P_{\textrm{d},i}$ and
$P_{\textrm{n},i}$ are the corresponding likelihoods of $y_{i}$ given
$\bm{\theta}$ and $\bm{x}_{i}$. The likelihood for the $i$-th galaxy to have
an Hi detection is modeled by

$P_{\textrm{d},i}(y_{i}|\bm{\theta},\bm{x}_{i})=\frac{1}{\sqrt{2\pi}\sigma_{i}}\exp\left\\{-\frac{(y_{i}-y_{0,i})^{2}}{2\sigma_{i}^{2}}\right\\},$
(8)

where $y_{0,i}=y_{0,i}(\bm{x}_{i})$ is the mean value predicted by Equation 2,
and $\sigma_{i}=\sigma_{i}(M_{\textsc{Hi},0,i})$ is the scatter given by
Equation 4. When calculating $\sigma_{i}$, we set the value of
$M_{\textsc{Hi},0,t}$ in Equation 4 to be
$M_{\textsc{Hi},0,t}=10^{8.5}\textrm{M}_{\odot}$ based on tests which show
that our results are insensitive to the exact value of $M_{\textsc{Hi},0,t}$.
For a non-detection with only an upper limit of $y_{i}$, the likelihood
$P_{\textrm{n},i}$ can be calculated by integrating $P_{\textrm{d},i}$ over
$y\leq y_{i}$, i.e.

$P_{\textrm{n},i}(y_{i}|\bm{\theta},\bm{x}_{i})=\int^{y_{i}}_{-\infty}P_{\textrm{d},i}(y^{\prime}|\bm{x}_{i},\bm{\theta})dy^{\prime}\,.$
(9)

We use the Python package $\tt{emcee}$ to perform Markov chain Monte Carlo
(MCMC) sampling over the parameter space (Foreman-Mackey et al., 2013), and
obtain the posterior distribution of the model parameters using Equation 7.
During this process we weight each galaxy in SampleX in such a way that the
stellar mass distribution estimated by the weighted galaxies matches the
stellar mass function (GSMF) in the local universe as measured by Chen et al.
(2019) (see their Table 1). The posterior distributions of all the model
parameters obtained this way are plotted in Figure 3 as blue contours.

Next, we use the ALFALFA HIMF $\Phi(\textrm{M}_{\textsc{Hi}})$ at $\log
M_{\rm{\textsc{Hi}}}/M_{\odot}>9.5$ to further constrain our model parameters,
adopting the likelihood derived above, $P(\bm{\theta}|\bm{D})$, as the prior.
Again, we write the last equation in Equation 5 in logarithmic form as

$\displaystyle\ln P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})=$ $\displaystyle\ln
P(\Phi_{\textsc{Hi}}|\bm{\theta})+\ln P(\bm{\theta}|\bm{D})+\textrm{const.}$
(10) $\displaystyle=$ $\displaystyle\sum_{j}\ln
P_{j}(\Phi_{\textsc{Hi},j}|\bm{\theta})+\ln
P(\bm{\theta}|\bm{D})+\textrm{const.},$

where $\Phi_{\textsc{Hi},j}$ is the HIMF measurement in the $j$-th mass bin,
and $P_{j}(\Phi_{\textsc{Hi},j}|\bm{\theta})$ is the likelihood of
$\Phi_{\textsc{Hi},j}$ given $\bm{\theta}$:

$P_{j}(\Phi_{\textsc{Hi},j}|\bm{\theta})=\frac{1}{\sqrt[]{2\pi}\Delta_{j}}\exp\left\\{-\frac{\left[\Phi_{\textsc{Hi},j}-\Phi^{\prime}_{\textsc{Hi},j}(\bm{\theta})\right]^{2}}{2\Delta_{j}^{2}}\right\\}$
(11)

Here $\Phi^{\prime}_{\textsc{Hi},j}(\bm{\theta})$ is the predicted HIMF in the
$j$-th mass bin derived by applying our HI estimator to SampleS using the $\rm
1/V_{max}$ weighting scheme (see subsection 4.1 for more details), and
$\Delta_{j}$ is the error of $\Phi_{\textsc{Hi},j}$.

Figure 3: Probability distribution of model parameters. Contours correspond to
1-$\sigma$ and 2-$\sigma$ from inside out. The blue contours show the
distribution of $P(\bm{\theta}|\bm{D})$, and the red contours show the
distribution of $P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})$.

Figure 4: Normalized distribution of $u-r$, $\log M_{*}$ and $\log
M_{\rm{\textsc{Hi}}}$ of galaxies. The cyan histogram represents the results
of real samples and the magenta histogram of mock samples. Both the grey solid
histogram in the top panel and the grey dashed histogram in the bottom panel
represent the SDSS volume-limited sample used to construct ALFALFA mock
catalog. The grey solid histogram in the bottom panel represents the mean
distribution of the 25 mock xGASS samples. The magenta shaded area shows the
$1-\sigma$ scatter of the 25 mocks.

Figure 5: Hi mass fraction as a function of galaxy properties for
ALFALFA(top)/xGASS(bottom) mock and real survey. From left to right, the
horizontal axis represents galaxy surface mass density, $u-r$ color index,
stellar mass and galaxy concentration respectively. The y-axis is the observed
Hi mass fraction for Hi-detected galaxies in real surveys and the estimated Hi
mass fraction for galaxies in mock surveys. The grey contours represent the
volume-limited sample used in the mock ALFALFA survey. The contours include
5%, 15%, 25%, …, 95% of the total sample. Galaxies outside these contours are
shown as grey dots. The boxes show the $25\%-75\%$ percentile range of
galaxies detected in mock surveys (magenta) or real surveys (cyan) in each
bin. The vertical error bars show the $5\%-95\%$ percentile range. The green
boxes show the distribution of the total xGASS sample (SampleX, detections +
non-detections)

The same MCMC sampling as used in the first step is applied here to derive
$P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})$ in Equation 10. The results are
plotted in Figure 3 as red contours. As expected, $P(\bm{\theta}|\bm{D})$ and
$P(\bm{\theta}|\Phi_{\textsc{Hi}},\bm{D})$ are consistent with each other and
the latter are better constrained thanks to the additional information
provided by the HIMF. In particular, the posterior distributions of $c$ and
$e$ are narrowed significantly. We adopt as our best model with maximum
likelihood: $a=-0.42\pm 0.20,\ b=-0.82\pm 0.19,\ c=0.15\pm 0.11,\ d=-2.57\pm
1.30,\ e=3.57\pm 0.95,\ c_{a}=0.16\pm 0.10$ and $c_{b}=-1.20\pm 0.76$. The
average Hi-to-stellar mass ratio of galaxies at fixed stellar mass is
predicted to have negative correlations with $\mu_{\ast}$, $u-r$ and
$R_{90}/R_{50}$, which is consistent with direct observations (e.g. Kannappan,
2004; Zhang et al., 2009; Catinella et al., 2010). There is some degeneracy
among model parameters. For instance, $e$ shows negative correlations with $a$
and a positive correlation with $b$. We note that $c_{a}$ and $c_{b}$ are
strongly degenerate, implying that there is still room to improve the model of
the scatter.

###  3.3 Tests with mock Hi catalogs

Our Hi mass estimator is not simply a predictor of the mean value for galaxies
of a given set of optical properties (as usually the case in previous
studies), but rather it includes the variance of individual galaxies. For an
individual galaxy with real Hi measurement, the estimated Hi mass may be
different from the real value because of the scatter. Therefore, we cannot
directly compare the estimated and real values of the Hi mass for individual
galaxies in ALFALFA and xGASS. Instead, we can test our Hi estimator only
statistically, by comparing the predicted distribution of Hi mass with that
obtained from real observations. To facilitate a fair comparison, we also need
mock catalogs constructed to include the same selection effects as xGASS or
ALFALFA. If our Hi estimator provides an unbiased Hi prediction for a sample
of general-population galaxies (rather than individual galaxies) selected from
an optical survey, the mock catalog should reproduce the Hi mass distribution
obtained from the real sample.

From the ALFALFA$\alpha$.100 sample, we select a 100% complete sample of
$6198$ Hi-detected galaxies with $z<0.05$ and
$\log\textrm{M}_{\ast}/\textrm{M}_{\odot}>9$, and in the sky area of
$138^{\circ}<\alpha<232^{\circ}$ and $0^{\circ}<\delta<36.5^{\circ}$, the
overlapping sky coverage of the SDSS and ALFALFA $\alpha.100$ footprints. We
require the lower mass limit of $\log\textrm{M}_{\ast}/\textrm{M}_{\odot}>9$
in order to keep consistency with the mass limit of the calibration sample
SampleX. In addition, galaxy samples selected with $z<0.05$ and
$\log\textrm{M}_{\ast}/\textrm{M}_{\odot}>9$ are actually volume-limited, and
so we do not need to worry about sample incompleteness when comparing our mock
catalog with the ALFALFA sample. Using the ALFALFA 100% completeness limit is
for the same consideration. The 100% completeness limit is obtained by
extrapolating the $25\%,\ 50\%$ and $90\%$ completeness limits derived by
Haynes et al. (2011) using the distribution of ALFALFA extragalactic sources
in the plane defined by the Hi integrated flux density ($S_{21}$) and the line
profile width ($W_{50}$). To construct the corresponding mock catalog, we
start with a volume-limited sample of SDSS galaxies selected from SampleS that
covers the same sky area, stellar mass and redshift ranges as the ALFALFA 100%
complete sample. This volume-limited sample contains $38,298$ galaxies. For
each galaxy we then apply our
$\log(\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast})$ estimator to obtain an Hi
mass, and calculate $S_{21}$ according to the estimated Hi mass and the
redshift of the galaxy. We estimate the line width $W_{50}$ for the galaxy
using the stellar mass Tully-Fisher relation given in Bloom et al. (2017) with
inclination considered. The galaxy is included in the mock catalog if its Hi
flux density is above the 100% detection limit corresponding to its flux
density and line width. A total of $5553$ galaxies are selected into the mock
catalog.

For xGASS, we take all the 983 galaxies including 595 Hi-detected galaxies
with $\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast}$ above a detection limit of
4%, and 388 non-detections below the limit. The detection threshold of the Hi-
to-stellar mass ratio is chosen to be slightly higher than the usually quoted
limit of $\sim 1.5\%$ to ensure that the subsample of detections is truly
complete down to the limit. We form a set of 25 mock catalogs, each
constructed by randomly selecting 983 galaxies from SampleS in the same
stellar mass and redshift ranges either as the GASS survey
($10.2<\log\textrm{M}_{\ast}/\textrm{M}_{\odot}<11.5$ and $0.025<z<0.05$) or
as the extended xGASS survey
($9.0<\log\textrm{M}_{\ast}/\textrm{M}_{\odot}<10.2$ and $0.01<z<0.02$). In
addition, we require each of the mock catalogs to follow the same stellar mass
distribution as the xGASS sample. The mock galaxies with
$\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast}>4\%$ are taken as detections and
others are non-detections. On average, the mock catalogs contain 588
detections, very close to the real sample (595).

In Figure 4 we compare the distribution of $u-r$, $\log\textrm{M}_{\ast}$ and
$\log\textrm{M}_{\textsc{Hi}}$ for the mock and real samples, with the upper
and lower panels for ALFALFA and xGASS, respectively. The mock catalogs well
reproduce the distribution of $\textrm{M}_{\textsc{Hi}}$ for Hi-detected
galaxies in both surveys, indicating that our estimator is able to provide
unbiased Hi estimates for the general population of galaxies.

The mock catalog of ALFALFA also reproduces the distributions of $u-r$ and
$\log\textrm{M}_{\ast}$ in the real sample. The distribution of the parent
SDSS volume-limited sample used to construct the ALFALFA mock catalog is
plotted in the upper panels for comparison. As expected, the ALFALFA
detections are biased to relatively blue and Hi-rich galaxies. The full sample
of xGASS including both detections and non-detections show a similar
distribution in $u-r$ (although with a slightly higher fraction of blue
galaxies) and the same distribution in $\textrm{M}_{\textsc{Hi}}$ when
compared to the SDSS volume-limited sample. This indicates that xGASS is a
more representative sample of the general population, and thus more
appropriate for calibrating the Hi mass estimator, than ALFALFA.

In Figure 5 we plot the Hi-to-stellar mass ratio as functions of $\mu_{\ast}$,
$u-r$, $\textrm{M}_{\ast}$$/$$\textrm{M}_{\odot}$ and $R_{90}/R_{50}$ for the
real and mock catalogs. For the mock catalogs and the SDSS volume-limited
sample the Hi-to-stellar mass ratios are the values predicted by our
estimator, while the observed values are plotted for real galaxies in ALFALFA
and xGASS. The mock catalogs of Hi-detected galaxies agree well with real Hi
samples for both surveys. When compared to the model predictions based on the
SDSS volume-limited sample, the samples of detection from both surveys and the
corresponding mock catalogs are biased for gas-rich galaxies, and this bias
becomes more significant for galaxies with larger stellar mass, redder color,
higher surface density, and more concentrated light distribution. When
including both detections and non-detections (the green boxes), the xGASS mock
catalogs closely follow the trend of the volume-limited sample, suggesting
again that the xGASS survey represents well the general population of
galaxies.

##  4 Hi gas contents in galaxies and in dark matter halos

###  4.1 The Hi Mass Function of galaxies

Using the constrained Hi mass estimator given above, we can predict the Hi
mass function (HIMF) starting from an optical sample. For this purpose we use
the SDSS galaxy sample (SampleS) for which we have estimated an Hi mass for
each galaxy using the best model parameters and their dispersion. We use two
different methods to estimate the HIMF, as detailed below.

In the first method, we use the $1/V_{\rm max}$-weighting scheme following Li
& White (2009) and Chen et al. (2019) in their estimates of the stellar mass
function of galaxies. For each galaxy $i$ we determine $V_{\textrm{max},i}$,
the total comoving volume of the survey out to the maximum redshift
$z_{\textrm{max},i}$ at which the galaxy in question would meet the apparent
magnitude criterion of our sample $r\leq 17.6$. Luminosity evolution and
$K$-correction are included in the calculation of $z_{\textrm{max},i}$. The
HIMF is then estimated as

$\Phi(\textrm{M}_{\textsc{Hi}})\Delta\textrm{M}_{\textsc{Hi}}=\sum_{i}(f_{\textrm{sp},i}V_{\textrm{max},i})^{-1}\frac{\rho_{\rm
u}}{\rho(V_{\textrm{max},i})},$ (12)

where the sum runs over all galaxies with Hi mass in the range
$\textrm{M}_{\textsc{Hi}}\pm 0.5\Delta\textrm{M}_{\textsc{Hi}}$. Here
$f_{\textrm{sp},i}$ is the spectroscopic completeness which varies across the
survey area and is defined as the fraction of photometrically-selected targets
that are spectroscopically observed and included in our sample. The last term
in the equation $\rho_{\rm u}/\rho(V_{\textrm{max},i})$ is the ratio of the
average mass density within the whole survey volume to that within
$V_{\textrm{max},i}$, as determined in Chen et al. (2019) based on the ELUCID
simulation (Wang et al., 2014), a constrained simulation in the SDSS volume,
in order to account for the effect of cosmic variance caused by the limited
sample volume at low redshift.

In the second method, we first estimate the $r$-band luminosity function (LF)
from SampleS, $\phi(M_{r})$, using the same $1/V_{\rm max}$-weighting scheme
as described above. We then divide all galaxies in SampleS into a successive
sequence of non-overlapping $M_{r}$ bins, $M_{r,k}\leq
M_{r}<M_{r,k+1}~{}(k=1,...,N_{m})$, with a fixed logarithmic interval of
$\log(M_{r,k+1}/M_{r,k})=0.2$ dex. For the $k$-th interval we obtain the Hi
mass distribution of sample galaxies $\Psi_{k}(\textrm{M}_{\textsc{Hi}}\pm
0.5\Delta\textrm{M}_{\textsc{Hi}})$, normalized so as to have
$\sum\Psi_{k}\Delta\textrm{M}_{\textsc{Hi}}=1$. The HIMF of the whole sample
can then be estimated as

$\Phi(M_{\textsc{Hi}})=\sum_{k=1}^{N_{m}}\left[\Psi_{k}(M_{\textsc{Hi}})\int_{M_{r,k}}^{M_{r,k+1}}\phi(M_{r})dM_{r}\right],$
(13)

where the sum runs over all the absolute magnitude bins for a given Hi mass
range $\textrm{M}_{\textsc{Hi}}\pm 0.5\Delta\textrm{M}_{\textsc{Hi}}$. In
principle the two methods should lead to identical results. The advantage of
the second one is that the complex selection effects of the galaxy sample are
already accounted for when estimating the LF, so that one can obtain the HIMF
simply by counting the sample galaxies in bins of $\textrm{M}_{\textsc{Hi}}$.
This method is valid only when the sample galaxies at a given luminosity are
not biased in their Hi mass distribution. This should be true for an
optically-selected sample like the SDSS that does not involve any Hi-related
selections. We use the LF estimated from SampleS by ourselves, not only for
self-consistency but also to take advantage of the $1/V_{\textrm{max}}$
weights given in Chen et al. (2019) with corrections for the cosmic variance
effect in the local universe.

Figure 6: The Hi mass function of galaxies. Left panel: The Hi mass function
of all galaxies in the local universe. The green open circles with errorbars
are the results of this work using $\rm 1/V_{max}$ weighting method. The green
dashed line is the fit to the green circles. The dark shaded region shows the
HIMF derived from our luminosity function. The thick dark dashed line is the
HIMF derived from conditional luminosity function (see subsection 5.2). The
magenta dashed line and open circles are the results of ALFALFA and xGASS (see
subsection 3.2). The grey shaded region indicates the Hi mass regime where
SampleS is incomplete and so the HIMF (green symbols) can only be considered
as lower limits of the true HIMF. Right panel: the Hi mass function of red and
blue galaxies. The open circles are derived from $\rm 1/V_{max}$ method. The
colored areas are derived from luminosity function. The solid lines are
derived from conditional luminosity function. Table 1: Schechter function
parameters of CHIMFs

$\log(M_{200}/M_{\odot})$ | $\Phi^{\ast}_{\textsc{Hi}}$ | $M^{\ast}_{\textsc{Hi}}$ | $\alpha$
---|---|---|---
All galaxies:
$[12.00,12.34)$ | $0.34\pm 0.04$ | $10.13\pm 0.04$ | $-1.17\pm 0.05$
$[12.34,12.68)$ | $0.26\pm 0.02$ | $10.29\pm 0.03$ | $-1.41\pm 0.02$
$[12.68,13.03)$ | $0.31\pm 0.02$ | $10.34\pm 0.02$ | $-1.54\pm 0.02$
$[13.03,13.37)$ | $0.92\pm 0.11$ | $10.14\pm 0.04$ | $-1.42\pm 0.03$
$[13.37,13.71)$ | $2.55\pm 0.24$ | $9.99\pm 0.04$ | $-1.23\pm 0.03$
$[13.71,14.05)$ | $2.79\pm 0.32$ | $10.15\pm 0.03$ | $-1.50\pm 0.04$
$[14.05,14.39)$ | $8.66\pm 1.10$ | $9.99\pm 0.04$ | $-1.41\pm 0.03$
$[14.39,14.73)$ | $8.85\pm 1.09$ | $10.09\pm 0.04$ | $-1.56\pm 0.03$
$[14.73,15.08)$ | $26.99\pm 3.15$ | $10.09\pm 0.03$ | $-1.53\pm 0.03$
Red galaxies:
$[12.00,12.34)$ | $0.10\pm 0.02$ | $9.84\pm 0.06$ | $-1.45\pm 0.06$
$[12.34,12.68)$ | $0.17\pm 0.03$ | $10.01\pm 0.05$ | $-1.43\pm 0.05$
$[12.68,13.03)$ | $0.16\pm 0.02$ | $10.22\pm 0.03$ | $-1.59\pm 0.03$
$[13.03,13.37)$ | $0.31\pm 0.06$ | $10.19\pm 0.06$ | $-1.51\pm 0.04$
$[13.37,13.71)$ | $0.89\pm 0.13$ | $10.00\pm 0.06$ | $-1.28\pm 0.04$
$[13.71,14.05)$ | $0.55\pm 0.09$ | $10.29\pm 0.05$ | $-1.68\pm 0.04$
$[14.05,14.39)$ | $3.75\pm 0.81$ | $9.97\pm 0.06$ | $-1.46\pm 0.05$
$[14.39,14.73)$ | $3.56\pm 0.71$ | $10.11\pm 0.06$ | $-1.69\pm 0.05$
$[14.73,15.08)$ | $9.02\pm 2.05$ | $10.16\pm 0.06$ | $-1.71\pm 0.05$
Blue galaxies:
$[12.00,12.34)$ | $0.30\pm 0.05$ | $10.13\pm 0.06$ | $-1.02\pm 0.07$
$[12.34,12.68)$ | $0.22\pm 0.03$ | $10.21\pm 0.05$ | $-1.26\pm 0.05$
$[12.68,13.03)$ | $0.19\pm 0.01$ | $10.36\pm 0.02$ | $-1.45\pm 0.01$
$[13.03,13.37)$ | $0.66\pm 0.05$ | $10.09\pm 0.03$ | $-1.34\pm 0.02$
$[13.37,13.71)$ | $1.72\pm 0.12$ | $9.97\pm 0.02$ | $-1.19\pm 0.02$
$[13.71,14.05)$ | $2.76\pm 0.22$ | $10.03\pm 0.03$ | $-1.30\pm 0.03$
$[14.05,14.39)$ | $5.23\pm 0.40$ | $9.98\pm 0.02$ | $-1.35\pm 0.03$
$[14.39,14.73)$ | $6.33\pm 0.49$ | $10.03\pm 0.03$ | $-1.36\pm 0.02$
$[14.73,15.08)$ | $23.48\pm 1.42$ | $9.96\pm 0.02$ | $-1.21\pm 0.02$
Satellite galaxies:
$[12.00,12.34)$ | $0.06\pm 0.01$ | $10.02\pm 0.02$ | $-1.65\pm 0.03$
$[12.34,12.68)$ | $0.14\pm 0.02$ | $10.04\pm 0.03$ | $-1.63\pm 0.03$
$[12.68,13.03)$ | $0.24\pm 0.02$ | $10.10\pm 0.03$ | $-1.67\pm 0.03$
$[13.03,13.37)$ | $0.79\pm 0.08$ | $10.04\pm 0.03$ | $-1.49\pm 0.03$
$[13.37,13.71)$ | $2.23\pm 0.18$ | $9.95\pm 0.03$ | $-1.27\pm 0.03$
$[13.71,14.05)$ | $2.82\pm 0.27$ | $10.09\pm 0.03$ | $-1.51\pm 0.04$
$[14.05,14.39)$ | $8.37\pm 1.07$ | $9.97\pm 0.04$ | $-1.42\pm 0.04$
$[14.39,14.73)$ | $8.58\pm 1.00$ | $10.08\pm 0.03$ | $-1.58\pm 0.03$
$[14.73,15.08)$ | $27.34\pm 3.05$ | $10.08\pm 0.03$ | $-1.53\pm 0.03$

Note. — $\Phi^{\ast}_{\textsc{Hi}}$, $M^{\ast}_{\textsc{Hi}}$, and $\alpha$
are the characteristic amplitude, characteristic mass and the low-mass end
slope of a Schechter function. In the last two halo bins, there is nearly no
blue central galaxies so no fitting parameters are given.

Table 2: Gaussian function parameters of CHIMFs

$\log(M_{200}/M_{\odot})$ | $\Phi^{\ast}_{\textsc{Hi}}$ | $M^{\ast}_{\textsc{Hi}}$ | $\sigma$
---|---|---|---
Central galaxies:
$[12.00,12.34)$ | $0.606\pm 0.017$ | $0.479\pm 0.010$ | $8.877\pm 0.013$
| $0.369\pm 0.024$ | $0.391\pm 0.016$ | $9.714\pm 0.026$
$[12.34,12.68)$ | $0.742\pm 0.024$ | $0.460\pm 0.010$ | $9.072\pm 0.016$
| $0.217\pm 0.021$ | $0.360\pm 0.020$ | $9.810\pm 0.037$
$[12.68,13.03)$ | $0.794\pm 0.032$ | $0.455\pm 0.009$ | $9.271\pm 0.021$
| $0.165\pm 0.009$ | $0.442\pm 0.014$ | $9.932\pm 0.024$
$[13.03,13.37)$ | $0.851\pm 0.029$ | $0.444\pm 0.007$ | $9.371\pm 0.017$
| $0.107\pm 0.009$ | $0.467\pm 0.016$ | $9.940\pm 0.041$
$[13.37,13.71)$ | $0.867\pm 0.024$ | $0.451\pm 0.006$ | $9.445\pm 0.014$
| $0.099\pm 0.007$ | $0.443\pm 0.026$ | $10.013\pm 0.034$
$[13.71,14.05)$ | $0.895\pm 0.022$ | $0.457\pm 0.007$ | $9.559\pm 0.011$
| $0.072\pm 0.010$ | $0.573\pm 0.024$ | $10.012\pm 0.069$
$[14.05,14.39)$ | $0.921\pm 0.037$ | $0.453\pm 0.009$ | $9.615\pm 0.016$
| $0.016\pm 0.001$ | $0.419\pm 0.033$ | $10.051\pm 0.041$
$[14.39,14.73)$ | $0.930\pm 0.054$ | $0.412\pm 0.015$ | $9.686\pm 0.024$
| $-$ | $-$ | $-$
$[14.73,15.08)$ | $0.948\pm 0.063$ | $0.425\pm 0.019$ | $9.719\pm 0.028$
| $-$ | $-$ | $-$

Note. — $\Phi^{\ast}_{\textsc{Hi}}$, $M^{\ast}_{\textsc{Hi}}$, and $\sigma$
are the amplitude, center and width of a Gaussian profile. For each halo bin
the CHIMF of central galaxies are fitted with two Gaussians, corresponding to
the two rows of parameters in the table (The first row for red and the second
row for blue).

In the left panel of Figure 6, we show the HIMF derived for our sample using
the two methods described above, with green circles with error bars for the
first method and the grey-shaded band for the second method. The errors of the
HIMF from both methods are estimated using the bootstrap re-sampling method,
which accounts for uncertainties due to sampling, but does not include the
cosmic variance and uncertainties in the Hi mass estimate. As expected, the
HIMFs from the two methods are in good agreement with each other over the
entire mass range. The HIMF increases as one goes from high to low
$\textrm{M}_{\textsc{Hi}}$ but at $\textrm{M}_{\textsc{Hi}}\sim
10^{8}\textrm{M}_{\odot}$ the HIMF flattens and even decreases slightly
towards the low-mass end. This behavior at the low-mass end is not real, but
caused by the incompleteness of the SDSS sample at stellar masses below
$\textrm{M}_{\ast}$$\sim 10^{8}$$\textrm{M}_{\odot}$, corresponding to a mean
Hi mass of $\textrm{M}_{\textsc{Hi}}$$\sim 10^{8}$$\textrm{M}_{\odot}$
according to our Hi mass estimator. Thus, our HIMF at
$\textrm{M}_{\textsc{Hi}}<10^{8}\textrm{M}_{\odot}$ should be taken as a
‘lower limit’. We fit our HIMF at
$\textrm{M}_{\textsc{Hi}}>10^{8}\textrm{M}_{\odot}$ with a single Schechter
function (Schechter, 1976) and plot the result as the dashed green line in
Figure 6. The parameters of the fit are $\Phi_{\textsc{Hi}}^{\ast}=2.61\pm
0.2\times 10^{-3}{\rm Mpc^{-3}dex^{-1}}$, $\log
M_{\textsc{hi}}^{\ast}/M_{\odot}=10.07\pm 0.02$, and $\alpha=-1.51\pm 0.02$.
Integrating the best fit of our HIMF gives the cosmic density of the Hi gas
locked in local galaxies: $\Omega_{\rm{\textsc{Hi}}}=4.1\pm 0.1\times
10^{-4}$, which is 17% and 5% higher than that from HIPASS (Zwaan et al.,
2005) and ALFALFA (Jones et al., 2018), respectively. For comparison we plot
the Schechter function fit of the ALFALFA HIMF derived by Jones et al. (2018)
and the xGASS HIMF derived by us in subsection 3.2. The HIMF estimated from
the optical sample agrees well with both the ALFALFA and xGASS HIMFs at
$\log\textrm{M}_{\textsc{Hi}}$$>9.5$ where we have all the HIMF estimates
available. At lower masses our HIMF differs from that of ALFALFA in that it
has increasingly higher amplitudes and a steeper slope as the Hi mass
decreases from $\textrm{M}_{\textsc{Hi}}\sim 10^{9}\textrm{M}_{\odot}$ down to
$\textrm{M}_{\textsc{Hi}}\sim 10^{8}\textrm{M}_{\odot}$. The difference is a
factor of about 2 at $\textrm{M}_{\textsc{Hi}}\sim 10^{8}\textrm{M}_{\odot}$.

In the right panel of Figure 6 we show the HIMF derived in the same way, but
separately for red and blue galaxies. For this purpose we have classified each
galaxy in SampleS as red or blue according to the luminosity-dependent $u-r$
color demarcation given by Baldry et al. (2004). As can be seen, the two
methods provide very consistent results even when galaxies are divided into
red and blue subsamples. The HIMF of the blue galaxy population is higher than
that of the red population over the entire $\textrm{M}_{\textsc{Hi}}$ range
except for the largest masses where the amplitudes of the two HIMFs become
comparable. Clearly, the total HIMF of all galaxies is dominated by blue
galaxies.

Figure 7: Conditional Hi mass function. The orange and green open circles
represent the CHIMFs of satellites and centrals respectively. The black solid
lines represent the CHIMFs of all galaxies. For comparison, the CHIMF of all
galaxies in halos of $\log M_{h}\in[13.4,13.7)$ is plotted as the black dashed
line in each panel. The red and blue solid lines are the CHIMFs of red and
blue galaxies. The fitting result of central galaxies is shown as red and blue
dashed lines for red and blue central galaxies. The colored numbers in each
panel shows the contribution of
central(green)/satellite(orange)/red(red)/blue(blue) galaxies to the total Hi
mass, calculated using corresponding CHIMFs above $\rm
M_{HI}>10^{8}M_{\odot}$. The grey shaded region shows the Hi mass range where
the CHIMFs suffer from incompleteness.

###  4.2 The conditional Hi Mass Function

The second statistic we can predict is the conditional Hi mass function
(CHIMF), $\Phi(M_{\textsc{Hi}}|M_{h})$, which is the HIMF of galaxies hosted
by dark matter halos of a given mass $M_{h}$. This statistic has not been
obtained directly from observation due to the lack of data. Here we perform
the first attempt to estimate the CHIMFs using SampleS with Hi masses
estimated above. We adopt the second method described in subsection 4.1,
except that we use the conditional luminosity function (CLF;
$\phi(M_{r}|M_{h})$) instead of the total LF in Equation 13. The CLFs of
galaxies in the local Universe were estimated by Lan et al. (2016, hereafter
Lan16) based on the SDSS galaxy group catalog (Yang et al., 2007) and the SDSS
photometric sample down to a $r$-band absolute magnitude of $M_{r}=-12$. Lan16
considered a number of halo mass bins ranging from $M_{h}\sim
10^{12}\textrm{M}_{\odot}$ up to $\sim 10^{15}\textrm{M}_{\odot}$. For each of
the bins they obtained the CLFs separately for red and blue galaxies (divided
by a luminosity-dependent $u-r$ color demarcation), as well as for central and
satellite galaxies (with the central/satellite classification provided by the
SDSS group catalog). Correspondingly, we estimate the CHIMFs for red/blue and
central/satellite galaxies separately, applying the same color cut and
central/satellite classification as in Lan16 to divide our sample galaxies
into red/blue or central/satellite subsamples. We divide galaxies of a given
subsample into bins of $M_{r}$. The CHIMF for the subsample in each halo mass
bin, $\Phi(M_{\textsc{Hi}}|M_{h})$, is then given by Equation 13, where the LF
is replaced by the corresponding CLF $\phi(M_{r}|M_{h})$. We note that we have
adopted a luminosity limit of $M_{r}=-14.4$ when estimating the CHIMFs, to
avoid large uncertainties caused by the rather noisy distributions of
$\textrm{M}_{\textsc{Hi}}$ at lower luminosities.

In Figure 7 we show the CHIMFs of red/blue and central/satellite galaxies as
well as the total CHIMF, for different halo mass bins as indicated in each
panel. The total CHIMF of the halo mass bin of
$13.37\leq\log(\textrm{M}_{\textrm{h}}/\textrm{M}_{\odot})<13.71$ is repeated
as the dashed black line in each panel for reference. From the CHIMFs we
calculate the fraction of the total Hi mass in dark halos of given mass as
contributed by red, blue, central and satellite galaxies with
$\textrm{M}_{\textsc{Hi}}$$>10^{8}$$\textrm{M}_{\odot}$. These mass fractions
are indicated in each panel with different colors for different types of
galaxies. A Schechter function can well describe all the HIMFs except those
for central galaxies where a double Gaussian profile provides a better
description for low-mass halos and a single Gaussian profile works better for
high-mass halos. When divided into red and blue subsamples, central galaxies
at fixed halo mass can each be modeled by a single Gaussian, as is shown in
the figure where red and blue dashed lines are the Gaussian fits of red and
blue centrals, respectively. The Schechter function parameters for the total
sample and the subsamples of red, blue and satellite galaxies in each halo
mass bin are listed in Table 1, while the parameters of the two Gaussian
functions for sub-samples of central galaxies are listed in Table 2.

Figure 8: The left panel shows the average Hi mass of dark matter halos and
the right panel shows the average Hi mass of central galaxies. In both panels,
the black open circles are the result of Guo et al. (2020) in $N_{g}\geq 1$
case. For central galaxies, they have been corrected for the confusion effect
using our estimator. The black solid line is our result derived from group
catalog. The brown, orange and green solid lines represent the result of
L-GALAXIES semi-analytical model, Illustris simulation and IllustrisTNG
simulation, respectively. The black upward and downward triangles represent
the Milky Way and M31, respectively. In the left panel, the pink dashed line
is the result of our CHIMF and the pink upward arrow with error bars
represents Virgo Cluster. In the right panel, the black dashed line is the
average central Hi mass derived from our CHIMFs, and the grey shaded region
represents the 16% to 84% percentile of the central CHIMFs in each halo mass
bin. The red and blue dashed lines are the mean Hi mass of red and blue
central galaxies from our CHIMFs. The magenta downward arrows represent the
upper limit of Hi masses of central galaxies in xGASS sample.

As can be seen from Figure 7, the amplitude of the total HIMF increases
monotonically with increasing halo mass, indicating a positive correlation of
the total Hi mass of galaxies with the dark matter mass of their host halos.
When galaxies are divided by color, we find that the red and blue populations
present similar CHIMFs for halos of different masses, in terms of both
amplitude and shape, which results in a comparable Hi mass fraction for the
two populations in halos of different masses. Quantitatively, the Hi mass
fraction in red (blue) galaxies increases (decreases) steadily from 21% (79%)
at $\textrm{M}_{\textrm{h}}$$\sim 10^{12}$$\textrm{M}_{\odot}$ to 56% (44%) at
$\textrm{M}_{\textrm{h}}$$\sim 10^{15}$$\textrm{M}_{\odot}$. It has been well
established that blue galaxies dominate the central galaxy population in
intermediate-to-low mass halos, and red galaxies are either (massive) centrals
in massive halos or (less-massive) satellites in halos of different masses
(Yang et al., 2009; Lan et al., 2016). Thus, massive halos are expected to
contain more red galaxies than blue galaxies. The fact that red and blue
galaxies have comparable CHIMFs at given halo mass can be understood from the
color-dependent Hi-to-stellar mass ratio of galaxies, which is smaller for
redder galaxies.

For central galaxies, the CHIMFs appear to depend weakly on halo mass, with
the average $\textrm{M}_{\textsc{Hi}}$ increasing slightly and the width
decreasing slightly with increasing halo mass. The two Gaussian components
correspond to contributions by central galaxies of red and blue colors,
respectively. Blue centrals contribute significantly to the total Hi of the
central galaxy population only in low-mass halos, with a contribution that is
comparable to that of red centrals in Milky Way-like halos and declines with
halo mass. At $\textrm{M}_{\textrm{h}}$$\gtrsim 10^{14}$$\textrm{M}_{\odot}$,
the CHIMF of centrals is dominated by red galaxies, with a negligible
contribution from blue centrals. Central galaxies as a whole contribute a
significant amount of Hi mass only in low-mass halos, with a fraction of 75%
and 58% in the two lowest-mass halo bins, respectively, and the fraction
decreases rapidly at larger halo masses, becoming less than 10% when the halo
mass exceeds $\textrm{M}_{\textrm{h}}\sim 10^{14}\textrm{M}_{\odot}$. In
contrast to centrals, the satellite galaxy population has a Schechter CHIMF
that varies with halo mass in a way similar to the total CHIMF. Consequently,
the satellite population contributes a larger fraction of the total Hi mass as
the halo mass increases.

Figure 9: Left panel: the Hi mass function. The magenta dashed line is
ALFALFA’s result. The orange solid line is derived from GAMA LF. The black
shaded region and solid line are derived from our LF with and without CV
correction respectively. The grey shaded region shows the mass range suffering
from sample incompleteness. Right panel: the orange solid line is the GAMA LF.
The brown and cyan points with errorbars represent the LF of Blanton et al.
(2005) derived by SWML and $\rm 1/V_{max}$ method respectively. The dark
region and points show our LF with and without CV correction.

###  4.3 The Hi-halo mass relation

One can obtain the total Hi mass locked in member galaxies in dark matter
halos by integrating the CHIMFs derived above. The result is shown as dashed
lines in Figure 8, where we plot the total Hi mass in all galaxies (left
panel) and central galaxies (right panel) as a function of halo mass. The
total Hi mass of all galaxies increases monotonically with halo mass, from
$\textrm{M}_{\textsc{Hi}}$$\sim 8\times 10^{9}$$\textrm{M}_{\odot}$ for Milky
Way-like halos with $\textrm{M}_{\textrm{h}}$$\sim
10^{12}$$\textrm{M}_{\odot}$ up to $\textrm{M}_{\textsc{Hi}}$$\sim 7\times
10^{11}$$\textrm{M}_{\odot}$ for Coma-like halos with
$\textrm{M}_{\textrm{h}}$$\sim 5\times 10^{14}$$\textrm{M}_{\odot}$. The total
Hi mass in central galaxies shows a rather weak dependence on halo mass, with
$\textrm{M}_{\textsc{Hi}}$ slightly increasing from $\sim 4\times
10^{9}$$\textrm{M}_{\odot}$ to $\sim 8\times 10^{9}$$\textrm{M}_{\odot}$.

The Hi-to-halo mass (HIHM) relation can also be estimated from a galaxy group
catalog (Lim et al., 2017) for which we have estimated the Hi mass for each
galaxy using our estimator. To that end, we divide galaxy groups into a set of
halo mass bins, and for groups in each bin we sum up the Hi mass of all member
galaxies (or central galaxies) to obtain the HIHM relation for all galaxies
(or central galaxies). Each galaxy is weighted in the same way as in Equation
12 in order to correct for the sample incompleteness and cosmic variance
effect. The relations derived this way are plotted as the solid black lines in
Figure 8 over the range from $\textrm{M}_{\textrm{h}}$$\sim
10^{11}$$\textrm{M}_{\odot}$to $\textrm{M}_{\textrm{h}}$$\sim 5\times
10^{13}\textrm{M}_{\odot}$. The error bars are estimated from the scatter
among 20 bootstrap samples. For central galaxies, the HIHM relation derived
from the group catalog is consistent with the average relation derived from
the CHIMFs. For all galaxies as a whole, however, the relation derived from
the group catalog agrees with the CHIMF-based result only at
$\textrm{M}_{\textrm{h}}$$\sim 1-2\times 10^{12}$$\textrm{M}_{\odot}$. Above
this mass, the group catalog gives lower Hi mass, and the difference increases
with halo mass, reaching $\sim 0.3$ dex at $\textrm{M}_{\textrm{h}}$$\sim
5\times 10^{13}$$\textrm{M}_{\odot}$. Possible reasons for this difference
will be discussed in the next section.

The HIHM relations estimated from the group catalog are in good agreement with
those of Guo et al. (2020) obtained by stacking the ALFALFA datacubes, as
shown by open circles in Figure 8. We have corrected the confusion effect for
the average $\textrm{M}_{\textsc{Hi}}$ mass of central galaxies measured by
Guo et al. (2020) by subtracting the contribution of satellite galaxies that
reside near their central galaxy. Following Guo et al. (2020), we use an
angular aperture of max{200kpc/$D_{A},8^{\prime}$} and a maximum velocity
separation of 300 km $s^{-1}$ to define the nearness. For comparison, the left
panel of the same figure shows the total Hi mass for the Milky Way (Kalberla &
Kerp, 2009), M31 (Chemin et al., 2009) and the Virgo Cluster. The Hi mass of
the Virgo Cluster is estimated by ourselves by applying our Hi estimator to
the Extended Virgo Cluster Catalog (EVCC; Kim et al., 2014). We only use
galaxies in EVCC that are brighter than $M_{r}=-14.4$, corresponding to the
luminosity limit adopted above for our CHIMFs, and are classified to be member
galaxies of the Virgo Cluster based on the spherical symmetric infall model
(Praton & Schneider, 1994). Thus, the resulting Hi mass should perhaps be
considered as a lower limit of the total Hi mass in the cluster, as some Hi
may be contained in galaxies that are not included in our estimate. Using
central galaxies in SampleX we can also estimate the average HIHM relation,
which is shown in the right-hand panel 8. We have included both detections and
non-detections in this estimation, and so the Hi masses from xGASS should be
regarded as upper limits. As one can see, our HIHM relations are in broad
agreement with all these previous measurements based on real observations of
the Hi gas, indicating that our estimator is reliable also for member galaxies
of groups.

Figure 10: The Hi mass in dark matter halos under different luminosity limits.
In each panel, the black points represent the average Hi mass of all galaxies
in halos derived using group catalog. The horizontal black dotted line shows
the result of Guo et al. (2020). The red/blue/green/orange lines represent the
Hi mass in red/blue/central/satellite galaxies from our CHIMFs. The pink
upward arrows show the result in the Virgo Cluster.

##  5 Discussion

###  5.1 HIMF at the low mass end

Figure 6 shows that the HIMFs predicted by our estimator are higher and
steeper than the HIMF obtained from ALFALFA at masses below a few $\times
10^{9}$$\textrm{M}_{\odot}$. This discrepancy is predominantly caused by
cosmic variance, which is corrected in our HIMFs but not in the ALFALFA
result. The cosmic variance is taken into account through the last term of
Equation 12, which was carefully determined by Chen et al. (2019) based on the
ELUCID simulation, a constrained simulation with initial conditions that are
reconstructed to reproduce the dark matter density field in the SDSS volume
(Wang et al., 2014). The local volume within $z\sim 0.03$ has a density that
is lower than the mean density of the universe. Consequently, the SDSS under-
represents faint galaxies in the universe because it can detect these galaxies
only at the very low redshift. ALFALFA is also limited to very low redshift,
and so is also affected by the cosmic variance.

The effect of cosmic variance is shown clearly in Figure 9. In the right-hand
panel, we show the $r$-band luminosity functions (LFs), $\Phi(M_{r})$,
estimated from SampleS with and without the cosmic variance correction, and
compare them with the LFs obtained by Blanton et al. (2005) using an early
SDSS sample, and with the LF obtained by Loveday et al. (2015) using the GAMA
survey (Driver et al., 2011). The latter two studies did not correct for the
cosmic variance effect, and it is not surprising that their LFs are flatter
and lower than our measurements at the faint end. Our LF without including the
cosmic variance correction shows a similarly lower amplitude at the low-mass
end, confirming the importance of cosmic variance in the local universe. As
pointed out in Blanton et al. (2005), the dip at $M_{r}\sim-17$ and the upturn
at the faint end seen both in our uncorrected LF and the LF of Blanton et al.
(2005) (the one estimated with the $1/V_{\textrm{max}}$ method) are produced
by cosmic variance caused by the presence of large-scale structure in the
sample volume, as the $1/V_{\textrm{max}}$ method is not capable of correcting
such variance. The LF of Blanton et al. (2005) estimated using the step-wise
maximum likelihood does not show such fluctuations at the faint end and is
more consistent with the GAMA LF, but both are still lower than our LF
corrected for the cosmic variance. In the left-hand panel of Figure 9 we
compare the HIMFs derived with Equation 13, using our LFs with and without
cosmic variance correction, as well as using the GAMA LF. The HIMFs without
the correction are in good agreement with the ALFALFA result, demonstrating
that the difference in the low-mass end between our corrected HIMF and the
ALFALFA result is produced by the different treatments of the cosmic variance
effect.

###  5.2 Conditional Hi mass functions

As demonstrated in subsection 3.3, our Hi-to-stellar mass estimator reproduces
well the statistical properties of the mock catalogs that have the same
selection effects as ALFALFA and xGASS. The application of our estimator to
the SDSS sample also yields HIMFs that are in good agreement with the ALFALFA
HIMF. All these indicate that our predicted CHIMFs are reliable. As a further
test, we combine our CHIMFs for different halo masses to predict a total HIMF
and compare it with the HIMFs estimated using other methods. The black dashed
line in the left panel of Figure 6 shows the HIMF obtained from the sum of
CHIMFs of different halo mass bins, each weighted by the number density of
halos in the halo mass bin calculated using the halo mass function from Driver
et al. (2022). Since the CHIMFs can be estimated only down to Milky Way-size
halos with $\textrm{M}_{\textrm{h}}$$\sim 10^{12}$$\textrm{M}_{\odot}$ due to
the mass limit in the SDSS group catalog, we have included the contribution of
halos of lower masses following the method described in Lan et al. (2016). As
can be seen from the figure, the HIMF estimated from the CHIMFs agrees very
well with the HIMFs obtained with other methods.

As we have shown, the Schechter function provides a good description for the
CHIMFs of all/red/blue/satellite populations in halos of different masses,
while a double Gaussian function is suitable for central galaxies (see Table 1
and Table 2). These results have important implications for models of galaxy
formation. For instance, a recent study by Li et al. (2022) used the Hi-to-
halo mass relation and clustering of Hi galaxies obtained by Guo et al. (2017)
and Guo et al. (2020) to constrain empirical models of the Hi-halo relation by
assuming that the Hi mass of galaxies depends only on dark matter halo mass
plus a secondary halo property, such as halo formation time and concentration.
Their best-fit model predicts CHIMFs that have a log-normal form, quite
different from the Schechter function found here (Figure 7). This implies that
the model must have missed some key ingredients that regulate the Hi gas
content of galaxies. Our Hi estimator involves four parameters, among which
the stellar mass and color may be modeled by halo mass and formation time (see
Wechsler & Tinker, 2018, for a review). The other two parameters, surface
stellar mass density and concentration, characterize the structure of galaxies
and have been found to depend weakly on environments and properties of their
host halos (e.g. Li et al., 2006; Blanton & Moustakas, 2009). On the other
hand, the Hi mass fraction of galaxies may be stripped by environmental
processes operating in dense regions, such as tidal stripping and ram-pressure
stripping (e.g. Gunn & Gott, 1972; Cowie & Songaila, 1977; Nulsen, 1982;
Boselli & Gavazzi, 2006; Wang et al., 2021, 2022). These effects are expected
to depend on both the stellar mass and surface mass density of galaxies (e.g.
Li et al., 2012; Zhang et al., 2013). It is thus likely that physical
processes related to both galaxy structure and environments may also play a
role in determining Hi gas contents of galaxies.

When divided according to color, red and blue central galaxies follow a
Gaussian CHIMF, rather than a Schechter function (Figure 7). This is similar
to the conditional stellar mass (or luminosity) function of central galaxies
found in previous studies (e.g. Lan et al., 2016; Yang et al., 2008). However,
as shown in Figure 8, the average Hi mass of central galaxies is limited to a
small range, $\textrm{M}_{\textsc{Hi}}$$\sim 4-8\times
10^{9}$$\textrm{M}_{\odot}$ for halos of Milky Way masses or larger. The weak
dependence on the halo mass holds even when central galaxies are divided into
red and blue subsamples (see Figure 8). This result implies that, unlike the
stellar content of central galaxies, which closely follows the growth of their
host dark matter halos, the gas contents and gas-related evolution of central
galaxies may be driven by internal and local processes. In massive halos, the
central galaxies are predominantly red. Considering the large stellar mass and
relatively small star formation rate as inferred by the red color, it is
perhaps not surprising to see the nearly constant Hi mass of these galaxies.
However, the weak dependence on halo mass itself suggests that their gas-
related quenching are unlikely regulated by the host halos, but rather some
internal processes may be able to shut off star formation in these galaxies,
making them red but do not reduce their Hi gas contents significantly.

The total Hi mass of galaxies as a function of dark matter halo mass, obtained
using our CHIMFs and the SDSS galaxy group catalog, should be understood as
the total Hi mass locked in optically-selected galaxies down to a certain
limiting luminosity. Note that the optical samples used to obtain the CLFs and
to construct the group catalog have different depths. The group catalog of
Yang et al. (2007) was constructed using the SDSS spectroscopic galaxy sample
limited to $M_{r}\sim-19.5$. On the other hand, the CLFs of Lan et al. (2016)
used the same group catalog to identify groups/halos, but with member galaxies
identified from a much deeper photometric sample. This explains the increasing
difference in the total Hi mass between the prediction of the CLFs and that of
the group catalog, as seen in Figure 8. The difference is larger for more
massive halos where groups are observed to higher redshifts and thus may miss
more satellite galaxies below the luminosity limit. To better understand this
effect, we have repeated the estimation of the Hi-to-halo mass relation using
different luminosity limits, from $M_{r}=-20.4$ down to $M_{r}=-14.4$. This
result is shown in Figure 10, with the solid black line representing results
for all galaxies, and colored lines for red/blue and central/satellite
galaxies, respectively. The total Hi mass based on the group catalog is
plotted as solid circles with error bars in each panel, and the lower limit of
Hi mass predicted for the Virgo EVCC sample is plotted as pink squares with
upward arrows in the last panel. For comparison, the Hi mass from Guo et al.
(2020) is plotted as a horizontal dotted line in each panel.

As expected, the Hi mass increases as the magnitude limit goes deeper, but the
increase is significant only when the luminosity limit is brighter than
$M_{r}\sim-18$. For limits deeper than $M_{r}\sim-18$, the total Hi mass
estimated from the group catalog increases very little, while the mass
estimated from the CLFs increases slowly, with a total increase of
$\sim$0.1-0.2 dex. This result is true for all types of galaxies and for halos
of different masses, except for satellites in low-mass halos where the Hi mass
increases more rapidly with the limiting luminosity. This increase has little
effect on the total Hi mass, as it is dominated by central galaxies. The total
Hi mass from Guo et al. (2020) matches our results at a limiting luminosity
that increases from $M_{r}\sim-17$ at $\textrm{M}_{\textrm{h}}$$\sim
10^{12}$$\textrm{M}_{\odot}$ to $M_{r}\sim-20$ at
$\textrm{M}_{\textrm{h}}$$\sim 3-5\times 10^{13}$$\textrm{M}_{\odot}$,
indicating that the ALFALFA datacubes may have missed an increasingly
significant fraction of Hi emission from faint satellites which dominate the
total Hi mass in massive halos. However, it is also possible that our
estimator somehow overpredicts the Hi mass for satellite galaxies in massive
clusters, if the gas-to-stellar mass ratio is lower for satellites in massive
groups than for galaxies of similar optical properties in the general field
(see subsection 5.4 for more discussion). Clearly, deeper Hi surveys are
needed to resolve the Hi gas content in satellite galaxies of massive halos.
For Virgo Cluster, our estimated Hi mass (lower limits) increases with the
limiting luminosity in the same way as the SDSS-based estimates, thus with a
constant offset that should be attributed to the incompleteness of the EVCC
galaxy catalog.

###  5.3 The Hi-halo mass relation in theoretical studies

It is interesting to compare the Hi-to-halo mass relations obtained in this
paper with theoretical predictions in the literature. In Figure 8 we plot the
Hi-halo relation predicted by the L-GALAXIES semi-analytical model (Fu et al.,
2013), by the hydrodynamical simulation Illustris (Vogelsberger et al., 2014),
and by its successor IllustrisTNG (Springel et al., 2018; Marinacci et al.,
2018; Naiman et al., 2018; Nelson et al., 2018; Pillepich et al., 2018).
Although all the models predict a positive correlation between
$\textrm{M}_{\textsc{Hi}}$ and the halo mass, they over-predict the Hi mass to
various degrees. In particular, at fixed halo mass the Illustris simulation
predicts a total Hi mass that is a factor of 3-5 higher than the observed
value, while the L-Galaxies model over-predicts the Hi mass for Milky Way-mass
halos but a better match to the data at lower and higher halo masses. The TNG
result falls in between the other two models. For central galaxies, the TNG
over-predicts the Hi mass at all halo masses, while Illustris over-predicts
the central Hi mass for low-mass halos and under-predicts it for high-mass
halos. The L-Galaxies model predicts no Hi gas in centrals in halos with
masses above $\sim 10^{13}$$\textrm{M}_{\odot}$ but too much Hi gas in
centrals of Milky Way halos. These discrepancies must be related to the
various physical recipes implemented in the simulations and the semi-analytic
model to regulate gas-related processes. Our measurements of the Hi-to-halo
mass relation, as well as the CHIMFs, for both the total galaxy population and
for different sub-populations, can provide important constraints on
theoretical models.

###  5.4 Possible biases in the Hi estimator due to environmental dependence

The fact that our model prediction for group galaxies matches observational
results (see Figure 8) suggests that galaxies of similar optical properties
have similar Hi mass, irrespective of whether a galaxy resides in dense
regions or in the general field. This result supports the basic assumption
behind our Hi estimator that there is a tight relation between the Hi gas
content and other intrinsic properties of a galaxy. However, gas stripping
effects can make this relation complex. On the one hand, the Hi gas may be
removed without strongly affecting star formation on a timescale of $<1$ Gyr
(Brown et al., 2017). In this case $NUV-r$ should be better than $u-r$ for
estimating the Hi gas content, as the former is dominated by younger stellar
populations and so can more closely trace the change in star formation. On the
other hand, there are also evidences indicating that star formation can be
enhanced in cluster galaxies that experience gas stripping (Roberts et al.,
2021; Ramatsoku et al., 2019). This, in addition to the fact that the xGASS
sample is selected for the general population of galaxies and lacks very rich
groups/clusters as pointed out by Catinella et al. (2013), may lead to bias in
our Hi estimator, particularly for galaxies that have recently experienced gas
stripping. In other words, if galaxies in rich clusters like Virgo follow
different scaling relations of Hi content from the average population (field
galaxies), a bias can still be present in our Hi estimator which is calibrated
by galaxies mostly outside galaxy clusters.

In fact, Zhang et al. (2013) tested whether the Hi estimator of L12 exhibits
any significant biases in dense environment by examining the difference
between the estimated and observed
$\log(\textrm{M}_{\textsc{Hi}}/\textrm{M}_{\ast})$ for group/cluster galaxies
as a function of the environment overdensity $\ln(1+\delta)$ (see their Figure
3). They found the L12 estimator only slightly underestimates the Hi-to-
stellar mass ratio by at most 0.1 dex on average, and this bias was caused by
the Hi-rich galaxies rather than the Hi-poor galaxies in the group/clusters
(thus opposite to the expectation for an environment-induced bias). Since our
estimator is motivated from and quite similar to the L12 estimator, we expect
our estimator to also provide unbiased Hi content for galaxies in different
environments. Nevertheless, considering the different calibration samples
adopted in L12 and in our work, we should still keep in mind the possible bias
in our estimator, especially when applying it to galaxies in dense regions. A
complete and deep survey of Hi 21cm emission covering a full range of
environment is needed if one were to accurately quantify the possible
dependence of an Hi estimator on environment.

###  5.5 Limitations and Prospects

Our Hi gas estimator is calibrated with the xGASS sample which is limited to
nearby galaxies with $z<0.05$ and
$\textrm{M}_{\ast}>10^{9}\textrm{M}_{\odot}$. Therefore, in principle our
estimator can be applied only to galaxy samples selected to these redshift and
mass limits, although this estimator is shown to be able to reproduce the
ALFALFA HIMF down to masses as low as $\textrm{M}_{\textsc{Hi}}$$\sim
10^{7}$$\textrm{M}_{\odot}$ (see Figure 9). For this reason all the results
presented in this work for $\textrm{M}_{\ast}$$<10^{9}$$\textrm{M}_{\odot}$
should be taken in caution, and future studies with real Hi observations of
galaxies at lower masses would be needed to verify our predictions, e.g. the
conditional HIMFs at the low-mass end.

In the near future, next generation Hi surveys, such as WALLABY (Koribalski et
al., 2020, the Widefield ASKAP L-band Legacy All-sky Blind surveY) and LADUMA
(Holwerda et al., 2012, Looking at the distant universe with the MeerKAT
Array) will provide unprecedented Hi data in the local and distant universe
and improve our understanding of the Hi gas content of galaxies. However, the
gap between Hi and optical surveys will remain for many years, and this is
particularly true for low-mass and high-redshift galaxies. In the local
universe, galaxies with $\log M_{*}/M_{\odot}>9.0$ already have relatively
complete Hi surveys, such as xGASS and RESOLVE. Although Hi surveys, such as
ALFALFA, THINGS (Walter et al., 2008, The H I Nearby Galaxy Survey) and LITTLE
THINGS (Hunter et al., 2012, Local Irregulars That Trace Luminosity Extremes,
The HI Nearby Galaxy Survey), have provided useful information about the Hi
gas properties for lower mass galaxies, the low-mass galaxy samples obtained
are still strongly biased for Hi-rich systems. For high-redshift galaxies, so
far the only way to observe the Hi gas content is through damped Lyman
$\alpha$ systems. The estimator presented here, therefore, provides a way to
predict and study Hi gas contents for large samples of optically selected
galaxies that cover large ranges of galaxy properties and redshift. For
instance, Zhang et al. (2021) have recently attempted to estimate the HIMF of
galaxies at $z\sim 1$ by applying a two-parameter Hi mass estimator calibrated
with nearby galaxies to the DEEP2 galaxy sample. The results obtained in the
paper demonstrate that such an approach is powerful and provides an important
avenue to understand cold gas contents in galaxies and in galaxy systems.

##  6 Summary

In this paper, we introduce a new Hi mass estimator in which the logarithm of
the Hi-to-stellar mass ratio
($\log($$\textrm{M}_{\textsc{Hi}}$$/$$\textrm{M}_{\ast}$$)$) is modeled by a
linear combination of four galaxy properties: surface stellar mass density
($\log\mu_{\ast}$), color index ($u-r$), stellar mass ($\log M_{\ast}$) and
concentration index ($\log(R_{90}/R_{50})$). The variance of individual
galaxies around the mean Hi mass is modeled by a Gaussian distribution. We
calibrate the estimator with the xGASS sample, which is constructed to be
representative for the Hi content of the local galaxy population. We constrain
our model parameters using Bayesian inferences and including both detections
and non-detections from the xGASS sample. Applying our estimator to mock
catalogs that include the same selection effects as the ALFALFA survey, we are
able to statistically reproduce the distributions of both the Hi mass and
other galaxy properties of the real sample. We then apply our estimator to the
SDSS spectroscopic galaxy sample to predict the Hi mass function (HIMF) of
local galaxies, as well as the conditional Hi mass function (CHIMF) and the
Hi-halo mass relation (HIHM) for dark matter halos of different masses.

Our main results can be summarized as follows.

  * •

The HIMF derived from the optical sample can be fitted with a single Schechter
function, but with a higher amplitude and a steeper slope at Hi masses below a
few $\times 10^{9}$$\textrm{M}_{\odot}$ when compared to the ALFALFA HIMF. We
show that this difference is essentially caused by the cosmic variance of the
local volume, which is corrected for our HIMF but not for the ALFALFA HIMF.

  * •

We use our estimator to predict the CHIMFs for dark matter halos with
$\textrm{M}_{\textrm{h}}$$\gtrsim 10^{12}$$\textrm{M}_{\odot}$, for both the
total galaxy population and sub-populations of red/blue and central/satellite
galaxies. A single Schechter function can well describe all the CHIMFs at all
halo masses, except for central galaxies where a double Gaussian is preferred.

  * •

Although red galaxies dominate the galaxy population in massive halos, red and
blue galaxies have similar CHIMFs and thus contribute similarly to the total
Hi mass at given halo mass. This result is true for all halo masses.

  * •

Central galaxies contribute significantly to the total Hi mass only in low-
mass halos, while in massive halos the majority of Hi gas is associated with
satellite galaxies.

  * •

The Hi mass of central galaxies shows very weak dependence on halo mass, even
when divided into red and blue sub-populations. This is quite different from
the total Hi mass of all galaxies which increases monotonically with
increasing halo mass.

  * •

Current hydrodynamical simulations and semi-analytic models of galaxy
formation fail to reproduce the Hi-halo mass relation for both the central
population and the total population.

## Acknowledgments

This work is supported by the National Key R&D Program of China (grant No.
2018YFA0404502), and the National Natural Science Foundation of China (grant
Nos. 11821303, 11733002, 11973030, 11673015, 11733004, 11761131004,
11761141012). This work has made use of the following softwares: Numpy (Harris
et al., 2020), Scipy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), emcee
(Foreman-Mackey et al., 2013), Astropy (Astropy Collaboration et al., 2013),
GetDist (Lewis, 2019), seaborn (Waskom, 2021), and h5py (Collette et al.,
2021).

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan
Foundation, the Participating Institutions, the National Science Foundation,
the U.S. Department of Energy, the National Aeronautics and Space
Administration, the Japanese Monbukagakusho, the Max Planck Society, and the
Higher Education Funding Council for England. The SDSS Web Site is
http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the
Participating Institutions. The Participating Institutions are the American
Museum of Natural History, Astrophysical Institute Potsdam, University of
Basel, University of Cambridge, Case Western Reserve University, University of
Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the
Japan Participation Group, Johns Hopkins University, the Joint Institute for
Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and
Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences
(LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for
Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico
State University, Ohio State University, University of Pittsburgh, University
of Portsmouth, Princeton University, the United States Naval Observatory, and
the University of Washington.

## References

  * Adelman-McCarthy et al. (2008) Adelman-McCarthy, J. K., Agüeros, M. A., Allam, S. S., et al. 2008, The Astrophysical Journal Supplement Series, 175, 297, doi: 10.1086/524984
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Baldry et al. (2004) Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681, doi: 10.1086/380092
  * Barnes & Haehnelt (2014) Barnes, L. A., & Haehnelt, M. G. 2014, MNRAS, 440, 2313, doi: 10.1093/mnras/stu445
  * Blanton et al. (2011) Blanton, M. R., Kazin, E., Muna, D., Weaver, B. A., & Price-Whelan, A. 2011, AJ, 142, 31, doi: 10.1088/0004-6256/142/1/31
  * Blanton et al. (2005) Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, ApJ, 631, 208, doi: 10.1086/431416
  * Blanton & Moustakas (2009) Blanton, M. R., & Moustakas, J. 2009, ARA&A, 47, 159, doi: 10.1146/annurev-astro-082708-101734
  * Blanton & Roweis (2007) Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734, doi: 10.1086/510127
  * Blanton et al. (2005) Blanton, M. R., Schlegel, D. J., Strauss, M. A., et al. 2005, The Astronomical Journal, 129, 2562, doi: 10.1086/429803
  * Bloom et al. (2017) Bloom, J. V., Croom, S. M., Bryant, J. J., et al. 2017, Monthly Notices of the Royal Astronomical Society, 472, 1809, doi: 10.1093/mnras/stx1701
  * Boselli & Gavazzi (2006) Boselli, A., & Gavazzi, G. 2006, PASP, 118, 517, doi: 10.1086/500691
  * Brown et al. (2017) Brown, T., Catinella, B., Cortese, L., et al. 2017, MNRAS, 466, 1275, doi: 10.1093/mnras/stw2991
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * 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. (2013) Catinella, B., Schiminovich, D., Cortese, L., et al. 2013, MNRAS, 436, 34, doi: 10.1093/mnras/stt1417
  * Catinella et al. (2018) Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, Monthly Notices of the Royal Astronomical Society, 476, 875, doi: 10.1093/mnras/sty089
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Chauhan et al. (2021) Chauhan, G., Lagos, C. d. P., Stevens, A. R. H., et al. 2021, MNRAS, 506, 4893, doi: 10.1093/mnras/stab1925
  * Chauhan et al. (2020) —. 2020, MNRAS, 498, 44, doi: 10.1093/mnras/staa2251
  * Chemin et al. (2009) Chemin, L., Carignan, C., & Foster, T. 2009, ApJ, 705, 1395, doi: 10.1088/0004-637X/705/2/1395
  * Chen et al. (2019) Chen, Y., Mo, H. J., Li, C., et al. 2019, ApJ, 872, 180, doi: 10.3847/1538-4357/ab0208
  * Collette et al. (2021) Collette, A., Kluyver, T., Caswell, T. A., et al. 2021, h5py/h5py: 3.5.0, 3.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.5585380
  * Cowie & Songaila (1977) Cowie, L. L., & Songaila, A. 1977, Nature, 266, 501, doi: 10.1038/266501a0
  * Davé et al. (2020) Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146, doi: 10.1093/mnras/staa1894
  * Driver et al. (2011) Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971, doi: 10.1111/j.1365-2966.2010.18188.x
  * Driver et al. (2022) Driver, S. P., Robotham, A. S. G., Obreschkow, D., et al. 2022, MNRAS, 515, 2138, doi: 10.1093/mnras/stac581
  * Eckert et al. (2015) Eckert, K. D., Kannappan, S. J., Stark, D. V., et al. 2015, ApJ, 810, 166, doi: 10.1088/0004-637X/810/2/166
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * 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, The Astronomical Journal, 130, 2598, doi: 10.1086/497431
  * Gunn & Gott (1972) Gunn, J. E., & Gott, J. Richard, I. 1972, ApJ, 176, 1, doi: 10.1086/151605
  * 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. (2017) Guo, H., Li, C., Zheng, Z., et al. 2017, ApJ, 846, 61, doi: 10.3847/1538-4357/aa85e7
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Haynes et al. (2011) Haynes, M. P., Giovanelli, R., Martin, A. M., et al. 2011, The Astronomical Journal, 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
  * Holwerda et al. (2012) Holwerda, B. W., Blyth, S. L., & Baker, A. J. 2012, in The Spectral Energy Distribution of Galaxies - SED 2011, ed. R. J. Tuffs & C. C. Popescu, Vol. 284, 496–499, doi: 10.1017/S1743921312009702
  * Hunter et al. (2012) Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, AJ, 144, 134, doi: 10.1088/0004-6256/144/5/134
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jones et al. (2018) Jones, M. G., Haynes, M. P., Giovanelli, R., & Moorman, C. 2018, MNRAS, 477, 2, doi: 10.1093/mnras/sty521
  * Kalberla & Kerp (2009) Kalberla, P. M. W., & Kerp, J. 2009, ARA&A, 47, 27, doi: 10.1146/annurev-astro-082708-101823
  * Kannappan (2004) Kannappan, S. J. 2004, ApJ, 611, L89, doi: 10.1086/423785
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121, doi: 10.1086/321545
  * Kim et al. (2014) Kim, S., Rey, S.-C., Jerjen, H., et al. 2014, ApJS, 215, 22, doi: 10.1088/0067-0049/215/2/22
  * Koribalski et al. (2020) Koribalski, B. S., Staveley-Smith, L., Westmeier, T., et al. 2020, Ap&SS, 365, 118, doi: 10.1007/s10509-020-03831-4
  * Lan et al. (2016) Lan, T.-W., Ménard, B., & Mo, H. 2016, MNRAS, 459, 3998, doi: 10.1093/mnras/stw898
  * Lewis (2019) Lewis, A. 2019, arXiv e-prints, arXiv:1910.13970.  https://arxiv.org/abs/1910.13970
  * 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. (2006) Li, C., Kauffmann, G., Jing, Y. P., et al. 2006, MNRAS, 368, 21, doi: 10.1111/j.1365-2966.2006.10066.x
  * Li & White (2009) Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177, doi: 10.1111/j.1365-2966.2009.15268.x
  * Li et al. (2022) Li, Z., Guo, H., & Mao, Y. 2022, arXiv e-prints, arXiv:2207.10414.  https://arxiv.org/abs/2207.10414
  * 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
  * Loveday et al. (2015) Loveday, J., Norberg, P., Baldry, I. K., et al. 2015, MNRAS, 451, 1540, doi: 10.1093/mnras/stv1013
  * Lu et al. (2020) Lu, Y., Yang, X., Liu, C., et al. 2020, arXiv e-prints, arXiv:2008.09804.  https://arxiv.org/abs/2008.09804
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113, doi: 10.1093/mnras/sty2206
  * Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, The Astrophysical Journal, 619, L1, doi: 10.1086/426387
  * 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
  * Mo et al. (2010) Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206, doi: 10.1093/mnras/sty618
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624, doi: 10.1093/mnras/stx3040
  * Nulsen (1982) Nulsen, P. E. J. 1982, MNRAS, 198, 1007, doi: 10.1093/mnras/198.4.1007
  * Obuljen et al. (2019) Obuljen, A., Alonso, D., Villaescusa-Navarro, F., Yoon, I., & Jones, M. 2019, MNRAS, 486, 5124, doi: 10.1093/mnras/stz1118
  * Paul et al. (2018) Paul, N., Choudhury, T. R., & Paranjape, A. 2018, MNRAS, 479, 1627, doi: 10.1093/mnras/sty1539
  * Pillepich et al. (2018) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648, doi: 10.1093/mnras/stx3112
  * Praton & Schneider (1994) Praton, E. A., & Schneider, S. E. 1994, ApJ, 422, 46, doi: 10.1086/173702
  * Rafieferantsoa et al. (2018) Rafieferantsoa, M., Andrianomena, S., & Davé, R. 2018, MNRAS, 479, 4509, doi: 10.1093/mnras/sty1777
  * Ramatsoku et al. (2019) Ramatsoku, M., Serra, P., Poggianti, B. M., et al. 2019, MNRAS, 487, 4580, doi: 10.1093/mnras/stz1609
  * Roberts et al. (2021) Roberts, I. D., van Weeren, R. J., McGee, S. L., et al. 2021, A&A, 652, A153, doi: 10.1051/0004-6361/202141118
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297, doi: 10.1086/154079
  * Schmidt (1963) Schmidt, M. 1963, ApJ, 137, 758, doi: 10.1086/147553
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Teimoorinia et al. (2017) Teimoorinia, H., Ellison, S. L., & Patton, D. R. 2017, MNRAS, 464, 3796, doi: 10.1093/mnras/stw2606
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Genel, S., Castorina, E., et al. 2018, ApJ, 866, 135, doi: 10.3847/1538-4357/aadba0
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * 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
  * Wang et al. (2014) Wang, H., Mo, H. J., Yang, X., Jing, Y. P., & Lin, W. P. 2014, ApJ, 794, 94, doi: 10.1088/0004-637X/794/1/94
  * Wang et al. (2021) Wang, J., Staveley-Smith, L., Westmeier, T., et al. 2021, ApJ, 915, 70, doi: 10.3847/1538-4357/abfc52
  * Wang et al. (2022) Wang, S., Wang, J., For, B.-Q., et al. 2022, ApJ, 927, 66, doi: 10.3847/1538-4357/ac4270
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435, doi: 10.1146/annurev-astro-081817-051756
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341
  * Wong et al. (2006) Wong, O. I., Ryan-Weber, E. V., Garcia-Appadoo, D. A., et al. 2006, MNRAS, 371, 1855, doi: 10.1111/j.1365-2966.2006.10846.x
  * Yang et al. (2008) Yang, X., Mo, H. J., & van den Bosch, F. C. 2008, ApJ, 676, 248, doi: 10.1086/528954
  * Yang et al. (2009) —. 2009, ApJ, 695, 900, doi: 10.1088/0004-637X/695/2/900
  * Yang et al. (2007) Yang, X., Mo, H. J., van den Bosch, F. C., et al. 2007, ApJ, 671, 153, doi: 10.1086/522027
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Zhang et al. (2021) Zhang, W., Kauffmann, G., Wang, J., et al. 2021, A&A, 648, A25, doi: 10.1051/0004-6361/202039878
  * Zhang et al. (2013) Zhang, W., Li, C., Kauffmann, G., & Xiao, T. 2013, MNRAS, 429, 2191, doi: 10.1093/mnras/sts490
  * Zhang et al. (2009) Zhang, W., Li, C., Kauffmann, G., et al. 2009, MNRAS, 397, 1243, doi: 10.1111/j.1365-2966.2009.15050.x
  * Zu (2020) Zu, Y. 2020, MNRAS, 496, 111, doi: 10.1093/mnras/staa1457
  * 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