THE CLUSTERING CHARACTERISTICS OF H i-SELECTED GALAXIES FROM THE 40% ALFALFA SURVEY

The ongoing ALFALFA survey is completing a census of galaxies in the local universe, out to z ∼ 0.06, using the seven-pixel ALFA receiver at the Arecibo Observatory to detect the 21 cm line of neutral hydrogen. Compared to previous blind neutral hydrogen surveys (e.g., HIPASS), ALFALFA's enhanced sensitivity, detection centroiding, volume, and sample size, resulting in a cosmologically representative sample, make it ideally suited for an accurate measurement of the correlation function of gas-rich galaxies.

The sample used here includes the sky coverage of the α.40 sample recently presented in Haynes et al. (2011), referred to as α.40 because it includes the data extracted from coverage of 40% of ALFALFA's skyprint. The statistical completeness and noise characteristics of the ALFALFA source catalog are well understood and have been discussed extensively elsewhere. Further details may be found in Saintonge (2007), Martin et al. (2010), and Haynes et al. (2011), which include discussions of the characteristics of the α.40 sample and the sensitivity of the ALFALFA survey. In particular, Haynes et al. (2011) discuss impacts of various volume restrictions on the derivation of the H i mass function. Here, we summarize the salient points.

Confidently detected sources in ALFALFA are assigned one of three object codes, where Code 1 refers to a reliable extragalactic detection with a high signal-to-noise ratio (S/N; > 6.5). For the sample used here, we neglect the other objects, Code 2 and Code 9 sources; Code 9 sources are high velocity clouds (HVCs) of hydrogen gas in the vicinity of the Milky Way and are thus not extragalactic, whereas Code 2 extragalactic sources have lower S/N and are only included in the catalog because they are corroborated by a known optical source at the same position and redshift. Furthermore, ALFALFA's ability to detect extragalactic signal near its redshift limit is degraded due to a strong source of terrestrial radio frequency interference, the FAA radar at the San Juan airport. We therefore include only objects within 15,000 km s⁻¹, which results in only a modest loss of source counts.

ALFALFA's sensitivity depends not only on the integrated flux, but also on the 21 cm spectrum's profile width W₅₀ (km s⁻¹). Because the mass of an H i source is a function of its distance and integrated flux, integrated flux can be thought of as a proxy for mass. The survey therefore is not volume-, flux-, or mass-limited, and the reconstructed selection function must take this complex sensitivity into account. Thus, when the sample is viewed as the distribution of galaxy masses as a function of distance, as in Figure 3 of Haynes et al. (2011), it is clear that α.40 is sensitive to very low H i mass galaxies nearby but only to significant masses at greater distances

We refer the reader to Figure 1 of Martin et al. (2010), which displays the dependence of ALFALFA's sensitivity on both integrated flux and profile width. The HIPASS survey recovered sources with the same dependence on these two parameters. ALFALFA 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 5.6 Jy km s⁻¹  for the same source in HIPASS. In Martin et al. (2010), we fit a linear relationship between integrated flux and profile width, with a break at 400 km s⁻¹, which describes the sensitivity of the survey. We use that same relationship and the selection function derived in that work throughout this paper.

The trimmed sample includes only Code 1 objects from the α.40 catalog within 15,000 km s⁻¹, for a total of ∼10,150 galaxies used in measuring the two-point correlation function.

The distance dependence of the selection function of α.40 was determined using the two-dimensional stepwise maximum likelihood (2DSWML) method, described fully in Appendix B of Martin et al. (2010). 2DSWML is related to the stepwise maximum likelihood (SWML) approach, but modified to account for the survey sensitivity's two-dimensional dependence on both integrated flux and profile width. 2DSWML splits the distribution of galaxy masses and profile widths in the α.40 sample into logarithmic bins, and then calculates the best-fit H i mass function (analogous to a luminosity function) which maximizes the joint likelihood of detecting all galaxies in the sample. This approach simultaneously measures the selection function for each detected galaxy in the sample. In this work, we will use that selection function S(D_i) for each galaxy i with a known distance D_i in Mpc.

