THE ARECIBO LEGACY FAST ALFA SURVEY. X. THE H i MASS FUNCTION AND $\Omega _{\rm H\,{\mathsc{i}}}$ FROM THE 40% ALFALFA SURVEY

The ongoing ALFALFA survey takes advantage of the new multipixel ALFA receiver at the Arecibo Observatory. When complete, the survey will have measured >30, 000 galaxies in the 21 cm line out to z ∼ 0.06 with a median redshift of ∼8000 km s⁻¹. The survey is more sensitive than HIPASS, with a 5σ detection limit of 0.72 Jy km s⁻¹ for a source with profile width 200 km s⁻¹ in ALFALFA compared to a 5σ sensitivity 5.6 Jy km s⁻¹ for the same source in HIPASS (Giovanelli et al. 2005). Narrow profile widths, down to ∼15 km s⁻¹, allow us to probe extremely small objects. ALFALFA detects objects with neutral hydrogen masses $M_{\rm H\,{\mathsc{i}}} \sim 3 \times 10^{7}$ M_☉ out to the distance of the Virgo Cluster. In addition to greater sensitivity, ALFALFA probes gas-rich galaxies in the local universe with greater velocity resolution (11 km s⁻¹ after Hanning smoothing versus 18 km s⁻¹) and a deeper limiting redshift (18,000 km s⁻¹ versus 12,700 km s⁻¹) than HIPASS. Our significantly improved survey depth for low-mass objects allows the ALFALFA survey to better constrain the low-mass slope of the H i mass function.

ALFALFA survey data are acquired in a minimally invasive drift scanning mode, in two passes ideally separated by several months, and individual 600 s drift scans are combined into three-dimensional data grids covering 24 in both R.A. and decl.; it therefore takes many nights of observations to complete a grid from which extragalactic sources can be extracted.

Confidently detected sources are assigned one of the three object codes, where Code 1 refers to a reliable extragalactic detection with a high signal-to-noise ratio (S/N >6.5), Code 2 refers to extragalactic sources with marginal S/N (4.5 < S/N < 6.5) confirmed by an optical counterpart with known optical redshift matching the H i measurement, and Code 9 refers to high velocity clouds (HVCs). For this analysis, we consider only objects designated Code 1, since we are interested in extragalactic objects with well-known selection criteria. Code 1 objects have a reliable S/N, a good match between the two polarizations that are independently observed by ALFALFA, a clean spectral profile and, in almost every case, a confident match with an optical counterpart. The signal detection pipeline, discussed at length in Saintonge (2007), combines a matched-filtering technique for identifying source candidates with an interactive process for source confirmation and parameter measurement. This technique is estimated to result in a reliability of candidate detections of ∼95% for Code 1 objects, with a completeness better than 90% for the narrowest galaxies above the prescribed S/N threshold. The subsample of Code 1 objects provides a robust sample for the HIMF.

Published ALFALFA catalogs contain a set of measured parameters (including coordinates, heliocentric velocity, line profile velocity width W50 measured at the 50% level of two profile peaks, integrated flux density S_int, S/N, and noise figure σ_rms) in addition to a distance estimate and a derived H i mass in solar units, obtained from the expression $M_{\rm H\,{\mathsc{i}}} = 2.356 \times 10^{5} D^{2}_{\rm Mpc} S_{\rm int}.$ Our distance estimates are subject to errors due to each galaxy's unknown peculiar velocity, which translate into mass errors. The fractional distance error due to peculiar velocity decreases with increasing distance (the so-called Eddington effect), so the lowest-mass galaxies which are only found nearby are most prone to this error, our treatment of which is discussed in detail in Section 3.2.

ALFALFA's ability to detect a signal depends not only on the integrated flux, but also on the profile width W₅₀ (km s⁻¹). Figure 1 displays the distribution of sources detected by ALFALFA. Rather than a single flux limit, the ALFALFA detection threshold is dependent on both S_int and profile width W₅₀, and we find that this relationship changes above W₅₀∼ 400 km s⁻¹. We fit the S_int,th relationship empirically to the data, rather than using the assumed expression above. Due to differences in the two methods we employ to calculate the HIMF, we consider two different threshold cuts, described separately in Sections 3.3 and 3.4.

Zoom In Zoom Out Reset image size

ALFALFA catalogs have been extracted for a large contiguous region in the southern Galactic hemisphere (i.e., anti-Virgo direction; 22^h <α< 03^h, 24° < δ<32°) and two regions in the northern Galactic hemisphere (i.e., Virgo direction; 16^h30^m <α< 07^h30^m, 4° < δ< 16°, and 24° < δ< 28°), with coverage totaling 2607 deg² or ∼40% of the final ALFALFA volume. This includes the previously published catalogs with a total of 2706 extragalactic source measurements (Martin et al. 2009; Stierwalt et al. 2009; Kent et al. 2008; Saintonge et al. 2008; Giovanelli et al. 2007) in addition to an upcoming large online data release (M. P. Haynes et al. 2010, in preparation).⁶ This primary data set includes both Code 1 (n = 10452) and Code 2 (n = 2750) galaxies in addition to Code 9 (n = 629) HVCs, where this figure includes measured subcomponents of larger cloud complexes.

From the primary data set, we have selected the 40% ALFALFA survey sample, hereafter α.40. This sample has been selected to include only Code 1 objects, and the total sample size is further reduced by the exclusion of galaxies found beyond 15,000 km s⁻¹, where radio frequency interference from the Federal Aviation Administration (FAA) radar makes ALFALFA blind to cosmic emission in a spherical shell ∼10 Mpc wide. The final α.40 sample contains 10,119 Code 1 galaxies, for a detection rate of 3.9 deg⁻² compared with the HIPASS detection rate of ∼0.2 deg⁻² (5317 extragalactic sources over 29,000 deg²; Meyer et al. 2004; Wong et al. 2006). While rich in absolute number, HIPASS does not extend deep enough in redshift to sample a cosmologically fair volume.

