The bivariate gas–stellar mass distributions and the mass functions of early- and late-type galaxies at $z\,{\sim}\,0$

Abstract

We report the bivariate $\rm HI$ - and $\rm H_{2}$ -stellar mass distributions of local galaxies in addition of an inventory of galaxy mass functions, MFs, for $\rm HI$ , $\rm H_{2}$ , cold gas, and baryonic mass, separately into early- and late-type galaxies. The MFs are determined using the $\rm HI$ and $\rm H_{2}$ conditional distributions and the galaxy stellar mass function (GSMF). For the conditional distributions we use the results from the compilation presented in Calette et al. [(2018) RMxAA, 54, 443.]. For determining the GSMF from $M_{*}\sim3\times10^{7}$ to $3\times10^{12}\ \text{M}_{\odot}$ , we combine two spectroscopic samples from the Sloan Digital Sky Survey at the redshift range $0.0033<z<0.2$ . We find that the low-mass end slope of the GSMF, after correcting from surface brightness incompleteness, is $\alpha\approx-1.4$ , consistent with previous determinations. The obtained $\rm HI\,$ MFs agree with radio blind surveys. Similarly, the $\rm H_{2}\,$ MFs are consistent with CO follow-up optically-selected samples. We estimate the impact of systematics due to mass-to-light ratios and find that our MFs are robust against systematic errors. We deconvolve our MFs from random errors to obtain the intrinsic MFs. Using the MFs, we calculate cosmic density parameters of all the baryonic components. Baryons locked inside galaxies represent 5.4% of the universal baryon content, while $\sim\! 96\%$ of the $\rm HI$ and $\rm H_{2}$ mass inside galaxies reside in late-type morphologies. Our results imply cosmic depletion times of $\rm H_{2}$ and total neutral H in late-type galaxies of $\sim\!1.3$ and 7.2 Gyr, respectively, which shows that late type galaxies are on average inefficient in converting $\rm H_{2}$ into stars and in transforming $\rm HI$ gas into $\rm H_{2}$ . Our results provide a fully self-consistent empirical description of galaxy demographics in terms of the bivariate gas–stellar mass distribution and their projections, the MFs. This description is ideal to compare and/or to constrain galaxy formation models.

1. Introduction

The determination of the matter–energy content of the Universe is one of the most important achievements from the recent advances in observational cosmology (e.g., Planck Collaboration et al. 2016; 2018). Current determinations are fully consistent with the spatially flat $\Lambda$ Cold Dark Matter ( $\Lambda\mbox{CDM}$ ) cosmology, with a present-day matter–energy content dominated by the cosmological constant, $\Omega_{\Lambda} = 0.689$ , and contributions of cold dark matter and baryon matter of $\Omega_{\rm cdm} \approx 0.262$ and $\Omega_{\rm bar} = 0.049$ , respectively (for a value of the normalised Hubble constant of $h=0.674$ , Aver et al. 2015; Cooke, Pettini, & Steidel 2018; Planck Collaboration et al. 2018). Therefore, the universal baryon mass density fraction is $f_{\rm bar,U} \equiv \Omega_{\rm bar} / \Omega_{\rm m} = 0.158$ , where $\Omega_{\rm m}= \Omega_{\rm cdm} + \Omega_{\rm bar}$ . How much of these baryons, and their different components, are locked inside galaxies? This paper addresses this question by quantifying the contribution from stars, atomic, and molecular gas in galaxies of different masses and morphological types.

According to the current paradigm of structure formation, non-baryonic dark matter played a major role in the evolution of the non-linear structures that we see today. Particularly, galaxies are believed to form and evolve within extended dark matter haloes, where multiple physical mechanisms are responsible for self-regulating star formation and thus setting up their observed properties (for reviews see, Mo, van den Bosch, & White 2010; Frenk & White 2012; Somerville & Davé 2015). As dark matter structures and galaxies evolve, baryons are redistributed from an initial smooth distribution to a more complex variety of structures. Of primordial importance for galaxy evolution is the amount of neutral hydrogen available for the formation of stars. Gas radiative cooling within the haloes regulates the inflow of cold gas to galaxies. The subsequent formation of stars is regulated by a complex interaction between cold gas inflows and the gas heating/outflows produced by the stars, a process that depends on halo mass. In low-mass halos, the stellar feedback, mostly form Supernova (SN) explosions, is able not only to heat the interstellar medium (ISM) but also to expel large gas fractions from the galaxy. In high-mass haloes, the long cooling time of shock-heated gas and the powerful feedback from rapidly accreting supermassive black holes that heats and/or expels the gas tend to suppress the star formation. Thus, it is not surprising that the expected fraction of baryons inside galaxies will differ from the universal baryon fraction, $f_{\rm bar,U}$ . Therefore, constraining the fraction of baryons and their different components in galaxies (mainly stars, atomic, and molecular gas) is essential to constrain the processes that have taken place during the evolution of the galaxies.

One of the main properties of galaxies are their stellar masses $M_{*}$ . Indeed, the abundance of galaxies as a function of $M_{*}$ provides important clues regarding the evolution of the galaxy population (e.g., Peng et al. 2010, 2012; Yang et al. 2012; Rodríguez-Puebla et al. 2017). Over the last two decades, there has been a remarkable progress in assembling large galaxy samples from multi-wavelength sky surveys that have led to robust determinations of the galaxy stellar mass function (GSMF; for recent discussions, and compilations of observations up to high redshifts see, Conselice et al. 2016; Rodríguez-Puebla et al. 2017). While there have been similar efforts in assembling galaxy samples for atomic gas mass, $M_{\rm HI}$ , based on radio blind observations (e.g., Zwaan et al. 2003; Meyer et al. 2004; Koribalski et al. 2004; Kovac et al. 2005; Martin et al. 2010; Haynes et al. 2011; Hoppmann et al. 2015; Haynes et al. 2018) or from follow-up subsamples based on optical/infrared surveys (e.g., Springob et al. 2005; van Driel et al. 2016), these are relatively shallow and/or in small volumes compared to the optical/infrared sky surveys, as well as strongly affected by selection effects. Therefore, the demographical analysis of $M_{\rm HI}$ is challenging especially when determining the low- and high-mass ends of the $\rm HI$ mass function ( $\rm HI$ MF) (for a more detailed discussion, see Jones et al. 2018)^a, as well as other statistics like the $\rm HI$ two-point correlation functions^b. The situation is not that different and even more challenging for the molecular gas as there are not blind galaxy samples in $\rm H_{2}$ . Nonetheless, there are some notable efforts to use optically selected samples combined with small and shallow CO surveys to indirectly derive, from the (uncertain) CO-to- $\rm H_{2}$ mass conversion factor, the galaxy MF in $\rm H_{2}$ , $\rm H_{2}$ MF (e.g., Keres, Yun, & Young 2003; Lagos et al. 2014; Saintonge et al. 2017; Andreani et al. 2018). Unfortunately, these CO surveys are also subject to incompleteness and selection effects or subject to a large fraction of galaxies with upper limits reported due to flux detection limits.

As mentioned above, galaxy formation is a non-linear and complex process. Remarkably, well-defined correlations (usually power laws) are, however, found from the observations. Among these are the correlations between the $\rm HI$ and stellar mass, $M_{\rm HI}$ – $M_{*}$ , and the $\rm H_{2}$ and stellar mass, $M_{\rm H2}$ – $M_{*}$ . While both correlations present large scatters, when divided into early- and late-type galaxies they tend to show different and tighter correlations (e.g., Calette et al. 2018, and more references therein). This is not surprising given that the formation histories of early- and late-type galaxies were different. Thus, understanding the contribution of these two populations to the abundance of galaxies traced by $\rm HI$ and $\rm H_{2}$ provides further key constrains to galaxy formation theory models.

In a recent work, (Calette et al. 2018, hereafter Paper I) were able to determine empirically not only the mean $M_{\rm HI}$ – $M_{*}$ and $M_{H_{2}}$ – $M_{*}$ relations and their scatters for early- and late-type galaxies but also the full conditional probability distribution functions (CPDFs) of $M_{\rm HI}$ and $M_{H_{2}}$ given $M_{*}$ , hereafter $\rm HI$ -CPDF and $\rm H_{2}$ -CPDF, respectively. In this paper, we combine the empirical CPDFs with the $z=0$ GSMF to derive the bivariate gas-to-stellar mass distributions and the MFs for the $\rm HI$ , $\rm H_{2}$ , cold gas, and baryon components, for all galaxies as well as for early and late types. Thus, the present paper represent a natural continuation of Paper I with some updates. These updates include new constraints on the best-fitting parameters to the observed CPDFs from Paper I.

In this paper, we compute the GSMF and its decomposition into early- and late-type galaxies. While there are many studies that have determined the GSMF in the past, they do not typically report systematic errors or do not deconvolve it from random errors (with a few exceptions, e.g., Bernardi et al. 2010, 2017; Obreschkow et al. 2018) or they are limited in the dynamical range of $M_{*}$ due to the limited depth of the sample and/or the cosmic variance in the galaxy sample (but see Wright et al. 2017). We combine here two large galaxy samples, the low-redshift sample, low-z, from the NYU Sloan Digital Sky Survey (SDSS) DR4 (Blanton et al. 2005a; b), and the new photometry pipeline for the SDSS DR7 from Meert et al. (2015, 2016). The low-z sample suffers from surface brightness incompleteness, but here we estimate and correct for the fraction of missing galaxies due to this selection effect. As for the SDSS DR7, the new photometry from Meert et al. (2015) shows that galaxy magnitudes were previously underestimated due to sky subtraction problems (see also, Simard et al. 2011); the impact of these new determinations has been studied previously in Bernardi et al. (2017). We extend the Bernardi et al. (2017) analysis using not only different definitions of galaxy stellar masses but by dividing into two morphology groups, early- and late-type galaxies.

The results reported in this paper integrate the $\rm HI$ - and $\rm H_{2}$ -CPDFs with new determinations for the GSMF to offer a full statistical description of the local galaxy demographics traced by the stellar, $\rm HI$ , $\rm H_{2}$ , total cold gas, and baryon mass components. This statistical description of the local galaxy demographics is much more complete than the typically employed GSMF for constraining models and simulations of galaxy formation. The new generation of semi-analytic models (e.g., Croton et al. 2016; Lagos et al. 2018; Henriques et al. 2019; Yung et al. 2019) and cosmological hydrodynamics simulations (e.g., Hirschmann et al. 2014; Vogelsberger et al. 2014; Schaye et al. 2015; Pillepich et al. 2018; Davé et al. 2019), and their post-processing outcomes, are now able to predict stellar, $\rm HI$ , and $\rm H_{2}$ masses for large galaxy populations in cosmological boxes (see e.g., Lagos et al. 2015; Diemer et al. 2018, 2019; Popping et al. 2019). The empirically based results presented here are optimal for comparing with these predictions as well as for calibrating theoretical models of galaxy evolution (see e.g. Romeo 2020). The results to be presented in this paper are the basis for further studies as the inference of the galaxy–halo connection extended to $\rm HI$ , $\rm H_{2}$ , cold gas, and baryon masses.

The present paper is organised as follows. In Section 2, we describe our method of using CPDFs in order to derive galaxy MFs traced by atomic, molecular, and cold gas masses as well as by the baryonic mass. In Section 3, we describe the samples we use to derive our local GSMF divided into early- and late-type galaxies. In Section 4, we present the results for our inventory of galaxy MFs and compare them with direct observational results. We also present our estimates for the cosmic density parameters related to the different baryonic components in galaxies. Section 5 discusses the impact of systematics and random errors. In Section 6, we present a summary and our main conclusions.

In this paper, we adopt cosmological parameter values that are close to the Planck mission: $\Omega_{\Lambda} = 0.693$ , $\Omega_{\rm m} = 0.307$ , $\Omega_{\rm bar} = 0.048$ , and $h = 0.678$ . All stellar masses are normalised to a Chabrier (2003) initial mass function (IMF).

2. Modelling the bivariate distributions and MFS from the conditional distribution functions

In this section, we describe the statistical method for deriving the $\rm HI$ and $\rm H_{2}\,$ MFs (as well as the total cold gas and baryon MFs), from the GSMF and the respective correlations of $M_{\rm HI}$ and $M_{H_{2}}$ with $M_{*}$ , or more generally, the respective full mass conditional distribution functions, CPDFs. In general, our approach allows to calculate bivariate distribution functions of the $\rm HI$ or $\rm H_{2}$ mass and the stellar mass. One can imagine that our methodology is equivalent to an optically selected volume-limited sample that it is complete in stellar mass, with HI and ${\rm H_{2}}$ gas masses determined for every galaxy in the sample, and for which any MF can be determined. When information about morphology is available, the CPDFs are useful for deriving the corresponding MFs into different morphological components. Here, we consider that the galaxy population is divided into two main morphological groups: early- and late-type galaxies. Following Paper I, our definition of early-type galaxies includes morphological types that comprise E and S0 galaxies or equivalently $T\leq 0$ from the Nair & Abraham (2010) morphology classification. Late-type galaxies are just the complement, from Sa to Irr. Below, we briefly describe the basic ingredients for calculating the MFs:

• Conditional Distribution Functions: For a fixed morphology, a galaxy of mass $M_{*}$ has the chance of having either a $\rm HI$ or $\rm H_{2}$ mass described by their corresponding CPDFs. We denote the CPDFs of early and late types by $P_{E} (M_j| M_{*}$ ) and $P_{L} (M_j| M_{*}$ ), respectively, where $ j = {\rm HI}$ or $\rm H_{2}$ . The HI-CPDF and H2-CPDF contain information about all the moments of the HI- and ${\rm H_{2}}$ -to-stellar mass correlations. We use the observed HI-CPDF and H2-CPDF from Paper I. In Section 2.2, we describe the functional forms for the CPDFs proposed in Paper I.

• Galaxy Stellar Mass Function: The GSMF is an important input since it allows us to project the CPDFs into their corresponding MFs. We derive the GSMF for all galaxies, as well for the early and late type, based on the SDSS. Section 3 describes our methodology to compute the observed GSMF over $\sim\!5$ decades in $M_{*}$ , as well as its decomposition into early- and late-type galaxies.

The reader interested in our resulting MFs and bivariate distributions may skip to Section 4.

2.1. Generalities

As discussed above, a CPDF, $P_j\text{(M}_j| M_{*}\text{)}$ , determines the chances that a galaxy of mass $M_{*}$ possess a specific galaxy property $M_j$ , with $ j = {\rm HI}$ , $\rm H_{2}$ , cold gas, or baryonic mass. Note that the units of $P_j$ is per $\text{M}_{\odot}$ . The relation between the distribution $P_j$ in bins per $\text{M}_{\odot}$ to dex-1, $\mathcal{P}_{j} $ , is given by

(1) \begin{equation}\mathcal{P}_{j} \text{(}M_j| M_{*}\text{)} = P_{j}\text{(}M_j| M_{*}\text{)} \times \frac{M_j}{\log e}.\end{equation}