For application to the correlation function, $S_{D_i}$ is calculated for every galaxy in the sample and then S_D is smoothed. This smoothed selection function can additionally be combined with an H i mass function to make predictions regarding the number of galaxies of a given mass that are expected to be found in the survey, or used to predict the redshift distribution of the survey under an assumption of homogeneity. Figure 1, previously published in Haynes et al. (2011), shows a histogram of the α.40 redshift distribution, with peaks and dips representing clusters and voids, respectively, along with an overplotted prediction based on the selection function and a non-clustered universe.

Zoom In Zoom Out Reset image size

Disagreements between the prediction and the observations are due both to the existence of large-scale structure in the survey volume and to the loss of survey sensitivity at certain velocities due to radio frequency interference. This contamination is quantified as a percentage of survey coverage at a given heliocentric velocity, or spectral weighting, as discussed in Martin et al. (2010) and earlier publications from the ALFALFA survey. For the purposes described here, the weights map has been translated into the cosmic microwave background (CMB) reference frame in order to most accurately model the predicted ALFALFA galaxy distribution (see Section 3.2). The selection function is used both for the creation of the random samples for estimation of the correlation function, and for the weighting of pair counts in that estimate.

We adopt the Landy–Szalay pairwise estimator (Landy & Szalay 1993) for the correlation function. The Landy–Szalay normalization of data–data (DD), random–random (RR), and data–random (DR) pair counts allows us to construct a random catalog that contains many more objects than the observed data catalog, thereby reducing the introduction of shot noise from the random set. The Landy–Szalay estimator is constructed from these normalized counts:

$$\begin{equation} \hat{\xi }_{{\rm LS}} = \frac{D_{{\rm DD}} - 2(D_{{\rm DR}}) + D_{{\rm RR}}}{D_{{\rm RR}}}. \end{equation} \tag{ 3 }$$

Because α.40 is not volume-limited, the pair counts must be weighted so that the measurement is not dominated by galaxies at the peak of the selection function. Following Meyer et al. (2007) and Hawkins et al. (2003), we apply a weighting w_ij = w_i × w_j for the contribution of each pair i, j to the Landy–Szalay estimator, given by

$$\begin{equation} w_i = \frac{1.0}{1.0 \: + \: 4\pi \: N_D \: S(r_i) \: J_3(s)}, \end{equation} \tag{ 4 }$$

where N_D is the number of galaxies in α.40, S(r_i) is the selection function measured for α.40 at r_i = cz_CMB, i/H₀, and

$$\begin{equation} J_3(s) = \int ^{s}_{0} s^{\prime 2} \: \xi (s^{\prime }) \: ds^{\prime } \end{equation} \tag{ 5 }$$

defined in terms of the redshift–space coordinate $s = \sqrt{\sigma ^2 + \pi ^2}$.

This expression for J₃ requires an assumed model for ξ(s), but the final measurement of the correlation function is not sensitive to this assumed input for object weighting; we assume a power-law form

$$\begin{equation} \xi (s) = \left(\frac{s}{s_0} \right)^{-\gamma } \end{equation} \tag{ 6 }$$

and we test our robustness by first assuming a fiducial value found for optically selected samples, s₀ = 5.0 h⁻¹ Mpc and γ = 1.8 and, after that, iterating to the value s₀ and γ measured for α.40. No statistically significant difference is observed through this iterative process, and we therefore proceed as other authors have, using the fiducial optical values reported here in our J₃ weighting. Following Fisher et al. (1994a, 1994b), we apply an artificial cutoff with a maximum value of s = 30 h⁻¹ Mpc in the expression for J₃.

We construct random samples that contain 20 times the number of objects in the α.40 data set. These random samples are carefully designed to include survey selection effects while excluding correlations due to large-scale structure. This is accomplished by predicting the distribution of cz_CMB from the survey selection and H i mass functions (see Figure 1) and then folding in the loss of volume as a function of velocity due to radio frequency interference, measured from the spectral weights map in Martin et al. (2010). Objects in the random set are randomly assigned a sky position within the right ascension and declination boundaries of α.40 and are then assigned a redshift from this predicted distribution. The resulting redshift distribution for one example instance of the random sample procedure is shown in Figure 2.

Zoom In Zoom Out Reset image size