In Figures 2 and 3, we present the redshift distribution of the 10,119 Code 1 objects in α.40 as a set of cone diagrams by region in the survey. The two most obvious features in Figure 2 are the prominent void in the foreground of the Pisces–Perseus supercluster, leading to the dearth of detections out to about 3000 km s⁻¹, and the portion of the main ridge of that supercluster that cuts across the diagram. In the top panel of Figure 3, the nearby Virgo Cluster is prominent, as is the Coma supercluster. ALFALFA probes a wide variety of environments in the local universe and will soon study the overall properties of H i-selected galaxies as a function of environment (A. Saintonge et al. 2010, in preparation).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 displays histograms of the statistical properties of the α.40 sample. From panels (a)–(d), these histograms represent the heliocentric velocity, velocity width W₅₀, integrated flux S_int, and S/N properties. In particular, note that the S/N is high for all detections, since Code 2 objects have been excluded from this analysis. For clarity, the histogram of the H i masses of galaxies in the sample is plotted separately, in Figure 5. On the low-mass end, where ALFALFA can place strong constraints on the faint-end slope of the HIMF, the α.40 sample contains ∼340 galaxies with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) < 8.0$ and ∼114 galaxies with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) < 7.5$; on the high-mass end, which is best probed by surveys with deep redshift limits, there are ∼50 galaxies with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) > 10.5$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The large sample size of ALFALFA, extending over a range of H i masses, is one of its key strengths in relation to the problem of characterizing the density of neutral gas in the present-day universe. With such a large number of galaxies, we can approach our calculation of the HIMF in two distinct ways. First, using the entire sample and a well-known characterization of our sensitivity, we can apply corrections and obtain the overall function without excluding sources. Second, however, we can make stringent integrated flux cuts and use only those galaxies bright enough to be detectable irrespective of other properties (e.g., profile width). The sample contains ∼3500 galaxies with an integrated flux >1.8 Jy km s⁻¹, which provides a strict cut above which our objects are detected regardless of profile width. This subsample size is comparable to the full sample size for previously published HIMFs such as HIPASS, but samples a fair cosmological volume. This subsample, referred to hereafter as α.40_1.8, provides a test case for analyzing the quality of the HIMF measurement for the full α.40 sample. The precise details of the calculation, of ALFALFA's sensitivity, and of the corrections applied to the HIMF calculated from α.40, make up the bulk of the following sections and of Appendices A and B.

The H i mass function, like galaxy luminosity functions, is usually parameterized as a Schechter function of the form

$$\begin{equation} \phi (M_{\rm H\,{\mathsc{i}}}) = \frac{dn}{d \log M_{\rm H\,{\mathsc{i}}}} = \ln 10 \ \phi _{*} \ \left(\frac{M_{\rm H\,{\mathsc{i}}}}{M_{*}} \right)^{\alpha + 1} \ e^{-\frac{M_{\rm H\,{\mathsc{i}}}}{M_{*}}}. \end{equation} \tag{ 1 }$$

The parameters of interest are the faint-end slope α, the characteristic mass log M_*, and the scaling factor ϕ_*.

$\phi (M_{\rm H\,{\mathsc{i}}})$ has historically been calculated in one of two ways. The Σ1/V_max method (Schmidt 1968) can be understood by analogy to a purely volume-limited sample, in which case the HIMF would be obtained by the galaxy counts divided by the total volume of the survey. The Σ1/V_max method treats each individual galaxy in this way, by weighting the galaxy counts by the maximum volume V_max,i within which a given source could have been detected. This weighting strategy allows the inclusion of low-mass galaxies, visible only in the nearby universe, in the same sample as rare high-mass galaxies, found only in larger volumes. Additionally, the weights may be adjusted in order to correct for a variety of selection effects, large-scale structure effects, and missing volume within the data set, so that a well-characterized survey can robustly measure the HIMF.

An alternative method, the 2DSWML approach, was applied to the HIPASS measurements of the HIMF (Zwaan et al. 2003, 2005). This method is designed to make the calculation of the HIMF less sensitive to local large-scale structure, since shallow blind H i catalogs are contaminated by the richness of the Local Supercluster. If 1/V_max is used without correction for this overdensity, the resulting HIMF will overestimate the contribution by low-mass galaxies and steepen the faint-end slope α. Stepwise maximum likelihood (SWML) methods, by contrast, are designed to reduce this bias, by assuming that the shape of the HIMF is the same everywhere and then obtaining the $\phi (M_{\rm H\,{\mathsc{i}}})$ that maximizes the probability of the observed distribution (Efstathiou et al. 1988). Given the dependence of the ALFALFA survey's sensitivity on both mass and profile width (Section 2.3), a 2DSWML approach is necessary to calculate the HIMF for the full sample (Loveday 2000). 2DSWML maintains the main advantages of the SWML method, which are its robustness against density fluctuations in the survey volume and its model-independent approach.

In this work, we apply both the 1/V_max and the 2DSWML methods for various reasons. Given our knowledge of our sample's characteristics and sensitivity, the 1/V_max method is simple to apply and straightforward to assess for potential bias. We can account for large-scale structure and other selection effects by applying well-motivated corrections (discussed in Section 3.3). Perhaps more significantly, this method also allows us to quantify and understand those effects on the ALFALFA survey. In particular, a goal of ALFALFA is to further probe the differences between H i mass functions in different environments; the 2DSWML assumption that the shape of the function is the same throughout a sample may not be valid. By contrast, the 2DSWML method is designed to be more resistant to effects from large-scale structure, and also results in a calculation of the selection function which can be used in future analysis of the sample via, for example, the two-point correlation function. A comparison of the 1/V_max and 2DSWML methods as applied to α.40 is considered in Section 6.

In both the 1/V_max and 2DSWML analyses, we have used five mass bins per dex and have found that the HIMF is not strongly affected by the choice of bin size. In the case of 2DSWML, we also bin by profile velocity width and find no significant difference for bin sizes between 2 and 20 bins per dex. The two main sources of error are counting statistics within the bins and mass errors.

Minimizing and taking into account distance errors are key to robust estimation of luminosity and mass functions, in particular at the faint end. Masters et al. (2004) considered how strongly distance uncertainties will tend to affect a given local volume survey's estimate of the faint-end slope of the mass function. In that work, the authors accounted for distance errors by constructing a mock catalog, with masses assigned from a chosen HIMF and with the spatial distribution determined from the density field of the IRAS Point Source Catalog Redshift survey (PSCz; Branchini et al. 1999). They concluded that a survey toward the Virgo Cluster, like a portion of the sample considered here, will overestimate distances to those galaxies if pure Hubble flow is used, since objects in that field are falling into Virgo. Since the H i mass depends on distance as D², this has serious consequences for the faint-end slope of the HIMF. Therefore, work in this region relies both on the development of well-constrained local velocity models from primary and secondary distance catalogs and on a careful consideration of the effects of distance uncertainties. We consider the Virgo Cluster as a special case of this general problem in Section 6.1.

These difficulties arise precisely because the lowest-mass objects can be detected only at small distances, so that the fractional distance errors due to deviations from Hubble flow most strongly affect the most interesting bins of the mass function. The best distance estimates, primary distances based on, e.g., Cepheids or the tip of the red giant branch, can only estimate distances to within ∼10% error, so beyond cz ∼ 6000 km s⁻¹ the uncertainties on distances obtained via a primary method and those obtained assuming pure Hubble flow become comparable, and the latter is typically used for simplicity. Within that distance, however, the distance uncertainties can have a very strong influence, up to 100% in the case of the Virgo Cluster.