The advantage of using $P_j\text{(}M_j|M_{*}\text{)}$ is that it contains information about all the moments of the distribution, in particular the mean $M_j-M_{*}$ relation and its standard deviation.

The joint distribution function of $M_{*}$ and $M_j$ , hereafter referred as the bivariate distribution function, is defined as:

(2) \begin{equation}\Phi(M_{j},M_{*}) = \frac{d^2 N\text{(}M_{j}| M_{*}\text{)}} {V d\log M_{j} d\log {M_{*}}}=\mathcal{P}_{j}\text{(}M_{j}| M_{*}\text{)} \phi_* \text{(}M_{*}\text{)},\end{equation}

where $d^2 N$ is the bivariate number of galaxies within the mass range $\log {M_{*}}\pm d\log{M_{*}}/2$ and $\log M_{j}\pm d\log M_{j}/2$ in a given volume V, and $\phi_* \text{(}M_{*}\text{)}$ is the GSMF in units of ${\rm Mpc}^{-3} {\rm dex}^{-1}$ . The integration (marginalisation) of $\Phi(M_{j},{M_{*}})$ over $M_{*}$ results in the total MF for $M_j$ , $\phi_j(M_j)$ , that is,

(3) \begin{align}\phi_{j} (M_j) &=\int \Phi(M_{j},{M_{*}})d\log{M_{*}}\nonumber\\ &= \int \mathcal {P}_{j}(M_j|{M_{*}}) \phi_* \text{(}M_{*}\text{)} d\log{M_{*}}.\end{align}

The above equation shows how the CPDFs are projected into a number density function via the GSMF. Note that integration of $\Phi(M_{j},{M_{*}})$ over $M_j$ gives the total GSMF^c.

As discussed previously, when studying the properties of galaxies, it is useful to separate them into, at least, two morphological components such as early types, or spheroid-dominated galaxies, and late types, or disc-dominated galaxies. Thus, the total GSMF can be formally represented as the contribution of these two types

(4) \begin{equation}\phi_*\text{(}M_{*}\text{)} = \phi_{*,E}\text{(}M_{*}\text{)} + \phi_{*,L}\text{(}M_{*}\text{)},\end{equation}

denoted respectively by $\phi_{*,E}$ , and $\phi_{*,L}$ . In terms of the fraction of early- and late-type galaxies ( $f_E$ and $f_L$ ), their corresponding galaxy stellar MFs are given respectively by $\phi_{*,E} = f_E\times\phi_*$ , and $\phi_{*,L} = f_L\times\phi_*$ , with $f_E+f_L = 1$ .

Early- and late-type galaxies are different in their HI- and ${\rm H_{2}}$ -to-stellar mass distributions. Thereby, Equation (3) can be generalised in terms of the distribution $\mathcal {P}_{i,j}(M_j|{M_{*}})$ , where the subscripts indicate $i = $ early or late type, and $ j = $ $\rm HI$ , $\rm H_{2}$ , cold gas, or baryonic mass. Then, the generalisation of Equation (3) to galaxies with morphological type i and mass component j is

(5) \begin{equation}\phi_{j, i} (M_j) = \int f_i\text{(}M_{*}\text{)} \mathcal {P}_{i,j}(M_j|{M_{*}}) \phi_* \text{(}M_{*}\text{)} d\log{M_{*}}.\end{equation}

Finally, the total CPDFs are calculated from the respective conditional distributions of early- and late-type galaxies as:

(6) \begin{align}\mathcal {P}_{j}\text{(}M_{j}| M_{*}\text{)} &= f_{E}\text{(}M_{*}\text{)}\times \mathcal {P}_{E,j} \text{(}M_{j}| M_{*}\text{)}\nonumber\\& \quad +f_{L}\text{(}M_{*}\text{)}\times \mathcal {P}_{L,j} \text{(}M_J| M_{*}\text{)},\end{align}

with $ j = {\rm HI}, {\rm H_{2}},$ cold gas, or baryonic mass.

2.2. The $\text{\em HI}$ and $\text{\em H}_{\text{\em 2}}$ conditional distribution functions

As shown in Equations (2) and (5), the conditional or bivariate distribution functions are useful to statistically determine the MFs. Evidently, in the case of atomic and molecular gas, we are assuming that for every galaxy that is optically selected, there must exist $\rm HI$ and $\rm H_{2}$ counterparts. The discussion on the possible existence of pure $\rm HI$ or $\rm H_{2}$ galaxies, those that will not be observed in optically selected samples but rather in radio blind surveys, is out of the scope of this paper. Note that if they exist, the chance of observing those galaxies is very low over the mass ranges that we will derive the MFs. For example, in the case of pure $\rm HI$ galaxies, the Arecibo Legacy Fast ALFA (ALFALFA) survey has found $\sim\!1.5\%$ of $\rm HI$ sources that were not clearly associated with an optical counterpart. Of those, $\sim\!75\%$ are likely tidal in origin (Haynes et al. 2011). Thus, $\sim\!0.4\%$ of HI source observed in the ALFALFA survey are purely gaseous galaxies candidates, most of them at the mass range $10^{7}<M_{\rm HI}/\text{M}_{\odot}< 10^{10}$ (Cannon et al. 2015). As we will show, our completeness limit for the $\rm HI$ MF is $M_{\rm HI}\sim 10^{8}\text{M}_{\odot}$ . The above fraction could be considered as an upper limit as some of these sources have already detected optical counterparts revealing unusual high HI mass-to-light ratios (Cannon et al. 2015). Thus, we conclude that our results are unlikely to be affected by excluding pure gas galaxies in our analysis.

2.2.1. The Calette et al. (2018) $\text{\em HI}$ and $\text{\em H}_{\text{\em 2}}$ conditional distribution functions

Here, we use the results from Paper I (Calette et al. 2018) who determined the $\rm HI$ - and $\rm H_{2}$ -to-stellar mass ratio distributions (CPDFs) as a function of $M_{*}$ from a large compilation of optically selected samples with radio observations. Next, we briefly describe the steps taken in Paper I to derive the $\rm HI$ and $\rm H_{2}$ CPDFs. The reader is referred to that paper for details.

The compiled data described in Paper I consist of a set of incomplete and inhomogeneous samples. We first homogenised all these samples to a common IMF, cosmology, radio telescope configuration and sensitivity, and CO-to-luminosity conversion factor. Then, we selected only those samples without obvious biases due to selection effects such as environment. Radio non-detections, reported in the literature as upper limits, are an important source of uncertainty when deriving distributions or correlations. In Paper I, we included non-detections to derive the $\rm HI$ and $\rm H_{2}$ CPDFs. Below, we briefly describe the treatment that we employed for radio non-detections.

In our compiled samples, most of radio non-detections are early-type galaxies representing a non-negligible fraction of intermediate and massive galaxies, which are (typically) gas poor. An important fraction of those galaxies are from the GASS (Catinella et al. 2013) and the CO Legacy Database for the GALEX Arecibo SDSS survey (COLD-GASS) (Saintonge et al. 2011) surveys at distance of $109<D/{\rm Mpc}<222$ . Compared to other more nearby samples of intermediate and massive early-type galaxies with measurements of $\rm HI$ and $\rm H_{2}$ mass, such as the ATLAS 3D (Serra et al. 2012) at $\bar{D}\,{\sim}\,25\, {\rm Mpc}$ , we noted that the upper limits of the GASS/COLD-GASS samples are $\sim\!1-2$ orders of magnitude larger than nearby samples (Paper I). The above lead us to first introduce a correction for the upper limits of the GASS/COLD-GASS surveys by a distance effect. Recall that radio non-detections or upper limits depend not only on the sensitivity of the radio telescope or integration time but also on the distance to the object. In Paper I, we corrected the upper limits of the GASS/COLD-GASS samples by a distance effect using nearby samples such as the ATLAS 3D survey. Briefly, our correction consists in using the distances and upper limits from nearby samples to estimate the upper limits in the GASS and COLD-GASS as if these two samples were at the same distance as the nearby ones. We validated our procedure using a mock galaxy survey by applying similardistance sensitivity effects as GASS/COLD-GASS surveys, for details see Paper I. For late-type galaxies, notice that most of them are detected in radio due to their large fractions of gas and it is not necessary to introduce the above corrections. Next, we describe the treatment of the upper limits to derive the $\rm HI$ and $\rm H_{2}$ CPDFs.

In our analysis from Paper I, we included upper limits, or left-censored data, using the Kaplan & Meier (1958) non-parametric estimator. This estimator provides a reconstruction of information lost by censoring. Feigelson & Nelson (1985) adapted this estimator for astronomical samples. We used the Astronomical SURVival statistics (ASURV) package based on Feigelson & Nelson (1985) to derive the $\rm HI$ and $\rm H_{2}$ CPDFs from our compiled samples. We have also applied the censoring Buckley & James (1979) regression method to derive the relationship and standard deviations between the $\rm HI$ - and $\rm H_{2}$ -to-stellar mass ratio and $M_{*}$ . We note that the regression results are consistent with the (logarithmic) mean and standard deviation values obtained from the CPDFs based on the Kaplan & Meier (1958) estimator.

2.2.2. The functional forms of the $\text{\em HI}$ and $\text{\em H}_{\text{\em 2}}$ conditional distribution functions

For the $\rm HI$ and $\rm H_{2}$ CPDFs of late-type galaxies, in Paper I, we found that they are described by a Schechter function. In the case of early-type galaxies, the CPDFs are better described by a (broken) Schechter function plus a uniform distribution at the low $-\mathcal{R}_j$ values. Following, we describe in more detail these functional forms.

We begin by introducing the following Schechter-type probability distribution function for the $\rm HI$ - or $\rm H_{2}$ -to-stellar mass ratios, $ \mathcal{R}_j=M_j/M_{*}$ , in the range $\log \mathcal{R}_j \pm d \log \mathcal{R}_j / 2$ :

(7) \begin{equation}\mathcal{S}_{i,j} (\mathcal{R}_j) = \frac{\ln(10)}{\mathcal{N}_{i,j}}\left(\frac{\mathcal{R}_{j}}{\mathcal{R}^*_{i,j}}\right)^{\alpha_{i,j} + 1}\exp\left({-}\frac{\mathcal{R}_{j}}{\mathcal{R}^*_{i,j}}\right),\end{equation}

where the morphology is represented with $i = $ early or late type, and the galaxy property is represented with $ j = $ $\rm HI$ or $\rm H_{2}$ . The parameters are the characteristic gas-to-stellar mass ratio, $\mathcal{R}^*_{i,j}$ , the normalisation parameter, $\mathcal{N}_{i,j}$ , which constrains the probability to be between zero and one,^d and the power law slope $\alpha_{i,j}$ for the part of the distribution of galaxies with low gas-to-stellar mass ratio.

• Late-type galaxies:

For late-type galaxies, that is $i = L$ , in Paper I, we found that the HI-CPDF and H2-CPDF are described by the Schechter-type distribution function given by (7) with the parameters $\alpha_{L,j}$ and $\mathcal{R}^*_{L,j}$ functions of $M_{*}$ as follows:

(8) \begin{equation}\alpha_{L,j} = \alpha_{0;L,j} \log {M_{*}} + \alpha_{1;L,j} ,\end{equation}

and

(9) \begin{equation}\mathcal{R}^*_{L,j} = \frac{\mathcal{R}^*_{0;L,j}}{\left(\frac{{M_{*}}}{\mathcal{M}^*_{L,j}} \right)^{\beta_{L,j}} + \left(\frac{{M_{*}}}{\mathcal{M}^*_{L,j}} \right)^{\gamma_{L,j}}}.\end{equation}

Consider that $\mathcal{S}_{L,j} (\log\mathcal{R}_j) d\log \mathcal{R}_j = \mathcal{S}_{L,j} (\log M_j - \log{M_{*}}) d\log (M_j / {M_{*}})$ . By definition $M_{*}$ is fixed, thus the $\rm HI$ and $\rm H_{2}$ CPDFs of late-type galaxies are given by:

(10) \begin{equation}\mathcal {P}_{L,j}(M_j | {M_{*}}) d\log M_j = \mathcal{S}_{L,j} (\log M_j - \log{M_{*}}) d\log M_j.\end{equation}

The above explicitly shows that the integration over conditional distribution functions can also be interpret as convolutions in Equation (3).

• Early-type galaxies:

In the case of early-type galaxies, $i = E $ , we showed in Paper I that both for the HI-CPDF and H2-CPDF are described as the sum of two distribution functions; the Schechter-type distribution function, $\mathcal{S}_{E,j}$ , and a uniform function, $\mathcal {U}_{0,j}$ , are

