{'fraction_non_alphanumeric': 0.03999779322520137, 'fraction_numerical': 0.026885872963330757, 'mean_word_length': 4.242865406864635, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Population of Galaxies that Contribute to The HI Mass Function

We look at the contribution of different galaxy populations to the atomic hydrogen (HI) mass function (HIMF) and the HI density parameter, ΩHI, in the local Universe. Our analysis is based on a sample of 7857 HI-selected galaxies selected from a volume common to the SDSS and ALFALFA surveys (40% catalog – α.40). We define different populations of galaxies in the color(u-r)-magnitude(Mr) plane and compute the HIMF for each of them. Additionally we compute the HIMF for dark galaxies; these are undetected in SDSS and represent ∼ 2% of the total sample. We find that the bright red population dominates the total HIMF for log10(MHIh70/M ) ≥ 10.4. The full red population – bright and faint – represents about ∼ 17% of the ΩHI budget, while that of the dark population is ∼ 3%. The HIMF about the knee, log10(MHIh70/M ) ∈ [8, 10.4], is dominated by the faint and bright blue populations, the latter dominating at larger masses in this interval. Their total contribution to ΩHI is ∼ 55 − 70%, the variation depending on the definition of population. The dominant populations at the low mass end, log10(MHIh70/M ) ≤ 8.0 are the faint blue and faint bluer populations, the latter’s dominance being sensitive to its definition. The full blue (blue–bluer bright and faint) population represents ∼ 80% of ΩHI. A bimodal HIMF suggested by our results is however not seen since the amplitude of the HIMF of the bright red population is small compared to that of the bright blue population.


INTRODUCTION
The relationship between gas, metals, feedback and stars in galaxies is crucial for our understanding of galaxy formation and evolution. We need a clearer picture of how dark matter halos, which harbor galaxies, are supplied with cold gas, the fuel for star formation; how local and global conditions in the galaxy are responsible in processing them into stars; and finally how the gas is polluted with metals and recycled back to the intergalactic medium due to feedback processes within the galaxy.
Three basic properties that describe galaxies are: (i) the star formation rate (SFR), (ii) the stellar mass (Mstar) and (iii) the cold neutral hydrogen gas mass both in molecular (MH 2 ) and atomic (MHI) phases. The amount of neutral gas tells us the amount of fuel currently available for future star formation. The SFR is the current rate of forming stars from the supply of cold gas that is available, while Mstar depends on the integrated star formation history. Within each galaxy, at sub-kpc scales, the observed correlation between surface density of molecular hydrogen, ΣH 2 , and the SFR surface density, ΣSFR, is stronger (Bigiel, et al. 2008;Leroy, et al. 2008) as compared to the correlation between the HI surface density, ΣHI, and ΣSFR. The HI is often distributed beyond the optical radius of the galaxy and is more diffuse whereas star formation and H2 are mostly concentrated within the optical radius and occur in clumpier regions (Leroy, et al. 2008). However when the HI gas cools and becomes denser it transitions to molecular gas, which cools further and becomes denser leading eventually to gravitational collapse to form stars. It is therefore common to correlate the total gas (HI+H2) to the SFR, the so-called Kennicutt-Schmidt law (Schmidt 1959(Schmidt , 1963Kennicutt 1998Kennicutt , 1989 for star formation, where ΣSFR ∝ Σ 1.4 gas . With the aid of multiwavelength observations in the optical, ultraviolet (UV) and infrared (IR) bands followed up by spectroscopic measurements, we are able to infer scaling relations between various properties of galaxies, derived and/or observed. These scaling relations are very important since any theoretical model of galaxy formation should be able to reproduce them. In practice, these relations are used as parameters in theoretical models of galaxy formation when studied in cosmological volumes, since the physics on small, subparsec scales is never resolved (Springel & Hernquist 2003;Di Matteo, Springel & Hernquist 2005). A multipronged approach is used to study the distribution of HI in the post reionization Universe; these include cosmological hydrodynamical simulations of galaxy formation (Davé, et al. 2017), semi-analytical models (SAM) of galaxy formation (Kim, et al. 2017) and halo-occupation distribution (HOD) (Paul, Choudhury & Paranjape 2018). All of these approaches invoke in some form or other some of the observed scaling relations and model predictions are tested against observations which were not used as their input. Since these scaling relations depend on the observed sample they may be biased, i.e. the relations may depend on how the sample is chosen. E.g. the MHI − Mstar relation differs if a sample is chosen by HI mass (Huang, et al. 2012) or stellar mass (Catinella, et al. 2010). On the other hand quantities like the luminosity function, mass function, correlation function (to name a few) are corrected for the survey selection and tell us about the underlying abundance and distribution of different galaxy types in the survey volume.
By analyzing data from optical, UV and IR surveys over the past decade we have formed a clearer picture of how galaxies, on average, have formed and evolved over the past ∼ 12.5 billion years from redshift z = 6 to today. The cosmic stellar mass density, ρ * (units of M Mpc −3 3), has increased monotonically by nearly 2.5 decades from log(ρ * ) 6.3 at z = 6 to log(ρ * ) 8.8 today (Madau & Dickinson 2014). These observations also tell us how the cosmic SFR density (SFRD denoted by ψ with units M yr −1 Mpc −3 ) has changed during this time. It increases from log(ψ) = −1.7 at z = 6 to a peak value of log(ψ) = −0.9 at z 2 and finally dropping by a decade to log(ψ) = −1.8 today (Le Floc'h, et al. 2005;Madau & Dickinson 2014;Casey, Narayanan & Cooray 2014). In contrast, surveys targeting gas content of galaxies have lagged behind in depth and number. The 21cm line of HI being a weak line its detection in emission is limited only to the local Universe.
Blind HI surveys like the HI Parkes All Sky Survey (HIPASS, Meyer, et al. 2004) and the Arecibo Fast Legacy ALFA survey (ALFALFA, Giovanelli, et al. 2005) have been used to accurately measure the HIMF (Zwaan, et al. 2003(Zwaan, et al. , 2005Martin, et al. 2010;Haynes, et al. 2011;Jones, et al. 2018) in the local Universe (z ≤ 0.05). The HIMF can then be integrated to obtain the HI density parameter ΩHI. At higher redshifts the 21cm flux gets further diluted and direct detection becomes difficult with existing instruments. Stacking the HI data on known optical counterparts or alternately cross-correlating the HI intensity maps with the optical catalog in a common volume, then becomes a useful tool in making detections. The stacking method has been applied for star-forming galaxies at z = 0.24 (Lah, et al. 2007), galaxies in a cluster environment at z = 0.37 (Lah, et al. 2009), for field galaxies at z = 0.1 − 0.2 (Rhee, et al. 2013), in the zCOSMOS 1 field at z = 0.37 (Rhee, et al. 2016) and the VVDS 2 field at z = 0.32 (Rhee, et al. 2018). Most of the HI observations using the stacking method were 1 http://cesam.lam.fr/zCosmos/ 2 https://cesam.lam.fr/vvds/ done with the Giant Metrewave Radio Telescope (GMRT). More recently stacking has been applied to estimate the HI content in filaments (Tramonte, et al. 2019). However at higher redshifts (z ∼ 1.3) the stacking method did not result in a detection of star forming galaxies in the DEEP2 3 field (Kanekar, Sethi & Dwarakanath 2016) whereas positive detections have been reported by cross-correlating the HI Intensity map with an optical survey at z ∼ 0.8 (Chang, et al. 2010;Masui, et al. 2013). In the stacking method one can estimate ΩHI after correcting for the optical survey's completeness limit, whereas the cross-correlation method constrains ΩHIbHIr. Here bHI and r are the HI bias parameter and galaxy-HI cross correlation coefficient. In summary these HI surveys constrain ΩHI = 4 ± 1.6 × 10 −4 out to z 0.4 (Rhee, et al. 2018).
Beyond z = 0.4 and out to z 5 the HI content of the Universe is derived from studying damped Lyman-α systems (DLA) seen, in absorption, in the spectra of quasars (Prochaska, Herbert-Fort & Wolfe 2005;Noterdaeme, et al. 2012;Neeleman, et al. 2016;Rao, et al. 2017). Combining both these approaches at low redshift (emission) and high redshift (absorption) one sees that ΩHI increases monotonically as (1 + z) 0.56 from ΩHI(z = 0) = 4 × 10 −4 to ΩHI(z = 5) = 1.1 × 10 −3 (for a full compilation of all observational results see Rhee, et al. 2018). The HI density at z = 0 is only 1.5× smaller than at z = 2 with the data consistent with a no-evolution picture. On the other hand the cosmic SFRD has decreased 10 fold in this interval. Clearly our picture of galaxy formation is incomplete at these redshifts since the decrease in SFRD is not commensurate with a depletion of HI. Upcoming HI surveys like the Square Kilometre Array (SKA) will help us in this direction.
However even with the existing data it is important to understand the dependence of the abundance and distribution of HI selected galaxies on different galaxy properties and environments. The clustering of HI selected galaxies in the ALFALFA survey has been measured (Martin, et al. 2012;Papastergis, et al. 2013;Guo, et al. 2017) suggesting that ALFALFA galaxies cluster weakly. The dependence of HIMF has also been explored on the environment with Jones, et al. (2018) reporting a decrease in the low mass slope with increasing density in ALFALFA whereas Said, Kraan-Korteweg & Staveley-Smith (2019) find that the low mass slope increases with increasing density in the HI Zone of Avoidance survey with the Parkes telescope. Similarly the the HI velocity width function has been studied for wall and void galaxies as well as red and blue galaxies (Moorman, et al. 2014), indicating a strong dependence on both environment and galaxy color. Zwaan, et al. (2003) have looked at the dependence of the HIMF on galaxy morphology as well as on the early-late type classification. In this paper we look at the dependence of the HIMF on the different populations classified from the color-magnitude plane. This is similar in spirit to the analysis of Zwaan, et al. (2003); Moorman, et al. (2014) carried out in HI surveys and also bears resemblance to similar analysis carried out for optical surveys where the contribution of the red and blue galaxies to the galaxy stellar mass function and luminosity functions are explored (Baldry, et al. 2004;Drory, et al. 2009).
Our paper is organized as follows. In section 2 we describe a sample common to ALFALFA and SDSS (Sloan Digital Sky Survey), in section 2 we report our measurements of the HIMF for different populations within this sample, in section 4 we discuss our results and their implication and finally we summarize and conclude in section 5. We assume a flat ΛCDM cosmology with Ωm = 0.3, ΩΛ = 0.7 and a value for the dimensionless Hubble constant h70 = 0.7.
This catalog has been listed with the following observed quantities (i) an unique entry number from the Arecibo General Catalog(AGC); (ii) right ascension and declination of the HI source and the most probable optical counterpart (OC) from SDSS DR7; (iii) heliocentric velocity (cz helio ) which is the midpoint of the HI flux density profile; (iv) velocity profile width (W50), measured as the full width at 50% of the peak HI flux density; (v) the integrated flux density of the HI source (S21). The derived quantities, contained in this catalog, are -(i) distance to the object in Mpc (D). For sources with cz helio > 6000 km s −1 this quantity is cz cmb /H0, where cz cmb is the velocity in cosmic microwave background (cmb) reference frame and H0 is the Hubble constant (taken to be 70 km s −1 Mpc −1 ). And sources with cz helio < 6000 km s −1 this quantity has been estimated using a local flow model (Masters 2005) The α.40 catalog also provides another property which is the Code number, with a value of 1,2 or 9. The sources with SNR > 6.5 are referred as Code 1 objects, Code 2 objects are the detections with SNR < 6.5 and Code 9 refers to the high velocity clouds (HVC). For our analysis we have considered only Code 1 galaxies. The number of Code 1 galaxies is 11941. We have also considered a cut for cz cmb < 15000km.s −1 to avoid radio frequency interference (RFI) generated by Federal Aviation Administration (FAA) radar at the San Juan airport (Martin, et al. 2010;Haynes, et al. 2011). This restricts the sample to redshift, z ≤ 0.05 and reduces the sample to 10785 galaxies.
Since ALFALFA is not totally overlapping with the SDSS footprint, we have defined a common boundary for both the surveys in this work in figure 1. For this work our final area of analysis has four subregions whose vertices right ascension-declination (RA,dec) are given as (  (17,16), (-6,16); and (iv) (-2.6,24), (13.6,24), (13.6,26), (-2.6,26). The total volume of these regions is 2.02 × 10 6 Mpc 3 and it corresponds to an angular area of ∼ 2093 deg 2 . In this region there are 8344 galaxies in ALFALFA. This includes removal of 4 OH high-redshift impostors (Suess et al. 2016).
In fig. 2 we look at the distribution of Code 1 objects in ALFALFA in the S21−W50 plane. We also show the distributions in three mass bins which correspond to the small-mass end or faint-end (thick solid line), the knee (dot-dashed line) and the high-mass end (thin-solid line). As seen from these distributions on the average velocity width increases with increasing mass. This is expected and is also seen in the MHI − W50 relation (see fig. 7 of Moorman, et al. 2014). One also sees that it is very unlikely to have a low(high) mass galaxy with a large(small) velocity width. However intermediate mass objects can have the full range of velocity widths.
The broken solid line in fig.2 is the sensitivity limit and is given by a 50% completeness relation of Code 1 objects in eq. 2 (Haynes, et al. 2011). This tells us that the detection of objects not only depends on the integrated flux, but also on the observed velocity width. At fixed S21 the detection is more likely for narrower HI profile widths.
log S21 = 0.5 log W50 − 1.207 : log W50 < 2.5 log W50 − 2.457 : log W50 ≥ 2.5 (2) After applying this completeness cut we are left with a sample of 7857 galaxies. Among these, 6076 galaxies have spectroscopic as well as photometric measurements in SDSS DR7 and 1633 galaxies have only photometric measurements. As for the remaining 148, we loosely refer to them as dark galaxies. These objects are not being identified in Figure 2. The distribution of galaxies in the S 21 − W 50 plane is shown for all Code 1 ALFALFA galaxies. The broken solid line is the 50% completeness relation (Haynes, et al. 2011). It is given in eq. 2. We have further divided the sample into three mass bins: log 10 [M HI /M ] ∈ [6. 0, 8.5[, [8.5, 10.2[, [10.2, 11.0[. The 1 − σ contours and the peaks of the distributions for these three populations are given by thick solid line (plus-circle), the dotdashed line (dot-circle) and the thin solid line (cross-circle).
the SDSS pipeline as potential galaxies. Although follow up observations have been made on some of the dark galaxies (Cannon, et al. 2015;Janowieki, et al. 2015;Leisman, et al. 2017), in this work we refer to all of them as dark. This translates to around ∼ 2% of galaxies which are dark or have no optical counterparts in SDSS.
The ALFALFA catalog also lists the SDSS objectIDS of the OC. Using these we extracted the photometric properties like the ugriz values (model magnitudes), then corrected for extinction (due to our own galaxy) for these 7709 (these are non-dark galaxies) galaxies. We also kcorrect (Blanton & Roweis 2007) the magnitudes to finally obtain rest frame magnitudes. For objects which do not have spectroscopic redshifts we have supplied the HI redshifts to kcorrect them. The kcorrect code also estimates additional properties like galaxy stellar mass, integrated star formation history and metallicity for these objects.
As of writing this paper the ALFALFA team has released the 100% catalog (α.100) ) which also include the RA and dec of the optical counterparts. However we find that there are many galaxies which have bright foreground stars due to which SDSS has masked the region covering the galaxy and photometric values are not provided. We therefore restrict ourselves with α.40 catalog and will revisit the α.100 sample in the future.

Subsamples and Populations of Galaxies
From our sample of 7857 galaxies we identify subsamples which define different populations of galaxies. These populations are disjoint sets and their union (including the dark galaxies) forms the full sample. The populations are based on dividing the color(u − r)-magnitude(Mr) plane of the HI selected galaxies into six disjoint regions. This is shown in figure 3.
In figure 3 the upper solid curve demarcates the red (above curve) from the blue (below curve) population as in Baldry, et al. (2004). This optimal divider is given by The vertical solid line divides the bright (leftward of line) and faint populations. The line has been chosen so that the fraction of bright red objects over all red objects is 0.87 which we refer to as the 1.5σ line. Similarly the lower solid curve is chosen so that the fraction of blue (above curve) galaxies over all blue galaxies is 0.87 (or the 1.5σ cut in color). This curve has been chosen to be parallel to the curve which demarcates the red and blue populations. We refer to the objects below(above) the curve as bluer(blue) galaxies. Similarly we define the 1σ sample (dashed lines). This breaks the sample of HI selected galaxies (which have optical counterparts) into 6 disjoint sets in the color-magnitude plane; we call each set a single population. The number for each population is quoted. The numbers in brackets are for the 1σ sample. The data points are: i. filled (open) triangles for bright (faint) red galaxies ii. filled (open) squares for bright (faint) blue galaxies iii. filled (open) circles for bright (faint) bluer galaxies and have been marked for the 1.5σ sample. Our definition of a population of galaxies is broadly a subsample of the total sample of 7857 HI selected galaxies. The dark population has no optical counterparts. The bright(faint) population is a subsample of galaxies with absolute magnitudes M i r , such that M i r < Mr (M i r ≥ Mr) for some reference value Mr (the vertical line in figure 3). Similarly we have defined three populations based on color (red, blue and bluer), based on equation 3. The six populations are formed by further splitting these three into bright and faint populations. As discussed earlier, each population is dependent on the boundaries that define them which we have called a 1σ or 1.5σ line. We caution the reader that our definition of sample may be confusing since both these samples represent the same set of 7709 galaxies which have optical counterparts. However the number of galaxies for any given population is different for the 1σ and 1.5σ samples.
We have chosen to define two samples to demonstrate that the qualitative results do not depend very strongly on sample definition. However as we will show in the next section, the definition of the faint bluer population (i.e. the 1.5σ or 1σ sample) determines whether it contributes significantly or not to some part of the HIMF. Looking at how the observed counts change, both in relative and absolute terms, when going from the 1.5σ to the 1σ sample we expect that the HIMF of the bright blue population to be the affected the most. In section 4 (figure 6) we will discuss the break in the MHI − Mstar, at Mstar 9, where Mstar is the logarithm of the galaxy stellar mass in units of M /h 2 70 . One can then use the mean Mr − Mstar relation to convert Mstar = 9 to Mr = −19. This value is sandwiched between the vertical lines (solid and dashed) in figure 3. . The distribution of HI selected galaxies in this work in the color-magnitude plane. The upper solid curve demarcates the red (above curve) from the blue (below curve) population as in Baldry, et al. (2004). The vertical solid (dashed) line divides the bright and faint populations. The line has been chosen so that the fraction of bright red objects over all red objects is 0.87 which we refer as the 1.5σ line. Similarly the lower solid curve is chosen so that the fraction of blue (above curve) galaxies over all blue galaxies is 0.87 (or the 1.5σ cut in color). This curve has been chosen to be parallel to the curve which demarcates the red and blue populations. We refer to the objects below(above) the curve as bluer(blue) galaxies. Similarly we define the 1σ sample with the help of dashed lines. This breaks the sample of HI selected galaxies (which have optical counterparts) into 6 disjoint sets in the color-magnitude plane. The number for each population is quoted. The numbers in brackets are for the 1σ sample.

RESULTS
The mass function φ(M ) is defined as the number density of objects in the mass range [M, M + dM ]. The HI mass function φ(MHI) can be expressed as where, dN is the total number of galaxies having HI mass between MHI and MHI+dMHI and V is the survey volume of interest. The HI mass function can be parameterized as a Schechter Function Here, α is the faint end slope, φ * is the amplitude and M * is the characteristic HI mass. A simple and intuitive way of calculating the HIMF is by the 1/Vmax method (Schmidt 1968). The underlying assumption in this method is that the sample of galaxies detected by a survey is a representative sample of galaxies in the Universe (in the same redshift range), or in other words it assumes homogeneity. For each detected galaxy 'i' a maximum detectable distance D i max is calculated based on inverting eq. 1 by using the mass M i HI of the galaxy and the limiting flux S i lim which is the property of the survey. In the case of ALFALFA the completeness relation (eq. 2) determines this limiting flux for the velocity width W i 50 of the galaxy. D i max is then converted to a volume V i max which is the volume in which the galaxy could be in, and still be detected by the survey. Finally the galaxies are binned in mass with relative weights 1/V i max to obtain the mass function. For the brighter or more massive objects V i max is larger than the survey volume. In such a case the relative weight is set to unity. This method has the advantage that it is non-parametric and does not require any prior knowledge to estimate the HIMF. However since galaxies cluster, the estimate of the HIMF will be sensitive to large scale structure in the local volume. Additionally some volumes of the survey may be inaccessible due to RFI and needs further correction. An estimate of these effects can be used to recalibrate the weights V i max (Martin, et al. 2010). Maximum likelihood (Sandage, Tammann, & Yahil 1979) and Step-Wise Maximum Likelihood (SWML) (Efstathiou, Ellis, & Peterson 1988) methods on the other hand are designed to be insensitive to large-scale structure. In the former the assumption is that the galaxy sample is drawn from an underlying distribution function, e.g. the HIMF for this work φ(MHI), and the likelihood method determines the parameters of this function. Although in most cases a Schechter function is the chosen function one has no way of testing whether it is the optimal function to describe data (Efstathiou, Ellis, & Peterson 1988). Baldry, et al. (2012) and Drory, et al. (2009) find that a single Schechter function does not describe the galaxy stellar mass function. In the latter, i.e. SWML, the underlying distribution φ(MHI) has no functional form but is discretized in steps or bins of mass and a uniform distribution is assumed in each bin. Hence the value φj which is the value of φ in the j th mass bin becomes the parameter. The joint likelihood of detecting all galaxies in the sample is maximized, with respect to the parameters φj, thus determining their values. This method works when the sample is flux-limited.
When the selection function depends on other properties of the galaxies one needs to consider an underlying bivariate or multivariate distribution for φ. One has to then generalize the SWML method to higher dimensions. Loveday (2000) estimated the bivariate luminosity function φ(MK, MB) and then marginalized over MB to obtain the K-Band luminosity function, φ(MK), starting with a bJselected sample in the Stromlo-APM Redshift Survey.
For a blind HI survey like ALFALFA the limiting flux (figure 2 and eq. 2) depends on the velocity width W50. A two-dimensional SWML (2DSWML) method similar to Loveday (2000) was applied by Zwaan, et al. (2003) to estimate the HIMF for HIPASS Galaxies. The bivariate distribution in this case is φ(MHI, W50) which can be marginalized over W50 to obtain the HIMF (Zwaan, et al. 2003;Martin, et al. 2010;Haynes, et al. 2011;Jones, et al. 2018) or marginalized over MHI to obtain the HI velocity width function (Zwaan, Meyer, & Staveley-Smith 2010;Moorman, et al. 2014). The details of our implementation are given in appendix A.
In figure 4 we show our estimate of the HIMF for all the populations including the total (crosses) and the dark (filled diamonds) populations. The columns are for the 1.5σ (left) and the 1σ (right) populations. We have broken our results into a dominant (top row) and a subdominant (bottom row) population to better illustrate our results. We call a dominant population one which dominates the HIMF over the rest of the populations in some mass range and also contributes greater than 10% to ΩHI (see table 2). The symbols for the six populations are the same as in figure 3. The curves are Schechter function fits for the total (solid), red (dotted), blue (dot-dashed), bluer (dashed) and dark (dot-dot-dotdashed) populations. Our Schechter function fits are summarized in table 1. In the rest of the paper we will quote the values of the characteristic mass M * and the amplitude of the Schechter function φ * in the units log(M * /M ) + 2 log h70 and (10 −3 h 3 70 M pc −3 dex −1 ) respectively. We will also quote the values of MHI in the same units as M * . In practice we bin the mass and the velocity width in logarithmic bins therefore the faint slope in figure 4 (α ) differs from α in equation 5 by 1, or α = α + 1.
The goodness of fits, χ 2 reduced , is given in the last column of table 1. For the dark, faint red, faint bluer populations the χ 2 reduced is of order unity. In the 1σ sample of the faint red population the biggest contribution of χ 2 reduced comes from the lowest mass bin. If we remove that point as an outlier then χ 2 reduced = 0.57 and α flattens to −1.06 from −1.16, φ * increases by about 22% from 0.31. The change in M * is negligible. The bright red and bright bluer populations have a low χ 2 reduced , irrespective of sample definition, and looking at the data points relative to their fitted curve we find that there is little variation between them. This means that the error bars are larger than the variation between data and model. The bright blue and faint blue populations on the other hand have larger χ 2 reduced . For the 1.5σ sample of the faint blue population χ 2 reduced = 4.10 and improves to χ 2 reduced = 0.96 for the 1σ sample. The 1σ sample for the faint blue population has a larger number of brighter galaxies as compared to the 1.5σ sample (figure 3) and because MHI is correlated to MR, (see figure 5) the high mass end of the HIMF is better represented in the 1σ sample, leading to a smoother change in data and a better fit. This can also be seen in the top row of figure 4. The bright blue population has χ 2 reduced ∼ 2 irrespective of sample definition. We start by looking at the bright population. Due to the monotonic relation between MHI and MR (figure 5), across populations, we expect the bright population to dominate the massive end of the HIMF and be subdominant at the low mass end. This is seen in figure 4. The characteristic mass, M * increases systematically from the bright bluer to the bright red population. There is little change in M * for both the bright red and bright blue populations with respect to sample definition. We also see little change in the HIMF with respect to sample definition for MHI ≥ 10.3 for both these populations. The bright bluer population on the other hand has M * = 9.84 for the 1.5σ sample and reduces to M * = 9.52 for the 1σ sample. At the low mass end the bright populations have shallower slopes α + 1 ≥ 0. This is expected as mentioned earlier since at the low mass end we expect the faint population to dominate. At this end, the sample definition affects the bright bluer population the most, where α = −0.92 (1.5σ sample) and increases to α = +0.86 (1σ sample), the change being the smallest for the bright red population. On the other hand the amplitude is most affected for the bright blue population, it changes from φ * = 4.85 (1.5σ sample) to φ * = 2.54 (1σ sample). The change is negligible for the bright red and about ∼ 72% for the bright bluer population. We point out that changes in the HIMF with respect to sample definition can For the 1.5σ sample the bright blue is the dominant population at the knee of the HIMF, however it contributes nearly equally with the faint blue population when we consider the 1σ sample. This is because the faint blue population has a net increase in observed galaxies from n gal = 1409 (1.5σ sample) to n gal = 2422 (1σ sample), the net change coming from the intersection of bright blue (1.5σ) and faint blue (1σ) (see figure 3); this is the primary reason for the increase of φ * by ∼ 4.8×. An interesting result is that the bright red population is the dominant population at MHI ≥ 10.4. This result is insensitive to sample definition. We find that for MHI ≥ 10.4 the bright red population represents ∼ 60% of total detections with ∼ 40% coming from the bright blue population which also translates to similar fractions in total HI mass at this end.
We now move to the faint population. All of them have steeper slopes as compared to their bright counterparts and do not dominate the HIMF at the high mass end. The faint red population is the most subdominant population. The observed counts of galaxies of the faint bluer population see the largest relative change with sample definition increasing from n gal = 452 (1.5σ sample) to n gal = 1369 (1σ sample), a factor ∼ 3×. This results in a small change in slope from α = −1.76 to α = −1.62 but a large, ∼ 10×, increase in amplitude from φ * = 0.10 to φ * = 1.09. In the 1.5σ sample the faint bluer population is a subdominant population, however it becomes the dominant population for masses MHI ≤ 8.1. The faint blue population is on the other hand a dominant population below the knee of the mass function. It dominates the HIMF for MHI ≤ 8.7 in the 1.5σ sample and for 8.1 ≤ MHI ≤ 9.7 in the 1σ sample. The dark population is characterized by a very steep slope α = −1.92, large characteristic mass M * = 10.03 and a very small amplitude φ * = 3.25 × 10 −2 and is a subdominant population. Extrapolating to masses below MHI ≤ 6 our results suggest the dark population will be the dominant population. However it is unclear how far down we can extrapolate since it is unlikely that there will be too many low mass, gas rich galaxies devoid of stars which will be able to self shield themselves from the photoionizing background.

The contribution of different galaxy populations to Ω HI
We can analytically integrate the HIMF, fitted to a Schechter function to obtain the cosmic HI density parameter, Alternately we can sum the binned measurements of the HIMF. Similar to Haynes, et al. (2011) we find that both methods give similar results, with the exception of the dark sample which has a very steep slope. As we argued in the previous section it is not physical to extrapolate the HIMF to very small masses. Hence we choose to quote our results by the summed method. We note that the results do not change if we integrate the Schechter function from MHI = 6.1 to ∞. We summarize our results in table 2. Column 2 is the estimate of ΩHI from each population and column 3 is the fractional contribution to Ω tot HI , the values in brackets are for the 1σ sample.
For the total sample we obtain ΩHI = (4.894 ± 0.469) × 10 −4 h −1 70 which is consistent at 1σ with ΩHI = (4.4 ± 0.1) × 10 −4 h −1 70 (summed) of Martin, et al. (2010). With the addition of Code 2 objects in the α.40 sample we see that it is only consistent at the 2σ level with Haynes, et al. (2011) who obtain ΩHI = (4.2 ± 0.1) × 10 −4 h −1 70 (summed). Since our M * is comparable with the α.40 results, and α is only a bit steeper, the main reason for this discrepancy can be traced to φ * (equation 6). Our value of φ * is ∼ 10% higher than Martin, et al. (2010) which translates to a 10% higher estimate of ΩHI at fixed M * and α. However the relative ratios should not be sensitive to this change.
From table 2 we see that the red population (bright and faint) have a non-negligible contribution of ∼ 17% to Ω tot HI . When combined with the dark population (3% of Ω tot HI ), this adds up to a non-negligible fraction of 20% (rounded). The full blue population (faint and bright blue and bluer) then contributes ∼ 80% (rounded) of Ω tot HI . We will discuss the implications of these numbers in the next section. The dominant sample about the knee of the HIMF are the bright blue and the faint blue populations (as M * galaxies) and together they contribute most to Ω tot HI , ∼ 73% (1.5σ sample) and ∼ 55% (1σ sample)

DISCUSSION
Since we are looking at the contribution of different galaxy populations to the total HIMF we would like to see whether these populations have different scaling relations, e.g. in the MHI − Mstar plane. Such relations have been explored for galaxies in the ALFALFA sample (Catinella, et al. 2010;Huang, et al. 2012;Maddox, et al. 2015) and the HI Parkes All-Sky Survey Catalog (HICAT) (Parkash, et al. 2018). In this work the stellar masses are estimated by kcorrect which uses the population synthesis code of Bruzual & Charlot (2003). Our estimates on Mstar differ from the estimates of these authors. The MHI − Mstar scaling relations are shown in figure 6.
In order to avoid crowding figure 6 we choose to compare our results for the total sample with Huang, et al. (2012) (crossed-circle, thin solid line) only. The 1σ scatter on the data points (total) is ∼ 0.5 dex. For the total sample we find that our results compare well (within the scatter) with Huang, et al. (2012) in figure 6. The scaling relations that we find are also consistent with Maddox, et al. (2015);Parkash, et al. (2018) (MHI selected sample) but differ from Catinella, et al. (2010) (Mstar selected sample). However our stellar masses are underestimated at lower masses. This difference can be attributed to the choice of sample but more so due to attenuation by dust, affecting the redder sample, which these authors have considered. In this work we have not attempted to correct for reddening due to dust while Huang, et al. (2012) have used the additional two UV bands in GALEX to correct for it. Not correcting for it should therefore change the average scaling relations. This is also evident when looking at the scaling relations for the three faint populations. The faint blue and faint bluer populations have similar slopes but these are steeper compared to their corresponding bright populations. The faint red population, on the other hand, has a shallower slope with respect to the bright red population as well as the faint blue and bluer populations. We also find that the HI fraction, fHI = MHI/Mstar, increases with decreasing (u-r ) color.
We see a clear transition in the scaling relations when going from the low mass to the high mass end. fHI gets suppressed for the total sample at about Mstar ∼ 9 consistent with Huang, et al. (2012) ;Maddox, et al. (2015). The transition scale also depends on the galaxy population. For the blue, bluer and red populations it occurs at Mstar = 10.1, 9.4, 8 respectively. The transition scale of Mstar ∼ 9 corresponds to a change in the dominant morphology of galaxy populations (Maddox, et al. 2015) and also a transition between hot and cold mode accretion seen in cosmological hydrodynamical simulations (Kereš, et al. 2009).
One interesting result that we have quantified in the last section is the non-negligible HI content of red galaxies. The red galaxies dominate the HIMF at the high mass end MHI ≥ 10.4 and ∼ 17% of the HI content, ΩHI, is locked up in them. Using the HOD framework for HI, Paul, Choudhury & Paranjape (2018) also find that the red population is the dominant population at higher masses. Since the AL-FALFA sample is an HI selected sample, with the majority of the galaxies belonging to the blue cloud, one may ask: why do the rarest, gas rich galaxies, predominantly belong to the red cloud? Looking at the morphology of these gas rich red galaxies we find that these are predominantly spirals and lenticular galaxies, but there also exist some elliptical galaxies. A number of spirals have prominent bulges which would classify them as early-type spirals and there are a number of galaxies which harbor dust streaks visible on their disk plane. A significant number of galaxies are edge-on or somewhere in between edge-on and face-on. Indeed HI has been detected in early type galaxies (Morganti, et al. 2006;Oosterloo, et al. 2007;Grossi, et al. 2009), but these do not go beyond MHI ∼ 10. Schiminovich, et al. (2010) find that the 47% of the total local SFR density is found for Mstar > 10 in the GASS sample. Given that there is little correlation between (u-r ) color and Mstar at these masses and the fact that GASS (and ALFALFA) detections are predominantly in the blue cloud (Catinella, et al. 2010) we would expect these red gas rich galaxies to contribute a negligible fraction to the local SFR density. This does not mean that individually all bright red, gas rich galaxies have low star formation rates but, rather, their numbers are so small that their total contribution is negligible. A fraction of these red galaxies would then be the ones transitioning from a blue star-forming phase to a red passive phase with little star formation and another fraction will be dusty star forming galaxies, while the rest would be passively evolving. The amount of reddening would be enhanced if they are edge-on and would redden the color of disky galaxies which are either transitioning to the red phase or contain considerable dust on their disks (Graham & Worley 2008;Tempel, et al. 2011).
Although the bright red galaxies are the dominant population at the high mass end (MHI ≥ 10.4) the observed counts (60% of total) are only 50% more than those of the bright blue (40% of total) population. When plotted on a logarithmic scale the differences between these two mass functions are not very large (see top row of figure 4). If the inclination and reddening are important and the magnitudes are corrected for them, then a fraction of gas rich red galaxies would move to the bright blue population and would bring the mass functions of these two populations closer to each other at the high mass end. Conversely if the reddening is increased due to inclination we would expect that the effects would be more dramatic in HI velocity width function. This is indeed seen for ALFALFA galaxies (Moorman, et al. 2014), where the HI velocity width function for the red and blue are well segregated at the high velocity end. We therefore believe that reddening due to dust and inclination can partially explain why the red sample is the dominant population at the high mass end of the HIMF.
The results of section 3 are essentially conditional HIMF integrated over a range in color and magnitude, which we have called the HIMF for different populations. We can repeat the exercise and compute the HIMF in finer intervals of color and magnitude to obtain a conditional HIMF (conditioned on luminosity and color). This will then tell us about the distribution of HI in the color-magnitude plane. As an application one can then make better estimates of the corrections applied to ΩHI with the stacking methods (Rhee, et al. 2013(Rhee, et al. , 2016(Rhee, et al. , 2018 at higher redshifts. A second application would be to inform a proposed HI survey, which galaxy populations to stack on to make a tentative detection. In the survey volume considered in this work we find that only 12% of the red population in SDSS have HI detections in ALFALFA. This number is 38% for blue galaxies. On the other hand 98% of ALFALFA galaxies have optical counterparts. The detections are due to a combination of total HI signal and observed HI velocity widths. The question we wish to ask is are the non-detections in the red cloud due to insufficient HI gas or due to large velocity widths or both? We argue that in the bright red sample the nondetections are due to insufficient HI gas as well as large widths and HI masses should decrease with either increasing stellar mass or halo mass. Although figure 6 suggests that the most massive galaxies (large stellar mass, Mstar) are also the richest in terms of their gas content, this relation is biased since it is from an HI selected sample. The appropriate sample is the GASS sample which is selected on stellar mass. We look at the MHI − Mstar relation in the final data release of GASS (DR5), which is summarized on table 1 of Catinella, et al. (2013). The relation has nearly a flat slope for Mstar ∈ [10.76, 11.30], being slightly positive if all non-detections have been assigned the limiting HI mass and slightly negative if all non-detections are assigned zero HI mass. The last bin however has only a handful of objects which are dominated by non-detections. However the GASS and ALFALFA catalogs are relatively shallow as compared to optical surveys like the SDSS and would miss a significant number of massive galaxies. The tail of the stellar mass function is dominated by the red galaxies and at Mstar = 11.3 the number of red galaxies is ∼ 10× the blue galaxies (Baldry, et al. 2012). These galaxies are probably central red galaxies (Drory, et al. 2009) and would be in halos of log 10 (M halo /M ) 14 − 14.5 (Behroozi, Conroy & Wechsler 2010) with virial temperatures Tvir ∼ 10 7 K and circular velocities Vcirc ∼ 6 − 9 × 10 2 km.s −1 . Most of these galaxies will then be the central galaxies of large groups of galaxies or clusters of galaxies. Given the large virial temperatures and large circular velocities it would be very unlikely to detect a considerable amount of neutral hydrogen in these systems. It is then very likely that the MHI − Mstar relation will not asymptote to a constant as indicated in Catinella, et al. (2013) but rather decrease with increasing stellar mass. This is suggested in the results of Kim, et al. (2017).
If the average MHI − Mstar becomes a non-monotonic function of stellar mass and therefore halo mass, HI abundance matching techniques, used to obtain MHI −M halo relation (Khandai, et al. 2011;Padmanabhan & Kulkarni 2017), will break down. The HI HOD models which also assume a step-like function (with the help of the error function) (Guo, et al. 2017;Paul, Choudhury & Paranjape 2018) for the average occupation of centrals, may need to be revised. A log normal form for the mean occupation function for centrals was compared to the step-like parameterization in the context of describing quasar clustering (Shen, et al. 2013), but it was found that the HOD parameters were not well constrained. Only more direct observations will shed light on the HI content of these massive galaxies and hopefully provide better inputs for the HOD parameterization.

SUMMARY
In this work we have measured the HIMF of different galaxy populations picked from the color-magnitude plane. The galaxies considered were from a local volume common to ALFALFA and the SDSS surveys in the redshift range z ∈ [0, 0.05]. After putting the relevant cuts in quality of detection, volume and completeness the final sample analyzed consists of 7857 galaxies. We divided the total sample first into bright and faint populations (cut in magnitude) and these were further split into three colors: red, blue and bluer. This forms a disjoint set of six populations. A seventh population which we call a dark population is one which does not have any optical counterparts in SDSS but has a detection in ALFALFA. The union of these seven populations is the total sample of 7857 galaxies. We have further considered a second sample which redefines the six populations by shifting the boundaries defining them. We have called them the 1σ and 1.5σ samples (section 2.1). The reason for doing this is to illustrate that our sample definition does not change the qualitative results that we report. The reason for splitting the full blue cloud into four (faint/bright for blue/bluer) instead of just two (faint/bright for all blue) was because we wanted to see the systematic effect of the tail of the blue population (especially the faint end) on our results. We indeed find that based on how we define our sample the faint bluer population becomes the dominant population at the low mass end of the HIMF (figure 4 and tables 1&2). We summarize our results below: • At fixed color (e.g. red, blue or bluer) the HIMF of the bright population dominates over their corresponding faint counterparts at the knee and the large mass end whereas the faint populations dominate at the low mass end.
• At fixed magnitude (e.g. bright or faint) there is no systematic trend at the low mass end and the knee of the HIMF, with decreasing color, for the faint population. However for the bright population, we see that the HIMF at the high mass end increases with increasing color. The bright red population is the dominant population at this end.
• The bright red sample dominates the total HIMF at MHI ≥ 10.4. When combined with the faint red sample it contributes about ∼ 16 − 17% of the ΩHI budget. The dark population contributes ∼ 3% to ΩHI. This has implications for detections done with stacking at higher redshifts which would target the blue star forming cloud for a first detection.
• The total blue cloud (blue and bluer) represents about ∼ 80% of the ΩHI budget.
• In the mass range about the knee, MHI ∈ [8, 10.4], the dominant populations are the faint and bright blue populations with the latter dominating at larger masses in this interval. Their total contribution to ΩHI is ∼ 55 − 70% depending on sample definitions.
• The dominant populations contributing to the low mass end of the HIMF the faint blue and faint bluer populations, the latter being dominant only for the 1σ sample definition.
• The bright bluer and faint red populations are subdominant populations contributing a total of ∼ 6 − 10% to ΩHI.
It would be interesting to see if simulations (Davé, et al. 2017), SAMs (Kim, et al. 2017) and HOD models are able to reproduce our results which provide additional constraints on the population of HI selected galaxies. In section 4 we argued that the effect of dust and inclination are responsible for the bright red population dominating the total HIMF at large masses. This was based on the results of Catinella, et al. (2013); Moorman, et al. (2014) but needs to be confirmed with a more detailed followup. In a forthcoming paper we will report on a more detailed analysis of the HI velocity width function. We are also working on measuring the conditional (conditioned on color or magnitude or both) HIMF which will put additional constraints on the properties of gas bearing galaxies. [h] No.
Conditions that the completeness relation eq. 2 has a change in slope at log 10 W50/ km.s −1 = 2.5, which coincides with the edge of the bin in W50 for our choice. Table A1 assumes this so that no further cases are considered.
In the 2DSWML one wishes to find φ jk that maximizes the joint probability or likelihood of finding all the galaxies in the sample simultaneously Using eq. A1 the joint likelihood is To obtain φ jk , we maximize rather the log-likelihood This gives us where, n jk = i V ijk is the number of galaxies in mass bin j and profile width bin k. φ jk are iteratively determined from eq. A8. To start the iteration we set the initial value of on the RHS of eq. A8 to be φ in mn = n jk / [Vsurv∆M ∆W ]. We set a relative tolerance for convergence of 10 −3 for all φ jk . We find that φ jk converges within 20 iterations. Finally the HIMF is obtained by summing over the velocity width bins k.

A1 Normalization of HIMF
One drawback for likelihood methods, as opposed to the 1/Vmax method, is that the normalization of the HIMF is not fixed. This is obvious from eq. A1. There are a number of ways to fix the normalization (Davis & Huchra 1982;Willmer 1997;Zwaan, et al. 2003;Martin, et al. 2010) which involve computing the selection function. Here we try a slightly different method. We assume that the high mass end of the HIMF is not affected by the selection function. One can test this assumption by comparing the ratio of normalized mass function from 2DSWML to that of the observed mass function which is related to the observed counts. In the top panel of figure A2 the observed HIMF is shown as filled pentagons (solid line) and the un-normalized and converged HIMF from the 2DSWML is shown as filled triangles (dashed line). This is done for the α.100 sample (Haynes, et al. 2018). The ratio of these is shown for the last 7 mass bins in the lower panel of figure A2. If the selection function affects the high mass end, the ratio at this end would not have a flat feature. We indeed find that in this particular example of α.100 the last 3 mass bins are unaffected at the level ranging from 0−0.003%; whereas the last but third bin is relatively suppressed by around 0.4%. For this method to work Figure A2. Upper panel: Red dashed line is the mass function estimated using 2DSWML method, which is not normalized. The black solid line is the number density per mass bin calculated for the same sample of HI galaxies. Lower panel: The blue solid line shows the ratio of the un-normalized mass function and the number density, multiplied by 100.
we need to test the flatness of this ratio which means that at least the last two points at the high mass end should be unaffected by the selection function. Starting from the high mass end we search sequentially at smaller bins which are unaffected by the selection function. The search is stopped when the selection function affects the particular bin at the level of 0.1% or greater. Finally we fix the normalization by matching the integrated counts over these bins to that of the observed HIMF.
We have compared this method to the one which normalizes the mass function to match the average observed counts as in Davis & Huchra (1982) ;Martin, et al. (2010). We find that they match at the relative level of ∼ 0.4%. Finally we compare in figure A3 our result with that of α.100 HIMF (Jones, et al. 2018). As we can see our implementation with some minor modifications reproduces the HIMF of ) extremely well.

A2 Error Analysis of HIMF
(i) Mass Errors: Since M HI ∝ S 21 D 2 the uncertainties on both integrated flux and distances lead to uncertainties on mass. Peculiar velocities of galaxies can affect distance estimates. This effect is larger in the local volume and smaller at higher redshifts where the Hubble flow dominates over peculiar velocities. The α.40 catalog also includes radial distances (Haynes, et al. 2011), using a local volume flow model (Masters 2005) for galaxies with czCMB < 6000km.s −1 . This model has an estimated local velocity dispersion of σv = 163km.s −1 . For these galaxies we take the error in the distance to be the maximum of σv and 10% of the distance. For galaxies cz CMB > 6000km.s −1 distances are estimated using the Hubble and we take the error on distances to be 10%. Using the errors on distance and fluxes and their observed values we generate 300 Gaussian realizations and recompute M HI for every object. We apply the 2DSWML method to all the re- Figure A3. Comparison of the HIMF in the α.100 sample. The open circles is the HIMF by Jones, et al. (2018). The open squares is the HIMF estimated by our implementation of the 2DSWML method.
alizations and find out the width of the distribution for every φ j which we quote as an error.
(ii)Poisson Errors: The observed count in some of the bins is very low which means that it is important to consider errors due to Poisson counting.
(iii)Sample Variance: We estimate this error by splitting the area into 26 contiguous regions of approximately equal area each. We compute the HIMF by removing one region at a time. Finally the jackknife uncertainty for φ j is computed as σ whereφ j is the jackknife mean and φ j i is the value for the i th jackknife sample.
(iv) Other Errors: There are many other sources of errors which are discussed in Jones, et al. (2018). E.g. the error associated with the 2DSWML method which one can either estimate using the information matrix (Efstathiou, Ellis, & Peterson 1988) or by making further mocks (Jones, et al. 2018). One needs to understand how these errors are correlated with a finite sample or Poisson errors. We also expect Poisson errors to be correlated to mass errors especially in the lowest and highest mass bins where the observed counts are low. Another source of uncertainty is the error in velocity width W 50 . Since the peak flux, S peak , is inversely proportional to W 50 we expect their errors to be correlated. Since we integrate over W 50 to obtain the HIMF, we do not consider errors in W 50 . In order to properly account for errors one will need their covariances. We have attempted to add further sources of errors in quadrature but we find that the error bars become progressively larger and the Schechter function fits have a χ 2 red < 1 which means that we may be overestimating the errors. For this work we stick to the errors (i) − (iii) and add them in quadrature. These errors are comparable to Martin, et al. (2010) ;Haynes, et al. (2011).