To minimize distance uncertainties, the ALFALFA survey has adopted a distance estimation scheme that makes use of a peculiar velocity flow model for the local universe (Masters 2005). This parametric multiattractor model, based on the SFI++ catalog of galaxies with Tully–Fisher distances (Springob et al. 2007), includes two attractors (Virgo and a Great Attractor) along with a dipole and quadrupole component. Distances to almost all α.40 galaxies within 6000 km s⁻¹ are estimated from the flow model. Beyond cz_CMB = 6000 km s⁻¹, the model is not well constrained, so distances are estimated from Hubble flow (H₀ = 70 km s⁻¹ Mpc⁻¹). Within 6000 km s⁻¹, some galaxies have measured primary distances, which are applied in our scheme, and other galaxies are known to belong to a group, in which case the group's mean velocity is used for distance estimation. The Masters (2005) flow model also provides error estimates, constrained by the fit of the model to the observed velocity field and with a minimal error based on the local velocity dispersion 163 km s⁻¹. When distances are estimated using pure Hubble flow, the error is estimated to be ∼10% via the assumption that peculiar velocities are ∼ a few hundred km s⁻¹.

Mass errors for individual galaxies in our sample are calculated from the measured error on the integrated flux and an estimated error on the distance, which is the larger of the local velocity dispersion 163 km s⁻¹, the distance error estimate of the Masters (2005) flow model, or 10% of the distance. Because the mass error shifts galaxies into different bins of the HIMF, the relationship between these errors and the final HIMF parameter errors is complex. We deal with these errors by calculating several hundred realizations of the HIMF after randomly assigning flux and distance errors to each galaxy to find the spread in each mass bin.

There is a complication on the high-mass end of the sample as well. Arecibo's relatively large beam size at 21 cm (∼3.5 arcmin) can cause source confusion at large distances, where we also find our largest-mass objects. When this occurs, ALFALFA may be detecting more than one individual gas-rich galaxy as a single source, but in cases of interaction it is also possible that the galaxies involved are part of a single, large H i envelope. While higher-resolution follow-up would be required to fully resolve this issue, we have investigated optical images and redshift catalogs for the highest mass ($\log M_{\rm H\,{\mathsc{i}}} >$ 10.5) ALFALFA detections and have found that the majority of these objects are not likely to be blends of H i emission from an interacting system and some others are close pairs that are likely to share a single gas envelope.

For each galaxy in α.40, V_max,i is calculated based on that galaxy's H i mass M_i, the minimum integrated flux S_min,i at which such a galaxy is detected in ALFALFA, and finally the distance D_max,i corresponding to that limit. The calculated V_max,i, corresponding to the effective search volume for that galaxy, excludes volume that is not covered by ALFALFA, including volumes where detection ability, and therefore effective search volume, is reduced by the appearance of radio frequency interference at the corresponding frequency. Galaxies are binned by mass and $\phi (M_{\rm H\,{\mathsc{i}}})$ is calculated by summing the reciprocals of V_max.

By weighting the count for each galaxy, the 1/V_max method can be corrected for a variety of known systematic effects. The major corrections applied to the HIMF for this sample address (1) missing volume, (2) the profile width-dependent sensitivity of the survey, and (3) the known large-scale structure in the local volume.

Sources of radio frequency interference contaminate the signal in regions of frequency space corresponding to spherical shells in the survey volume. This effectively reduces the search volume of the overall survey. Figure 6 shows the average relative weight, compared to 100% coverage, within the α.40 survey volume as a function of velocity. The large dip between 15,000 and 16,000 km s⁻¹ is due to the FAA radar at the San Juan airport, and because of this extreme loss of volume at large distances we restrict the α.40 sample to only those galaxies within 15,000 km s⁻¹. Given our knowledge that the relative weight is less than 1.0 at specific distances, the V_max value calculated for a specific galaxy is reduced to reflect the loss of effective search volume. This correction is not significant for the lowest-mass galaxies, but more generally, the correction is very small. The effect on the final Schechter parameters for the HIMF is on the order of 2%.

Zoom In Zoom Out Reset image size

As discussed in Section 2.3, ALFALFA's detection ability is dependent on the profile velocity width of the signal, W₅₀, in km s⁻¹, rather than strictly on the integrated flux of the signal. To obtain an expression for this detection limit, we used the data itself, as displayed in Figure 1. The dependence of ALFALFA sensitivity on both flux and profile width, described in Section 2.3, has the further complication of affecting the survey's completeness, and this must be accounted for in order to extract the underlying HIMF. The distribution in Figure 1 indicates that ALFALFA finds many galaxies with low fluxes and narrow widths, but there is a deficiency of galaxies with low fluxes and large widths. Because we have no knowledge of the true distribution below ALFALFA's detection capability, we have developed a completeness correction that takes advantage only of the data, making no assumptions about the potentially intrinsically small unobserved population. The profile width completeness correction most strongly affects galaxies with ${\sim} 9.0 < \log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) < 10.0$ and has a very small influence (<2%) on both the faint-end slope α, since low-mass (i.e., narrow velocity width) galaxies are not affected, and log(M_*), since the counts in the high-mass bins are large enough to robustly constrain this. This is essentially a galaxy counting correction, so its primary influence is on ϕ_*, increasing that parameter by a factor of 20%. The full details of this completeness correction are described in Appendix A. The validity of this completeness correction, which we have applied to the full sample, is tested in Sections 4.1 and 5.1 by calculating the HIMF using an integrated flux cut, which allows us to neglect the biased sensitivity dependence on width. By comparing the resulting HIMF in both cases, we assess the impact of this correction.

The most significant bias in the 1/V_max calculation of the HIMF is that due to the large-scale structure of the galaxy distribution. Blind H i surveys tend to be relatively shallow and are thus biased by the overdensity of the local volume, which particularly affects the lowest-mass H i-rich galaxies that are only found nearby. If a correction for large-scale structure is not applied, we overestimate the impact of low-mass galaxies on the overall HIMF, therefore boosting the faint-end slope α artificially. We discuss this correction in Appendix A. The large-scale structure volume correction has only a very weak effect on log(M_*), but the effects on α (∼10%) and ϕ_* (∼30%) are large. Since this correction is so significant, it is sensitive to the details of the density reconstruction used. Agreement between the 1/V_max and 2DSWML results provides the best indication of the quality of this correction.

However, large-scale structure introduces the further bias of selectively reducing counts in mass bins that are primarily detectable in void volumes, and the weighting scheme correction cannot account for that. The voids in the Pisces–Perseus region between 3000 and 8000 km s⁻¹, visible in Figure 2, in particular, bias that portion of the α.40 sample against galaxies with $8.5 < \log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) < 9.0$, leading to a systematic undercounting in those bins. Because the 1/V_max method is sensitive to large-scale structure, this undercounting introduces a spurious "bump" feature into the HIMF, described in detail in Section 4.1.

As discussed in Sections 3.1 and 3.3, the main disadvantage of the 1/V_max method is its potential sensitivity to large-scale structure. If large-scale structure corrections were not adopted, the density of low H i-mass galaxies would be systematically overestimated, since most of these galaxies are detectable only in the very local, substantially overdense universe, including the Virgo Cluster and the Local Supercluster. This would bias the low-mass slope of the Schechter fit to the HIMF (α), weakening one of the major strengths of the ALFALFA data set, which is its ability to probe the population of extremely low H i-mass galaxies over a wide solid angle for the first time.