The correlation function is measured in bins of separation. While the correlation function is expressed as a function of several different coordinates while iterating toward the real-space correlation function ξ(r), the bin counts and thus the measured correlation functions are correlated with one another in every such coordinate system. Because structures, such as clusters, will contribute an overabundance of pairs to a set of several bins, the measurement in each bin is not independent of the others. In plots of the correlation function shown here, we display the on-diagonal elements of the covariance matrix (i.e., the standard deviations) as uncertainties on each point. However, in order to work with our measurement to estimate the power-law shape of the correlation function of gas-rich galaxies, we must construct a full covariance matrix and take off-diagonal elements into account.

To construct the covariance matrix C, we carry out our pair-counting routine on more than 500 bootstrap resamplings of the data, and a single catalog of random objects is reused in each case. Each of the bootstrap measurements of ξ(σ, π) contains N_g galaxies selected at random from α.40, with replacement. From this set of realizations, we construct the covariance matrices for ξ(s), Ξ(σ) and ξ(r). The covariance between two correlation function bins b_l and b_m is given by

$$\begin{equation} C(l, \: m) = \sum _{i=1}^{N_{{\rm realizations}}} \frac{(b_{l, \: i} - \bar{b_l})-(b_{m, \: i} - \bar{b_m})}{N_{{\rm realizations}} - 1}. \end{equation} \tag{ 7 }$$

The significant off-diagonal elements of the covariance matrix make it difficult to obtain a power-law model fit by minimizing the χ² values weighted by the variance. The covariance matrix, however, is not an inescapable quality of the data, but is actually dependent on the basis in which the data are projected. In this case, we have some number of bins N_b representing a set of variables b (bin centers in h⁻¹ Mpc), and can choose to work in an orthonormal basis with N_b coordinate axes in which the covariance matrix C is diagonalized. This basis is defined by the principal component eigenvectors of the measurement, and we borrow elements of principal component analysis in order to obtain model parameter fits and uncertainty estimates.

The principal components are linear combinations of the original N_b variables arranged such that the first principal component corresponds to an orthonormal axis through N_b-dimensional space that explains the largest proportion of variance in the data set. These principal component vectors are defined by the eigenvectors of the covariance matrix of the original data set.

Following Fisher et al. (1994a), we calculate the principal eigenvectors and construct a diagonalizing matrix, R, the columns of which are these eigenvectors, and a new covariance matrix, $\tilde{C}$, projected in the new basis set. Since all of the covariance has been accounted for in the definition of the principal components, $\tilde{C}$ has no off-diagonal elements, and the variance is captured in the on-diagonal elements $\tilde{\sigma }$.

Given $\tilde{C}$ and R, a set of models with varying values for s₀ and γ can be projected into the principal component basis via $\tilde{b}_{{\rm model}}$ = R^T b_model, for comparison to the measured $\tilde{b}$. We find the value of each parameter that minimizes the expression

$$\begin{equation} \chi ^2 = \frac{1}{N_b - 2} \sum ^{N_b}_{1} \frac{\tilde{b}_{l} -\tilde{b}_{l,{\rm model}}}{\tilde{\sigma }_{l}}. \end{equation} \tag{ 8 }$$

Finally, we construct error ellipses to fully describe the likely parameter space of the power-law model for the correlation function (Press et al. 1992).

From the three-dimensional galaxy coordinates available to us, we can construct ξ(σ, π). This calculation of ξ(σ, π) is the fundamental measurement upon which the results presented in the rest of this work are based. The resulting image, shown in Figure 3 with contours overplotted, clearly reveals the redshift–space distortions that lead to the difficulty in estimating the real-space correlation function. The radial coordinate, π, appears weakly stretched at small angular separation σ because of the Eddington effect in clusters, though because H i-selected galaxies are known to avoid dense cluster regions, this effect is less prominent for α.40 than for optically selected samples. In the other dimension, π is flattened on large scales because of the coherent motion of galaxies toward attractors. This "squeezing" effect is determined by the clustering bias (with respect to dark matter) of the sample and by the underlying matter density fluctuations. In a future work, we will explore the shape of ξ(σ, π) to model the matter density field.

Zoom In Zoom Out Reset image size