(11) \begin{equation}\mathcal {E}_{j} (\mathcal{R}_j) = \left\{\begin{array}{l@{\quad}l}\mathcal {U}_{0,j} & \mathcal{R}_{0,j} \le\mathcal{R}_j < \mathcal{R}_{1,j} \\A \times \mathcal{S}_{E,j} (\mathcal{R}_j) & \mathcal{R}_{1,j} \le \mathcal{R}_j\end{array},\right.\end{equation}

where $\mathcal{R}_{0,j} = \mathcal{R}_{1,j} / 10,$ ^e and $ \log \mathcal{R}_{1,j} = r_{0,j} \log{M_{*}} + r_{1,j}$ , while the uniform distribution is given by

(12) \begin{equation}\mathcal {U}_{0,j} \text{(}M_{*}\text{)} = \frac{p_{0,j} \log{M_{*}} + p_{1,j}}{\Delta},\end{equation}

and

(13) \begin{equation}A = \left( 1 - \mathcal {U}_{0,j} \times \Delta\right) \times \frac{\mathcal{N}_{i,j}}{\eta_{i,j}(\mathcal{R}_{1,j})},\end{equation}

where in Paper I we assumed that $\Delta = \log_{10} = 1$ dex, the symbol $\eta_{i,j}(\mathcal{R}_{1,j})$ takes into account the fraction of galaxies in the Schechter-type mode for galaxies with gas ratio above $\mathcal{R}_{1,j}$ .^f The HI-CPDF and $\rm H_{2}$ -CPDF of early-type galaxies are

(14) \begin{equation}\mathcal {P}_{E,j}(M_j | {M_{*}}) d\log M_j = \mathcal{E}_{j}(\log M_j - \log{M_{*}}) d\log M_j.\end{equation}

2.2.3. Constraints on the best-fitting parameters

In Paper I, the best-fit parameters for late-type galaxies, Equations (8)–(9), and for early-type galaxies, Equations (11)–(13), were constrained using the observed $\rm HI$ - and $\rm H_{2}$ -CPDFs on various stellar mass bins. Computing CPDFs over $M_{*}$ bins requires the GSMF in addition of the fraction of early- and late-type galaxies (see Section 4.1 for more details). Since we are using slightly different inputs, namely the GSMF and the fractions of early- and late-type galaxies, than in Paper I, we prefer to perform our own fits to the same data, for consistency. The results are presented in Section 4.1. The differences with the parameters reported in Paper I are actually small.

2.3. The cold gas and baryonic conditional distribution functions

Once we have constructed the HI-CPDF and $\rm H_{2}$ -CPDF, we can now define the conditional distributions for the cold gas and baryon masses, $M_{\rm gas}$ and $M_{\rm bar}$ .

The total cold gas content in a galaxy is composed of $\rm HI$ , $\rm H_{2}$ , helium, and metals; helium and metals account for roughly 30% of the cold gas, $M_{\rm He} + M_{\rm Z}\approx0.3 M_{\rm gas}$ Therefore, $M_{\rm gas} = M_{\rm HI} + M_{H_{2}} + M_{\rm He} + M_{\rm Z} = 1.4 \times (M_{\rm HI} + M_{H_{2}}\text{)}$ . For simplicity, let $M_{\rm HI}$ and $M_{H_{2}}$ be two independent random variables. Section 5 discusses the validity of this assumption. Then, $M_{\rm gas}$ is a random variable with the conditional distribution function:

(15) \begin{align}P_{\rm gas}(M_{\rm gas}|{M_{*}}) &=\frac{1}{1.4} \int P_{\rm HI}\left(0.71 M_{\rm gas} - M_{H_{2}}|{M_{*}}\right) \nonumber\\[2pt] & \quad \times P_{\rm H_2}(M_{H_{2}}|{M_{*}}) d M_{H_{2}}, \nonumber\\[2pt] &=\frac{1}{1.4} \int P_{\rm HI}(M_{\rm HI}|{M_{*}}) \nonumber\\[2pt] & \quad \times P_{\rm H_2}(0.71 M_{\rm gas} - M_{\rm HI}|{M_{*}}) d M_{\rm HI},\end{align}

or after some algebra, the same distribution function but per bin in log space is:

(16) \begin{align}\mathcal{P}_{\rm gas}({M_{\rm gas}}|{M_{*}}) &= \int \frac{\mathcal{P}_{\rm HI}(0.71{M_{\rm gas}} - {M_{H_{2}}}|{M_{*}})}{1 - 1.4\ {M_{H_{2}}}/{M_{\rm gas}}} \nonumber\\[2pt] & \quad \times \mathcal{P}_{\rm H_2}({M_{H_{2}}}|{M_{*}}) d \log {M_{H_{2}}}, \nonumber\\[3pt] &= \int \mathcal{P}_{\rm HI}({M_{\rm HI}}|{M_{*}}) \nonumber\\[3pt] &\quad\times \frac{\mathcal{P}_{\rm H_2}(0.71 {M_{\rm gas}} -{M_{\rm HI}}|{M_{*}})}{1 - 1.4\ {M_{\rm HI}}/{M_{\rm gas}}} d \log {M_{\rm HI}}.\end{align}

For the baryonic conditional distribution functions, we again assume that $M_{\rm gas}$ and $M_{*}$ are two independent random variables. Thus, $M_{\rm bar} = M_{\rm gas} + M_{*}$ is a random variable with a distribution function given by

(17) \begin{align}P_{\rm bar}(M_{\rm bar}|{M_{*}}) &= \int P_{\rm gas} (M_{\rm bar} - \mathcal{M}_*|{M_{*}}) \nonumber \\[2pt] &\quad\times \delta(\mathcal{M}_* - {M_{*}}) d \mathcal{M}_*\\[3pt] &= P_{\rm gas} (M_{\rm bar} - {M_{*}}|{M_{*}}),\nonumber\end{align}

where $P_{\rm gas}$ is the conditional distribution function for gas, Equation (15), and the Dirac- $\delta$ function appears explicitly for the $M_{*}$ term. Similarly as above, we find that

(18) \begin{equation}\mathcal{P}_{\rm bar}(M_{\rm bar}|{M_{*}}) = \frac{\mathcal{P}_{\rm gas}(M_{\rm bar} - {M_{*}}|{M_{*}})}{1 - {M_{*}} / M_{\rm bar}}.\end{equation}

Finally, we derive the gas and baryon MFs using Equations (5), (16), and (18), the last two valid for early- and late-type galaxies.

3. The GSMF of all, early-, and late-type galaxies

The preceding section described a methodology to use the GSMF as an interphase that transforms galaxy CPDFs into MFs, see Equation (5). In this section, we briefly describe how we determine the local GSMF for masses above $M_{*}\,{\sim}\,3\times10^{7} \text{M}_{\odot}$ , as well as the GSMF s for early- and late-type galaxies. For a more detailed description of the galaxy samples utilised here and the different corrections we apply, the reader is referred to Appendices A–C.

3.1. The galaxy samples and the GSMF

To estimate the GSMF over a large dynamical range, we use two galaxy samples. Next, we shortly describe the procedure and our determinations.

(1) For masses above $M_{*} = 10^9\ \text{M}_{\odot}$ , we use the SDSS DR7 based on the photometric catalogue from Meert et al. (2015) and Meert et al. (2016)^g at the redshift interval $0.005<z<0.2$ . Previous studies have concluded that the measurements of the apparent brightnesses based on the standard SDSS pipeline photometry are underestimated due to sky subtraction problems, particularly, in crowded fields (Bernardi et al. 2010; Blanton et al. 2011; Simard et al. 2011; Bernardi et al. 2013; He et al. 2013; Mendel et al. 2014; Kravtsov, Vikhlinin, & Meshscheryakov 2014; Meert et al. 2015; D’Souza, Vegetti, & Kauffmann 2015; Bernardi et al. 2016; Meert et al. 2016. New determinations of the GSMF based on the new algorithms for obtaining more precise measurements of the sky subtraction, and thus to improve the photometry, have concluded that the bright end of the luminosity/MF has been systematically underestimated (Bernardi et al. 2017). While there are various groups working in improving the determination of galaxy apparent brightnesses, see references above, Bernardi et al. (2017) showed that all those studies agreed up to 0.1 dex in the GSMF. In this paper, we use the apparent Sérsic r, g, and i band luminosities reported in Meert et al. (2015) and (2016) derived for the SDSS DR7 based on the PyMorph software pipeline (Vikram et al. 2010; Meert, Vikram, & Bernardi 2013). This software has been extensively tested in Meert et al. (2013) and shows that it does not suffer from sky subtraction problems. All magnitudes and colours are K+E corrected at a redshift rest frame $z=0$ , see Appendixes A and B. As described in Appendix A, for every galaxy, we estimate $M_{*}$ from five colour-dependent mass-to-light ratios but we define as our fiducial $M_{*}$ the geometric mean of all the determinations. Using the $1/{V}_{\rm max}$ method, we derive six $\mbox{GSMF}$ s based on the mass definitions described above. Consistent with Bernardi et al. (2017), we find that the differences in mass-to-light ratios introduce large discrepancies in the GSMF, especially at the high-mass end. In Figure 11 from Appendix A, we find that a shift of $\sim\pm0.15$ dex in the $M_{*}$ axis recovers systematic errors in the GSMF due to different mass-to-light ratios.

(2) For masses below $M_{*} =10^9\ \text{M}_{\odot}$ , we use the SDSS DR4 NUY-VAGC low-z sample,^h at the redshift interval $0.0033<z<0.005$ , and ideal to study the low mass/luminosity galaxies (Blanton et al. 2005a; Blanton et al. 2005b). As before, all absolute magnitudes and colours were K+E corrected at a redshift rest frame $z=0$ . Also, we derive $M_{*}$ from five colour-dependent mass-to-light ratios and, again, we define our fiducial $M_{*}$ as the geometric mean of all the determinations. We construct the GSMF using the $1/V_{\rm max}$ method and include missing galaxies due to surface brightness incompleteness, as described in Appendix C. For surface brightness incompleteness, we follow closely the methodology described in Blanton et al. (2005b). The latter correction is relevant for the low-mass end. Based on the conclusions from Baldry et al. (2012), we use a simple correction for the low-mass end in order to correct for the local flow model distances from Willick et al. (1997) to the one by Tonry et al. (2000).

Our final GSMF is the result of combining the SDSS NUY-VAGC low-z sample, for galaxies with masses $M_{*}\approx 3\times 10^{7}\text{M}_{\odot}$ to $M_{*}\approx10^{9}\text{M}_{\odot}$ , and the SDSS DR7 sample for galaxies with $M_{*}\gtrsim 10^{9}\text{M}_{\odot}$ , based on our fiducial $M_{*}$ determination. Figure 1 presents our final GSMF with the black solid circles and error bars. The black solid line shows the best fit to the data (described below), and the grey shaded area shows a shift in the $M_{*}$ axis of $\pm 0.15$ dex. As discussed above, in Appendix A, we find that this is a good approximation to the systematic errors in the GSMF due to differences in the mass-to-light ratios. In the same figure, we include comparisons to previous works. In order to account for differences in cosmologies, we scale previous studies to our cosmology using the following relations:

(19) \begin{equation}\phi_{*, {\rm us}} = \phi_{*, {\rm lit}}\left(\frac{h_{\rm us}}{h_{\rm lit}}\right)^3,\end{equation}

and

(20) \begin{equation}M_{*, {\rm us}} = M_{*, {\rm lit}}\left(\frac{h_{\rm lit}}{h_{\rm us}}\right)^2,\end{equation}

where $h_{\rm us} = 0.678$ and $h_{\rm lit}$ is the respective value reported in the literature. Nonetheless, the impact of accounting for different cosmologies is small.

Figure 1. Observed GSMF when combining the SDSS NYU-VAGC low-redshift sample and the SDSS DR7 sample, black filled circles with error bars. We reproduce our results in the upper and the middle panels. The best-fit model composed of a Schechter function with a sub-exponential slope and a double power law function is shown as the black solid line. The shaded area shows an estimate of the systematic errors with respect to the best-fitting model. The bottom panel shows the residuals for our best-fitting model as a function of $M_{*}$ . We include comparisons to some previous observational determinations of the GSMF: in the upper panel we show determinations that are complete down to $\sim\! 10^9\ \text{M}_{\odot}$ , mostly based on the SDSS DR7, while in the middle panel we show determinations based on the GAMA survey, which are complete down to $\sim\! 3{-}5\times 10^7\ \text{M}_{\odot}$ , but suffer from cosmic variance at high masses due to the small volume.

In the upper panel of Figure 1, we reproduce the $\mbox{GSMF}$ s from previous determinations with stellar mass completeness above $M_{*}\,{\sim}\,10^9\ \text{M}_{\odot}$ . The violet triangles with error bars are the determinations from Moustakas et al. (2013), who used a spectroscopic sample of SDSS DR7 galaxies from the New York University Value Added Galaxy Catalog (NYU-VAGC) with redshifts $0.01<z<0.2$ combined with observations from GALEX. The red squares with error bars are the estimation obtained in Bernardi et al. (2013) from a sample of SDSS DR7 galaxies with photometry based on the PyMorph software pipeline at $z\,{\sim}\,0.1$ . Here, we reproduce their result based on Sérsic luminosities. Additionally, we compute the GSMF using the stellar mass estimates from Sérsic photometry from Mendel et al. (2014) who used the Simard et al. (2011) SDSS DR7 sample of g and r band photometry and extended to u, i, and z bands, blue filled circles with error bars. We show the best-fitting model from D’Souza et al. (2015), who estimated the GSMF by stacking images of galaxies with similar stellar masses and concentrations to correct Model magnitudes from the SDSS DR7, dark green solid line. Finally, we compare our result to Thanjavur et al. (2016), who derived the GSMF using the analysis from Mendel et al. (2014).

Our GSMF agrees well with previous determinations at the $\sim 10^{9.3}-10^{11}\ \text{M}_{\odot}$ range. At the high-mass end, it is shallower than previous determinations (e.g., Moustakas et al. 2013) except to Bernardi et al. (2013), who use Sérsic photometry from the SDSS DR7. As extensively discussed in Bernardi et al. (2017), there are two systematic effects that could lead to differences when comparing to previous determinations from the literature; assumptions on mass-to-light ratios and estimations of galaxy surface brightness. In the case of Moustakas et al. (2013) and D’Souza et al. (2015), who used cmodel and Model magnitudes, the comparison is not obvious due to systematic effects in both mass-to-light ratios and photometry (Bernardi et al. 2017). In the case of Mendel et al. (2014) and Bernardi et al. (2013), effects on photometry are not the dominant ones but mass-to-light ratios. Nonetheless, those differences are within the expected systematic effect, especially at the massive end, (Bernardi et al. 2017, see also Figure 11). We therefore conclude that when comparing to other previous determinations, the differences that we observe are consistent with the differences expected from systematic effects. Indeed, Figure 1 shows that most of the previous determinations are within our region of systematic errors. Thus, hereafter we will assume that our shift of $\pm0.15$ dex in the $M_{*}$ axis approximately captures systematics not only from stellar population models but also from photometry.

The middle panel of Figure 1 presents comparisons to some previous determinations from deep but small-volume samples. The purple dots with error bars are from Baldry et al. (2008), who used the SDSS NYU-VAGC low-z sample but did not include missing galaxies due to surface brightness incompleteness. In addition, we compare to Baldry et al. (2012), who used the Galaxy And Mass Assembly (GAMA) survey for galaxies at $z<0.06$ , and complete down to $r=19.4$ mag for two-thirds of the galaxy sample and to $r=19.8$ for one-third of the sample. Finally, we reproduce the observed GSMF from Wright et al. (2017), who also used the GAMA survey to estimate the GSMF.

At low masses, our results are in excellent agreement with the GAMA $\mbox{GSMF}$ s. This is encouraging since the GAMA survey does not suffer from surface brightness incompleteness, at least within the stellar mass range that we are comparing our results. This is an indication that the surface brightness corrections described in Appendix C are able to recover the slope of the GSMF at low masses. Consistent with the values reported in Baldry et al. (2012) and Wright et al. (2017), we find that the faint-end slope of the GSMF is $\alpha \approx -1.4$ , below we describe in more detail the fitting model for the GSMF. The above is also in good agreement with Sedgwick et al. (2019) who recently determined the low-mass end of the GSMF by identifying low surface brightness galaxies based on data of core-collapse supernovae. The authors used the IAC Stripe 82 legacy project (Fliri & Trujillo 2016) and the SDSS-II Supernovae Survey (Frieman et al. 2008).

At the massive end we notice, however, some apparent tension between our and the GAMA results. Effects due to cosmic variance (due to the small redshift and angular coverage of the GAMA sample) could explain those differences as well as systematics in the mass-to-light ratios. Indeed, we see that some of the data are within the systematic errors. In addition, note that Figure C5 from Appendix C shows that using the mass-to-light ratios from Taylor et al. (2011), utilised in the Baldry et al. (2012) GSMF, tend to underestimate the high-mass end of the GSMF.

3.2. Best-fitting model to the GSMF

To provide an analytic form to our GSMF, we choose to use a function composed of a Schechter function with a sub-exponential decreasing slope and a double power law function. Note that the resulting high-mass end of our GSMF is shallower than an exponential function, and, thus, better fitted to a power law (see also Tempel et al. 2014). The Schechter sub-exponential function is given by:

(21) \begin{equation}\phi_{*, {\rm S}}(M_{*})=\phi^{*}_S\ln 10 \left(\frac{M_{*}}{\mathcal {M}_S}\right)^{1+\alpha_S} \exp\left[{-\left(\frac{M_{*}}{\mathcal {M}_S}\right)^{\beta}}\right],\end{equation}

where $\phi^{*}_S$ is the normalisation parameter in units of Mpc^–3 dex^–1, $\alpha$ is the slope at the low-mass end, $\mathcal {M}_S$ is the characteristic mass, and $\beta$ is the parameter that controls the slope at the massive end; note that $\beta=1$ corresponds to a Schechter function. The double power law function is given by:

(22) \begin{equation}\phi_{*, {\rm D}}(M_{*}) =\phi^{*}_D\ln 10\left(\frac{M_{*}}{\mathcal {M}_D}\right)^{1+\alpha_D}\left[1 + \left(\frac{M_{*}}{\mathcal {M}_D}\right)^{\gamma}\right]^{\frac{\delta-\alpha_D}{\gamma}},\end{equation}

where $\phi^{*}_D$ is the normalisation parameter in units of Mpc^–3 dex^–1, $\alpha$ and $\delta$ control the slope at low and high masses, respectively, while $\gamma$ determines the speed of the transition between the low- and high-mass regimes; and $\mathcal {M}_D$ is the characteristic mass of the transition. Finally, the analytic form for fitting the observed GSMF is given by

(23) \begin{equation}\phi_{*, {\rm model}}(M_{*}) = \phi_{*, {\rm S}}(M_{*}) + \phi_{*, {\rm D}}(M_{*}),\end{equation}

where we assumed that $\mathcal {M}_S =\mathcal {M}_D $ .

We find the best-fit parameters $\vec{p}_{\rm GSMF} = (\phi_{S}^{*}, \alpha_{S}, \mathcal{M}_{S}, \beta, \phi_{D}^{*}, \alpha_{D}, \delta,\gamma)$ that maximise the likelihood function $\mathcal{L}\propto \exp({-\chi^2/2})$ using the Markov chain Monte Carlo (MCMC) method algorithm described in Rodríguez-Puebla et al. (2013). Here

(24) \begin{equation}\chi^2 = \sum_{i=1}^{N_{\rm obs}}\left(\frac{\phi^{i}_{*,{\rm SDSS}} - \phi_{*, {\rm model}}^i }{\sigma^{i}_{\rm SDSS}}\right)^2,\end{equation}

with $N_{\rm obs} $ as the number of observational data points of the GSMF each with an ith value of $\phi^{i}_{*,{\rm SDSS}}$ and an error of $\sigma^{i}_{\rm SDSS}$ . The ith value of our model is given by $\phi_{*, {\rm model}}^i$ .

We sample the best-fit parameters by running a set of 10 chains with $1\times 10^{5}$ MCMC models each. Table 1 lists the best-fit parameters. For our best-fitting model, we find that $\chi^2=85.42$ from a number of $N_{\rm obs} = 50$ observational data points. Our model consist of $N_{\rm p} = 8$ free parameters, thus the reduced $\chi^2$ is $\chi^2 / {\rm d.o.f.} = 2.03$ . The upper and middle panels of Figure 1 show our best-fitting model as the black solid line and the bottom panel shows the residuals as a function of $M_{*}$ . Our best-fitting model has an error of $\sim\!2\%$ at the range $M_{*}\sim\!2\times10^{9}-5\times10^{11}\text{M}_{\odot}$ and an error lower than $\sim\!10\%$ at the mass range $M_{*}\sim\!7\times10^{8}-1\times10^{12}\text{M}_{\odot}$ . For lower masses, errors can be up to $\sim\!20\%$ .

Table 1. Best-fitting parameters for the GSMF Equations (21)–(23)

A valid question is how much we improve the analytic prescription when using a Schechter sub-exponential plus a double power law function model confronted to a double Schechter function model, commonly employed by previous authors (see e.g., Baldry et al. 2012; Wright et al. 2017). We have explored this possibility but assuming a Schechter function, $\beta=1$ in Equation (21), and a Schechter sub-exponential function, that is, we are adding a extra degree of freedom due to the shallow decay at the high-mass end. Based on this alternative, we repeat our fitting procedure but this time finding that $\chi^2=662.817$ from a number of $N_{\rm obs} = 50$ observational data points. Now, our model consists of $N_{\rm p} = 6$ free parameters resulting in a reduced $\chi^2$ of $\chi^2 / {\rm d.o.f.} = 15.06$ . This is considerably worse when combining Schechter sub-exponential and double power law functions. Thus, hereafter, we will consider only the latter model.

3.3. The GSMFs of early- and late-type galaxies

Our main goal for this paper is to construct bivariate distributions as well as MF based on the observed gas mass CPDFs and the $\mbox{GSMF}$ s of early- and late-type galaxies. In this section, we determine the GSMF of early- and late-type galaxies from the SDSS DR7 spectroscopic sample with the public automated morphological galaxy classification by Huertas-Company et al. (2011).ⁱ The morphological classification in Huertas-Company et al. (2011) was determined based on support vector machine algorithms. Here we use their tabulated probabilities for each SDSS galaxy as being classified as an early type, P(E). For masses below the completeness of the SDSS DR7 sample, we use an extrapolation of the observed fraction of early-type galaxies. We will come back to this point later in this section.

From a catalogue of galaxies with visual morphological classification (UNAM-KIAS; Hernández-Toledo et al. 2010), we find that galaxies with types $T\le0$ are mostly those with $P(E)>0.65$ , and those with $P(E)\le0.65$ correspond mostly to $T>0$ ; here T is the Fukugita et al. (2007) notation.^j Based on the above, we consider as early-type galaxies those with a probability $P(E)>0.65$ , while late-type galaxies those with $P(E)\le0.65$ . We checked that our morphology definition between early- and late-type galaxies is consistent with the morphological classification based on the concentration parameter $c=R_{90}/R_{50}$ . That is, the division between early- and late-types is approximately at $c=2.85$ (see below and also, Hyde & Bernardi 2009; Bernardi et al. 2010).

We calculate the SDSS DR7 GSMF of early- and late-type galaxies using the $1/\text{V}_{\rm max}$ method described in Appendix A. Figure 2 shows the corresponding $\mbox{GSMF}$ s of early and late types in the upper left and right panels, respectively. For comparison, we show the $\mbox{GSMF}$ s for red and blue galaxies based on a colour cut limit in the ${(g-r)^{0.0}-M_{*}}$ diagram. In this diagram, we find that a rough division criteria from blue to red galaxies is given by the colour limit of ${(g-r)^{0.0}=0.66}$ .^k In the same figure, we compare our results to different determinations from the literature as we describe below. All the data in this figure have been renormalised to our cosmology.

Figure 2. SDSS DR7 GSMFs for early- and late-type galaxies, left and right upper panels, respectively. Early- and late-type galaxies are defined as those with $P(E)>0.65$ ( $P(E)\le0.65$ ) from the tabulated probabilities of Huertas-Company et al. (2011). This is equivalent to morphological types that comprises E and S0 galaxies or $T\le0$ (Sa to Irr galaxies or $T>0$ ). We compare to various previous determinations from the literature as indicated by the legends, see also the text for details. Our determinations are in general in good agreement with previous determinations from SDSS spectroscopic samples, while a tension is evident with determinations from the GAMA survey. We also present our resulting GSMFs for blue and red galaxies. These GSMFs follow closely those by morphology from the GAMA survey. The bottom panel shows our number density-weighted fractions of early-type and red galaxies as a function of $M_{*}$ . Their corresponding best-fit models Equation (25) are shown with solid and dashed lines, respectively.

Recently, Moffett et al. (2016a) visually classified morphologies in the GAMA survey and reported the GSMF for different morphologies. Here, we reproduce their GSMF from E to Sa galaxies as early types and the complement as late types. Contrary to our definition, Sa galaxies are included in the early-type group; this is because the authors report S0 and Sa galaxies as one morphology group. As shown in Figure 2, the GSMF of early-type galaxies from Moffett et al. (2016a) results in an overabundance of low-mass galaxies compared to other studies. We reproduce the results from Thanjavur et al. (2016) with bulge-to-total ratios of $B/T>0.8$ as early types and $B/T\leq0.8$ as late types. Thanjavur et al. (2016) used the bulge-to-disc decomposition from Simard et al. (2011) SDSS DR7 spectroscopic sample, and stellar masses derived from Mendel et al. (2014). We also include results from Kelvin et al. (2014). Similarly to Moffett et al. (2016a), Kelvin et al. (2014) visually classified morphologies in the GAMA survey. We again use their GSMF from E to Sa galaxies for early types since the authors combined S0–Sa galaxies as in Moffett et al. (2016a). The filled triangles with error bars show the GSMF from Bernardi et al. (2010) for galaxies with concentration parameter $c>2.86$ for early types and $c\leq2.86$ for late types.^l

Finally, using the Nair & Abraham (2010) morphology catalogue, who visually classified 14,034 spectroscopic galaxies from the SDSS DR4, we derive the GSMF for early-type galaxies.^m We utilise their morphological notation and define early-type galaxies as those objects with $-5\le T\le0$ (E-S0s), equivalent to $T\le1$ in the Fukugita et al. (2007) notation. We additionally derive the GSMF with morphologies between $-5\le T\le1$ in the Nair & Abraham (2010) notation, which include Sa galaxies.

In general, our results agree with previous determinations, especially with those from the SDSS spectroscopic samples. In contrast, the GSMF of early-type galaxies from the visual classification of the GAMA survey are systematically above our results at the low-mass end, $M_{*}\lesssim2\times10^{10}\text{M}_{\odot}$ , but closer to our classification based on galaxy colour. While it is not clear the reason of the differences outlined above (the inclusion of Sa galaxies as early types, environment, etc.), in Appendix C.4, we will discuss the impact of using galaxy colour instead of morphology when deriving the $\rm HI$ , $\rm H_{2}$ , cold gas, and baryonic MFs separated into two main galaxy populations.

Finally, the bottom panel of Figure 2 shows the resulting fraction of early-type galaxies as a function of stellar mass, $f_E\text{(}M_{*}\text{)}$ . In addition, we show the fraction of red galaxies when using our $g-r$ colour cut limit, $f_r\text{(}M_{*}\text{)}$ . We find that the fraction at which early-type galaxies is 50% is at $M_{*}\sim\!10^{11}\text{M}_{\odot}$ , while at $M_{*}\sim\!1.6\times10^{10}\text{M}_{\odot}$ and $M_{*}\sim\!8\times10^{11}\text{M}_{\odot}$ the fractions are 16% and 84%, respectively. For red galaxies, the fraction of 50% is at $M_{*}\sim\! 10^{10}\text{M}_{\odot}$ , while at $M_{*}\sim\!3\times10^{9}\text{M}_{\odot}$ and $M_{*}\sim\!10^{11}\text{M}_{\odot}$ the fractions are 16% and 84%, respectively. Note that the characteristic mass at which the fraction of early-type galaxies is 50% is a factor of $\sim\!10$ larger than for red galaxies. In general, $f_E\text{(}M_{*}\text{)}$ rises slower than the fraction $f_r\text{(}M_{*}\text{)}$ . In the same figure, we present the best-fit model to the data. After exploring different functions, we find that two sigmoid functions accurately describe the functionality of $f_E\text{(}M_{*}\text{)}$ or $f_r\text{(}M_{*}\text{)}$ :

(25) \begin{equation}f_{k}\text{(}M_{*}\text{)} = \frac{1-A}{1+e^{-\gamma_1(x_{C,1}+x_{0,1})}}+ \frac{A}{1+e^{-\gamma_2(x_{C,2}-x_{0,2})}},\end{equation}

where $k=E$ or r, $x_{C,i} = {M_{*}} / \mathcal{M}_{C,i}$ , with $i = 1,2$ . The best-fit parameters for the two fractions are listed in Table 2.

To derive the analytic model for the GSMF of early- and late-type galaxies, we use the best-fit model to our GSMF, Section 3.2, and the best-fit model for $f_E\text{(}M_{*}\text{)}$ . For masses below $5\times 10^8\ \text{M}_{\odot}$ we extrapolate $f_E\text{(}M_{*}\text{)}$ . This is an acceptable approximation since as seen in Figure 2, the fraction $f_E\text{(}M_{*}\text{)}$ tends to $\sim\!0$ below $M_{*}=10^9\ \text{M}_{\odot}$ . Recall that our main goal in this paper is to derive the MFs for $\rm HI$ , $\rm H_{2}$ , cold gas, and baryons by combining the observed $\rm HI$ and $\rm H_{2}$ CPDFs with the GSMF, both for early- and late-type galaxies, over a large mass range. Thus, at this point, we are in a position to determine these MFs.

4. Results

In this section, we present our fits to the $\rm HI$ - and $\rm H_{2}$ -CPDFs for early- and late-type galaxies from Paper I, the corresponding correlations (first and second moments), the bivariate $\rm HI$ - and $\rm H_{2}$ -stellar mass distributions, and the $\rm HI$ and $\rm H_{2}\,$ MFs. We also present the total cold gas and baryonic MFs. We will show that our empirically inferred $\rm HI$ and $\rm H_{2}\,$ MFs agree with direct determinations from blind or optically/infrared (selected) radio galaxy surveys. Previous works related to our approach are Obreschkow & Rawlings (2009), Lemonias et al. (2013), and Butcher et al. (2018).

Table 2. Best-fit parameters to the fraction of early-type and red galaxies

Table 3. Best-fit parameters of the $\rm HI$ and $\rm H_{2}$ mass CPDFs for late- and early-type galaxies

For those interested in using our results, we provide a Python code containing all the necessary information to reproduce the results presented here, for details see Section 6.

4.1. The $\text{\em HI}$ and $\text{\em H}_{\text{\em 2}}$ conditional distribution functions

Section 2.2 describes the functional forms for the $\rm HI$ - and $\rm H_{2}$ -CPDFs of early- and late-type galaxies proposed in Paper I. Here, using the determinations of the CPDFs for early- and late-type galaxies from Paper I, we find the best-fit parameters of the proposed functional forms: a Schechter-type function and a Schechter-type + Uniform function, respectively (see Section 2.2). While Paper I reported their corresponding best-fit parameters, here we opt for an update based on our own determinations of the GSMFs, for consistency. There are two reasons for doing this: (i) when fitting a CPDF that is determined over stellar mass bins, one should take into account contributions to this CPDF from the different masses. Weighting the conditional distributions by the GSMF takes care on that, see Equations (6) and (8) of Section 5 from Paper I; and (ii) Paper I used the fraction of bulge-dominated galaxies from Moffett et al. (2016b) as a proxy to the fraction of early-type galaxies. As discussed in Section 3.3, the results from Moffett et al. (2016a), and thus Moffett et al. (2016b), overestimate the fraction of early-type galaxies compared to the SDSS morphological catalogues. The above could be due to the inclusion of Sa galaxies into the group of early types. We used the above to argue in favour of our derived fraction of early-type galaxies based on the automated morphological classification from Huertas-Company et al. (2011).

Following Paper I, we use the Bayesian approach described previously through a MCMC method applied jointly to all the data (the CPDFs in different $M_{*}$ bins) to find the best-fit parameters of the proposed functions. These are listed in Table 3. Figure 3 shows our best-fitting models for late-type galaxies compared to the CPDFs from Paper I. Figure 4 shows the same but for early-type galaxies. We notice that our best-fit parameters are very similar to those determined in Paper I.

Figure 3. $\text{HI}$ and $\text{H}_{\text{2}}$ mass CPDFs for late-type galaxies. The results for the compilation sample from Paper I are shown as filled circles with error bars. Note that the above results include non-detections since the authors used the Kaplan & Meier (1958) estimator for uncensored data in their analysis. Our best-fitting models are shown as the solid lines.

Figure 4. Same as Figure 3 but for early-type galaxies. Note that the CPDFs of early-type galaxies reported in Paper I account for upper limits corrected by distance selection effects when necessary and the treated with the Kaplan & Meier (1958) estimator, see Section 2.2.1.

Figure 5. Logarithmic, left panels, and arithmetic, right panels, averaged mass ratios $\mathcal{R}_j$ as a function of $M_{*}$ from our analysis, with $j={\text{HI}}$ , ${\text{H}_{\text{2}}}$ . Blue and red lines are for early- and late-type galaxies, respectively, while the black lines correspond to all galaxies. The shaded areas show the respective standard deviations. Notice that $\log \langle \mathcal{R}_j (M_{*})\rangle \geq \langle \log \mathcal{R}_j (M_{*})\rangle$ and the dispersion reduces for the arithmetic mean. The open circles with error bars in the upper left panel correspond to the data from ALFALFA galaxies with SDSS spectral and stellar mass counterparts according to Maddox et al. (2015).

4.2. The $\text{\em HI}$ - and $\text{\em H}_{\text{\em 2}}$ -to-stellar mass correlations

Next, we explore the resulting first and second moments from our best-fitting models to the observed $\rm HI$ - and $\rm H_{2}$ -CPDFs, shown in Figure 5. The left panels of the figure present the logarithmic mean $\langle \log \mathcal{R}_j\rangle$ and its corresponding standard deviation, $\sigma_{\log \mathcal{R}_j}$ , $j=$ $\rm HI$ or $\rm H_{2}$ , as a function of $M_{*}$ for early- and late-type galaxies as well as for all the galaxies. At low masses, the correlation of all galaxies approaches the one of the late-type galaxies, while at high-mass end it approaches early types. The above trends are just the consequence of the observed fraction of early and late types. Figure 5 shows that early- and late-type galaxies follow different $\langle \log \mathcal{R}_j\rangle-M_{*}$ correlations. Therefore, due to the strong bimodality of these correlations, conclusions based on some subset of galaxies as representative of all galaxies would lead to incorrect results.

In the literature, sometimes the gas-to-stellar mass relations are reported using the arithmetic mean (though the results are plotted in logarithmic diagrams). The right panel of Figure 5 shows $\log \langle \mathcal{R}_j\rangle$ versus $M_{*}$ from our empirical CPDFs. As is clearly seen, there are notable differences when computing different ways of averaging the distributions: (i) $\log \langle \mathcal{R}_j\text{(}M_{*}\text{)}\rangle > \langle \log\mathcal{R}_j\text{(}M_{*}\text{)}\rangle$ , being larger the difference for the early-type galaxies;ⁿ and (ii) the standard deviations from the arithmetic mean is smaller than from the logarithmic mean.

In the left upper panel of Figure 4, we reproduce the results from (Maddox et al. 2015) for the ALFALFA galaxies with SDSS spectral and stellar mass counterparts. It is clear that the ALFALFA survey is biased towards galaxies with high $\rm HI$ -to– $M_{*}$ ratios. In other words, the ALFALFA survey mainly detects galaxies in the upper envelope of the full distribution of $\mathcal{R}_{{\rm HI}}$ (see also Maddox et al. 2015) and is dominated mostly by late types (see also e.g., Haynes et al. 2011).

4.3. The bivariate mass distribution functions

Figure 6 shows the resulting bivariate stellar- $\rm HI$ mass distribution function, $\Phi(M_{\rm HI},{M_{*}})$ , see Equation (2). The colour code shows various number density levels for $\Phi(M_{\rm HI},{M_{*}})$ . Notice that $\Phi(M_{\rm HI},{M_{*}})$ is for all galaxies, that is, it includes the contribution from early- and late-type galaxies. The discontinuity in the number density isocontours at the bottom right of the diagrams is related to contributions from the non-detections from early types. Recall that for the CPDFs of early types, we assumed an uniform function (or top-hat) for the lowest values of the gas-to-stellar mass ratios $\mathcal{R}_{\rm HI}$ , where the non-detection piled up,^o see Section 2.2.2 and Figure 4. In the bottom right and upper panels of the same figure, we present respectively our measurements of the $\rm HI$ MF and GSMF with the solid black lines. We compare the $\rm HI$ MF with blind $\rm HI$ galaxy surveys based on ALFALFA (Jones et al. 2018; Papastergis et al. 2012; Martin et al. 2010) and on $\rm HI$ Parkes All Sky Survey (HIPASS) (Zwaan et al. 2005). While in the next subsection, we discuss in more detail the comparison with previous works, for the moment we note that our total $\rm HI$ MF is in good agreement with the above direct observations. In the case of ALFALFA, this is a revealing result given the strong selection bias of this survey towards $\rm HI$ -rich and late-type galaxies as seen in Figure 6 (open circles reproduce the results from Maddox et al. 2015, see also the discussion of the previous subsection and Figure 4). As we will discuss in the next section, the above reflects that the total $\rm HI$ MF is dominated by late-type galaxies.

Figure 6. Atomic gas–stellar mass bivariate (joint) distribution function. The colour code shows various number density levels as indicated by the legends. Due to the rising slope of the MFs at low masses, most of the galaxies are located at small $\text{HI}$ and stellar masses. Note that the discontinuity seen at the low- $\text{HI}$ and high-stellar masses is due to the assumption of an uniform function for the lowest values of gas-to-stellar mass ratios of early-type galaxies where non-detections piled up. Recall that in our analysis, non-detections (upper limits) are included using the non-parametric estimator Kaplan & Meier (1958) for censored data in Paper I. The solid lines show the mean $\langle\log M_{\rm HI}\rangle$ as a function of $M_{*}$ , both for early- and late-type galaxies. The upper panel shows the GSMF which is the result of integrating the bivariate distribution function along the $M_{\rm HI}$ axis, while the bottom right panel shows the same but for the $\text{HI}$ MF which results from integrating along the $M_{*}$ axis. We compare to some previous observational determinations of the MFs and the relationship between $M_{\text{HI}}$ and $M_{*}$ derived in Maddox et al. (2015) for the ALFALFA survey with SDSS spectral and stellar mass counterparts.

Figure 6 explicitly shows the contribution of galaxies of different stellar masses to the $\rm HI$ MF. Particularly, we observe that the low-mass end of the $\rm HI$ MF is composed mainly by low $M_{*}$ galaxies but there is also a non-negligible contribution from a population of high $M_{*}$ galaxies. Most of these high $M_{*}$ galaxies are early-type (quenched) galaxies for which there is a significant fraction of non-detections ( $\sim\! 55\%$ ). In Paper I, we included non-detections for the determination of the HI-CPDF based on methods of censored data (Kaplan & Meier 1958). Nonetheless, the contribution of non-detections is only marginal because the fraction of early-type galaxies at those masses is low, see the bottom panel of Figure 2. In addition, Figure 6 shows that the completeness limit in the $\rm HI$ MF, due to our stellar mass limit of $M_{*}=10^{7}\text{M}_{\odot}$ , is at $M_{\rm HI}\sim\! 10^{8}\text{M}_{\odot}$ (see below), which excludes a large region of galaxies with non-detections.

Figure 7. Same as Figure 6 but for molecular gas. Note that while there are more non-detections for $\text{H}_{\text{2}}$ observations, these are mostly from early-type galaxies that represent only a small fraction overall in the $\text{H}_{\text{2}}$ mass bivariate distribution function. We also compared to previous determinations from Keres et al. (2003).

Similarly to Figure 6, Figure 7 presents the resulting bivariate stellar– $\rm H_{2}$ mass distribution function for all galaxies and the mean $\langle\log M_{\rm H2}\rangle$ for early- and late-type galaxies. The resulting total $\rm H_{2}$ MF is shown with the solid line in the bottom right panel and compared to the Keres et al. (2003) $\rm H_{2}$ MF based on the CO luminosity function. At the low-mass end, there is a substantial population of galaxies with non-detections, roughly $\sim\! 78\%$ , which are mostly early-type galaxies. As above, non-detections have been included in the statistical analysis of the H2-CPDFs, for details see Paper I. Nevertheless, from the contour density level, their contribution is marginal. Finally, we can conclude that our $\rm H_{2}$ MF is complete for $M_{\rm H2}\gtrsim10^{7}\text{M}_{\odot}$ .

4.4. The mass functions

Next, we discuss in detail the MFs presented above. In particular, we focus on the determinations separately for early- and late-type galaxies based on the morphology classification described in Section 3.3. The various panels in Figure 8 present the MFs for atomic and molecular gas, total cold gas, and baryons, as indicated by the labels. In all the panels, the MFs for late-type galaxies are shown as blue filled circles with error bars, while for early-type galaxies are shown as red circles with error bars. These symbols are for a stellar mass limit of $M_{*}>10^{7}\ \text{M}_{\odot}$ . Instead, we use blue/red open circles when the MFs is incomplete. We also calculate the MFs in the hypothetical case of a mass limit of $M_{*}=0$ shown with the blue/red dashed lines. Note that the total MFs for $\rm HI$ and cold gas are not plotted. This is because in these cases the MFs of late-type galaxies are hardly distinguishable from the total one at all masses. We also omit the total $\rm H_{2}$ MF. This is because it is hardly distinguishable from the MF of late-type galaxies at $M_{\rm H2}\lesssim 2\times 10^{10}\ \text{M}_{\odot}$ while for larger masses it is indistinguishable from the MF of the early-type galaxies. In the case of the baryon masses, the total MF is plotted with black filled circles.

Figure 8. Results on the galaxy MFs of early- and late-type galaxies for atomic gas, left upper panel, molecular gas, right upper panel, cold gas, bottom left, and baryons, bottom right panel. In all the panels, late-type galaxies are shown as the blue circles with error bars, while early-type galaxies are shown as the red circles with error bars, when using a stellar mass limit of $M_{*}=10^{7}\text{M}_{\odot}$ . Filled blue/red circles indicate when the MFs are complete, while open circles clearly show that the MFs became incomplete. The dashed lines are for MFs when using a stellar mass limit of $M_{*}=0$ . The total MFs for $\text{HI}$ and cold gas are not shown because they are practically indistinguishable from the respective MFs of late-type galaxies. Our results are in good agreement with observational determinations of the total MFs. For only early-type galaxies, we compare our results with those from the ATLAS 3D sample (red triangles). While we observe some tension, we suspect that selection effects are more likely to artificially increase the amplitude of their MF at low masses.

We note that our determinations for the MFs are the result of the convolution between random errors and the intrinsic MFs, similarly as it happens with the direct observational determinations of MFs. In Section 5.1, we discuss the impact from random errors and present the intrinsic MFs, after deconvolving by these errors.

4.4.1. The $\text{\em HI}$ MF

We compare our results with previous direct observational determinations of the total $\rm HI\,$ MFs. According to our results, late-type galaxies dominate the $\rm HI$ MF, even at the highest masses, so that it is practically equivalent to the total $\rm HI$ MF. In Figure 8, we reproduce the best fit to the $\rm HI$ MF from Jones et al. (2018) based on the final catalogue of the blind $\rm HI$ ALFALFA Survey, dark grey solid line. The violet solid line shows the best-fit model reported in Butcher et al. (2018) to the $\rm HI$ MF from the Nançay Interstellar Baryons Legacy Extragalactic Survey, which is a project that complemented recent and/or ongoing large blind $\rm HI$ surveys. Results from Papastergis et al. (2012) and Martin et al. (2010) based on the 40% sample of the ALFALFA Survey are shown respectively with black open circles and squares. The skeletal symbols with error bars show the results from Zwaan et al. (2005), who used the HIPASS. Note that our $\rm HI$ MF for late-type galaxies, which dominates the total MF, is in good agreement with direct inferences from $\rm HI$ blind surveys, particularly those based on ALFALFA as discussed in Figure 6. As for the $\rm HI$ MF of early-type galaxies, we plot the determinations from the ATLAS 3D (Serra et al. 2012) and HIPASS (Lagos et al. 2014) samples shown as the red solid triangles and skeletal symbols, respectively.

Our resulting $\rm HI\,$ MFs are in good agreement with direct determinations from radio observations. This is particularly true for $\rm HI$ masses above the completeness limit corresponding to our $M_{*}$ limit for the GSMF. These masses are $M_{\rm HI} \sim\! 10^{8}\text{M}_{\odot}$ for late-type galaxies and $M_{\rm HI} \sim 10^{7}\text{M}_{\odot}$ for early-type galaxies. Even when extrapolating to a limit mass of $M_{*}=0$ , at the low-mass end, we find a good agreement with direct determinations, though the early-type galaxies from the ATLAS 3D sample present a higher amplitude for masses below $10^{7}\text{M}_{\odot}$ . However, those extrapolations should be taken with caution as it is not clear whether we can extrapolate our bivariate distribution functions to such low masses. In conclusion, we consider that the remarkable consistency between our $\rm HI$ MFs and radio blind surveys above $M_{\rm HI}\,{\sim}\,10^{8}\text{M}_{\odot}$ is reassuring and validates the HI-CPDFs determined in Paper I. Recall that the observational data used in that paper were derived from various heterogeneous samples, affected by many selection effects, including those related to the non-radio detections. Therefore, the agreement between the $\rm HI$ MF with that of the blind radio observations is far from trivial, unless adequate corrections are introduced and the data are adequately analysed from the statistical point of view.

4.4.2. The $\text{\em H}_{\text{\em 2}}$ MF

In the upper right panel of Figure 8, we present the results for the $\rm H_{2}$ MF. The $\rm H_{2}$ MF is largely dominated by late-type galaxies below $M_{H_{2}}\,{\sim}\,2\times10^{10}\text{M}_{\odot}$ , but for larger masses, early-type galaxies are more abundant. In the same panel, we reproduce the total $\rm H_{2}$ MF from Keres et al. (2003), who used a CO luminosity function from FIR and $B-$ band limited galaxy samples and adopting a constant CO-to- $\rm H_{2}$ conversion factor of $\alpha_{\rm CO} = 4.76$ , open black circles with error bars. The dashed line shows the best fit to a Schechter function derived in Obreschkow & Rawlings (2009). Additionally, we show the results from the ATLAS 3D sample for early-type galaxies (Lagos et al. 2014) with the filled triangles. The magenta solid line shows the results from Obreschkow & Rawlings (2009), who derived the $\rm H_{2}$ MF by introducing a phenomenological model for the $\rm H_{2}$ -to- $\rm HI$ mass ratio that depends on the galaxy morphological type and its total cold gas mass.

When comparing to the $\rm H_{2}$ MF from Keres et al. (2003), we observe a good agreement with our results. At the low-mass end, though, the Keres et al. (2003) MF seems to be slightly shallower than ours. It is not clear the origin for this discrepancy. One possibility is due the constant $\alpha_{\rm CO}$ factor employed by the authors. Based on previous empirical studies, Paper I showed that ignoring the dependence of $\alpha_{\rm CO}$ with $M_{*}$ ^p flattens the low-mass end of the $\rm H_{2}$ MF, consistent with the results from Keres et al. (2003). Another possibility is an effect of the incompleteness of the CO luminosity function. As for Obreschkow & Rawlings (2009), our results are consistent for masses below $M_{H_{2}}\,{\sim}\,3\times10^{9}\text{M}_{\odot}$ . For larger masses, the Obreschkow & Rawlings (2009) MF falls much stepper than ours. Similar to our analysis, Obreschkow & Rawlings (2009) used the relationships between galaxy properties to derive their MF. As mentioned above, their phenomenological model employed the dependence of the $\rm H_{2}$ -to- $\rm HI$ mass ratio with morphology and cold gas mass. While the above differences could be arguably referred to the use of different estimators for the $\rm H_{2}$ gas masses, it could be also an indication that random errors are larger when using only one galaxy parameter. Recall that in this paper, we are using $M_{*}$ . In Section 5.1, we will show that after deconvolving from random errors, our intrinsic $\rm H_{2}$ MF becomes steeper at the high-mass end and more consistent with the derivation from Obreschkow & Rawlings (2009). Nonetheless, it is difficult to conclude the origin of the above differences given the different nature of the models employed in both studies.

As for early-type galaxies, our results are consistent with those from the ATLAS 3D (Lagos et al. 2014) at the high-mass end but they lie slightly below at the low-mass end. It is unclear the reason of the above discrepancy for low- $M_{\rm H_2}$ early-type galaxies, though large-scale and environmental selection effects could boost the inferences of the MF when using the 1/V $_{\rm max}$ estimator, see for example, Appendix C and Baldry et al. (2008). Recall that in the case of $\rm HI$ , the ATLAS 3D also presents an slightly enhancement at the low-mass end of the $\rm HI$ MF. Thus, selection effects are more likely to artificially increase the amplitude of the MF for low-mass galaxies in the ATLAS 3D sample.

4.4.3. The gas and baryonic MFs

Similarly to the $\rm HI$ MF, the cold gas MF is completely dominated by late-type galaxies, even at the high-mass end. In Figure 8, we compare our results with the phenomenological determination from Obreschkow & Rawlings (2009, pink solid line). These authors combined their inference of the $\rm H_{2}$ MF with the $\rm HI$ MF from Zwaan et al. (2005) to derive the gas MF. Despite the differences mentioned above for the Obreschkow & Rawlings (2009) $\rm H_{2}$ MF, their total cold gas MF is in excellent agreement with our one. This is not surprising as it is just reflecting that $\rm HI$ is much more abundant than $\rm H_{2}$ .

Finally, we show our results for the baryonic MFs in the bottom right panel of Figure 8. The baryonic MF is very similar to the GSMF at the high-mass end but at the low-mass end, it is steeper as the contribution of cold gas becomes more relevant. On the other hand, late-type galaxies dominate the baryonic MF for $M_{\rm bar}<10^{11}\ \text{M}_{\odot}$ , while at the high-mass end, early-type galaxies are more abundant than the late-type ones. We reproduce with the green solid line the baryonic MF from Baldry et al. (2008). These authors combined the GSMF from the low-z survey of the SDSS DR4, the same galaxy survey used here for low masses, with a closed-box model and the mass–metallicity relation to derive cold gas masses for their baryonic MF. The open black circles show the MF from Papastergis et al. (2012), who defined baryonic mass as $M_{\rm bar}=1.4 \times M_{\rm HI}+M_{*}$ . We notice that these previous baryonic MF determinations are in good agreement with our results at the mass range $\sim\! 2\times 10^9-2\times 10^{11}\ \text{M}_{\odot}$ , while for lower and higher masses there are some differences, which are easy to understand.

The MF from Baldry et al. (2008) is steeper than our MF at low masses. This might be a consequence of the fact that the Baldry et al. (2008) GSMF itself is steeper compared to other determinations, for example, Baldry et al. (2012). As these authors discuss, the disagreement between the Baldry et al. (2008) and Baldry et al. (2012) $\mbox{GSMF}$ s is just the result of different flow models for distances, which affect significantly to nearby low-mass galaxies. Recall that our GSMF has been corrected to be consistent with the flow model of Tonry et al. (2000). Additionally, the gas masses in Baldry et al. (2008) were actually obtained from a close-box model constrained with the empirical mass–metallicity relation. The combination of these two assumptions are likely the result of a steep baryonic MF at low masses, which differs from our results and those of Papastergis et al. (2012).

Regarding the high-mass end, our baryonic MF falls shallower than those of Papastergis et al. (2012) and Baldry et al. (2008). This is because, our GSMF is shallower. As discussed in Section 3, there are two main systematic effects that could lead to different GSMFs, mass-to-light ratios, and the determination of galaxy surface brightness (especially due to sky subtraction problems). Both effects are likely to affect the high-mass end of the baryonic MF. In addition, due to the small volumes of the surveys used in Baldry et al. (2008) and Papastergis et al. (2012), cosmic variance enhances the differences.

Table 4. Cosmic density of $\text{HI}$ , $\text{H}_{\text{2}}$ , gas, stars, and baryons for all, LTGs and ETGs. The fraction of each component is denoted as $f_j = \Omega_j / \Omega_{\rm bar,U}$ with $\Omega_{\rm bar,U} = 0.048$ .

LTG = late-type galaxies; ETG = Early-type galaxies

Figure 9. Density parameter $\Omega$ of $\rm HI$ , $\rm H_{2}$ , cold gas, and baryonic mass locked in all galaxies as well as in early- and late-type galaxies (coloured filled circles; the errors are smaller than the circle size). The $\Omega$ parameter values are reported as fractions in per cent of the universal matter (left axis) and baryonic (right axis) densities. The grey boxes show the range of values from previous determinations and the horizontal lines correspond to the mean of these values.

4.5. Cosmic density parameters and relevant timescales

4.5.1. Cosmic density parameters

The cosmic density parameter measures the average mass density of some matter species with respect to the critical density, $\rho_c$ . Here, we determine the mass density in stars, $\rm HI$ , $\rm H_{2}$ , cold gas, and baryons that are locked inside galaxies using their MFs. The differential mass density function $d\rho_j(M_j)$ for some galaxy mass component $M_{j}$ in the mass range $\log M_j\pm d\log M_j /2$ is $d\rho_j(M_j ) = M_{j} \times \phi_j (M_j ) d\log M_j$ , where $\phi_{j}$ is in units of Mpc-3 dex-1. Thus, the cosmic mass density is given by:

(26) \begin{equation}\rho_j = \int_{-\infty}^{\infty} d\rho_j(M_j ),\end{equation}

with the cosmic density parameter

(27) \begin{equation}\Omega_j = \frac{\rho_j}{\rho_c},\end{equation}

where the critical density is $\rho_c = 2.775\times 10^{11} h^{-1}\text{M}_{\odot}/ (h^{-1}{\rm Mpc})^3 = 1.2756 \times 10^{11} \text{M}_{\odot}/ {\rm Mpc}^3 h_{67.8}^{2}$ .^q The limits of integration in Equation (26) reflect that we are considering all the spectrum of masses from galaxies. In reality, this is not possible, due to completeness limits in galaxy samples. Here, we consider the following mass limits for all our components: $M_{\rm inf} = 10^7\ \text{M}_{\odot}$ and $M_{\rm upper} = 10^{12.6}\ \text{M}_{\odot}$ . We notice that using values smaller than $M_{\rm inf}$ and/or larger than $M_{\rm upper}$ do not substantially change our results. This is because the multiplicity functions, $\propto M_{j} \times \phi_j (M_j )$ , have a maximum around the knee of the MFs.

Figure 9 shows the different values of each $\Omega_{j}$ corresponding to all galaxies and separately for early- and late-type galaxies, listed in Table 4. The $\Omega_j$ values are presented as the fractions in per cent with respect to the universal matter density ( $\Omega_{\rm m}=0.307$ , left axis) and the universal baryonic density ( $\Omega_{\rm bar,U}=0.048$ , right axis). To estimate errors in our cosmic density parameters, we use all our MCMC models for the HI-CPDF and H2-CPDF in addition of all our MCMC models to the fit of the GSMFs. We found that the largest uncertainties arise from the uncertainties in the CPDFs.

In the past, there have been some assessments of the cosmic density parameters at $z\,{\sim}\,0$ . Usually, these studies do not report cosmic density parameters for different populations and for different components at the same time. As mentioned in the Introduction, it is important to distinguish between different populations given that late- and early-type galaxies follow different formation histories. Studies close to ours are the ones by Fukugita, Hogan, & Peebles (1998), Fukugita & Peebles (2004), and Read & Trentham (2005). Below, we present and compare our results with many previous determinations from the literature.

• $\rm HI$ cosmic density: The atomic hydrogen in late-type galaxies is $\sim\! 27$ times larger than in early-type galaxies, which means that $\sim\!96\%$ of $\rm HI$ mass is in late-type galaxies. Previously, Zwaan et al. (2005), Read & Trentham (2005), Martin et al. (2010), Braun (2012), Delhaize et al. (2013), Hoppmann et al. (2015), Butcher et al. (2018), and Jones et al. (2018) have derived the $\rm HI$ cosmic density parameter using either blind $\rm HI$ galaxy surveys (HIPASS and ALFALFA) or indirect techniques. The mean value from these determinations is $\Omega_{\rm HI} = 4.2\times10^{-4}$ with a lower bound of $\Omega_{\rm HI}^- = 3\times10^{-4}$ and an upper bound of $\Omega_{\rm HI}^+ = 6.2\times10^{-4}$ , shown as the grey box in Figure 9^r. Recently, using a spectral stacking technique and from Westerbrook Synthesis Radio Telescope observations of 1895 galaxies crossed with the SDSS, Hu et al. (2019) found $\Omega_{\rm HI} = (4.15\pm 0.26)\times10^{-4}$ . Our determined value is in good agreement with these previous determinations, in particular with the latter one.

• $\rm H_{2}$ cosmic density: The molecular hydrogen cosmic density in late-type galaxies is $\sim\!21$ times larger than in early-type galaxies. This implies that $\sim\!95\%$ of $\rm H_{2}$ mass is in late-type galaxies. Using the CO surveys from Young et al. (1995), Keres et al. (2003) determined that $\Omega_{\rm H2}=(1.64\pm0.63)\times10^{-4}$ , while from the observations in CO from Maeda, Ohta, & Seko (2017) they report $\Omega_{\rm H2}=0.51\times10^{-4}$ . Obreschkow & Rawlings (2009) used a phenomenological model to derive $\Omega_{\rm H2}= (1.01\pm0.39)\times10^{-4}$ . Read & Trentham (2005) find that $\Omega_{\rm H2}=2.68\times10^{-4}$ . The above ranges of values are shown with grey box in Figure 9. As can be seen, our results are consistent with the range of determinations described above, especially with the results from Obreschkow & Rawlings (2009).

• Cold gas cosmic density: Most of the cold gas is located in late-type galaxies, $\sim\!96\%$ . Keres et al. (2003) found that $\Omega_{\rm gas}=(6.34\pm1.62)\times 10^{-4}$ , which includes the resulting abundance of $\rm HI$ from Zwaan et al. (1997). Obreschkow & Rawlings (2009) used their best phenomenological model to the $\rm H_{2}$ -to- $\rm HI$ ratio with the HIPASS results from the Zwaan et al. (2005) sample to derive $\Omega_{\rm gas}=(6.49\pm 1.18)\times 10^{-4}$ , while using the values for $\rm HI$ and $\rm H_{2}$ masses from Read & Trentham (2005), we calculate that $\Omega_{\rm gas}=7.95\times 10^{-4}$ after correcting from helium and heavier metals. Our value of $\Omega_{\rm gas}=(6.85\pm 0.92)\times 10^{-4}$ is consistent with the above results.

• Stellar cosmic density: The stellar cosmic density in late-type galaxies that we derive from the SDSS is approximately $\sim\!1.8$ larger than in early-type galaxies. Thus, $\sim\!64\%$ of the mass in stars at $z\,{\sim}\,0.1$ is in late-type galaxies. From the compilation by Madau & Dickinson (2014), the stellar cosmic density lies within $\Omega_{\ast}=(28.06-17.71)\times 10^{-4}$ , while the derivations from Wright et al. (2017) and Baldry et al. (2008) are respectively $\Omega_{\ast}=17.14\times 10^{-4}$ and $\Omega_{\ast}=29.73\times 10^{-4}$ . Our result for the cosmic density for all galaxies, $\Omega_{\ast}=(20.20\pm0.08)\times 10^{-4}$ , is consistent with the above values.

• Baryonic cosmic density: Finally, we find that there is a factor of $\sim\!2.4$ more baryons in late-type galaxies than in early types, and thus, $\sim\!71\%$ of the baryons are in late-type galaxies. Read & Trentham (2005) found that $\Omega_{\rm bar} = 35\times10^{-4}$ which is a factor of $\sim\!1.3$ larger than our results. We find a cosmic density parameters ratio of $\Omega_{\ast} / \Omega_{\rm bar} \approx 1.3$ . Finally, our baryon density parameter is $\approx 5.4\%$ of the universal baryon fraction, $f_{\rm bar,U}=0.156$ , or equivalently $\sim\! 18$ times lower than $f_{\rm bar,U}$ . Most of the baryons are definitively not locked inside galaxies.

4.5.2. Cosmic timescales

We are now in a position to derive some relevant cosmic timescales, such as the mean galaxy depletion times. We focus on late-type galaxies because most of the star formation occurs in those galaxies. To do so, we use the observed cosmic star formation rate (CSFR) at $z\,{\sim}\,0.1$ . According to Madau & Dickinson (2014), who derived the CSFR from far-UV and IR rest frame luminosities, the CSFR is $\dot{\rho}_\ast \sim\! 90\times10^{-4}\ \text{M}_{\odot}$ $\textrm{yr}^{-1}\ \textrm{Mpc}^{-3}$ after correcting to a Chabrier (2003) IMF. Unfortunately, the authors report the CSFR for all galaxies but not divided into different groups. The recent study by Sánchez et al. (2019), based on the fossil record analysis of a sample of more than $4\times10^{4}$ galaxies from the SDSS MaNGA survey, report similar values to the the total CSFR of $\dot{\rho}_\ast = 114.82\pm67.61\times10^{-4}\ \text{M}_{\odot}$ $\textrm{yr}^{-1}\ \textrm{Mpc}^{-3}$ or $\dot{\Omega}_\ast = \dot{\rho}_\ast/\rho_{\rm crit} = (9\pm 5)\times10^{-14}$ yr^–1 corrected to a Chabrier (2003) IMF. The authors also derived the CSFRs for star-forming galaxies; $\dot{\Omega}_{\ast, \rm SFG} = (6.5\pm 3.8)\times10^{-14}$ yr^–1. In the following, we use their value for star-forming galaxies as a representative determination for late-type morphologies, that is, $\dot{\Omega}_{\ast, \rm L}\approx \dot{\Omega}_{\ast, \rm SFG}$ .

We begin our discussion by estimating the mean molecular hydrogen depletion time of late-type galaxies, $\bar{t}_{\rm dep,L}({\rm H_{2}})=\Omega_{\rm H_2, L}/\dot{\Omega}_{\rm \ast,L}$ . The $\rm H_{2}$ depletion timescale is defined as the time at which a galaxy (or a molecular cloud) would consume its $\rm H_{2}$ gas reservoir by forming stars at the current SFR. From our cosmic density parameters, we find that $\bar{t}_{\rm dep,L}\text{(}{\rm H_{2}}\text{)}\approx 1.3$ Gyrs. This is consistent with the mean depletion time $\bar{t}_{\rm dep}\text{(}{\rm H_{2}}\text{)}=0.96$ Gyr reported in Saintonge et al. (2017) for star-forming galaxies in a volume complete sample. Note, however, that for local individual galaxies the molecular depletion time could vary from $\sim\! 0.9$ to 3 Gyrs (e.g., Kennicutt 1998; Bigiel et al. 2008; Leroy et al. 2008; 2013). Also, we estimate the mean total cold gas depletion time of late-type galaxies, $\bar{t}_{\rm dep,L}({\rm gas})=\Omega_{\rm gas,L} / \dot{\Omega}_{\rm \ast,L}$ , and find $\bar{t}_{\rm dep,L}({\rm gas})\approx 10.14$ Gyrs, that is, $\sim\!8$ times larger than for the molecular gas component. The values we find for these two timescales are consistent with the proposal that, on average, for local late-type galaxies, (i) the global conversion of molecular gas into stars is inefficient (recall that the $\rm H_{2}$ depletion times of observed local star-forming regions are actually 40–500 Myr, for example, Lada, Lombardi, & Alves 2010; Lada et al. 2012); and (ii) the global conversion of atomic to molecular hydrogen gas is also inefficient, or equivalently, the molecular cloud formation efficiency is low. Thus, the mean star formation efficiency, SFE, of local late-type/star-forming galaxies is low despite the fact that they contain a considerable amount of interstellar gas; according to Table 4, on average approximately 36% of the baryons in these galaxies are in form of cold gas.

According to Leroy et al. (2008), the SFE of a galaxy is the inverse of the neutral H gas depletion time, that is, the time required for current star formation to consume the neutral H reservoir. The SFE can be estimated as the product of the net efficiency of converting $\rm H_{2}$ gas into stars, ${\rm SFE}\text{(}{\rm H_{2}}\text{)} =1 / {t}_{\rm dep}\text{(}{\rm H_{2}}\text{)}$ , and the net efficiency of molecular cloud formation given by the mass fraction of $\rm H_{2}$ with respect to the total neutral H mass in the galaxy, that is, $M_{H_{2}}/\text{(}M_{\rm HI} + M_{H_{2}}\text{)}$ . Thus, using our estimations of the cosmic parameters for late-type galaxies, we calculate the cosmic (mean) SFE of late-type galaxies as:

(28) \begin{align}{\rm SFE}_L({\rm H}) &= {\rm SFE}_L\text{(}{\rm H_{2}}\text{)} \times \frac{\Omega_{\rm H_2,L}}{\Omega_{\rm HI,L}+\Omega_{\rm H_2,L}}\nonumber\\[3pt] &= \frac{1}{\bar{t}_{\rm dep,L}\text{(}{\rm H_{2}}\text{)}} \times \frac{\Omega_{\rm H_2,L}}{\Omega_{\rm gas,L}/1.4} = 1.4 \times \frac{\dot{\Omega}_{\ast, \rm SF}}{\Omega_{\rm gas,L}} \nonumber\\[5pt] & =\frac{1.4}{\bar{t}_{\rm dep,L}({\rm gas})} = 1.38\times 10^{-10} {\rm yr^{-1}}. \end{align}

The inverse of this efficiency is the cosmic neutral H depletion time, $\bar{t}_{\rm dep,L}({\rm H})\approx7.25$ Gyrs. Note that the relationship between the SFE based on the neutral H gas reservoir and the SFE based on the total cold gas reservoir is ${\rm SFE}_L({\rm H}) = 1.4\times {\rm SFE}_L({\rm gas})$ , or equivalently, $\bar{t}_{\rm dep,L}({\rm H})= \bar{t}_{\rm dep,L}({\rm gas})/1.4$ . The factor 1.4 takes into account He and metals.

We calculate also the cosmic SF timescale of late-type galaxies, which is given by the inverse of the cosmic-specific SFR, $\bar{t}_{\rm SF,L}= \Omega_{\rm \ast,L}/\dot{\Omega}_{\rm \ast,L}\approx 20.3$ Gyrs; this is a factor of $\sim\!1.5$ larger than the present age of the Universe. The cosmic SF timescale can be understood as the time required for the current cosmic SFR density to double the current cosmic stellar mass content. Interestingly enough, the ratio $\bar{t}_{\rm dep,L}({\rm H})/\bar{t}_{\rm SF,L} = (\Omega_{\rm H2} + \Omega_{\rm HI}) / \Omega_{\ast} = 0.36$ , that is, the gas reservoir of late-type galaxies has not yet been dramatically consumed by star formation. Including He and metals in the gas reservoir, this ratio increases to $\sim\! 0.5$ .

5. Discussion

In this paper, we employed a statistical approach that allows to project the observed $\rm HI$ - and $\rm H_{2}$ -CPDFs into their corresponding MFs, when using the GSMF as an interface or pivotal function, Section 2. Additionally, the cold gas and baryon MFs are obtained from the above. Our empirical approach makes use of the following observational data as input:

1. 1. The local GSMF over a large dynamical range and separated into early- and late-type galaxies.

2. 2. The observed CPDFs of $\rm HI$ and $\rm H_{2}$ as a function of $M_{*}$ , both for early- and late-type galaxies.

As a result, our approach provides a fully self-consistent and complete empirical description of the demographics of the local population of early- and late-type galaxies for a broad mass range. Furthermore, by construction, our MFs are derived separately for early- and late-type galaxies. As discussed in Section 4.4, our $\rm HI$ and $\rm H_{2}\,$ MFs are actually consistent with several previous determinations from radio blind or optically/infrared (selected) galaxy samples. Actually, the above level of agreement is not trivial due to the chain of assumptions and corrections for the data sets we used here and in Paper I, and it reinforces the robustness of the observational information employed. Note, however, that the above agreement is only valid above our completeness limit for the GSMF of $M_{*}=10^{7}\text{M}_{\odot}$ , which corresponds to a completeness limit of $M_{\rm HI}\sim 10^{8}\text{M}_{\odot}$ and $M_{H_{2}}\,{\sim}\,10^{7}\text{M}_{\odot}$ , respectively, for the $\rm HI$ and $\rm H_{2}\,$ MFs. In that regards, we are unable to constraint the very low-mass end of the $\rm HI$ and $\rm H_{2}\,$ MFs.

Below, we highlight aspects that we consider are relevant for the success of our empirical approach:

• The $\rm HI$ and $\rm H_{2}$ CPDFs for early- and late-type galaxies. We used the CPDFs from Paper I, where we derived the CPDFs from a compilation of many incomplete and inhomogeneous samples, carefully homogenised to a common IMF, cosmology, CO-to-luminosity conversion factor, and accounting for selection biases.

• The effect from upper limits in radio surveys. In addition to the above-mentioned homogenisation and corrections, it was important to take into account the upper limits reported in the original sources, when radio detections were not achieved. The fraction of non-detections in the compilation from Paper I was non negligible, especially for early-type galaxies. Non-detections were corrected by distance-sensitivity effects. Instead of ignoring radio non-detection or using the upper limits as the true values, as is commonly done in the literature, we derived the CPDFs by including them in our statistical analysis based on the non-parametric Kaplan–Meier estimator for censored data.

Next, we briefly discuss below some potential caveats on our approach. Over the next subsections, we will discuss them in more detail and show that they do not affect our main conclusions.

• The assumption that the $\rm HI$ and $\rm H_{2}$ masses are two independent random variables. In reality, this is not true, for example, Obreschkow & Rawlings (2009) showed that the $\rm H_{2}$ -to- $\rm HI$ mass ratio depends on the morphological type. Note, however, that this was partially taken into account in our approach in a statistical sense. Recall that we use the observed mass CPDFs separately for early- and late-type galaxies, that is, the dependence with morphology is roughly included, as shown in Paper I, see Figure C3 from that paper and Figure 5.

• Differences on the mass-to-light ratios. Figure A1 in Appendix A shows that the different mass-to-light ratios used to estimate $M_{*}$ lead to different $\mbox{GSMF}$ s, with differences up to $0.5--1$ dex in number densities at the high-mass end (see also Bernardi et al. 2017). While we choose to use the geometric mean over the five mass-to-light ratios described in Appendix A, one could naturally question that the agreement of our $\rm HI$ and $\rm H_{2}\,$ MFs with the observed ones is relative because using a different GSMF could result in different MFs. In Section 5.1, we explore and quantify the impact of systematics from varying mass-to-light ratios and show that its effect is marginal in the obtained $\rm HI$ and $\rm H_{2}\,$ MFs.

• Random errors from stellar mass estimates. Inevitably, random errors propagate to our MFs resulting in a Eddington (1940) bias effect. Thus, the comparison with the results based on radio surveys is not trivial as they do not suffer of an Eddington (1940) bias effect due to $M_{*}$ errors. Nonetheless, measurements of the $\rm HI$ and $\rm H_{2}$ masses are also subject to random errors. In Subsection 5.1.1, we deconvolve our MFs with the random errors, not only as a method to compare with results from radio surveys but also for obtaining the intrinsic MFs to be used to constrain the predictions from galaxy formation models.

• The morphological classification from the SDSS DR7. To derive our GSMF separated into early- and late-type galaxies, we used the morphological classification based on the Huertas-Company et al. (2011) vector machine analysis of the SDSS DR7. As shown in Figure 2, we find that the obtained early- and late-type GSMFs using this classification are consistent with other determinations for the SDSS but disagree with those based on the visual classification from the GAMA survey. While we explore in detail this effect in Appendix C.4, we do not include it as one of the main sources of uncertainty.

We conclude this section by emphasising the robustness of the MFs derived when combining observational gas-to-stellar mass correlations from small data sets with the GSMF, (see also, Lemonias et al. 2013; Butcher et al. 2018). While this is an indirect method to study the demographics of the galaxy distribution, it is a valid and valuable approach that gives results that are comparable to direct observations and generalise them into a full bivariate distribution.

Figure 10. Impact of random and systematics errors in the baryonic, stellar, cold gas, atomic, and molecular gas MFs for all and separately for early- and late-type galaxies. The dashed lines show the ‘observational’ MF from Section 4.4, while the solid lines show the MF after deconvolving from random errors, i.e., the intrinsic MFs. Systematic errors are shown with the shaded areas. While the impact of random errors affects notably the total cold, atomic, and molecular gas MFs, the impact of systematic uncertainty on $M_{*}$ is apparently marginal on them. However, the systematic uncertainties on $M_{*}$ are noticeable in the stellar and baryonic MFs.

5.1. The impact of random and systematics errors

When deriving stellar, $\rm HI$ and $\rm H_{2}$ masses, there are two sources of errors that will inevitably propagate over the MFs: the random and systematic errors. In this Section, we discuss the impact of both sources of errors on our results.

5.1.1. Random errors

The estimation of masses from both photometric and radio observations is subject to random errors. Here, we determine their impact on our resulting MFs. For simplicity, we assume that random errors follow lognormal distributions with a constant dispersion and independent of galaxy morphology. For the stellar masses, we assume a dispersion of $\sigma = 0.1$ dex following Behroozi, Conroy, & Wechsler (2010), Mendel et al. (2014), and Rodríguez-Puebla et al. (2017). For HI masses, $\sigma = 0.14$ dex, and for $\rm H_{2}$ masses, $\sigma = 0.22$ dex, following Calette et al. (2018, and more references therein). As for the gas and baryonic masses, we assume errors of respectively $\sigma = 0.14$ dex and $\sigma = 0.1$ dex as they are dominated by HI and $M_{*}$ components, respectively, especially at high masses, where random errors have a larger impact. Thus, our ‘observational’ MFs^s are the result of the convolution of the distribution of random errors and the respective intrinsic MFs. That is, our ‘observational’ MFs are given by $\phi_{\rm obs} = \mathcal{G}\ast \phi_{\rm int}$ , where the symbol $\ast$ denotes the convolution operation, $\mathcal{G}$ is the distribution of random errors, and $\phi_{\rm int}$ is the intrinsic MF. For more details, the reader is referred to Appendix D. There we describe our numerical algorithm for deconvolving the intrinsic MF, $\phi_{\rm int}$ .

In Figure 10, we reproduce with dashed lines the ‘observational’ MFs derived in Section 4.4. Their corresponding intrinsic MFs are shown with solid lines. In the same figure, we present the ratios $\phi_{\rm obs}/\phi_{\rm int}$ to show the effect of the deconvolution. The effect of deconvolving from random errors is small at low-intermediate masses but it increases at the massive end since the MFs are steeper (Eddington 1940). This is simply because the convolution depends on the logarithmic slope of the intrinsic MF (e.g., Cattaneo et al. 2008); the steeper the slope, the larger the effect on the MFs. This is also the reason why we observe a lower impact in the baryonic and stellar MFs compared to the $\rm HI$ , $\rm H_{2}$ , and cold gas MFs; the latter fall steeply at the high-mass end. For example, the intrinsic $\rm HI$ MF is a factor of $\sim\!4$ lower than the ‘observational’ one at $M_{\rm HI}\,{\sim}\,6\times10^{10}\ \text{M}_{\odot}$ , while the intrinsic $\rm H_{2}$ MF is an order of magnitude lower than the ‘observational’ MF at $M_{\rm H2}\sim\!2\times10^{10}\ \text{M}_{\odot}$ . The intrinsic gas MF is an order of magnitude lower than the ‘observational’ MF at $M_{\rm gas}\sim\!10^{11}\ \text{M}_{\odot}$ . Note that for the $\rm HI$ , $\rm H_{2}$ , and cold gas MFs the impact of random errors is more noticeable in late-type galaxies than in early-type ones.

5.1.2. Systematic errors

In addition to random errors, systematic errors have an impact when determining the MFs. The IMF is one of the most important sources of systematic errors for the GSMF. In this paper, we assumed an universal IMF given by the Chabrier (2003) function. While there is much debate on the IMF (see e.g., Bastian, Covey, & Meyer 2010; Conroy et al. 2013; Bernardi et al. 2018), exploring the different alternatives is beyond the scope of this paper.

The stellar masses are calculated typically using colour-dependent mass-to-light ratios based on results from stellar population synthesis (SPS) models (for a recent review, see Conroy 2013). Thus, the calculated stellar masses depend on the used SPS model. This introduces a systematic uncertainty in $M_{*}$ . Indeed, systematics in $M_{*}$ from SPS can be as large as $\sim\!0.25$ dex, see for example, Pérez-González et al. (2008); Muzzin et al. (2009); Moustakas et al. (2013); Rodríguez-Puebla et al. (2017) and references therein. Recently, Bernardi et al. (2017) showed that systematics from SPS introduces errors that are as large as $\sim\!0.5$ dex in the normalisation of the GSMF at the high-mass end. In Appendix A, we found similar differences using various recipes of colour-dependent mass-to-light ratios. While in this paper, we calculate five different stellar masses for every galaxy and decided to use the geometric mean of the five as our fiducial definition of $M_{*}$ , the above inevitable introduces the question of which stellar mass definition shall we use when deriving our MFs. Additionally, Bernardi et al. (2017) determined that systematics in photometry are of $\sim\!0.1$ dex. In order to quantify the impact of stellar populations in our MFs, in Appendix A, we noted that a constant shift of $\sim\!\pm0.15$ dex in the stellar mass axis reproduces systematic errors in the GSMF. In addition, Figure 1 shows that the same shift in the stellar mass axis could also explain differences from photometry. Thus, hereafter, we will use a shift of $\pm0.15$ dex in the stellar mass axis as our fiducial model for systematic errors in the GSMF. Note that we are assuming that this shift will be independent of morphology and we are ignoring systematic errors in the atomic and molecular gas components.

Figure 10 shows the impact of systematic errors from SPS models and photometry as the shaded areas around their corresponding $\phi_{\rm int}$ (solid lines). The effects of systematics is non-negligible at the massive end of the stellar and baryonic MFs; we observe differences up to $\sim\! 0.6 $ dex in their normalisations. This is approximately the same both for early- and late-type galaxies. The impact of systematic errors in the gas, $\rm HI$ , and $\rm H_{2}\,{\rm{MFs}}$ is marginal; we notice a shift in their normalisations of $\sim\!0.07$ at their low-mass ends but they increases respectively to $\sim\!0.4$ , $\sim\!0.4$ , and $\sim\!0.3$ dex at their massive ends. The above is due to the steeper slopes observed at the high-mass end from these MFs. In conclusion, the impact of systematic uncertainties in $M_{*}$ is only marginal for the derived $\rm HI$ , $\rm H_{2}$ , and cold gas MFs, making our results robust against this source of uncertainties.

As mentioned in the Introduction, the value of deriving robust MFs is that they can be used as key tools for constraining the processes that govern the evolution of the galaxies. However, using direct measurements from observations to constrain galaxy formation models is not trivial due to random and systematic errors, as discussed here. We end this subsection by emphasising the importance of deconvolving from random errors and understanding the impact of systematic errors when reporting results on galaxy demographics.

6. Summary and conclusions

We present a self-consistent empirical approach that unifies local galaxy gas-to-stellar mass correlations and the MFs of galaxies traced by their different baryonic components. We make available a Python code that displays tables and figures with all the relevant statistical distributions and correlations discussed in this paper.^t Next, we summarise our main results which can be used for comparing with theoretical predictions or as input for modelling galaxy mock catalogues:

• Conditional probability distribution functions (CPDFs): Section 2.2 presents the functional forms for the $\rm HI$ and $\rm H_{2}$ mass conditional distributions given $M_{*}$ (the CPDFs), which are described by Equations (7)–(12). Our best-fit parameters to the empirical information presented in Paper I are listed in Table 3, while Figures 3 and 4 show the data with their corresponding best fits in various stellar mass bins. Theoretical predictions for the $\rm HI$ , $\rm H_{2}$ , and cold gas CPDFs can be confronted with our empirically constrained distributions, for all galaxies as well as for early and late types in case the morphological classifications are available. If these predictions are limited in stellar and/or gas masses, then our (analytical) $\rm HI$ -CPDFs and $\rm H_{2}$ -CPDFs and their moments can be easily calculated over the same mass ranges as the theoretical predictions for a comparison. The $\rm HI$ - and $\rm H_{2}$ -CPDFs combined with the GSMF allowed us to calculate the respective bivariate mass distributions for all the galaxy population as plotted in Figures 6 and 7.

• Moments of the CPDFs: The (analytical) $\rm HI$ - and $\rm H_{2}$ -CPDFs contain the information about any moment of the distributions. Figure 4 (see also Figures 6 and 7) shows the $\langle \log M_j \rangle+-\log M_{*}$ relationships, with $j={\rm HI}$ , $\rm H_{2}$ , for early- and late-type galaxies as well as for all galaxies. In addition, we present these relationships using the arithmetic mean, $\langle M_j\rangle$ . As expected, these relationships lie above from those calculated with the logarithmic mean, $\langle\log M_j\rangle$ . Moreover standard deviations can vary significantly if they are computed with respect to the arithmetic or logarithmic mean, which also depends on the shape of the distributions. Other statistical measures that can be used to characterise the population distributions are medians and percentiles, for example. As mentioned above, any statistical quantity can be computed with our CPDFs and confronted with both theoretical and/or observational results.

• The Galaxy Stellar Mass Functions: In Section 3.1, we determined the GSMF from the SDSS DR7 based on the photometric catalogue from Meert et al. (2015) and (2016) for masses above $M_{*}=10^9\ \text{M}_{\odot}$ . For masses down to $\sim\!3\times10^{7}\text{M}_{\odot}$ , we used the low-z SDSS DR4 (Blanton et al. 2005a; b), and corrected it from surface brightness incompleteness and fluctuations due to large-scale structures. We determined also the fractions of early- and late-type galaxies using the SDSS DR7 morphological classification of Huertas-Company et al. (2011). Stellar masses were derived from five colour-dependent mass-to-light ratios. We used as our fiducial definition the geometric mean of these five stellar masses derived for each galaxy. We also determined the impact of systematic errors in $M_{*}$ due to mass-to-light ratio uncertainties in our MFs.

• Calculated MFs: Section 4.4 presents the results of calculating with our approach the MFs for atomic, molecular, cold gas, and baryons for early- and late-type galaxies, as well as for all galaxies. As discussed in Section 5.1, random errors in mass determinations artificially decrease the slope of the ‘observational’ MFs, an effect that affects especially the high-mass end, and that would lead to incorrect conclusions when comparing to theoretical predictions. Figure 10 presents our MFs deconvolved from random errors, that is, the intrinsic MFs, for different baryon matter components, and separately for early- and late-type galaxies. In the same section, we studied the effects on the MFs from systematic errors in $M_{*}$ , also shown in Figure 10. In Appendix A and Figure 1, we showed explicitly that systematic errors in the GSMF due to mass-to-luminosity ratios and photometric uncertainties are well represented by a shift in the $M_{*}$ -axis of $\pm 0.15$ dex. The effect of random errors in the baryonic MF is of the same order, while for the gas MFs the propagated systematic errors in $M_{*}$ have a negligible effect. Note that our MFs are complete only above a given mass limit, $\sim\! 3\times 10^7\ \text{M}_{\odot}$ for the GSMF, $\sim\! 10^7\ \text{M}_{\odot}$ for the $\rm H_{2}$ MF, and $\sim\! 10^8\ \text{M}_{\odot}$ for the $\rm HI$ , cold gas, and baryonic MFs.

From the results summarised above, we highlight the following conclusions:

• The low-mass slope of our GSMF, corrected for surface brightness incompleteness, is $\alpha \approx -1.4$ , consistent with recent determinations based on the deeper surveys such as GAMA (Wright et al. 2017), and estimations based on the search of low surface brightness galaxies from core-collapse supernovae (Sedgwick et al. 2019). The slope for the high-mass end is shallower than previous determination most likely as the result of the new photometric catalogue employed in this paper (Meert et al. 2015). Similar results have been reported in Bernardi et al. (2017).

• The total GSMF is well fitted by a function composed of a sub-exponential Schechter function and a double power law function. This fitting model has an error of less than $\sim\! 2\%$ in the mass range $2\times 10^9-5\times 10^{11}\ \text{M}_{\odot}$ . At the smallest and largest masses, the deviations increase to values above $\sim\!20\%$ . In contrast, the commonly employed double Schechter function model performs considerably worse.

• Systematic errors due to SPS models, that affect results on mass-to-light ratios, introduce a systematic effect on the normalisation of the GSMF, especially at the massive end. We find differences between $\sim\! 0.5- 1$ dex, consistent with the result discussed in Bernardi et al. (2017).

• The $\rm HI$ , $\rm H_{2}$ , and cold gas MFs are mostly dominated by late-type galaxies. In general, we notice that our $\rm HI$ MF is in good agreement with previous determinations from blind surveys. Similarly, the $\rm H_{2}$ MF is consistent with previous determinations based on CO follow-up optically selected samples. When we compare to the HIPASS and ATLAS 3D surveys for early-type galaxies, our $\rm HI$ MF is consistent with those observations. However, our $\rm H_{2}$ MF for early-type galaxies is in tension at the low-mass side with the MF derived from the ATLAS 3D survey.

• Our ‘observational’ MFs were deconvolved from random errors to obtain the intrinsic MFs. The effect of random errors is small at the low-mass end but larger at the high-mass end of our MFs. This is because the convolution depends on the logarithmic slope of the intrinsic MFs. Because the baryonic and stellar MFs are shallower at the massive end, the effects are relatively small, but the atomic, molecular, and cold gas MFs have steeper slopes resulting in a larger effect.

• While for the stellar (and hence baryonic) MF systematic errors due to mass-to-light ratio uncertainties introduce a non-negligible effect, especially at the high-mass end, for the atomic, molecular, and gas MFs, the effects of systematics are small. We thus conclude that our determinations for the gas MFs are robust against systematic errors in the the $M_{*}$ determination.

• We determined the $z\,{\sim}\,0$ cosmic densities of $\rm HI$ , $\rm H_{2}$ , cold gas, stars, and baryons locked in galaxies calculated from the respective MFs. Our results are in good agreement with previous determinations from different local censuses. Most of the atomic and molecular H gas is in late-type galaxies, $\sim\!96\%$ of the mass density, while this fraction decreases to $\sim\!70\%$ and $\sim\! 65\%$ for baryons and stars. We find that the fraction of $\rm HI$ and $\rm H_{2}$ in galaxies with respect to the universal baryon fraction is respectively $\sim\! 1 \%$ and $\sim\! 0.2\%$ , while the respective fractions for mass in stars is $\sim\! 4\%$ . Baryons in galaxies (the ionised and hot gas were not included) are $\sim\! 5.4\%$ of the universal baryon fraction.

• Based on the values reported in the literature for the local CSFR of star-forming (late-type) galaxies, we estimated the cosmic $\rm H_{2}$ and total gas depletion times of late-type galaxies. These timescales, $t_{\rm dep}\text{(}{\rm H_{2}}\text{)}\approx 1.3\ \text{Gyr}$ and $\bar{t}_{\rm dep,L}({\rm gas})\approx 10.14\ \text{Gyr}$ , respectively, imply that galaxies, on average, are inefficient to convert their molecular gas into stars, and are inefficient to transform their atomic gas into molecular gas. The depletion time for the total neutral hydrogen is $\bar{t}_{\rm dep,L}({\rm H}) = \bar{t}_{\rm dep,L}({\rm gas}) / 1.4\sim\! 7.25$ Gyrs. On the other hand, the average cosmic SF timescale (the inverse of the cosmic sSFR) is $\bar{t}_{\rm SF,L}\approx 20.3$ Gyrs, which implies that the ratio $\bar{t}_{\rm dep,L}({\rm H})/\bar{t}_{\rm SF,L} = 0.38$ . This shows that the gas reservoir of late-type galaxies has not yet been dramatically consumed by star formation.

Here, we provided a statistical description for calculating any moment to characterise the gas-to-stellar mass correlations, the $\rm HI$ - and $\rm H_{2}$ -stellar mass bivariate distributions as well as all the respective MFs. One of our motivations for this paper is to provide the community with a full self-consistent phenomenological description of the local galaxy population for various properties and divided into the two main morphological types in order to be confronted with theoretical results, both from semi-analytical models and cosmological hydrodynamical simulations. The next generation of sensitive radio telescopes will be able to survey large samples of extragalactic sources in $\rm HI$ and $\rm H_{2}$ gas, something that is a common practice with current optical surveys. Thus, robust and unbiased bivariate distributions and MFs of $\rm HI$ and $\rm H_{2}$ gas over large mass ranges will be routinely derived in the future along with the relationships of the gas contents with their optical/IR properties. Preparatory to that, and to pave the road to these surveys, studies based on radio follow-up observations of (relative small) optically selected galaxy samples provide valuable information that can be used for the gas demographics of galaxies. In this work, we have exploited the results from many of these studies, and by means of the conditional (or bivariate) approach we were able to derive the abundances of local galaxies as traced by different baryonic components and separated into the two main groups of galaxies, early, and late types.

The present work is the second paper of a series. In Paper I, we derived the CPDFs of $\rm HI$ and $\rm H_{2}$ as a function of $M_{*}$ , separately for early- and late-type galaxies, for an extensive compilation and homogenisation of radio data from the literature. In the present work, we made extensive use of these data. In the future, we will use the MFs derived here to extend the galaxy–halo connection for different baryonic components, and we will show that not only the $\rm HI$ and $\rm H_{2}\,$ MFs derived here are in good agreement with radio blind or optically selected surveys but also with the observed galaxy spatial clustering as a function of $\rm HI$ gas mass.