The original SWML method is applicable to samples selected by integrated flux. It assumes that the observed galaxy sample is drawn from a common H i mass function throughout the survey volume, denoted by $\phi (M_{\rm H\,{\mathsc{i}}})$. Unlike most maximum likelihood methods, which assume a functional form for $\phi (M_{\rm H\,{\mathsc{i}}})$ (Sandage et al. 1979), SWML splits the distribution in bins of $m = \log (M_{\rm H\,{\mathsc{i}}}/M_\odot)$ and assumes a constant distribution within each logarithmic bin. In this way, the value of the distribution in each of the bins becomes a parameter, ϕ_j (j = 1, 2, ..., N_m), which is adjusted in order to maximize the joint likelihood of detecting all galaxies in the sample, hence yielding a maximum likelihood estimate of the mass distribution. Since the values of the parameters are free to vary independently, the procedure above is completely general and does not assume any functional form for the distribution a priori.

In the case where the sample is not integrated flux limited and the selection function depends on additional observables, the SWML technique has to be extended to take into account the underlying galaxy distribution in all the physical properties that enter the calculation of the selection function. In the case of α.40, the limiting integrated flux depends on the galaxy profile width W₅₀ and thus the method needs to consider the joint two-dimensional distribution of galaxies in both H i mass and observed velocity width, $\phi (M_{\rm H\,{\mathsc{i}}},W_{50})$. 2DSWML relies on the assumption that the sample is statistically complete. Since ALFALFA's sensitivity to a source is dependent on both its integrated flux and its velocity width W₅₀ (Section 2.3), we fit a strict completeness threshold to the observed relationship as seen in Figure 1 and exclude galaxies falling below this completeness cut.

The details of the 2DSWML method and its application to α.40 are given in Appendix B.

The simplest source of error in the estimate of the HIMF is from Poisson counting errors in the bins, which is added to the other sources of error considered next. The relationship between errors on corrections applied to individual galaxies and errors on the final HIMF points and measured parameters is complex. Mass errors, for example, may shift galaxies in the sample from one bin to another as discussed in Section 3.2, so it is not possible to analytically calculate the error on a particular bin. In order to treat these errors appropriately, we create >250 realizations of the HIMF for each of the results shown in Sections 4 and 5. The error on S_int is measured in the ALFALFA source extraction pipeline, and we have estimated errors on the distance for each galaxy in the sample. Each of these contributes to the mass error, and we apply a Gaussian random error to each galaxy's mass in each realization. The spread in the bin values across the ensemble of realizations contributes to the overall error in each point. We consider errors due to uncertain parameter estimation in the relationship between $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot })$ and the Gumbel distribution parameters μ and β in the same way. This results in an HIMF that has taken known sources of error into consideration.

Sources of systematic bias remain, particularly for the 1/V_max measurement which is sensitive to the large-scale structure in the galaxy distribution. The effects of large-scale structure and of cosmic variance will be reduced as the survey continues, increasing its volume and coverage of varied cosmological environments.

In order to account for errors that are more difficult to quantify, we follow the example of Zwaan et al. (2005) and jackknife resample 21 equal-area regions. The resampling technique will help account for residual large-scale structure beyond that which we have corrected, and also for any systematic survey effects that change spatially across the sky or temporally throughout the survey's observations.

3.5.1. 2DSWML Error Estimates

The 2DSWML approach introduces another source of error. We assign errors on the parameters ϕ_jk, introduced in Section 3.4, via the inverse of the information matrix following Loveday (2000) and Efstathiou et al. (1988). The general form of the information matrix for a likelihood function $\newmathcal {L}$ that depends on a set of parameters $\btheta$ is given by

$$\begin{eqnarray} \mathbf {I}(\btheta) &=& - \left[\begin{array}{@{}cc@{}} \frac{\partial ^2}{\partial \theta _m \partial \theta _n}\ln \newmathcal {L}+\frac{\partial }{\partial \theta _m} g \frac{\partial }{\partial \theta _n}g & \frac{\partial }{\partial \theta _n}g \\ \frac{\partial }{\partial \theta _m}g& 0\end{array} \right], \end{eqnarray} \tag{ 2 }$$

where g is a constraint of the form $g(\btheta)=0$. We choose to apply the constraint $g=\sum _j \sum _k \: (\frac{M_{{\rm H\,{\mathsc{i}}},j}}{M_{\rm H\,{\mathsc{i}},ref}})^{\beta _1}(\frac{W_{50,k}}{W_{50,\rm ref}})^{\beta _2} \: \phi _{jk} \: \Delta m \Delta w-1=0$, with β₁ = β₂ = 1 and reference values for the H i mass and W₅₀ equal to the α.40 sample mean. The result is an error estimate for the parameters ϕ_jk, i.e., the value of the HIMF in each mass bin, and is added in quadrature to the other sources of error described above.

The global H i mass function derived from the α.40 sample via the 1/V_max method is presented in the top panel of Figure 7. Overplotted error bars have been derived as described above; mass errors due to errors on flux and distance estimates are reflected in the errors on the HIMF points, rather than on the mass-axis bin positions, since these errors change the bin counts.

Zoom In Zoom Out Reset image size

The best-fit Schechter function describing this HIMF is overplotted as a dashed line. The derived parameters are ϕ_* (h³₇₀ Mpc⁻³ dex⁻¹) = 6.0 ±.3 × 10⁻³, log (M_*/M_☉) + 2 log  h₇₀ = 9.91 ± 0.01, and α = −1.25 ± 0.02. However, the large-scale structure in the ALFALFA survey regions has introduced a "bump" into this measurement of the HIMF. The feature visible in Figure 7 at $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) \sim 9.0$ does not appear to be intrinsic to the H i-rich galaxy population. However, the previous work on luminous galaxies has suggested that the shape of luminosity and mass functions may be more complex than single Schechter functions. Luminosity functions in clusters, such as Coma and Fornax, are inconsistent with single values of α; Trentham (1998) has recommended a "composite" luminosity function that steepens for both bright and faint objects and flattens out in between, which provides a "dip" feature. Single Schechter functions provided a poor fit to 2dFGRS luminosity functions (Madgwick et al. 2002), and results from the Sloan Digital Sky Survey also suggest that a second (Baldry et al. 2004) or third (Li & White 2009) Schechter function component best describes the underlying population of galaxies at low redshift. While it is possible that the feature in Figure 7 suggests a complex shape in the HIMF given these findings, it is more likely that the feature is spurious, as we discuss below.