In order to obtain the real-space correlation function ξ(r), we must take the intermediate step of projecting ξ(σ, π) along the π axis (in practice, using the discrete bins of size Δπ), resulting in what is known as the "projected correlation function" and symbolized as Ξ(σ)/(σ):

$$\begin{equation} \frac{\Xi (\sigma)}{\sigma } = \frac{2}{\sigma } \sum _{0}^{\pi _{{\rm max}}} \xi (\sigma, \: \pi) \: \Delta \pi. \end{equation} \tag{ 9 }$$

Following previous work, the maximum value of the integration, π_max, is selected so that the summation is convergent but is kept as low as possible to avoid the introduction of noise from poorly measured intermediate scales. From the original estimation of ξ(σ, π), which counted pairs up to distances of ∼60 h⁻¹ Mpc, we carry the sum along the π axis up to a scale of π_max 29.7 h⁻¹ Mpc. Further, we confirm that the resulting correlation function is not strongly sensitive to the chosen value of π_max, but extending the integration to scales that are too large for the sample to sufficiently measure introduces scatter and noise into the correlation function estimate.

Ξ(σ)/(σ) is closely related to the function in which we are truly interested, ξ(r), where r is the real-space distance, via

$$\begin{equation} \frac{\Xi (\sigma)}{\sigma } = \frac{2}{\sigma } \int _{\sigma }^{\infty } \xi (r) \: \frac{r \: dr}{(r^2 \: - \sigma ^2)^{1/2}}. \end{equation} \tag{ 10 }$$

In order to evaluate the real-space correlation function, some assumptions must be made about its form. Two options are usually explored in the literature: a power-law form, or a stepwise-function form which makes no assumptions about shape but does assume that the binning used well represents an underlying smooth correlation function (i.e., the "direct inversion" method). If we assume a power law of the form ξ(r) = (r/r₀)^−γ, we find

$$\begin{equation} \frac{\Xi (\sigma)}{\sigma } = \left(\frac{r_0}{\sigma }\right)^{\gamma } \: \frac{\Gamma (1/2) \: \Gamma ((\gamma - 1)/2)}{\Gamma (\gamma /2)}. \end{equation} \tag{ 11 }$$

In Equation (11), the function Γ is the well-known Gamma function. Equation (11) can be recast in terms of fitting parameters:

$$\begin{equation} \frac{\Xi (\sigma)}{\sigma } = \left(\frac{r_0}{\sigma } \right)^{\gamma } \: A(\gamma). \end{equation} \tag{ 12 }$$

Following Meyer et al. (2007), we rearrange Equation (12), obtain the best-fit power law of the form ξ(σ) = $a_1 \; \sigma ^{a_2}$ using the χ² minimization given by Equation (8), and then relate those parameters to r₀ and γ which represent the best-fit power law for ξ(r).

In the next section, we derive and discuss ξ(r) using the mechanisms described in this section.

The projected correlation function Ξ(σ)/σ (recast for the figure and the fitting as Ξ(σ)) is displayed in Figure 4, along with error bars reflecting the on-diagonal elements of the full covariance matrix. The dashed line is the best-fit model obtained by χ² minimization using the full covariance matrix. In Table 1, we list the parameters for the fit and their uncertainties, along with the fits obtained if only the on-diagonal elements (the standard deviations, σ) are used to carry out the standard least-squares fit. For comparison, we also include the clustering reported by the HIPASS team (Meyer et al. 2007; note that those authors ignored the off-diagonal elements in their error analysis), the clustering found by Basilakos et al. (2007) using the same HIPASS data set, and the clustering of several optically selected samples of interest. We also display the Passmoor et al. (2011) results, which used a small, publicly available early subset of the ALFALFA data.

Zoom In Zoom Out Reset image size

Error ellipses are displayed in Figure 5, with the dashed contour giving the 1σ single-parameter uncertainties listed in Table 1.

Zoom In Zoom Out Reset image size

While both the full covariance analysis and the assumption of bin independence give similar results, the larger uncertainties on the full covariance analysis give an indication of the need to be conservative. Parameter uncertainties previously reported in the literature (i.e., Meyer et al. 2007) significantly underestimate their reported statistical uncertainties. Even with the greater sensitivity, larger sample size, and deeper redshift range of α.40, the correlation function analysis allows for quite a large range of clustering scenarios.

