© ESO 2018

1. Introduction

The optical luminosity function (LF) and the HI mass function (HIMF) are two of the most important and fundamental population tracers of the volume density distribution of galaxies in the universe. They yield clues to both the baryonic and dark matter content as well as their evolutionary histories. Consequently, both the LF and HIMF are frequently used to constrain galaxy formation models (see, e.g., Benson et al. 2003; Lu et al. 2014).

The LF was the first tracer of volume density to be studied in detail as it was historically easier to detect galaxies at optical wavelengths than in the radio spectrum. The first attempts to fit an analytic form to the LF were carried out by, e.g., Zwicky (1957), Kiang (1961), and Abell (1965), but the most successful form was developed by Schechter (1976) in which the LF can be characterized by:

$$\begin{array}{c}\hfill \mathrm{\Phi }(L)\mathrm{d}L={\varphi }^{\ast }(\frac{L}{L^{\ast }}{)}^{\alpha }\text{exp}(-\frac{L}{L^{\ast }})d(\frac{L}{L^{\ast }}),\end{array}$$(1)

where ϕ ^*, L ^* and α represent the normalization constant, characteristic luminosity of the “knee” of the function and the faint end slope, respectively. This analytic form has been shown to be a good estimator of the LF across differing environments and many subsequent studies have attempted to analyze how the LF’s parameters change as a function of environment and redshift (see, e.g., Felten 1985; Efstathiou et al. 1988; Loveday et al. 1992, Blanton et al. 2001; 2003; Montero-Dorta & Prada 2009; McNaught-Roberts et al. 2014; Loveday et al. 2015).

The HIMF can be modeled successfully using the same functional form as the LF. It has been analyzed to a somewhat lesser extent due mainly to completeness limitations of past HI surveys. The HIMF has primarily been studied using several approaches: blind HI surveys (e.g., Zwaan et al. 1997; 2003; Kilborn et al. 1999; Kovac et al. 2005; Martin et al. 2010; Hoppmann et al. 2015) and targeted HI surveys (e.g., Springob et al. 2005, S05 hereafter), and extrapolations from optical surveys using proxies to determine HI content (e.g., Rao & Briggs 1993; Solanes et al. 1996).

Regardless of the methods used for determining the LF and HIMF, the end goal has always been to determine or analyze the functional form of these two volume tracers independently. One of the main goals of NIBLES, the Nançay Interstellar Baryons Legacy Extragalactic Survey (see van Driel et al. 2016, Paper I; Butcher et al. 2016, Paper II), is to determine the correlation between the optical LF and the HIMF – more specifically, we want to analyze the HIMF and other galaxy properties as a function of optical luminosity. NIBLES is a 21 cm HI line survey at the 100m class Nançay Radio Telescope (NRT) of 2610 galaxies selected from the Sloan Digital Sky Survey (SDSS; see, e.g., York et al. 2000). The sample was selected to have radial velocities 900 <  c z <  12 000 km s⁻¹ and to be distributed evenly over the absolute z-band magnitude range of galaxies in the local volume (∼ − 13.5 to −24), which was used as a proxy for total stellar mass – see Papers I and II for further details.

Here we present the derivation and analysis of the optical luminosity-HI mass bivariate luminosity function (BLF) based on the NIBLES data. This two dimensional distribution function will provide additional constraints to galaxy formation models and allow more detailed analyses of galaxy populations and evolutionary histories.

In Sect. 2 we summarize the selection criteria of the NIBLES galaxies, in Sect. 3 we explain our selection of galaxies from the NIBLES sample for the bivariate luminosity and HIMF analysis. Section 4 discusses the method we use to derive our bivariate function and results are discussed in Sect. 5. In Sect. 6 we compare our resulting HIMFs and Ω_HI to those of other surveys, both optical and blind HI, and compare our gas-to-light distributions to those of ALFALFA. In Sect. 7 we discuss future work and possible ways to improve the accuracy of our BLF and we present our conclusions in Sect. 8.

2. The NIBLES sample

The NIBLES galaxy selection criteria are:

• Must have both SDSS magnitudes and optical spectrum;

• Must lie within the local volume (900 <  c z <  12 000 km s⁻¹);

• Uniform sampling of each 0.5 magnitude wide bin in absolute z-band magnitude, M _z ;

• Preferentially observe the most nearby objects;

• No a priori selection on color.

Throughout this paper distances are based on a Hubble constant of H ₀ = 70 km s⁻¹ Mpc⁻¹. Our lower velocity limit of c z = 900 km s⁻¹ was chosen to avoid large distance uncertainties due to peculiar velocities as well as the problems with the SDSS deblender encountered when trying to fit objects with large angular diameters (see, e.g., West 2005).

We note that the initial targets for NIBLES were selected from SDSS DR5, but the optical data was subsequently updated to DR9 after our observations were concluded to take advantage of the improved SDSS photometric processing pipeline. NIBLES galaxies were initially selected at random from DR5 until the most densely populated 0.5 mag wide M _z bins contained at least 150 galaxies. The exceptions were the more sparsely populated highest and lowest magnitude bins, which contained all the SDSS galaxies available at the time.

In DR5, association between photometric and spectroscopic targets was made using a position-based match only. However, this led to (sometimes severe) underestimates of total luminosities for many galaxies due to improper association between the spectroscopic target and a photometric object that does not contain the most flux from the target galaxy. This problem has been remedied since DR8, through an additional flux-based association between spectroscopic targets and photometric objects. (This affected 465 cases marked with flag F in Paper I.) The effect of updating the photometry of our sample from DR5 to DR9 is that galaxies were systematically shifted to brighter magnitudes. After this reassessment of total galaxy magnitudes, we still have adequate sampling of the luminosity function over the range (−23.5 < M _z ≲ − 13.5) with more than 50 galaxies per 0.5 magnitude wide bin. In addition to the new flux-based association between photometry and spectroscopy, the post DR8 releases contain re-processed photometry, which affect how we normalize our bivariate luminosity function (BLF). This will be discussed further in Sect. 4.2.

NIBLES is aimed to complement other recent and/or ongoing large HI surveys in the local volume, in particular, blind surveys such as ALFALFA (e.g., Haynes et al. 2011, H11 hereafter, and references therein). One advantage of NIBLES, which targets selected objects only, over blind HI surveys is our increased on-source integration time allowing us to reduce the rms noise of our observations. Each NIBLES source was initially observed for about 40 min of telescope time, resulting in a mean rms noise of ∼3 mJy. Clear detections were not observed again, but for weak or non-detections, observations were repeated (as time allowed) resulting in a bimodal rms noise distribution (see Fig. 6 in Paper I) with a mean rms of 2.3 mJy for the latter sources.

Furthermore, we have shown (Paper I) that the NRT fluxes are on average about 15% lower than those of the HI flux calibrator galaxies of O’Neil (2004). For the purposes of this analysis, we therefore choose to increase our HI line flux values by 15% and re-calculate the HI masses accordingly. We also adjust other HI surveys examined here to the same flux scale (see Paper I for details on flux comparisons).