Such features occur because the 1/V_max method is sensitive to large-scale structure. Because the survey's H i mass sensitivity varies with distance (i.e., α.40 is not a volume-limited sample), each mass bin in the HIMF corresponds to some preferred distance at which ALFALFA is most sensitive to galaxies in that mass bin. Extended large-scale structures can therefore change the shape of the HIMF in bins corresponding to the distance of those features. Because of the large sample size of α.40, it is possible to investigate separately the three survey regions represented by the cone diagrams in Figures 2 and 3 and to isolate the structures that contribute to such features. Specifically, the "bump" feature in Figure 7 is due to a lack of sources in the foreground of the Great Wall and an overabundance within the Great Wall, clearly evident in Figure 3. The large-scale structure correction (Section 3.3 and Appendix A) reduces this feature, but cannot totally eliminate it, in part because density maps used to correct for large-scale structure are smoothed to ∼a few Mpc scales and can underestimate extremes in the density contrast. Features such as this one will be reduced as the ALFALFA survey continues and the sample grows. The 2DSWML method is not sensitive to large-scale structure and does not produce this feature (Section 5).

This feature appears significant in part because our statistical errors on the HIMF points are so small, but it leads to a poor fit and an underestimate for the faint-end slope α. This is clear in Figure 8, which displays the residual between the 1/V_max HIMF points and the derived best-fit Schechter function in the top panel and shows that the Schechter function systematically overestimates and underestimates the H i mass function due to this feature.

Zoom In Zoom Out Reset image size

While this feature is well understood, it has the undesirable effect of artificially reducing the faint-end slope α. In an effort to reduce the effect of this spurious feature and to better fit the points, we fit the sum of a Schechter function and a Gaussian; the Gaussian component serves to filter out the feature, leading to a better estimate of α. The results are shown as the solid line in Figure 7 with the residuals shown in the bottom panel of Figure 8. This fit significantly improves the reduced χ², and the residuals are small and, near $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) \sim 9.0$, more randomly scattered about 0 in contrast to the top panel of Figure 8. However, there is larger uncertainty in the parameters in this case, since each function is constrained by fewer points. The Schechter function parameters, displayed in Table 1, are log (M_*/M_☉) + 2 log  h₇₀ = 9.95 ± 0.04 and α = −1.33 ± 0.03. The Schechter function measurement of ϕ_* (h³₇₀ Mpc⁻³ dex⁻¹) = 3.7 ±.6 × 10⁻³, however, has been affected by the addition of the second component to the fit, and we therefore defer to the 2DSWML measurement of that parameter.

Table 1. H i Mass Function Fit Parameters

Sample and	α	ϕ*	log  (M*/M☉)	$\Omega _{\rm H\,{\mathsc{i}}}$, Fit	$\Omega _{\rm H\,{\mathsc{i}}}$, Points
Fitting Function		(10−3 h370 Mpc−3 dex−1)	+ 2 log  h70	(× 10−4 h−170)	(× 10−4 h−170)
1/Vmax	−1.25 (0.02)	6.0 (0.3)	9.91 (0.01)	4.4 (0.2)	4.4 (0.1)
Schechter + Gaussiana	−1.33 (0.03)	3.7 (0.6)b	9.95 (0.04)
1/Vmax, Non-Virgo	−1.20 (0.02)	6.1 (0.3)	9.90 (0.01)	4.1 (0.2)	4.3 (0.1)
Schechter + Gaussiana	−1.33 (0.04)	3.1 (0.6)b	9.95 (0.05)
2DSWML	−1.33 (0.02)	4.8 (0.3)	9.96 (0.02)	4.3 (0.3)	4.4 (0.1)
2DSWML, Non-Virgo	−1.34 (0.02)	4.7 (0.3)	9.96 (0.01)	4.3 (0.3)	4.4 (0.1)
1/Vmax, α.401.8	−1.30 (0.03)	4.6 (0.3)	9.96 (0.02)	4.0 (0.3)	4.0 (0.1)
1DSWML, α.401.8	−1.36 (0.06)	4.5 (0.9)	9.96 (0.04)	4.4 (0.9)	4.3 (0.3)
HIPASS (Zwaan et al. 2005)c	−1.37 (0.06)	5 (1)	9.86 (0.04)	3.7 (0.5)
Leo Group (Stierwalt et al. 2009)d	−1.41 (0.2)

Notes. ^aIn the 1/V_max case, pure Schechter functions provide a poor fit to the faint-end slope α, which explains the difference in α for two fitting functions. The Gaussian component parameters are not shown in the table, given that they are not expected to be physical. ^bWe defer to the 2DSWML measurement of ϕ_*, due to the spurious feature in the 1/V_max results. ^cReported statistical and systematic errors combined in quadrature. ^dThe excluded parameters ϕ_* and M_* in the Leo group are highly uncertain due to the lack of high-mass galaxies in its small volume.

Download table as:  ASCIITypeset image

The Gaussian parameters are not included in Table 1, since they are used to filter out the "bump" feature and are not expected to have physical meaning. The best-fit Gaussian has peak height (h₇₀ Mpc⁻³) 5 ± 1 × 10⁻³, mean log (M_μ/M_☉) + 2 log  h₇₀ 9.28 ± 0.06, and spread in log (M_μ/M_☉) + 2 log  σ = 0.41 ± 0.03.

We conclude that the proper values of α and log (M_*/M_☉) extracted from the 1/V_max method are −1.33 ± 0.03 and 9.95 ± 0.04, respectively. Table 1 lists both the spurious 1/V_max Schechter function parameters as well as the parameters found when a Gaussian is added to fit the spurious feature. The addition of the Gaussian brings the 1/V_max results for the parameters α and M_* into excellent agreement with the 2DSWML method and the flux-limited α.40_1.8 subsample results.

As an additional test of our corrections for profile width sensitivity, we have derived the 1/V_max HIMF from the integrated flux-limited subsample α.40_1.8 (described in Section 2.4). This mass function is corrected for large-scale structure and includes mass errors, but is not subject to the same bias against broad H i profiles. The α.40_1.8 HIMF is well fit by a pure Schechter function. The results are listed in Table 1. The α.40_1.8 HIMF is consistent with those derived from the full α.40 sample. We therefore conclude that our survey sensitivity is well characterized and that our measurements based on the full sample are complete and representative. However, since this limited sample does not probe the galaxies at the extremes of the mass function, it is subject to larger errors on the points and in the parameters.

4.1.1. Measurement of $\Omega _{\rm H\,{\mathsc{i}}}$

The density $\Omega _{\rm H\,{\mathsc{i}}}$ of neutral hydrogen in the local universe, expressed in units of the critical density, can be calculated in two ways from the derived H i mass function. Integrating analytically over the best-fit Schechter function gives $\Omega _{\rm H\,{\mathsc{i}}} = \phi _{*} \, M_{*} \, \Gamma (2+\alpha)$ = (4.4 ± 0.3) ×10⁻⁴ h⁻¹₇₀, slightly (16%) higher than the final HIPASS value 3.7 ×10⁻⁴ h⁻¹₇₀ (Zwaan et al. 2005). Using the binned points directly, we find the same result: $\Omega _{\rm H\,{\mathsc{i}}}=$ (4.4 ± 0.1) ×10⁻⁴ h⁻¹₇₀. This agreement is an indication that our findings are well represented in the high-mass bins by our Schechter function fit, despite the spurious feature. $\Omega _{\rm H\,{\mathsc{i}}}$ carries a small error since it is negligibly affected by the mass and distance errors on the faint end.