We confirm the HIPASS result that H i-selected galaxies are among the most weakly clustered known class of galaxies, most comparable to, but still less clustered than, the faint late-type subsampling in 2dFGRS and the IRAS galaxy redshift survey, which was also biased toward star-forming galaxies. Similar to the results we will present in Section 6 for α.40, the flux-limited sample of IRAS galaxies considered in Fisher et al. (1994a) was found to be antibiased relative to cold dark matter on small scales but unbiased on intermediate scales (∼10 h⁻¹ Mpc) and positively biased on the largest scales (beyond 10 h⁻¹ Mpc). Similarly, the QDOT sample (Moore et al. 1994), also taken from the IRAS parent catalog but based on a lower flux limit and employing a different sampling strategy, was found to be unbiased with respect to dark matter.

Guzzo et al. (1997) provide an interesting comparison to our findings, as it was also based on an analysis of 21 cm galaxy profiles observed with the Arecibo Observatory, although that work was based on an optically selected, magnitude-limited sample rather than a blindly H i-selected sample as in this case. This work also samples a region that partially overlaps with α.40. Guzzo et al. (1997) split their sample by morphological type and determined the variation of clustering strength between early- and late-type (spiral and irregular) galaxies in Zwicky's catalog within the Pisces-Perseus region. They found that the early types were significantly more clustered than the late types, as reflected more generally in Table 1, but their volume-limited sample is significantly more clustered than the H i-selected ALFALFA sample.

Our findings for the clustering of H i-selected galaxies are in agreement with previous results, particularly with our understanding that ALFALFA galaxies tend to be blue, spiral, and late-type galaxies which are already known to be weakly clustered (Norberg et al. 2002; Kauffmann et al. 2004; Zehavi et al. 2005). Apart from the estimates of uncertainties, the clustering of α.40 is in agreement with the HIPASS findings but not with Passmoor et al. (2011), which is not unexpected giving the weaknesses—in particular extremely small sample size—of the latter.

ξ(r) will only be tidily related to Ξ(σ)/σ if we assume that the underlying physics of galaxy formation dictates a power-law form for ξ(r). In the previous section, we calculated the correlation function of gas-rich galaxies under that assumption, but it is also possible to avoid that assumption and obtain ξ(r) by direct inversion of the projected correlation function Ξ. The inversion method tests the power-law assumption, though it is a noisy measurement and results in large scatter. As an independent test of the shape of the correlation function, ξ(r) determined via this inversion method will be especially useful at scales above ∼2.5 h⁻¹ Mpc, where the correlation function shows features inconsistent with the power-law assumption.

Following Meyer et al. (2007), Hawkins et al. (2003), and Saunders et al. (1992), we take our measurement of Ξ(σ) to represent an underlying step function form with values Ξ_l in intervals with centers σ_l, rearrange Equation (10), and interpolate between bins to give

$$\begin{eqnarray} &&\xi (r = \sigma _i)= {-} \frac{1}{\pi } \sum _{j \ge i} \frac{\Xi (\sigma _{j+1}) - \Xi (\sigma _j)}{\sigma _{j+1} - \sigma _j} \times \ln \frac{\sigma _{j+1} + \sqrt{\sigma _{j+1}^{2} {-} \sigma _{i}^{2}}}{\sigma _j + \sqrt{\sigma _{j}^{2} {-} \sigma _{i}^{2}}}. \nonumber \\ \end{eqnarray} \tag{ 13 }$$

The sum in Equation (13) is truncated so that σ_max = π_max for the value of π_max used in Equation (9).

The projected correlation function and ξ(r) obtained by the inversion method are in excellent agreement, as shown in Figure 6, where the points are the inversion ξ(r) with on-diagonal uncertainties and the overplotted dashed line is the best-fit power-law model for the projected correlation function. It is also clear that the inversion method, as remarked earlier, is more vulnerable to variance. Its use is motivated as a check on the assumed shape of the correlation function. Furthermore, because the power-law assumption is known to be useful only on small-to-intermediate scales, the correlation function obtained via the inversion method allows us to extend to large scales, ∼30 h⁻¹ Mpc.