2.1. Effects of H I self-absorption

Our total HI masses of individual galaxies are calculated using M_HI = 2.36 × 10⁵ · D² · F_HI where F_HI is the integrated measured line flux (in Jy km s⁻¹) the galaxy distance in Mpc. We make no corrections for possible HI self-absorption (HISA), for two reasons – lack of consensus on its magnitude, and consistency with other HI surveys.

In principle, HISA will occur in dense interstellar clouds, and it can decrease the measured integrated HI line flux of an entire galaxy and thus lead to an underestimation of its total HI mass, provided that such clouds are sufficiently abundant.

Using single-dish radio telescope 21 cm line profiles, e.g., for the HIPASS blind HI survey an HISA correction factor of 1.1 and 1.3 was estimated for the F_HI of an average galaxy by Zwaan et al. (1997) and Lang et al. (2003), respectively, whereas Zwaan et al. (2003) concluded that no measurement of the possible effects of HISA could be made as they could not be disentangled from optical dust extinction.

Using high-resolution interferometric HI images of three Local Group galaxies (M31, M33 and the LMC), Braun (2012) derived a HISA correction factor of on average 1.34 to be applied to their measured integrated F_HI values.

The HISA correction factor is expected to be the largest for edge-on, gas-rich galaxies, which are commonly found in blind HI surveys. However, the NIBLES sample was selected without a galaxy color bias and hence contains a larger fraction of relatively gas-poor early-type galaxies than a typical blind HI survey.

Furthermore, we want to maintain consistency with recent blind HI surveys such as HIPASS (e.g., Koribalski et al. 2004; Meyer et al. 2004) and ALFALFA (e.g., Martin et al. 2010; H11) which did not apply corrections for possible HISA effects.

3. Bivariate selection criteria

As part of our classification process of the NIBLES sample, galaxies were flagged according to various properties (see Paper I for details). For our bivariate analysis we exclude all galaxies with C1, C2, U and R flags. The HI confusion flags C1 and C2 indicate, respectively, that the HI line profile of a galaxy is either definitely or probably contaminated by another detection within the telescope beam. The uncertain SDSS flux flag U indicates galaxies we deemed to either be missing a large portion of the photometric flux due to improper Petrosian radii fitting in the SDSS or due to foreground star contamination. The resolved HI source flag R indicates sources which may be missing a significant amount of HI flux.

In Papers I and II we defined a marginal detection category which is also excluded here, and those galaxies are treated as non-detections in this paper.

After applying the above cuts, our final sample consists of 2256 galaxies.

For this analysis, we use the SDSS r-band Petrosian magnitudes rather than the z-band magnitudes used to select sources for the NIBLES sample. We use the r-band because the SDSS processing pipeline sets the photometric aperture for each galaxy in this band, which therefore contains the most accurate Petrosian magnitude for a given galaxy.

Additionally, we do not make corrections for Galactic foreground extinction in order to maintain consistency with the luminosity functions from Montero-Dorta & Prada (2009) and Blanton et al. (2003) which we use in our analysis. However, since the SDSS galaxies are at high Galactic latitudes, applying an extinction correction has a negligible impact on our results: the luminosities of the NIBLES galaxies are diminished by only 0.04 dex on average when we apply the extinction corrections from the SDSS catalog data (based on Schlegel et al. 1998).

Note that while our analysis here utilizes non-extinction-corrected Petrosian magnitudes, the magnitudes listed in Papers I and II are extinction-corrected model magnitudes.

4. Method

The basic method we use is to count the number of galaxies in each logarithmic luminosity and HI mass bin to determine the HI mass distribution of the NIBLES sample. Then we scale the resulting bivariate distribution to a luminosity function to derive the two-dimensional volume density. For clarity, we outline the process using the two dimensional stepwise maximum likelihood method (2DSWML) commonly used to determine mass and luminosity functions (see, e.g., Efstathiou et al. 1988; Zwaan et al. 2003).

4.1. NIBLES M_HI-L_r bivariate distribution

To determine the M_HI distribution of the NIBLES sample, we need to determine the probability that a given galaxy will have a particular HI mass. To do this, we divide the NIBLES sample into N_M logarithmic HI mass bins and N_L logarithmic luminosity bins of width ΔM and ΔL, respectively. The probability that galaxy i will have an HI mass M_HIi given its luminosity L_k − ΔL/2 ≤ L_i ≤ L_k + ΔL/2 is:

$$\begin{array}{c}\hfill p(M_{\mathrm{H}{i}_i}|L_k)=\frac{\theta (M_{\mathrm{H}{i}_i}|L_k)}{{\int }_0^{\infty }\theta (M_{\mathrm{H}i}|L_k)\mathrm{d}M},\end{array}$$(2)

where θ represents an unknown HI mass – luminosity distribution.

The likelihood of a given HI mass for a particular luminosity is simply the product of the probabilities:

$$\begin{array}{c}\hfill \mathcal{L}=\prod _{i=1}^{N_{\mathrm{g}}}{p}_i,\end{array}$$(3)

where N_g is the total number of galaxies in the sample. To maximize the likelihood, we take the natural logarithm:

$$\begin{array}{cc}\hfill ln\mathcal{L}=\sum _i^{N_{\mathrm{g}}}ln{p}_i=& \sum _i^{N_{\mathrm{g}}}\sum _j^{N_{\mathrm{M}}}{B}_{\mathit{ij}}(M_i-M_j|L_k)ln{\theta }_{\mathit{jk}}\hfill \\ \hfill & -\sum _i^{N_{\mathrm{g}}}ln(\sum _j^{N_{\mathrm{M}}}{H}_{\mathit{ik}}{\theta }_{\mathit{jk}}\mathrm{\Delta }M)\hfill \end{array}$$(4)

set the derivative equal to zero, and solve for θ_jk, i.e., the probability that a galaxy will have an HI mass in bin j and a luminosity in bin k. The solution to this equation is straight-forward and need not be discussed here. The important point to note is that it reduces to:

$$\begin{array}{c}\hfill {\theta }_{\mathit{jk}}\mathrm{\Delta }M\mathrm{\Delta }L=n_{\mathit{jk}}/N_k,\end{array}$$(5)

where n_jk is the total number of galaxies in the jk HI mass and luminosity bin, and N_K is the total number of galaxies in luminosity bin k.

Equation 4 is a summation over all galaxies in the sample for each M_j and L_k bin, counting only the galaxies that fall into a particular M_j bin for a given L_k. The B_ij denotes whether or not a galaxy is in the j HI mass and k luminosity bin and H_jk denotes whether a galaxy is in the k luminosity bin.

When assigning a particular galaxy to a mass bin, we take into account the uncertainty of its HI mass, in the sense that we assign an occupation number that is proportional to the fraction of the galaxy’s possible mass range that falls within the bin. This has the effect of smoothing out the bivariate distribution for the low sample size bins. Therefore, B_ij can take on the following values:

$$\begin{array}{c}\hfill B_{\mathit{ij}}=\{\begin{array}{cc}1,\hfill & \text{if}|{\chi }_{-}-M_j|\le \mathrm{\Delta }M/2\text{and}\hfill \\ \hfill & |{\chi }_{+}-M_j|\le \mathrm{\Delta }M/2\text{and}|L_i-L_k|\le \mathrm{\Delta }L/2\hfill \\ \gamma ,\hfill & \text{if}{\chi }_{-}+\mathrm{\Delta }M/2<M_j\text{and}\hfill \\ \hfill & {\chi }_{+}-\mathrm{\Delta }M/2>M_j\text{and}|L_i-L_k|\le \mathrm{\Delta }L/2\hfill \\ (\frac{1}{2}+{ϵ}_{-})·\gamma ,\hfill & \text{if}|{\chi }_{-}-M_j|\le \mathrm{\Delta }M/2\text{and}\hfill \\ \hfill & |{\chi }_{+}-M_j|\ge \mathrm{\Delta }M/2\text{and}|L_i-L_k|\le \mathrm{\Delta }L/2\hfill \\ (\frac{1}{2}-{ϵ}_{+})·\gamma ,\hfill & \text{if}|{\chi }_{+}-M_j|\le \mathrm{\Delta }M/2\text{and}\hfill \\ \hfill & |{\chi }_{-}-M_j|\ge \mathrm{\Delta }M/2\text{and}|L_i-L_k|\le \mathrm{\Delta }L/2\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}\end{array}$$(6)

where χ₋ and χ₊ are the minimum and maximum HI mass values, γ = ΔM/2δM_i and ϵ₊ and ϵ₋ are $(M_j-{\chi }_{+})/\mathrm{\Delta }M$ and $(M_j-{\chi }_{-})/\mathrm{\Delta }M$, respectively.

H_ik takes the following values:

$$\begin{array}{c}\hfill H_{\mathit{ik}}=\{\begin{array}{cc}1,\hfill & \text{if}|L_i-L_k|\le \mathrm{\Delta }L/2\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$(7)

By design, this density distribution has

$$\begin{array}{c}\hfill \sum _{j=1}^{N_{\mathrm{M}}}{\theta }_{\mathit{jk}}\mathrm{\Delta }M\mathrm{\Delta }L=\varphi (L_k)\mathrm{\Delta }L,\end{array}$$(8)

where ϕ represents the normalized density distribution of the NIBLES sample. To scale our distribution to a volume density we use the r-band luminosity function from (Montero-Dorta & Prada 2009, M09 hereafter) which was also derived from a pre-DR8 SDSS data release, like NIBLES. We also scaled the M09 luminosity function to match our Hubble constant of 70 km s⁻¹ Mpc⁻¹.

We can now calculate the luminosity – HI mass volume density distribution as follows:

$$\begin{array}{c}\hfill \mathrm{\Phi }(M_j,L_k)\mathrm{\Delta }M\mathrm{\Delta }L={\theta }_{\mathit{jk}}\mathrm{\Delta }M·{\int }_{L_k-\mathrm{\Delta }L/2}^{L_k+\mathrm{\Delta }L/2}\mathrm{\Phi }(L)\mathrm{d}L.\end{array}$$(9)

From this distribution, we can calculate the HIMF and LF by:

$$\begin{array}{c}\hfill {\displaystyle \sum _{k=1}^{N_L}}\mathrm{\Phi }(M_j,L_k)\mathrm{\Delta }M\mathrm{\Delta }L=\mathrm{\Phi }(M)\mathrm{\Delta }M\end{array}$$(10)

and

$$\begin{array}{c}\hfill {\displaystyle \sum _{j=1}^{N_{\mathrm{M}}}}\mathrm{\Phi }(M_j,L_k)\mathrm{\Delta }M\mathrm{\Delta }L=\mathrm{\Phi }(L)\mathrm{\Delta }L.\end{array}$$(11)

Since our HI masses are distributed according to luminosity (due to our selection criteria), we must treat the luminosity uncertainty differently than the HI mass uncertainty. Due to our normalization to a luminosity function, any luminosity uncertainties applied in the method described above would be normalized out in Eq. (9). We therefore apply the luminosity uncertainties via a bootstrapping method described in Sect. 4.2 and assign a final value to each θ_j,k bin from the mean values of the bins in that bootstrapped distribution.

4.2. Differences between SDSS DR5 and DR9 luminosities

We calculate uncertainties in the BLF using two separate methods: the standard deviation of the binomial distribution, to address sampling errors, and a bootstrapping method to address galaxy luminosity uncertainties. The standard deviation of the binomial distribution for each bin in our BLF is given by:

$$\begin{array}{c}\hfill \sigma =\sqrt{\frac{{\theta }_{\mathit{jk}}(1-{\theta }_{\mathit{jk}})}{n}}\end{array}$$(12)

where n is the total number of galaxies in each luminosity bin.

The bootstrap method addresses errors from the uncertainties in galaxy luminosities in the SDSS database. The listed Petrosian magnitude uncertainties in the SDSS data releases are on average a factor of four lower than the uncertainties in the HI masses. However, these values will be underestimated as the listed uncertainty is only a function of the flux uncertainty within the aperture fitted to each galaxy by the processing pipeline (see York et al. 2000), but does not account for uncertainties introduced by fitting errors. As the SDSS does not provide estimates for the accuracy of this fit, we estimate it ourselves by comparing the scatter of the difference between DR5 and DR9 Petrosian magnitudes. For each DR9 photometric source from the NIBLES catalog, we pulled the nearest corresponding primary photometric source from DR5 for comparison. Every source in the NIBLES sample was visually checked to ensure that the SDSS DR9 photometric object was centered on the galaxy and contained an adequate Petrosian radius to account for all, or most, of the flux (see Paper I for details). The difference between DR5 and DR9 apparent r-band Petrosian magnitudes as a function of DR9 r-band luminosity is shown in Fig. 2.

There is an increasing discrepancy between the DR5 and DR9 magnitudes with increasing r-band luminosity, as indicated by the magenta dots which represent the mean value of DR5–DR9 magnitudes in each 0.5 dex wide bin in L_r. The standard deviation within each bin is fairly consistent at about 0.1 mag, but slightly lower in the two lowest luminosity bins. The scatter in this relationship is about three to four times larger than the uncertainties listed in the SDSS database. The uncertainties in luminosity as indicated by this scatter would place them in the same range as the HI mass fractional uncertainties.

The relationship shown in Fig. 2 provides us with both a measure of typical total photometric uncertainties and the offsets between the SDSS DR5 and DR9 data. The latest derivation of an SDSS luminosity function is from M09 which uses photometry from DR6 (which uses the same photometric processing pipeline as DR5). Since the more luminous galaxies have a greater difference between DR9 and DR5 magnitudes, we must take this offset into account when normalizing NIBLES DR9 data to the M09 luminosity function. To accomplish this, we create a bootstrapped sample as follows: we convert the y-axis from Fig. 2 to differences in luminosity and bin the data in 0.02 dex wide bins. Then, for each galaxy in our bivariate sample, we randomly assign a luminosity offset which is drawn from this luminosity difference distribution within its DR9 luminosity bin. This gives us a new set of luminosities from which we derive a new BLF. We repeat this process 100 times, generating a new BLF each time. We then calculate the standard deviation in each θ_j,k bin across our bootstrapped sample. The final, total uncertainty for each θ_j,k bin is the quadrature sum of the standard deviations of the bootstrapped sample and the binomial distribution in each bin. Additionally, each θ_j,k bin itself is assigned the median value of our bootstrapped sample, as mentioned in the previous section.

5. Results

Figure 3 shows the 2D HI mass – r-band luminosity distribution of the NIBLES sample (θ_jk from Eq. 5) binned in 0.5 dex intervals in both log(L_r) (in L_⊙) and log(M_HI) (in M_⊙). The distribution is shown such that the sum over each luminosity bin is equal to one. The plot shows (see also Fig. 6) that the galaxies in the highest-density bins have log(M_HI/L_r) ratios that lie within about 0.5 dex of zero, with this ratio becoming more negative with increasing luminosity. The highest and lowest HI mass bins, being on the edges of the distribution, are not fully sampled. In particular the log(L_r) = 6.75 and 10.75 bins only contains four and six galaxies, respectively, with their means in log(M_HI) being 6.8 and 10.6. The density in the log(M_HI) = 6.25 bin is much lower than in the next highest bin, which is due entirely to fractional occupation numbers being assigned in this bin as a result of our method described in Sect. 4.

Normalizing the bivariate distribution to the M09 luminosity function yields the two-dimensional luminosity-HI mass distribution function shown in Fig. 4. Panel a shows our resulting bivariate M_HI–L_r distribution. The input LF from M09 is shown in panel b and the resulting HIMF obtained by summing over luminosity is shown in panel c. Similarly to Fig. 3, the lowest HI mass bin (log(M_HI) = 6.25) is not completely sampled. We therefore do not plot this bin in the corresponding HIMF of panel c of Fig. 4.

Values for the mass distribution function from panel a in Fig. 4 are listed in Table 1 as log(Φ(M_HI,L_r)ΔM ΔL) in units of Mpc⁻³ dex⁻³, together with their fractional uncertainties.

In Fig. 5 we show the HIMF from panel c of Fig. 4 broken down into its corresponding fully-sampled luminosity bins. The low mass end shows clear flattening, in contrast to blind HI surveys such as (Zwaan et al. 2003, Z03 hereafter) and H11. This flattening is the result of two effects: (1) HI detectability limits for low-luminosity galaxies at V > 900 km s⁻¹, and (2) low-luminosity but high HI-mass galaxies that were not sampled due to being at or below the detection limit of the SDSS. We explore these two effects in the following two subsections.

5.1. HI sensitivity limits

The first cause of the HIMF faint-end flattening is that the NIBLES sample is optically selected with a recession velocity of V_hel > 900 km s⁻¹, in contrast to the abovementioned HI-selected surveys, Z03 and H11, which do not have a lower velocity limit – the former is about five times less sensitive and the latter has about the same sensitivity as NIBLES (see Paper I). Our sample is therefore probing galaxies near the limit of the NRT HI detectability for SDSS galaxies below log(L_r) = 8.0. These galaxies tend to be more distant than galaxies of similar luminosity in other surveys without the NIBLES minimum velocity restriction. Consequently, our detected galaxies in this low luminosity range tend to be biased toward higher HI mass-to-light ratios (referred to as “gas-to-light” ratios henceforth). In fact, Fig. 5 shows that for luminosities below log(L_r) = 8.0, the peak in volume density occurs at HI masses corresponding to a log(M_HI/L_r) ratio of greater than zero. For example, the log(L_r) = 7.75 bin has a peak HI mass volume density in the log(M_HI) = 8.25 bin.

This can be more clearly seen in Fig. 6 where we show the densities in Fig. 5 in the form of log(M_HI/L_r) curves for the same luminosity bins as shown in Fig. 5. We have again excluded the lowest luminosity bin due to poor sampling. The Luminosity bins below log(L_r) = 9.0 show a trend of increasing slope with decreasing luminosity for bins above log(M_HI/L_r) = 0. This increasing slope indicates that as luminosity decreases, the detections are on average probing more HI-rich subsamples within each luminosity bin (as to be expected). The more HI-poor subsamples, on the other hand, remain undetected, which lowers the measured HI density of the corresponding HI mass bins. For the low-luminosity bins log(L_r) = 7.25 − 9.25, this manifests itself as an increased slope for densities below log(M_HI/L_r) = 0. The relatively low gas-to-light ratios of these galaxies imply that the densities may be suppressed due to the galaxies within these bins being undetected simply due to distance effects. Accounting for these non-detections would tend to increase the volume densities for the gas-to-light ratios below log(M_HI/L_r) = 0, thus increasing the observed slope of the low-mass end of the HIMF. We explore this further in Sect. 5.3.

An interesting aspect of these trends in gas-to-light ratio is that for luminosity bins between log(L_r) = 7.25 and 9.75, the volume density of galaxies increases as a function of gas-to-light ratio at the same rate regardless of luminosity. In other words, the slopes of the densities as a function of gas-to-light ratio curves are all very similar up to their respective maximum ratios, which depend on luminosity. Further, the slopes for ratios above the point of maximum density also appear to be relatively similar.

The increasing gas-to-light ratio with decreasing luminosity trend also manifests itself as a difference in detection fraction (∼10%) between the log(L_r) = 9.25 and the lowest luminosity bins as shown in Fig. 7. This is mainly due to the combination of distance effects and sensitivity combined with the NIBLES selection criteria discussed previously. These low-luminosity galaxies are near the detection limit of both the SDSS and NRT, therefore limiting the both number of sources in the NIBLES sample and the relative number of HI detections. This decrease in detection fraction is to be expected and is explored further in Appendix A. The steeply decreasing detection fraction with increasing luminosity above log(L_r) = 9.25 (see Fig. 7) is mainly due to the increasing HI gas deficiency of larger galaxies.

5.2. Missing galaxies

The second cause of the HIMF faint-end flattening is simply due to the lack of low luminosity galaxies in the sample. From Fig. 5 it is evident that any given HI mass bin has the largest contributions from two adjacent luminosity bins, with smaller contributions from additional luminosity bins. For example, the log(M_HI) = 8.75 bin has equal density contributions from galaxies in the log(L_r) = 8.25 and 8.75 bins. However, as we approach the low mass end of the HIMF, these contributions to the HI mass density from adjacent luminosity bins run out. If we exclude the poorly sampled log(L_r) = 6.75 bin, the HIMF’s slope becomes positive at the low mass end, tracing the rapidly decreasing volume density of the log(L_r) = 7.25 bin.

As a result of the bias toward high gas-to-light ratio galaxies combined with the lack of available low luminosity sources, we find a flattening of the faint end slope of the HIMF. However, as we show next, we can make plausible assumptions to correct for the HI mass distributions of the undetected galaxies, as well as the HI mass distributions for galaxies in low-luminosity bins not probed by NIBLES (see Sect. 5.3) and examine their resulting effects on the HIMF (see Sect. 5.4). To do this, we examine trends in the log(M_HI/L_r) ratios as a function of luminosity to extrapolate possible distributions in the poorly-sampled and non-sampled low-luminosity bins.

5.3. Gas-to-light ratio distributions

The distributions of gas-to-light ratios as a function of luminosity are well constrained in the NIBLES sample. In Fig. 8, top panel, and also in Fig. A.1 we show that their gas-to-light ratios as a function of r-band luminosity follow a fairly consistent trend across the entire luminosity range. A least-squares fit to the mean log(M_HI/L_r) values of each 0.5 dex wide luminosity bin gives the relation log(M_HI/L_r) = −0.33 · log(L_r) + 2.79.

While the mean of the gas-to-light ratios in each bin show a trend with luminosity, the same is not true for their standard deviation or skewness (middle and bottom panels of Fig. 8, respectively). The standard deviation values are all similar (0.39 ± 0.02 on average) and do not depend on luminosity, while the skewness values display much more scatter (–0.6 ± 0.3 on average); the given uncertainties are the standard deviations across all the luminosity bins.

The skewness values are all rather small, going from about −0.4 for bin 7.75 to −1.0 for bin 10.25 and could be an indication of a weak trend. We did not consider bins outside the log(L_r) range 7.75–10.25 range as the 10.75 bin is incompletely sampled and in the bins below 7.75 the galaxies lie only in the upper halves of the bins (see Fig. 8). The lack of uniformity in these bins affects their measured values for the mean, standard deviation and skewness when compared to the more uniformly sampled bins. We use the consistent traits across luminosity bins log(L_r) = 7.75–10.25 to construct a reasonable extrapolation of an expected BLF that extends to lower luminosities and compare the results with the HIMFs of Z03 and H11.

5.4. Reconstructed HIMF

Due to the detection limits of both the SDSS and NRT as discussed in Sect. 5.2, the NIBLES HIMF displays an approximately flat low-mass slope. However, since we understand the nature of the data, we can use the observed trend in mean gas-to-light ratio with luminosity together with reasonable estimates of skewness and standard deviation values to create extrapolated gas-to-light distributions in the luminosity bins that lie outside our sample.

The distribution of gas-to-light ratios in each luminosity bin is well described by a Gaussian with negative skewness and an offset corresponding to the slope of the log(M_HI/L_r) vs. log(L_r) relation in the previous section, i.e., ∼0.3 · log(L_r). As with the discussion in the previous section, Gaussian fits to the log(M_HI/L_r) ratios from the BLF yield consistent values of variance and skewness across all the luminosity bins from Fig. 8. We extrapolate these trends down to log(L_r) = 5.25 due to the increase in average gas-to-light ratio with decreasing luminosity. Specifically, the density of the log(M_HI)=7.25 bin contains a significant contribution from the log(L_r) = 5.25 bin. Additionally, these extrapolated bins are necessary to compensate for the missing low luminosity galaxies as mentioned in Sect. 5.2.

We construct the extrapolated gas-to-light distributions such that within each luminosity bin we have:

$$\begin{array}{c}\hfill {\mathfrak{R}}_j={\int }_{\mathrm{\Delta }M/2L}^{-\mathrm{\Delta }M/2L}\varphi \text{and}{\mathfrak{F}}_k={\int }_{minM/L}^{maxM/L}\varphi ,\end{array}$$(13)

where ${\mathfrak{R}}_{\mathfrak{j}}$ is the jth gas-to-light ratio bin, ϕ represents the skewed Gaussian function, and ${\mathfrak{F}}_{\mathfrak{k}}$ is the detection fraction in the kth luminosity bin. ΔM/L is the log(M_HI/L_r) bin size.

Uncertainties for the reconstructed bins are calculated from the quadrature sum of the uncertainties in the fit parameters of the skewed Gaussian function.

Once we have reconstructed a gas-to-light distribution for each luminosity bin, we use it to generate the HI mass distribution and recalculate the BLF as outlined in Sect. 4.

In the NIBLES sample, the detection fraction increases with luminosity up to the log(L_r) = 9.25 bin, which is consistent with sensitivity and distance effects (see Appendix A). To analyze the effects of this change in detection fraction with luminosity, we construct two extrapolated BLFs. The first, ignoring the change in detection fraction with luminosity (our optically corrected sample), and the second taking it into account (our corrected sample). For the optically corrected BLF, we set the detection fraction of all luminosity bins below 7.75 equal to that of the 7.75 bin, since this is our lowest luminosity, well sampled bin. We make no attempt to extrapolate a decreasing detection fraction with decreasing luminosity due to the scatter in the observed trend (see Appendix A). For the corrected BLF, we set the detection fraction for all bins below 9.25 equal to the detection fraction of the 9.25 bin because the observed fall-off in detection fractions below this bin are consistent with decreases caused by sensitivity and distance effects.

In Fig. 9 we show (in gray) the reconstructed HIMFs per luminosity bin derived from the optically corrected BLF, together with the reconstructed HIMF from the corrected BLF for comparison. The shaded uncertainty regions from Fig. 5 are omitted here for clarity. The low-luminosity sources shown here in gray have the effect of boosting the low-mass slope of the HIMF to a level agreeing with the observed slopes in recent blind HI surveys. However, even using the extrapolated bins down to log(L_r) = 5.25, there is an apparent turn-over beginning at the log(M_HI) = 7.25 bin due to the lack of sources at even lower luminosity. This is illustrated by the log(M_HI) = 7.25 point falling below the HIMF Schechter fit in Fig. 9. These extrapolations show the need to probe extremely low luminosity sources when attempting to construct an HIMF from optically selected sources.

The need to probe extremely low luminosity galaxies is not unique to our survey: it applies to all surveys used to derive an HIMF. However, as blind HI surveys are not attempting to recover a two-dimensional HI-optical BLF distribution using an HI survey of optically selected galaxies like NIBLES, their lack of data at the low luminosity end is not problematic in this regard.

The difference in low-mass slope between the optically corrected and corrected BLF HIMFs is evident. The approximately 10% difference in detection fraction between the log(L_r) = 9.25 and 6.75 bins manifests itself as an approximately 6% steeper low-mass slope in the HIMF. The Schechter fit parameters for the two HIMFs are:

optically corrected: Φ = 0.0085 ± 0.0015, log(M_⋆) = 9.72 ± 0.06, α = − 1.26 ± 0.04

corrected: Φ =0.0066 ± 0.0012, log(M_⋆) = 9.79 ± 0.07, α = − 1.37 ± 0.03

These results imply that the fall-off in detection fraction for luminosity bins below log(L_r) = 9.25 is due primarily to sensitivity and distance effects.

For our corrected BLF, volume density values for luminosity bins of 8.25 and below are listed in Table 2 as log(Φ(M_HI,L_r)ΔMΔL ) in units of Mpc⁻³ dex⁻¹, together with their fractional uncertainties.

6. Discussion

In addition to being useful tools for measuring the baryon density of the local universe, both the HIMF and optical LF have also been instrumental in providing constraints on galaxy formation models (see, e.g., Lagos et al. 2011; Duffy et al. 2012; Davé et al. 2013). Currently these constraints are used independently, but our bivariate M_HI-L_r distributions provide cross-constraints on these models that will provide additional insights into the physical processes necessary to realize these properties. Additionally, the small variation in the observed standard deviation and skewness values of the gas-to-light ratios in each luminosity bin can provide additional modeling constraints.

As ours is the first optical luminosity-HI mass bivariate distribution constructed, we have run checks on its consistency through comparisons with published HIMFs derived in other studies. Because each of these studies has its own design, there are a number of factors that need to be accounted for when making these comparisons. We begin by comparing our HIMFs with those of blind HI surveys in Sect. 6.1, followed by a comparison between expected gas-to-light ratio distributions from our BLF and recently studied almost dark galaxies in Sect. 6.2. We compare our uncorrected HIMF with others from optically selected samples in Sect. 6.3, followed by a discussion and comparison of our Ω_HI with those from other surveys.

6.1. Comparison with HIMFs from blind HI surveys

Blind HI survey HIMFs have become one of the standard litmus tests used in comparisons of cosmological simulations with observations. Here (see Fig. 10) we compare our results from the optically corrected BLF to those of Z03 and H11.

Three differences between the surveys have either negligible effects or require simple adjustments. First, although H11 uses peculiar velocity corrections to the Hubble flow, this difference with our method has a negligible impact on the comparison since the largest discrepancies will occur for recession velocities below the minimum NIBLES velocity of 900 km s⁻¹. Second, we have applied offsets to the Z03 and H11 HIMFs (+0.07 and +0.11 dex in log(M_HI), respectively) corresponding to the flux scale differences between those surveys and NIBLES to bring them into agreement with the O’Neil (2004) flux standard – see Paper I for details. Third, we have rescaled the Z03 data to the 70 km s⁻¹ Mpc⁻¹ Hubble constant adopted in this paper and Paper I.

There is good agreement in Fig. 10 between our optically corrected HIMF and those of Z03 and H11. Although the low-mass data (log(M_HI) = 7.5 − 9.5) show small systematic vertical offsets, the slope of our optically corrected HIMF is in very good agreement with those of the other two, only differing by about 0.01 (see also Sect. 5.4). The agreement between the HIMFs is brought about by our re-scaling of the Z03 and H11 data flux scales.

Despite the suggestion in Martin et al. (2010) and H11 that ALFALFA found a 1 dex greater volume density of high-mass sources than the HIPASS survey, after adjusting both to the same calibration standard, the high-mass data points appear to be in good agreement as shown in Fig. 10. HI surveys have become sufficiently large with statistically significant numbers of rare high-mass galaxies that calibration differences have a significant impact on their interpretation.

The NIBLES optically corrected HIMF was constructed from trends observed in the gas-to-light ratio distributions, utilizing an extrapolation for galaxies that lie outside the observable limits (see Sect. 5.3) of the NIBLES sample. This extrapolation appears justified by the fact that the reconstructed low-mass volume density agrees well with those of the Z03 and H11 blind surveys. To achieve this agreement we included extrapolations down to the log(L_r) = 5.25 bin (see Fig. 5). This corresponds to an r-band absolute magnitude of –8.5, i.e., in the same range as globular clusters. The need to extrapolate to such low-luminosity galaxies indicates that optically selected surveys would have to achieve impractically high sensitivity to accurately sample the low-mass end of the HIMF, as will be discussed further in Sect. 6.3.

We discussed in Sect. 5.2 how the sensitivity limits prohibiting detection of very low luminosity objects suppresses the low-mass slope of the HIMF. This is true even for the optically corrected BLF, which agrees well with the HIMFs of Z03 and H11. After fully correcting for the sensitivity effects discussed in Sects. 5.1 and Appendix A, the resulting corrected BLF produces an HIMF with an approximately 6% steeper low HI-mass slope, in agreement with that found in Zwaan et al. (2005). This steeper slope reflects the contribution of some extremely low luminosity HI-rich objects, examples of which have been encountered in several HI surveys.

6.2. Almost dark galaxies

Blind HI surveys have detected a number of low-luminosity but HI-rich galaxies, some of which initially appeared to be optically dark sources. However, the HI and optical properties of these sources appear to be consistent both with the gas-to-light ratios predicted by our BLF, as well as the upper-most limits of M_HI/L found in other surveys (e.g., Rosenberg et al. 2005; Doyle et al. 2005).

As noted before, the HI detectability of ALFALFA is comparable to that of NIBLES, but as a blind HI survey ALFALFA can be used to search for extragalactic sources without optically cataloged counterparts, and it can detect very nearby objects since it has no lower radial velocity limit. Follow-up studies in both the optical and with radio interferometers have been made of ALFALFA detections that initially appeared to be (nearly) without optical counterparts. These comprise ALFALFA (Almost) Dark Galaxies Project candidates including SHIELDS (Survey of HI in Extremely Low-mass Dwarfs) objects, and the Leo P dwarf. Deep optical images usually revealed a likely optical counterpart, but not always (see, e.g., Cannon et al. 2011, 2015; Giovanelli et al. 2013; Rhode et al. 2013; Janowiecki et al. 2015; McNichols et al. 2016).

Among the ALFALFA (Almost) Dark Galaxies pilot surveys, we compare the gas-to-light ratios of the two isolated objects from Cannon et al. (2015) and the sole optically detected galaxy from Janowiecki et al. (2015) with those from the corrected BLF. In our BLF, both of the Cannon et al. (2015) galaxies, with luminosities in the log(L_r) = 7.25 and 6.75 bins, fall in the log(M_HI/L_r) = 1 bin. According to our two-dimensional corrected BLF distribution, we would expect to find ∼30–40% of HI detected galaxies of similar luminosities in this gas-to-light ratio bin. The optical detection from Janowiecki et al. (2015) (in the log(L_r) = 7.25 and log(M_HI/L_r) = 1.5 bins), is where we would expect to find only ∼2% of galaxies of similar luminosity. However, it should be noted that the M_HI/L_r ratio of the latter source is likely to be overestimated, as it clearly contains optical flux outside the measurement apertures used in the paper.

Correlating the α.40 catalog to optical sources from the SDSS DR9 give results that are in reasonable qualitative agreement with our BLF. We only used sources with r-band petrosian radii larger than 3^″, as sources with smaller radii tend to have unreliable magnitudes. This optical correlation results in ∼25% and 35% of galaxies in the log(L_r) = 6.75 and 7.25 bins, respectively being contained in the log(M_HI/L_r) = 1 bin. For the log(L_r) = 7.25 bin, ∼6% of the galaxies are in the log(M_HI/L_r) = 1.5 bin.

An order of magnitude smaller gas-to-light ratio (log(M_HI/L_r) = 0.2) was found for the nearby Leo P dwarf, at a distance of only 1.6 Mpc. With a log(L_V)∼5.8 it is one of the faintest optical ALFALFA detections. It also has one of the lowest HI masses observed at log(M_HI) ∼ 6 (see, e.g., Giovanelli et al. 2013; Rhode et al. 2013; Bernstein-Cooper et al. 2014; McQuinn et al. 2015). The gas-to-light ratio of the Leo P dwarf is also consistent with our corrected BLF. Based on our extrapolations to such low luminosities, roughly 27% of the population of these sources would have gas-to-light ratios in the log(M_HI/L_r) = 0 bin or lower.

From the comparison with our BLF, it appears that the four (Almost) Dark Sources discussed above have fairly common gas-to-light ratios for their luminosities. The reason they are (almost) dark is simply that at such low luminosities the mean gas-to-light ratio grows quite large, making these sources relatively easier to detect in HI surveys.

To illustrate this point, a galaxy with a high gas-to-light ratio (log(M_HI/L_r) ∼ 1.5) with log(L_r) ∼ 5 would be detectable to only about 2 Mpc in the SDSS but out to 13 Mpc in HI. This comparison is based on the assumption of W₅₀ = 36 km s⁻¹, an rms of 2.3 mJy and a signal to noise ratio of three, and an estimated SDSS r-band galaxy detection limit of 17.77 mag (see Loveday 2002); for NIBLES, we also noted that galaxies fainter than this limit often have unreliable magnitudes (see Paper I for details).

These differences in detectability in optical vs. HI for these high gas-to-light ratio galaxies explains the absence of an optical counterpart in HI surveys. While these galaxies make interesting case studies for galaxy formation and evolution, they so far do not appear to contain exceptionally high gas-to-light ratios given their low luminosities.

6.3. Comparison with HIMFs from optical samples

NIBLES is the first HI survey to use optically selected sources that were subsequently followed up with a uniform set of HI observations. This is in contrast to previous works such as R93 and S05 which used a combination of archival optical and HI data.

R93 used published conversions from optical luminosities to HI masses for different classes of galaxies to derive an HIMF from existing literature sources, while S05 used 2771 spiral and irregular galaxies from the Uppsala General Catalog of Galaxies (UGC; Nilson 1973) with optical diameters greater than 1′. The closest equivalent in our study is the uncorrected HIMF, which does not attempt to correct for the missing low luminosity sources.

The HIMFs obtained by these differing approaches can be seen in Fig. 11. The S05 data correspond to the ∑1/V_max method results in that paper, and for the R93 fit we converted the HIMF from their Fig. 5 from HI mass density to the now commonly used number density and scaled their data to the Hubble constant of 70 km s⁻¹ Mpc⁻¹ used in NIBLES.

At low HI-masses (log(M_HI) < 8.5) the density values are relatively similar for all three studies. However, near the knee of the Schechter functions, above log(M_HI) ∼ 8.5, a discrepancy appears with the NIBLES uncorrected HIMF having a higher volume density than the other two. Considering that the NIBLES uncorrected HIMF has density values in good agreement with the blind HI surveys, for the high HI-mass galaxies (log(M_HI) ∼ 8.5−10), this comparison suggests that the morphology and color criteria of R93 and S05 introduced biases. The NIBLES approach, which did not select by either morphology or color, results in a significantly more complete sampling of the volume density of even luminous, high-HI-mass galaxies than selecting optical sources generally expected to contain HI.

For the low HI-mass galaxies (log(M_HI) < 8.5) however, the similarities between the low HI mass volume densities in all three studies combined with the analyses presented in Sects. 5.3 and 6.1 suggest that optically selected samples generally lack sufficiently low-luminosity sources to construct an HIMF with the same low-mass slope as derived from blind HI surveys.

6.4. Cosmic density of HI, Ω_HI

We can use our HI Mass Function to calculate the cosmic mass density of HI, Ω_HI. We find Ω_HI = 4.2 ± 0.2 from the corrected HIMFs, in units of ${10}^{-4}{h}_{70}^{-1}$. This differs by only 0.1 from the value obtained from our uncorrected HIMF, and it is also very similar to the values obtained from other surveys: 4.0 ± 0.5 derived from HIPASS detections by Z05 and 4.2 − 4.4 ± 0.1 derived from the summation of binned measurements by Martin et al. (2010) and H11. The only measured 21cm HI line Ω_HI value at a significantly higher redshift, z ∼ 0.32, is 5.0 $\pm 1.8\times {10}^{-4}{h}_{70}^{-1}$ (Rhee et al. 2016) is based on stacked GMRT HI spectra of 165 galaxies with known optical redshifts. So the 21 cm Ω_HI remains constant within the estimated uncertainties out to at least z ∼ 0.3, but at a level which is only about half of that measured from Damped Lyman Alpha (DLA) systems at z > 2 (see Rhee et al. 2018, and references therein).

Theoretical studies of the cosmic evolution of Ω_HI have yielded a range of results. Since our estimate, while agreeing with those from other studies (see above, e.g., Z05), has formally smaller uncertainties, it can thus be used to constrain theoretical studies. For example, the theoretical HIMF of Popping et al. (2015) over-predicts the number of low HI mass galaxies compared to our results, but this discrepancy is less compared to the HIMFs of Martin et al. (2010) or H11. As a consequence, Popping et al. (2015) over-predicts Ω_HI by about a factor of four. On the other hand, more recent models, such as those of Davé et al. (2017), under-predict Ω_HI by about 30% compared to our observed results. The latter difference between an observational and a theoretical estimation is significant at the ∼6σ level (we have a ≈5% uncertainty). The discrepancy likely arises due to the model of Davé et al. preferentially under-predicting the HI gas fraction of high mass galaxies (log(M_⋆) ≳ 10.5) compared to observations, like our NIBLES results (Paper I).

Underlying models predicting the HI content can be divided into two basic types of prescriptions – either pressure-based for regulating the phase relative fractions of atomic to molecular gas (Blitz & Rosolowsky 2004), or dependent on the intensity of the interstellar radiation field molecular cloud opacity (Krumholz et al. 2009). While both theoretical approaches explain certain aspects of the gas and star formation properties of galaxies, it is not clear that they provide appropriate descriptions of the integrated properties of the ensemble of galaxies or even, in phenomenological models of galaxy formation and evolution, that their discrepancies with observed HIMFs and Ω_HI values are mainly driven by problems with these theoretical prescriptions.

One thing is certain, determining a robust and accurate characterization of the local HI mass function and Ω_HI through the BLF, as we have done, does provide additional needed constraints on the gas cycle in galaxies – how gas is acquired, what processes regulate the gas phases within galaxies, and how gas is lost – and thus on the physics of galaxy evolution.

7. Future work

The Bivariate Luminosity Functions presented in this work are based on HI detection fractions of optically selected sources. These detection fractions drop off as a function of decreasing gas-to-light ratio, as one would expect. However, there is a currently unknown HI distribution lying below the detection threshold of NIBLES which will expand the plane of the BLF if probed. We are currently working on implementing a BLF incorporating four times more sensitive follow-up observations of the NIBLES undetected sample with the Arecibo L-band wideband receiver (see also the pilot survey results in Butcher et al. 2016).

Further investigation of the BLF presented here is possible using the ALFALFA survey data, even though it has about the same HI sensitivity as NIBLES. The low luminosity and low HI mass portion of the BLF plane could be explored using the large statistics provided by the ALFALFA sources. However, the optical properties of many of these low luminosity sources are not accurately measured in SDSS and would require the re-processing of all photometric sources with r-band Petrosian magnitudes greater than 17.77. Additionally, follow-up observations (utilizing the L-band wideband receiver) of low luminosity sources not detected in ALFALFA would also be required to accurately sample the BLF plane in this parameter space.

The NIBLES BLF can also be used to construct a trivariate luminosity–HI-mass–dynamical mass function of galaxies utilizing line widths to infer the dynamical mass. This trivariate function could be used to examine what role the dark matter halo plays in regulating galaxy evolution.

8. Conclusions

We have developed the first optical luminosity-HI mass BLF based on the NIBLES sample of SDSS galaxies in the local volume (900 km s⁻¹ < cz <12 000 km s⁻¹), selected to span their entire range of stellar masses. The HIMF constructed from our uncorrected BLF, which was derived using only HI detections, agrees very well with other HIMFs constructed from optically selected samples.

We observed a very consistent gas-to-light ratio distribution (with similar standard deviations and skewness) displaying a progressive increase in the mean log(M_HI/L_r) value with decreasing luminosity across the entire luminosity range of the NIBLES sample, which may suggest that galaxy populations in all luminosity bins undergo the same evolutionary processes.

Using the predictable offset in mean gas-to-light ratio and the consistency in gas-to-light distribution as a function of luminosity, we constructed an optically corrected BLF by extrapolating these properties in luminosity bins down to log(L_r) = 5.25. Estimating bivariate densities down to this low luminosity results in a corrected HIMF with the same low mass slope as those derived from the blind HI surveys Z03 and H11. This suggests that the low mass end of the HIMF may contain significant contributions from galaxies with very low luminosities.

We also observed a ∼10% decline in detection fraction from log(L_r) = 9.25 to the lowest luminosity sources due to distance effects and sensitivity. Correcting for this decline leads to a slightly more negative low mass slope of the HIMF, in agreement with the Zwaan et al. (2005) HIPASS results.

The combination of sensitivity and minimum radial velocity cut-offs makes the detection of very low HI mass sources difficult for both the NIBLES and HIPASS surveys. This is reflected in the Zwaan et al. (2005) HIPASS HIMF low mass slope, devoid of sources below log(M_HI) = 7. The NIBLES corrected HIMF, on the other hand, has a low mass slope obtained by accounting for the above mentioned trends in our data as well as observational limitations. We are therefore able to improve on HIMFs from previous optically selected HI surveys and produce HIMFs that are in agreement with those of other blind HI surveys.

We note that based on the differences in detectability in the optical vs. HI for relatively high gas-to-light ratio sources with low luminosity, the optical counterparts of such HI detections may be very difficult to detect. Indeed, this appears to be the case with the ALFALFA (almost) Dark Galaxies observed to date.

Our corrected BLF and resulting HIMF explain why HIMFs tend to have steeper low mass slopes than the faint end of optical LFs: the increasing mean gas-to-light ratios of the lower luminosity sources cause a source’s HI mass to be counted in higher HI mass bins relative to where it would be counted in an optical luminosity function. This increases the density relative to luminosity, resulting in a boosted low-mass slope of the HIMF.

selected samples like NIBLES appear to lack a significant percentage of low luminosity objects that are detectable in HI, whereas our optically corrected and corrected BLFs can reproduce blind HI survey HIMFs, and our Ω_HI agrees with other surveys. All this gives us confidence that our corrected BLF provides an accurate representation of the galaxy population in the local universe. Therefore, it can be used to provide further constraints on the gas cycle in galaxies and insights into galaxy evolution.

Investigation of the low luminosity (log(L_r) < 7) regime of the BLF with higher sensitivity optical surveys could further confirm the accuracy of our extrapolated gas-to-light ratio distributions.

Acknowledgments

The Nançay Radio Telescope is operated as part of the Paris Observatory, in association with the Centre National de la Recherche Scientifique (CNRS) and partially supported by the Région Centre in France. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. We thank the anonymous referee for their helpful comments and suggestions for improving this paper.

References

Appendix A: NIBLES detection fractions

The fraction of galaxies detected for NIBLES peaks at the log(L_r) = 9.25 bin and decreases from there with both increasing and decreasing luminosity (see Fig. 7). The decrease in detection fraction with increasing luminosity is to be expected since more luminous galaxies tend to be relatively gas poor compared to their lower luminosity counterparts. However, the decrease in detection fraction with decreasing luminosity is somewhat counter-intuitive since the average gas-to-light ratio increases with decreasing luminosity. However, as mentioned in Sect. 5, this can be explained by a combination of the NRT sensitivity and the NIBLES selection criteria. We give a more detailed description of this phenomenon here.

We test the effect of our selection criteria on our detection fraction as follows.

In Fig. A.1 we show the gas-to-light ratio (loglog(M_HI/L_r)) as a function of L_r for NIBLES detections. Our fit to the data is shown in black, along with upper and lower bounds around 92% of the data in red. For each galaxy, we randomly assign an HI mass based on its luminosity that falls within these two bounds. HI surveys in general do not have a set minimum detectable HI line flux due to the dependence of detectability on line width. NIBLES has the added complication that we do not have a uniform rms noise for each source either. Therefore, for illustration purposes, we assume a 0.2 Jy km s⁻¹ minimum line flux for detectability, which is consistent with the lowest fluxes detected in NIBLES. This gives us a minimum detectable HI mass as a function of distance.

We then compare this minimum detectable HI mass to the assigned artificial HI masses based on the observed log(M_HI/L_r) distribution shown in Fig. A.1. Figure A.2 shows an example of our artificially generated HI masses as a function of distance, with the blue line indicating our minimum detectable HI mass. The fraction of galaxies that lie above their corresponding minimum detectable mass in each luminosity bin gives us our artificial detection fractions. We ran 100 iterations of the artificial mass distributions and used the mean values of each luminosity bin to derive our artificial detection fractions as a function of luminosity. These detection fractions are shown in Fig. A.3. Note that these fractions are only meant to illustrate the decreasing detection fractions of low mass galaxies as a function of luminosity due to distance effects alone, not to accurately describe the effects of HI depletion in larger, more luminous galaxies. There is roughly a 10% decrease in detection fraction for the sources with log(L_r) ≤ 8.25 compared to those with log(L_r) ≥ 9.25. This is similar to the approximately 10% observed decrease in detection fraction discussed in Sect. 5.1 as a function of decreasing luminosity.

Since we can construct artificial HI mass distributions for the NIBLES sample that show a similar decrease in detection fraction with decreasing luminosity as the observed detection fractions, we can conclude that it is likely that the approximately 10% decrease in observed detection fractions is due, at least in part, to distance effects only. We can therefore correct for this in the BLF.

All Tables

All Figures