In Figure 9, we show the contribution of each 1/V_max mass bin to $\Omega _{\rm H\,{\mathsc{i}}}$ as filled circles. The total density of neutral hydrogen in the local universe is dominated by galaxies with $9.0 < \log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) < 10.0$, and in these bins we measure the HIMF to be larger than Zwaan et al. (2005) do, thus finding a larger value of $\Omega _{\rm H\,{\mathsc{i}}}$. The ALFALFA survey extends further in redshift than HIPASS, with a median redshift ∼8000 km s⁻¹ compared to ∼ 3000 km s⁻¹, allowing us to detect significantly more high-mass objects (Section 6.2).

Zoom In Zoom Out Reset image size

The HIMF derived from α.40 through the 2DSWML method is shown in Figure 10. The derived parameters are ϕ_* (h³₇₀ Mpc⁻³ dex⁻¹) = (4.8 ± 0.3) × 10⁻³, log (M_*/M_☉) + 2 log  h₇₀ = 9.96 ± 0.02 and α = −1.33 ± 0.02. To test the robustness of this HIMF estimate, we also applied a one-dimensional SWML approach to the flux-limited α.40_1.8 sample and found results consistent with the global, two-dimensional result(ϕ_* = (4.5 ± 0.9) × 10⁻³, log(M_*) = 9.96 ± 0.04 and α = −1.36 ± 0.06).

Zoom In Zoom Out Reset image size

5.1.1. Measurement of $\Omega _{\rm H\,{\mathsc{i}}}$

Measurements of the H i mass function can be sensitive to large-scale structure in the survey volume. As discussed above, we correct for large-scale structure in the 1/V_max method to ameliorate this effect, but our 2DSWML measurement could also be sensitive to this large nearby overdensity. To test the robustness of the 1/V_max correction and of our derived HIMF, we consider the result obtained when we exclude the portion of α.40 that crosses the Virgo Cluster. Many of our low-mass objects are contributed by this nearby overdensity, and our large-scale structure correction mechanism is the largest in this region; if we are correcting appropriately, we should obtain the same result regardless of the inclusion of the Virgo sources. This test is imperfect, given that the local volume generally is overdense. We exclude all galaxies lying within our adopted Virgo field, covering 12^h <α< 13^h and the full declination extent of the α.40 survey (Trentham & Hodgkin 2002), reducing the sample size to ∼ 9200 for 1/V_max and ∼ 8600 for 2DSWML. Errors are measured as described above, but in this case we jackknife resample over only 18 subregions.

Our results, within the errors, are the same whether or not we exclude the Virgo overdensity. This is true both for parameters and for our measurement of $\Omega _{\rm H\,{\mathsc{i}}}$. In the case of 1/V_max, we again find that a Schechter summed with a Gaussian provides a better fit to the data by accounting for features introduced by large-scale structure in the foreground of the Pisces–Perseus supercluster. In Table 1, we compare our findings for samples inclusive and exclusive of Virgo. Additionally, we list the HIPASS H i mass function and the Stierwalt et al. (2009) HIMF of ALFALFA sources in the Leo group. In the case of the α.40 and α.40_1.8 samples, we also list the value of $\Omega _{\rm H\,{\mathsc{i}}}$ found by integrating the Schechter function fit or using the HIMF bin points. Each table entry is accompanied by 1σ errors in parentheses.

We find a value of $\Omega _{\rm H\,{\mathsc{i}}}$ that is 16% higher than the complete HIPASS survey value (Zwaan et al. 2005). That HIPASS result is excluded by our 2σ errors, but the more preliminary HIPASS result (Zwaan et al. 2003) is in agreement with our result while carrying significantly larger error than we find. We also find log  (M_*/M_☉) = 9.96, so that the break in our HIMF occurs at masses 0.1 dex higher than was found in either of the HIPASS analyses. Since the high-mass end of the HIMF is sensitive to M_*, HIPASS significantly undercounts the highest-mass gas-rich galaxies. When our Schechter function is extrapolated to log (M_*/M_☉) = 11.0, we predict an order of magnitude more galaxies than HIPASS. At more modest values, log  (M_*/M_☉) = 10.75, this is reduced to a factor of ∼ 5.

In Figure 12, we show the mass of α.40 detections as a function of their distance in Mpc and compare that to the HIPASS completeness and detection limits. The dashed vertical line shows the 12,700 km s⁻¹ redshift cutoff of HIPASS assuming H₀ = 70 km s⁻¹ Mpc⁻¹, demonstrating the ALFALFA survey's ability to probe the rare highest-mass galaxies at large redshifts. While the α.40 sample extends only to 15,000 km s⁻¹ in order to avoid radio frequency interference, the full ALFALFA bandwidth allows us to probe to 18,000 km s⁻¹. Given that the survey was designed to be sensitive at those greater redshifts, we are still able to observe many galaxies at the limit of α.40, while the Zwaan et al. (2005) sample becomes very sparse near the survey's redshift limits.

Zoom In Zoom Out Reset image size

This improved measurement of the HIMF has implications for work that relied upon the HIPASS results. Present-day H i surveys are limited in their ability to probe redshift space, even when they are targeted (z < 0.5), so models of evolution of the H i mass function rely on the measurement at z = 0. Higher-precision measurements provide better constraints for evolutionary models. Numerical models of galaxy formation and evolution (Power et al. 2010) depend on the z = 0 HIMF to assess the success of the models and to extrapolate that result to predictions for future H i surveys. For example, Abdalla et al. (2010) predicted the ability of future H i line surveys with an instrument like the Square Kilometer Array (SKA) to constrain dark energy through measurements of the baryon acoustic oscillation scale. Those authors consider models of the HIMF evolution that are sensitive to the value M*. Typically, these galaxy models also depend on the assumed H₂/H i ratio to convert simulated cold gas into atomic and molecular components (e.g., Baugh et al. 2004), so updated estimates of either $\Omega _{\rm H\,{\mathsc{i}}}$ or $\Omega _{\rm H_{2}}$ affect our ability to produce realistic models of gas-rich galaxies.

We confirm previous findings that $\Omega _{\rm H\,{\mathsc{i}}}$ at z = 0 is inconsistent with the value inferred from damped Lyman absorber (DLA) systems at z ∼ 2 and that significant evolution is required to reconcile measurements in the two epochs (Prochaska & Wolfe 2009; Prochaska et al. 2005; Noterdaeme et al. 2009; Rao et al. 2006), while providing a tighter constraint on the present-day energy density of cold gas.

Obreschkow et al. (2009, hereafter O09) used the Millennium Simulation catalog, the De Lucia & Blaizot (2007) virtual catalog of galaxies, and a physically motivated prescription to assign realistic gas (H i, He, and H₂) masses at a range of redshifts. While this catalog has a limited ability to realistically trace detailed galaxy evolution and limited mass resolution—down to about 10^8.0 M_☉ of neutral hydrogen, which is comparable to the particle size in the Millennium run (Springel et al. 2005)—it serves as the best currently available comparison of observed gas-rich disks with the underlying theory of dark matter halos.