Zoom In Zoom Out Reset image size

By relaxing the power-law assumption, we can further examine divergence from a power-law shape, a well-known phenomenon found in clustering studies of other populations of galaxies. The ALFALFA correlation function shows a "shoulder" at scales of ∼ a few h⁻¹ Mpc, as observed by, e.g., Guzzo et al. (1991), Hawkins et al. (2003), and Zehavi et al. (2004). Under the assumption of an inflationary, cold dark matter universe and using a halo occupation distribution (HOD) model, Zehavi et al. (2004) infer that the well-known shoulder is due to two distinct regimes in which galaxy–galaxy pairs are counted. On large scales, pairs are counted from separate dark matter halos, while on small scales, pairs are counted in the same dark matter halo and are subject to nonlinearity.

Because the galaxies probed by ALFALFA are gas rich, and because an H i-selected sample is biased toward gas-dominated low-mass objects that would be classified as low surface brightness (LSB) dwarfs in an optical survey, we expect that the characteristics of the single-halo regime would differ from that observed in the case of an optically selected sample. For example, tidal interactions and stripping within dense halos, which would decrease the pair counts of H i-selected galaxies, would change the relative contributions of the two regimes. Watson et al. (2011) find that a resulting power-law correlation function, when the contributions from both the one-halo and two-halo regimes are included, is only found under conditions in a narrow mass and redshift range for the general population of galaxies. Given a different HOD model as a function of galaxy properties, the shoulder or break from the power law would be more prominent. The ALFALFA correlation function therefore provides yet another approach from which we can better understand the evolution of H i and the distribution of gas-rich galaxies in the present universe. In a future paper, as the ALFALFA data set continues to grow, we will present an HOD analysis of the ALFALFA correlation function as an extension of the present work.

Our estimation of ξ(σ, π) included two choices upon which our results could be dependent. We selected both the logarithmic binning intervals for pair counting as well as the value π_max for projecting ξ(σ, π) into Ξ(σ). Investigations of alternative schemes suggest that our results are not strongly dependent on either the binning or the choice of π_max, though π_max is specifically selected to lead to a stable solution for ξ without introducing the noise and scatter from scales that are poorly probed by the α.40 sample.

We consider whether extreme redshift distortions nearby, for small values of cz_CMB, could be contaminating our results. To test this possibility, we repeat the measurement, this time excluding galaxies within cz_CMB < 2000 km s⁻¹, as well as 3000 km s⁻¹. This reduced the sample size to ∼9300 and ∼8900, respectively, but we found no difference in the final fitting parameters r₀ and γ. We conclude that there is no advantage to be gained in eliminating nearby galaxies from α.40 for the correlation function analysis.

The expression for J₃ in Equation (5) requires an expression for the shape of the correlation function. In calculating ξ(σ, π)), there is therefore a presumably small dependence on the as-yet-unknown parameters s₀ and γ. As briefly mentioned in Section 3.1, one possible way to avoid any potential problems is to use the fiducial optical sample parameters s₀ = 5.0 and γ = −1.8 to calculate the parameters for an H i-selected sample, and then iterate toward a stable solution. In attempts to do this, we find that there is no significant difference between the parameters estimated via these two methods, and confirm that ξ is not dependent on the precise form of J₃. Such iteration does not provide any advantage. We demonstrate this in Figure 7, which displays the error contours on the power-law fit parameters for a sample limited to cz_CMB > 2000 km s⁻¹, and for J₃ using the approximate parameters estimated in Section 4.1. The results are very close to those in Figure 5 and the 1σ estimated parameters are identical.

Zoom In Zoom Out Reset image size

The correlation function of gas-rich galaxies has implications for the improvement of galaxy simulations, by providing an observational constraint for the results of simulations. This work will allow a better match between simulations and the observed relationship between gas mass and clustering properties. Simulations are just now progressing to the point where reasonable, realistic cold H i gas masses can be assigned to galaxies. In this section, we will compare the results of the correlation analysis of α.40 with presently available cold dark matter simulations.