6.3.1. Simulated H i Mass Function

O09 derive an H i mass function that is, in its gross properties, consistent with HIPASS (Zwaan et al. 2005), ignoring spurious features near the mass resolution limit of the simulation. The O09 gas masses are obtained by combining the cold particle masses from the Millennium Run with a model to split the cold gas into molecular hydrogen and atomic hydrogen and helium components. Figure 13 compares the O09 HIMF, including only galaxies with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) > 8.0$ and at redshift z = 0, with the 1/V_max and 2DSWML HIMFs derived in this work. $\Omega ^{\rm sim}_{\rm H\,{\mathsc{i}}}$ = 3.4 ×10⁻⁴ inferred from the O09 HIMF is in good agreement with this work and with HIPASS. While it is clear that the overall statistical distribution of the cold gas prescription generally recovers the overall density and the gross properties of the statistical distribution, the details of the O09 HIMF disagree with observations, particularly at the extreme low-mass end where the Millennium Run work suffers from poor resolution and inadequate merger histories.

Zoom In Zoom Out Reset image size

It is also worth noting that O09 report that they overpredict the number of high-mass sources in comparison to HIPASS and suggest that this may be due to opacity in observed disks at these masses. However, we find that they underpredict high-mass galaxies at z = 0, the opposite effect. This is likely due to the O09 analysis of the HIMF, which is not limited to the final galaxies evolved to z = 0; rather, their HIMF also includes galaxies at higher-redshift simulation snapshots which are presumably more gas-rich than their present-day counterparts. This would therefore overpredict the abundance of high-mass galaxies.

6.3.2. Faint-end Slope

As has been found in the previous work, the faint-end slope of the α.40 HIMF is significantly shallower than the Press–Schechter prediction of α ∼ −1.8 (Press & Schechter 1974). Potentially, this difference can be linked to baryon loss and the suppression of accretion via photoheating in the low-mass dark matter halos. Simulations suggest that dark matter halos with masses below ∼6.5 × 10⁹ h⁻¹ M_☉ result in baryon-poor galaxies in present-day voids and other environments (Hoeft et al. 2008, 2006; Hoeft & Gottloeber 2010). In principle, the discrepancy could be explained by an argument invoking the mass scale at which photoheating becomes important.

A fitting function has been proposed (Gnedin 2000) to describe the behavior of baryon fraction as a function of underlying halo mass:

$$\begin{equation} f_{b} = f_{b0} \left[ 1 + (2^{\gamma /3} - 1) \left(\frac{M_{c}}{M_{\rm tot}} \right)^{\gamma } \right]^{-3/\gamma }, \end{equation} \tag{ 3 }$$

where the parameters f_b0 and M_c are, respectively, the baryon fraction in large halos and the characteristic halo mass where f_b = f_b0/2.

If decreasing baryon fraction with decreasing halo mass explains the difference between low-mass slopes in baryonic (stellar and H i) and halo mass functions, then this fitting function should consistently predict baryonic and cold gas mass functions with values of α ∼ −1.3. In the low-mass limit, the first term of Equation (3) can be dropped and the total mass in a halo can be assumed to be dominated by the dark matter, M_tot ≈ M_D. Via the definition f_B = M_B/M_D, we have

$$\begin{equation} M_{D}=\frac{M_{B}}{f_{B}}=\frac{M_{B}}{f_{b0}} (2^{\gamma /3} - 1)^{3/\gamma }\left(\frac{M_{D}}{M_{c}} \right)^{-3}. \end{equation} \tag{ 4 }$$

Compressing all constants gives the relation M_D ∝ M^1/4_B, which can then be used to relate the low-mass ends of the baryonic and dark matter mass functions. On the faint end of the dark matter mass function, the exponential term of the Schechter function can be dropped. From $\frac{d \log M_{D}}{d \log M_{B}}=\frac{1}{4}$ we can, finally, conclude that

$$\begin{equation} \phi (M_{B}) = \frac{dn}{d \log M_{B}} \propto \ \phi _{*} \ \left(\frac{M_{B}}{M_{*,B}} \right)^{(\alpha _{D} + 1)/4}, \end{equation} \tag{ 5 }$$

where (α_D + 1)/4 = α_B + 1. Starting from the Press–Schechter prediction of a faint-end slope α_D ∼ −1.8, the consideration of baryon fraction leads to α_B ∼ −1.2, which is more consistent with H i and stellar mass functions. In principle, the discrepancy between dark matter simulations and observed baryon mass functions could be explained by the photoheating simulations of Hoeft et al. (2006) and Hoeft & Gottloeber (2010).

The baryon fraction of O09's simulated galaxies loosely follows this descriptive baryon fraction function (Equation (3)). However, the halo mass scale at which the baryon loss starts to drop steeply is about 2 orders of magnitude larger than the scale found by the detailed hydrodynamical simulations of Hoeft et al. (2008). Additionally, there is large scatter in the mass interval of interest, since the simulation's resolution is poor for the halo masses where baryon loss becomes important. The level of agreement between O09 and Hoeft et al. (2008) is therefore difficult to quantify, and we use the latter's determination of f_b in what follows.

Equation (3) suggests that the baryonic content of low-mass galaxies in α.40 may be severely biased with respect to the underlying halo mass distribution. If simulations accurately predict the relationship between initial halo masses and resulting baryon fractions after reionization and photoheating, then the application of f_b should provide an estimate of the resulting baryon mass function at z = 0. This depends on the extremely naive assumption that the cold H i gas content is depleted in the same fraction as the baryons overall.

The publicly available GENMF code⁷ produces halo mass function fits to the Reed et al. (2007) N-body simulations at high resolution, from 10⁵ to 10¹² h⁻¹ M_☉. We adopt their mass function at z = 0, with their suggested parameters Ω_M≈ 0.238, Ω_Λ≈ 0.762, and σ₈ = 0.74 (at z = 0), and apply Equation (3) to extract the predicted baryon mass function and fit the faint-end slope. The results are displayed in Table 2 for an exemplary set of values for f_b,0, M_c, and γ.

Table 2. Faint-end Slopes of Modeled Baryon Mass Functions

fb,0	Mc	γ	α
0.20	9.0	1.0	−1.30
0.20	9.5	1.0	−1.27
0.16	9.0	1.0	−1.31
0.16	9.5	1.0	−1.28
0.16	9.0	1.5	−1.24
0.16	9.5	1.5	−1.22
0.16	9.0	2.0	−1.21
0.16	9.5	2.0	−1.19
0.15	9.0	1.0	−1.31
0.15	9.5	1.0	−1.28
0.15	9.0	2.0	−1.21
0.15	9.5	2.0	−1.19

Download table as:  ASCIITypeset image

Through this approach, it is possible to modify the underlying halo mass function (α≈ −1.8) to meet our observations (α≈ −1.3). The suggestion that low-mass halos may re-accrete cold gas at late times (Ricotti 2009), if substantiated, could further change the shape of the resulting baryon mass function. While this approach indicates that we may be close to resolving the missing satellites problem and the discrepancy between predicted and observed faint-end mass function slopes, the precise requirements of baryon depletion mechanisms are not well constrained by available simulations.

Giovanelli et al. (2005) predicted, from the precursor survey observations, that ALFALFA in full two-drift mode could expect an approximate integrated flux detection threshold, S_int,th in Jy km s⁻¹, dependent upon profile width as follows:

$$\begin{equation} S_{\rm int, th} = \left\lbrace \begin{array}{@{}ll} 0.15 \, {\rm S/N} \, (W_{50}/200)^{1/2}, & W_{50} < 200 \\ \ms 0.15 \, {\rm S/N} \, (W_{50}/200), & W_{50} \ge 200. \end{array} \right. \end{equation} \tag{ A1 }$$

In practice, however, ALFALFA outperforms this detection threshold, and we therefore use the data itself to fit a detection limit as described in Section 3.3.

The width-dependent sensitivity correction is based on the distribution of observed profile widths. We also assume that the distribution of observed galaxies gives an indication of the true underlying distribution. We are therefore interested in working with as many sample galaxies as possible, and thus we consider a detection threshold S_int,th as a function of W₅₀ that indicates the limits of ALFALFA's detection ability, rather than a strict completeness limit as in the 2DSWML case (Section 3.4).

The completeness correction is based on the relationship of galaxy mass to the distribution of profile widths W₅₀. It is known that H i profile widths and masses are correlated, and we observe a mass-dependent spread in the distribution of profile width. We determine the profile width distribution as a function of mass by binning α.40 galaxies by $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot })$ and fitting to each histogram a Gumbel (or Extreme Value Type 1) distribution:

$$\begin{equation} f(x) = \frac{1}{\beta } \, e^{\frac{x-\mu }{\beta }} \, e^{-e^{\frac{x-\mu }{\beta }}}, \end{equation} \tag{ A2 }$$

where μ parameterizes the center of the distribution and β its breadth. The profile width distributions feature narrow central peaks and extended skewed tails, which the Gumbel distribution is designed specifically to model.

We find that the center of the profile width distribution increases linearly with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot })$, and the breadth decreases linearly with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot })$. We derive a relationship between $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot })$ and the parameters μ and β in order to extrapolate to any mass and infer the underlying distribution of W₅₀ to which a given galaxy belongs, $P(W_{50}, M_{\rm H\,{\mathsc{i}}})$. The probability of detecting a galaxy in a given mass bin depends on the profile width distribution for that bin, as well as the limiting profile width W$_{50,\rm lim}$ beyond which that galaxy would not be detectable by ALFALFA. We are seeking a correction factor C that will account for the profile width-integrated flux bias and that satisfies the relationship

$$\begin{equation} N_{\rm galaxies} (M_{\rm H\,{\mathsc{i}}}) = C \, N_{\rm obs}(M_{\rm H\,{\mathsc{i}}}), \end{equation} \tag{ A3 }$$

where N_galaxies is the corrected galaxy count to be input for the calculation of the HIMF and N_obs is the observed galaxy count. In terms of the derived distribution $P(W_{50}, M_{\rm H\,{\mathsc{i}}})$, we have

$$\begin{equation} C = \frac{\int _{-\infty }^{+\inf } P (W_{50}, M_{\rm H\,{\mathsc{i}}}) \, dW_{50}}{\int _{-\infty }^{W_{50,{\rm lim}}} P (W_{50}, M_{\rm H\,{\mathsc{i}}}) \, dW_{50}}. \end{equation} \tag{ A4 }$$

Since a bin is made up of galaxies with varying W$_{50,\rm lim}$, we apply this correction to each individual galaxy, rather than on a mass bin-by-bin basis. The sum over effective search volume, Σ1/V_max, therefore becomes ΣC/V_max.

To be conservative, we have included the errors on our derived linear relationships between $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot })$ and the Gumbel distribution parameters μ and β in our final error analysis for the H i mass function.

The 1/V_max method would be biased by large-scale structure if we counted galaxies in overdense regions with the same weight as their counterparts in voids. Instead, we want to consider the effective search volume V_max,eff in such a way that overdense regions are counted as contributing more effective volume to the overall survey.

We modify Σ1/V_max to include weighting by the average density n(V_max) interior to D_max, normalized to the average density of the universe. The expression for measuring the HIMF then becomes Σ1/n(V_max)V_max (Springob et al. 2005). We obtain n(V_max) from the PSCz density reconstruction of Branchini et al. (1999), using their Cartesian map of evenly spaced grid points out to 240 Mpc h⁻¹ smoothed to 3.2 Mpc h⁻¹ and using our assumed value h = 0.7. For values D_max> ∼ 85 Mpc, the average density interior to D_max becomes equal to the average density in the PSCz map, so no correction is needed. The large-scale structure correction is therefore small compared to the Poisson counting error for galaxies with $\log (M_{\rm H\,{\mathsc{i}}}/M_{\odot }) > 9.0$, which are found at large distances.

This weighting scheme for galaxy counts in overabundant and underabundant regions corrects the relative counts between different environments, so that clusters and superclusters do not dominate the shape of the measured HIMF.

To transform the calculated probability density function into an H i mass function (e.g., transform the unitless {ϕ_kΔm} into space densities), we evaluate the amplitude of the HIMF by matching the integral of the distribution to the inferred average density of galaxies in the survey volume $\bar{n}$, as in Zwaan et al. (2003). Davis & Huchra (1982) discuss various estimators for $\bar{n}$ that strike different balances between stability against poor knowledge of the selection function of the survey and immunity to large-scale structure. Since we believe we have a good understanding of the selection function out to cz= 15,000 km s⁻¹, we choose to adopt the estimator that is least prone to bias, denoted by n₁, defined as

$$\begin{equation} n_1 = V_{\rm survey}^{-1} \int \frac{n(D) \: dD}{S(D)}, \end{equation} \tag{ B8 }$$

where n(D) dD is the number of galaxies in a spherical shell of thickness dD and radius D, and V_survey is the total survey volume. The selection function S(D) is the fraction of galaxies detectable at distance D and is given by

$$\begin{equation} S(D) = \frac{\int _{w_{\rm min}}^{w_{\rm max}} \int _{m_{\rm lim}(w,D)}^{m_{\rm max}} \: \phi (m,w) \: dm \ dw}{\int \int \: \phi (m,w) \: dm \ dw}. \end{equation} \tag{ B9 }$$

In the case of the 2DSWML method, we evaluate n₁ by the expression

$$\begin{equation} n_1 = V_{\rm survey}^{-1} \sum _i \frac{1}{\sum _j \sum _k H_{ijk} \phi _{jk} \Delta m \Delta w}. \end{equation} \tag{ B10 }$$

Equation (B10) corresponds to weighing each detected galaxy in the survey by the inverse of the selection function at the galaxy's distance, effectively correcting each detection by the fraction of galaxies that cannot be detected at distance D_i.