We are limited in our ability to compare our observations to simulations by what is available publicly. Martin et al. (2010) took advantage of the Obreschkow et al. (2009, hereafter O09) simulation, which assigned cold gas to galaxies from the De Lucia & Blaizot (2007) catalog of Millennium Simulation galaxies. In that work, we found that O09's simulation provided a reasonable fit to the observed H i mass function. However, this catalog may not be adequate for comparison with an observed correlation function. In particular, O09 caution that the mass resolution of the simulation prevents them from applying their findings to faint, LSB, or low-mass galaxies. Given the known correlations between galaxy type and clustering, and between H i mass and luminosity, it would be difficult to use this catalog to explore the relationship between current simulations and ALFALFA's observations. Furthermore, O09 did not themselves carry out this analysis.

Kim et al. (2011) have provided another option for comparison, using a set of four GALFORM semi-analytical models that treats a range of processes which influence gas reservoirs, including cooling, ram-pressure stripping, mergers, star formation, and supernova feedback. They report results over a range of redshifts, but for this work only their results at z = 0 are of interest. They find that the galaxy–galaxy correlation function of the simulations is consistent with those found for HIPASS, and confirm that their simulation shows gas-rich galaxies as being significantly less clustered than dark matter. Differences between the models and the scales at which those differences are important can be used to highlight potential problems in the assumptions, such as models that overpredict the gas richness of satellites.

In Figure 8, we compare the models presented in Kim et al. (2011) with the observed α.40 correlation function for H i-selected galaxies. The models include that of Bower et al. (2006), labeled as Bow06; a modified version of the same, labeled MHIBow06; a version that uses a slightly different background cosmology that is in better agreement with the Wilkinson Microwave Anisotropy Probe (WMAP) parameters, labeled GpcBow06; and, finally, the model of Font et al. (2008), labeled as Font08. In all models, only galaxies with M_cold > 10^9.5 h⁻² M_☉, where M_H i = 0.76 M_cold/(1 + .04), are included, which matches the HIPASS galaxy selection but may be more massive than would be ideal for matching α.40, which probes very small gas masses (Martin et al. 2010; Haynes et al. 2011). The models are described in detail in Kim et al. (2011), and here we only discuss the main differences that may be relevant for a comparison to α.40.

Zoom In Zoom Out Reset image size

Bow06, MHIBow06, and Font08 all use the Millennium Simulation to track galaxies and halos, while GpcBow06 uses a different method involving merger trees and a large box size. Bow06 and Font08 are able to match optical luminosity functions, but both overpredict the abundance of H i in low-mass galaxies. MHIBow06 was created by adjusting the star formation timescale in Bow06, thereby fixing this excess while maintaining the agreement with optical properties of galaxies. GpcBow06, finally, also has a modified star formation prescription which better fits the H i mass function compared to Bow06.

In Figure 8, it is clear that MHIBow06 and Bow06 fit the observed H i correlation function on small scales, while Font08 drastically overpredicts and GpcBow06 drastically underpredicts the strength of clustering for gas-rich galaxies on those scales. At large scales, both GpcBow06 and Bow06 underpredict the clustering strength, while Font08 and MHIBow06 follow it closely. Although not a perfect match, the MHIBow06 model appears to be most consistent with the clustering of gas-rich galaxies over the full range of accessible scales.

Part of these differences may be due to the mass resolution of the models, given that α.40 probes to significantly lower masses than the HIPASS survey for which these models were designed. What is clear, however, is that α.40 already provides constraints that can begin to differentiate between successful and unsuccessful models, and the full ALFALFA sample should be able to provide very robust constraints for testing simulations that take H i into account.

The bias between any two classes of objects indicates their relative clustering strength. For cosmological purposes, we are interested in comparing the clustering of types of galaxies with the underlying dark matter halo distribution, in order to understand how well future surveys would probe the true (baryonic + dark) mass distribution. The comparison is achieved through the linear bias parameter at z = 0 b₀:

$$\begin{equation} \xi _{{\rm gal}}(r) = b_{0}^2 \: \xi _{{\rm DM}}(r). \end{equation} \tag{ 14 }$$

In general, b₀ = b₀(r), based on the linear theory prediction that b₀ is independent of scale. In real galaxies, however, this is expected to become true only above intermediate scales of ∼10 h⁻¹ Mpc. If b₀ > 1.0, as is the case for red galaxies which tend to be found in clusters (e.g., Guzzo et al. 1997; Norberg et al. 2002; Zehavi et al. 2005; Li et al. 2006; Swanson et al. 2008), then the distribution is positively biased with respect to the dark matter. For galaxies like H i-selected populations, b₀ < 1.0 and they are said to be antibiased.

As a proxy for the underlying dark matter distribution at z = 0, we use the correlation function of dark matter halos from the Millennium Simulation, given in Springel et al. (2005) as a function over the same scales that we are interested in. The Millennium Simulation, however, used an early WMAP estimated value for the parameter σ₈ of ∼0.9 that is generally recognized to have been high. We have adjusted our calculations to use the recommended value σ₈ = 0.8 from the seven-year WMAP (WMAP7) results (Larson et al. 2011).

In Figure 9, we compare that correlation function to the α.40 observation of ξ(r) for H i-selected galaxies using the inversion method. The dark matter correlation function deviates strongly from a power law on small scales, an indication of the well-known fact that bias is scale dependent. Dark matter is, as expected, significantly more strongly clustered than this particular population of galaxies, but the bias becomes less significant at intermediate and large scales.

Zoom In Zoom Out Reset image size

In Figure 10, we display the bias parameter as a function of scale. The error bars are based only on the α.40 uncertainties and they assume that there is no uncertainty in the Millennium Simulation's measurement of the correlation function. This figure reflects what we already understand about the clustering properties of H i selected galaxies: on small scales, the clustering of gas-rich galaxies is weaker, and on ever-larger scales the distribution of gas-rich galaxies begins to more closely reflect the underlying matter distribution. Basilakos et al. (2007) measured the linear bias parameter on large scales for the HIPASS sample using a different technique. Exploring modeled dark matter power spectra for different values of b₀, including bias and assuming a concordance cosmology, they identified the most likely bias parameter. They found b₀ = 0.7 ± 0.1, in general agreement with the findings for α.40 though we have measured the bias parameter as a function of scale. The preliminary work of Passmoor et al. (2011) is generally consistent with this result, though that work does not capture our finding that the sample becomes unbiased on large scales.

Zoom In Zoom Out Reset image size

Because of the very limited redshift extent of ALFALFA and the α.40 sample, our work cannot comment on the evolution of clustering or bias for H i-selected objects. However, as a robust measurement of these properties at z = 0, we do provide a benchmark observational constraint with which theoretical models will need to agree. The earliest work comparing the bias of H i at z = 0 to the anticipated evolution, by Basilakos et al. (2007), determined b₀ ∼ 0.68 today and predicted b₄ would range from ∼2 to 4 by z = 4. Recently, Marín et al. (2010) have used a different, simple bias model, incorporating observational constraints, which relates M_H i to M_DM to estimate H i masses of Millennium Simulation halos, and then investigated the bias of H i with respect to the halo distribution. At z = 0, they estimate that the overall linear bias parameter on large scales is ∼0.8. Their Figure 6 is more comparable to our Figure 10, and shows the same overall rise of the bias with increasing scale found here. Their models also predict that the bias will rise sharply with redshift, with the linear bias parameter reaching b₄ ∼ 2 by z ∼ 4.

The α.40 observations and Marín et al. (2010) predictions have implications for large-scale 21 cm galaxy surveys and intensity mapping projects with such instruments as the SKA. If the theoretical results reflect the true evolution of H i gas in the universe, then these projects can expect strong 21 cm signals at a range of redshifts. Perhaps more importantly, the α.40 observation of the correlation function at low redshift provides a robust baseline constraint for the development of SKA model predictions, and future simulations will need to match both the H i mass function and the correlation function at z ∼ 0.

Because the H i-selected galaxy bias is likely to be strongly dependent on H i mass, with low-mass objects severely antibiased with respect to the underlying dark matter distribution, high-redshift surveys which are sensitive only to the high-mass end of the HIMF should expect to be mildly antibiased at low redshifts and increasingly positively biased at intermediate to high redshifts. We will explore the mass, color, and luminosity dependence of the correlation function as the focus of a future work.